【1次元熱伝導方程式】Fortranからgnuplotでグラフの１次元アニメーション作成

こんにちは（@t_kun_kamakiri）(‘◇’)ゞ

Fortranを学生時代に使っていたのですが、最近になって使い始めています。
完全に趣味です(^^)/

それはそうと、Fortranは数値計算する上では直感的でとてもわかりやすい言語だと個人的には思っているのですが、めんどくさい点があるのです。

それは、Fortranってグラフ作成などの可視化ツールがない・・・・(_;

他のプログラミング言語に詳しくないから何とも言えないですが、Pythonとか対話型のプログラミング言語なら可視化のライブラリが用意されているのに・・・

Fortranでよく使うグラフ作成ツールと言えば、gnuplot（ニュープロット）かなと思っています。

カマキリ

gnuplotはタダだし・・・

そこで次の疑問が・・・

疑問

1次元グラフをどうやってgnuplotでアニメーションするのか？

正直、↑これは学生の時からやり方わからなくて放置していたのですが、社会人になって改めて考えるとあっさりできたので記事にしておこうと思います。

※Fortranのコードとスクリプトは本記事にアップされています。
※FortranのコードはFortran90/95の書き方にしているつもりですが正確に書けているかは不明です。
※FortranはWSLを使ってインストールしています。
※gnuplotはこのサイトを参考にすれば良いかと思います。

僕の使用環境はWindows10でWindows Subsystem for Linux（WSL）を使って、UbuntuからFortranやgnuplotをインストールして使用しています。

カマキリ

説明がくどいな～って思った方は「まとめ」へGo！！

1次元熱伝導方程式をFortranで解く
全体の流れ
アニメーションを作成してみよう(^^)/
まとめ

1次元熱伝導方程式をFortranで解く

まずは解くべき方程式を明らかにしておきましょう。

1次元熱伝導方程式

\begin{align*}\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial^2 x}\tag{1}\end{align*}

熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。

今回は初期状態を以下のようにしておきます。

中央だけ中央だけ温度が高く、両端は温度０(K)の条件となっています。

ですので、この初期状態から時間発展させると高い温度は低い温度へ熱移動が生じるので、下記のような変化を示します。

これをgnuplotを使ってアニメーションを作ろうということになります。

全体の流れ

アニメーション作成までの全体の流れを押さえておくと少し気持ちのハードルは下がるかもしれないので、先に全体の流れを示しておきます。

1次元熱伝達方程式は偏微分方程式であるため、それを数値的に解くためには離散化してあげる必要があります。
今回は、「時間に対しては陽解法」「空間に対しては2次精度」で解くことを考えています。
次にForranを使って離散化した方程式を数値的に解くことを考えます。
（※コードは後ほど）
Fortranで出力したデータをgnuplotで順番に読み込んでいくためのスクリプトをFortranを使って作成します。

これらのことができれば、あとはgnuplotで「3.」で作成されたスクリプトを実行するだけアニメーションを作成することができます。

※1次元熱伝導方程式を題材にしたのは、簡単な方程式の方がアニメーション作成の方法を理解しやすいと思ったからです。
ゆえに、離散化も特にこだわってはいません。

1次元熱伝導方程式の離散化

1次元熱伝導方程式

\begin{align*}\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial^2 x}\tag{1}\end{align*}

この偏微分方程式を離散化するにあたって、

時間微分の項（左辺）
空間微分の項（右辺）

をどのように扱うのかが問題です。

いっちばん手っ取り早い離散化で進めていきましょう。

ある地点(x_{i})での時間微分の項

\begin{align*}\frac{\partial T}{\partial t}\simeq \frac{T^{n+1}(x_{i})-T^{n}(x_{i})}{\Delta t}\tag{2}\end{align*}

微分を差分に置き換えただけです。

数値計算では時間をステップ数で数えるため、時間変数の代わりにステップ\(n\)として文字の頭に添えています。

ある時間\(t\)での空間微分の項

