とらりもんHOME  Index  Search  Changes  Login

とらりもん - 微分方程式をfortran で解く Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

! 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 の詳細については[[こちら|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

http://ryuiki.envr.tsukuba.ac.jp/~tsurumi/torari/bibun1.png

解析解と数値解には若干のずれがありますね。

!1次元移流方程式の数値計算
さきほどは、1変数の微分方程式を解きました。今度は時間的、空間的に変化する、こんな微分方程式を解いてみましょう。

http://pen.envr.tsukuba.ac.jp/~tsurumi/torari/iryu.png

風上差分を使います。

http://ryuiki.envr.tsukuba.ac.jp/~tsurumi/torari/sabun1.png

これを式(1)に代入するとこうなって

http://ryuiki.envr.tsukuba.ac.jp/~tsurumi/torari/sabun2.png

t=0のときのf(x,0)を与えます。
今回は適当に三角波で与えました。

http://ryuiki.envr.tsukuba.ac.jp/~tsurumi/torari/triangle.png

さっそくプログラムを組んでみましょう。
プログラム例はこちら:

$ 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

http://ryuiki.envr.tsukuba.ac.jp/~tsurumi/torari/adv_anime.gif