微分方程式をfortran で解く
1変数の微分方程式を解く
- du/dz=z (I.C. u(0)=3)
簡単に前進差分で差分化しちゃいます。 差分については各々で勉強しましょう。。
$ vi bibun1.f90 ! bibun1.f90 ! du/dz=z を初期条件u(0)=3 で解く ! 2008/11/19 tsurumi parameter (NX=10) dimension u(0:NX) dz=0.1 u(0)=3 do i=0,NX-1 z=i*dz du=z*dz u(i+1)=u(i)+du end do do i=0,NX write(*,*)i*dz,u(i) end do end
計算できましたか?この計算結果を.txt ファイルにして、解析解u(z)=z^2/2+3 と比べてみます。 グラフはgnuplot を使って描きます。gnuplot の詳細についてはこちら。
$ ./bibun1>>bibun1.txt $ gnuplot gnuplot> set title 'comparison between Suuchikai and Kaisekikai' gnuplot> set xlabel 'z' gnuplot> set ylabel 'u(z)' gnuplot> set size 4,4 gnuplot> set key left top gnuplot> plot 'bibun1.txt' ps 5 pt 7 title 'suuchikai' gnuplot> replot 'kaiseki1.txt' ps 5 pt 7 title 'kaisekikai' gnuplot> set terminal postscript eps enhanced color 90 'bold' gnuplot> set output 'bibun1.eps' gnuplot> replot gnuplot> q $ convert bibun1.eps bibun1.png $ eog bibun1.png
解析解と数値解には若干のずれがありますね。
1次元移流方程式の数値計算
さきほどは、1変数の微分方程式を解きました。今度は時間的、空間的に変化する、こんな微分方程式を解いてみましょう。
風上差分を使います。
これを式(1)に代入するとこうなって
t=0のときのf(x,0)を与えます。 今回は適当に三角波で与えました。
さっそくプログラムを組んでみましょう。 プログラム例はこちら:
$ vi adv.f90 ! adv.f90 ! 1次元移流方程式を風上差分で解く ! compile gfortran adv.f90 -o adv ! output ./adv ! 2007/11/14 TSURUMI Yuki parameter (N=99) dimension f(0:N),fn(0:N) data cp,xp,dt,dx,u /0.5,5,0.1,1.0,0.5/ do i=2,N-1 x=dx*(i-1) if(i.lt.xp) then f(i)=cp/xp*x else if(i.lt.2*xp) then f(i)=-cp/xp*x+2*cp else f(i)=0 end if end do f(1)=0 do i=1,N write(10,*)i,f(i) end do do kk=11,100 do k=1,20 fn(1)=f(1) do i=2,N-1 fn(i)=(1-u*dt/dx)*f(i)+u*dt/dx*f(i-1) end do do i=1,N f(i)=fn(i) end do end do do i=1,N write(kk,*)i,f(i) end do end do end
さらにこれを可視化しちゃいます。
$ vi adv.sh
#!/bin/sh i=10 while test $i -le 99;do gnuplot<<EOF set xrange [0:80] set yrange [0:1] plot 'fort.${i}' w l lw 10 set terminal png medium set output 'adv_${i}.png' replot EOF i=`expr $i + 1` done convert -delay 10 adv_*.png adv_anm.gif firefox adv_anm.gif
シェルを組んだら実行権限を与えてから実行しましょう!! 実行権限?
$ chmod +x adv.sh $ ./adv.sh
Keyword(s):
References:[fortranを用いた数値計算]