\begin{align*}\kappa\frac{\partial^2T}{\partial^2 x}\simeq \kappa\frac{T^{n}(x_{i}+\Delta x)-2T^{n}(x_{i})-T^{n}(x_{i}-\Delta x)}{\Delta x^2}\tag{3}\end{align*}

これは空間に対して2次の精度となります。

2次の精度とは関数をテーラー展開してときに何字の項まで残したかというものです。

【補足説明です】

テーラー展開については下記の記事をご参考ください。

前回の記事はこちら

テーラー展開とその利用方法

\begin{align*}T(x_{i}+\Delta x)=T(x_{i})+\frac{\partial T}{\partial x}|_{(x_{i})}\Delta x+\frac{1}{2}\frac{\partial^2T}{\partial x^2}|_{(x_{i})}\Delta x^2+o(\Delta ^3)\end{align*}

\begin{align*}T(x_{i}-\Delta x)=T(x_{i})-\frac{\partial T}{\partial x}|_{(x_{i})}\Delta x+\frac{1}{2}\frac{\partial^2 T}{\partial x^2}|_{(x_{i})}\Delta x^2+o(\Delta ^3)\end{align*}

この両方の差から、\(\frac{\partial^2 T}{\partial x^2}\)の2次精度(3)式を得ることができます。

補足説明終了です。

(2)(3)式より、数値計算で解くべき離散化の式がわかりました（↓以下です。）

\begin{align*}T^{n+1}(x_{i})=T^{n}(x_{i})+\kappa dt\frac{T^{n}(x_{i}+\Delta x)-2T^{n}(x_{i})-T^{n}(x_{i}-\Delta x)}{\Delta x^2}\tag{4}\end{align*}

では、これをFortranで解きましょう(^^)/

Fortranで解く

Fortranで解くと言ってもプログラムを書かないことには何も始まらないので、プログラムを書きます。

以下がプログラムです。

!ファイル名　1D_Thermal_Eq.f90
module constants
  implicit none
  integer::i,nt
  integer, parameter :: imax = 100,ntmax = 2500
  real(8):: T(0:imax),xmax,x,dx,dt,kp
  integer,parameter :: SP = kind(1.0)
  integer,parameter :: DP = selected_real_kind(2*precision(1.0_SP))
end module constants

module print
  use constants
  implicit none
 
 contains
  subroutine print__profile(f)
    real(DP), dimension(0:imax), intent(in) :: f
    integer :: counter = 0
    character(len=*), parameter :: base = "./data/temp.j=middle."
    character(len=4) :: serial_num
    character(len=40) filename
  
    write(serial_num,'(i4.4)') counter

    write(filename,"(a,i5.5,a)") base//serial_num
    open(10,file=filename)
    do i = 0 , imax
       write(10,*) i*dx, f(i)
    end do
    close(10)
    counter = counter + 1
    
  end subroutine print__profile
end module print

program thermal_explicit
  use constants
  use print

  implicit none

  dt   = 0.004
  kp   = 1
  xmax = 10
  dx   = xmax/imax
  
  open(10,file="initial.txt")
  open(11,file="final.txt")

   !initial state
   do i = 0, imax
      x = dx*i
      if (i < imax/2) then
         T(i) = x
       else
         T(i) = xmax - x
      endif
      write(10,*) i,x,T(i)
   enddo
   call print__profile(T(:))

   
   !Time step
   do nt = 1, ntmax
      do i = 1, imax-1
         T(i) = T(i) + kp*dt/(dx*dx)*(T(i+1)- 2*T(i) +T(i-1))
      end do
      if (mod(nt,100)==0) call print__profile(T(:))
      write(*,*)nt
   end do
   
   !Final state
   do i = 0, imax
      x = dx*i
      write(11,*) i,x,T(i)
   end do
end program thermal_explicit

!ファイル名　1D_Thermal_Eq.f90

module constants

implicit none

integer::i,nt

integer, parameter :: imax = 100,ntmax = 2500

