とらりもんHOME  Index  Search  Changes  Login

微分方程式を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

bibun1.png

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

1次元移流方程式の数値計算

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

iryu.png

風上差分を使います。

sabun1.png

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

sabun2.png

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

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

adv_anime.gif

Last modified:2020/03/14 17:00:31
Keyword(s):
References:[fortranを用いた数値計算]