とらりもん - 微分方程式を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
* 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