real(8):: T(0:imax),xmax,x,dx,dt,kp

integer,parameter :: SP = kind(1.0)

integer,parameter :: DP = selected_real_kind(2*precision(1.0_SP))

end module constants

module print

use constants

implicit none

contains

subroutine print__profile(f)

real(DP), dimension(0:imax), intent(in) :: f

integer :: counter = 0

character(len=*), parameter :: base = "./data/temp.j=middle."

character(len=4) :: serial_num

character(len=40) filename

write(serial_num,'(i4.4)') counter

write(filename,"(a,i5.5,a)") base//serial_num

open(10,file=filename)

do i = 0 , imax

write(10,*) i*dx, f(i)

end do

close(10)

counter = counter + 1

end subroutine print__profile

end module print

program thermal_explicit

use constants

use print

implicit none

dt = 0.004

kp = 1

xmax = 10

dx = xmax/imax

open(10,file="initial.txt")

open(11,file="final.txt")

!initial state

do i = 0, imax

x = dx*i

if (i < imax/2) then

T(i) = x

else

T(i) = xmax - x

endif

write(10,*) i,x,T(i)

enddo

call print__profile(T(:))

!Time step

do nt = 1, ntmax

do i = 1, imax-1

T(i) = T(i) + kp*dt/(dx*dx)*(T(i+1)- 2*T(i) +T(i-1))

end do

if (mod(nt,100)==0) call print__profile(T(:))

write(*,*)nt

end do

!Final state

do i = 0, imax

x = dx*i

write(11,*) i,x,T(i)

end do

end program thermal_explicit

プログラムだけ載せられてもよくわからんという人もいるかもしれないので、説明を加えておきます(‘ω’)

↑変数の定義をモジュールとしてまとめておきます。

↑ここが今回の最も重要な部分で、アニメーションを作るためには各ステップ（出力したいステップ）のデータファイルを順に出力する必要があります。

そのためのサブルーチンを書いています。

↑ここからがメインのプログラム部分です。

ここではメインプログラムで使うための変数を指定しています。

↑ここで初期状態を作成しています。

↑ここで(4)式の離散化した式を解いています。

100ステップごとに温度の空間分布を出力するようにif文で制御しています。

↑ちなみに「data」という名前のディレクトリを自分で作っておく必要があります。

以上がプログラムの中身です。

ここまでできたら、以下のコマンドで計算を実行させます。

gfortran 1D_Thermal_Eq.f90
./a.out

1 2	gfortran 1D_Thermal_Eq.f90 ./a.out

実行ができれば「data」という名前のディレクトリの中に、以下のように大量のデータが保存されていることが確認できます。

これらのデータが指定したタイムステップで出力された温度の空間分布のデータです。

あとは、これらのデータを順番にgnuplotで読み込めばアニメーションが完成ということになります。

そのためには、gnuplotのスクリプトが必要になります。

スクリプトは書くのが少々面倒なので、スクリプト自体をFortranで書くことにします。

アニメーション用のgnuplotのスクリプトをFortranで作成

gnuplotのスクリプトを作成するためのFortranのプログラムがこちらです。

program animation_script
  implicit none
  integer, parameter :: imax = 100,xmax = 10
  integer :: counter, counter_end = 24,dx
  character(len=*), parameter :: base='./data/temp.j=middle.'
  character(len=4)            :: serial_num
  
  dx   = xmax/imax

  print *, "#"
  print *, "# gnuplot script generated by movie_script.f90"
  print *, "#"
  print *, " "
  print *, "set xlabel 'x'              # x-axis"
  print *, "set ylabel 'temperature'    # y-axis"   
  print *, "set xrange[0:", xmax, "]     # j-grid min & max"
  print *, "set yrange[0.0:10.0]         # temperature min & max"
  print *, "set grid"
  print *, "set xlabel font 'Arial,20' "
  print *, "set ylabel font 'Arial,20' "
  print *, "set tics font 'Arial,20'"

  do counter = 0 , counter_end
     write(serial_num,'(i4.4)') counter
     print *, "plot '"//base//serial_num//"' w lp"
     
     if ( counter==0) then 
        print *, "pause 2"
     else
        print *, "pause 1"
     end if
  end do
  print *, "pause -1"
end program animation_script

program animation_script

implicit none

integer, parameter :: imax = 100,xmax = 10

integer :: counter, counter_end = 24,dx

character(len=*), parameter :: base='./data/temp.j=middle.'

character(len=4) :: serial_num

dx = xmax/imax

print *, "#"

print *, "# gnuplot script generated by movie_script.f90"

print *, "#"

print *, " "

print *, "set xlabel 'x' # x-axis"

print *, "set ylabel 'temperature' # y-axis"

print *, "set xrange[0:", xmax, "] # j-grid min & max"

print *, "set yrange[0.0:10.0] # temperature min & max"

print *, "set grid"

print *, "set xlabel font 'Arial,20' "

print *, "set ylabel font 'Arial,20' "

print *, "set tics font 'Arial,20'"

do counter = 0 , counter_end

write(serial_num,'(i4.4)') counter

print *, "plot '"//base//serial_num//"' w lp"

if ( counter==0) then

print *, "pause 2"

else

print *, "pause 1"

end if

end do

print *, "pause -1"

end program animation_script

これはあまり説明を必要としないのですが、簡単に説明すると「do文」のところで先ほどの「data」という名前のディレクトリ内の大量のデータを順番に読み込む作業を行っています。

では、以下のコマンドを実行してみます。

gfortran movie_script.f90
./a.out >movie

1 2	gfortran movie_script.f90 ./a.out >movie

「./a.out >movie」とするとmovieというファイルが出来上がり、その中にFortranの「print文」の内容を全て書き出してくれます。

この「movie」というファイルがgnuplotのスクリプトとなります。

以上のことが終えれば、あとはgnuplotで「movie」という名前のスクリプトを走らせるだけです。

アニメーションを作成してみよう(^^)/

「movie」というスクリプトは2通りの方法で実行できます。

方法１

↑このようにコマンドを打っても良いですし、

方法2

↑このように「gnuplot」と打って、「load “movie”」とスクリプトファイルをloadしても良いです。

アニメーションはこんな感じ

カマキリ

やったー

まとめ

はじめての人にもできるように、わかりやすくまとめてみました(‘◇’)ゞ

が、説明がくどく長い記事になってしまったことをお許しください(‘ω’)

まだまだFortranもgnuplotも勉強不足ですが、何とかできたかなっていう感じです。

ごちゃごちゃ説明してしまいましたが、下記のzipファイルを解凍して下記のコマンドを打ってお試しください。

20190428_thermal1D

gfortran 1D_Thermal_Eq.f90
./a.out
gfortran movie_script.f90
./a.out >movie
gnuplot movie

gfortran 1D_Thermal_Eq.f90

./a.out

gfortran movie_script.f90

./a.out >movie

gnuplot movie

1次元熱伝導方程式をFortranで解く

全体の流れ

1次元熱伝導方程式の離散化

Fortranで解く

アニメーション用のgnuplotのスクリプトをFortranで作成

アニメーションを作成してみよう(^^)/

まとめ

【2次元熱伝導方程式】Fortranからgnuplotでグラフの３次元アニメーション作成

【マンデルブロー集合で学ぶFortran入門】繰り返し(do文)と制御構文(if文)

【バーンズリーのシダで学ぶFortran入門】乱数生成

Fortranで解いた問題をParaviewを使ってvtkでアニメーション作成

【マンデルブロー集合で学ぶFortran入門】変数とデータ

1次元の移流方程式を前進差分、後退差分、中心差分で解く。【Fortranコードあり】

【グロス・ピタエフスキー方程式（Fortranコードあり）】2次元ダークソリトンの崩壊

COMMENT コメントをキャンセル