スライド 1

Fortran と有限差分法の
入門の入門の…
高橋芳幸
目次
•
•
•
•
Fortran 入門 1-5
偏微分方程式の数値解法 1-2
差分の諸々
まとめ
Fortran 入門 1
• 表示
-------------------------------------program test1
implicit none
write( 6, * ) “Hello world!”
end program test1
--------------------------------------
“おまじない”
画面に出力
Fortran 入門 2
• 計算
-------------------------------------program test2
implicit none
real(8) :: num1, num2, num3, num4
num1 = 3.0
num2 = 5.0
num3 = 6.0
num4 = num3 * (num1 + num2)
write( 6, * ) num4
end program test2
--------------------------------------
“おまじない”
変数宣言
値の設定
計算 ( 6×( 3 + 5 ) )
画面に出力
Fortran 入門 3
• 繰り返し
-------------------------------------program test3
implicit none
“おまじない”
real(8) :: num
integer :: i
変数宣言
num = 0
do i = 1, 10
num = num + i
end do
（繰り返し）計算
write( 6, * ) num
画面に出力
end program test3
--------------------------------------
( 1 + 2 + 3 + … + 10 )
Fortran 入門 4
• 配列
-------------------------------------program test4
implicit none
real(8) :: num(5)
integer :: i
“おまじない”
変数（配列）宣言
do i = 1, 5
num( i ) = i * 5
end do
配列に値の設定
do i = 1, 5
write( 6, * ) i, num( i )
end do
画面に出力
end program test4
--------------------------------------
Fortran 入門 5
• 条件分岐（と入力）
-------------------------------------program test5
implicit none
“おまじない”
real(8) :: num
変数宣言
write( 6, * ) ‘Input number’
read( 5, * ) num
キーボードから読み込み
if( num > 10 ) then
write( 6, * ) ‘input number is greater than 10’
else
write( 6, * ) ‘input number is less than or equal to 10’
end if
条件分岐（画面に出力）
end program test5
--------------------------------------
偏微分方程式の数値解法 1
題材とする方程式
  x  

 2 x
t
 ( x,0)  5  2 sin( x)
初期条件
は有限差分法（のとある方法）を用いると次
のように離散化される.
 nj   ( x j , tn ), x j  jx, tn  nt
 nj 1  nj
 2 x j
t
 nj1   nj  2 x j t
 0n
 n1
 02
 21
 11

x1  x
0
1
0
 1n
 12
 11
 2n
 22
 21
x1  x x2  2x
 nj
…
…
…
 2j
 1j
… x j  j x
tn  nt
t2  2t
t1  t
偏微分方程式の数値解法 1
プログラム例
program partdiff1
implicit none
real(8) :: pi, x(100), dx, dt, psi(100), psi_n(100)
integer :: t, i
pi = acos( -1.0d0 )
dt = 0.2
dx = 2.0 * pi / 100
do i = 1, 100
x( i ) = -pi + dx * ( i – 1 )
end do
do i = 1, 100
psi( i ) = 5.0 + 2.0 * sin( x( i ) )
end do
do i = 1, 100
write( 6, * ) 0, x( i ), psi( i )
end do
write( 6, * )
write( 6, * )
演算はどうすれば良いでしょう?
end program partdiff1
πの設定
座標の設定
初期値の設定
初期値の出力
偏微分方程式の数値解法 1
プログラム例
program partdiff1
implicit none
real(8) :: pi, x(100), dx, dt, psi(100), psi_n(100)
integer :: t, i
pi = acos( -1.0d0 )
dt = 0.2
dx = 2.0 * pi / 100
do i = 1, 100
x( i ) = -pi + dx * ( i – 1 )
end do
do i = 1, 100
psi( i ) = 5.0 + 2.0 * sin( x( i ) )
end do
do i = 1, 100
write( 6, * ) 0, x( i ), psi( i )
end do
write( 6, *)
write( 6, *)
演算はどうすれば良いでしょう?
end program partdiff1
do t = 1, 10
do i = 1, 100
psi_n( i ) = psi( i ) - 2 * x(i) * dt
end do
do i = 1, 100
psi( i ) = psi_n( i )
end do
do i = 1, 100
write( 6, * ) t, x( i ), psi( i )
end do
write( 6, * )
write( 6, * )
end do
計算
 nj1   nj  2 x j t
変数の入れ替え
結果の出力
偏微分方程式の数値解法 1
実行例
プログラム名を main.f90 とする.
% ifort main.f90
% a.out > result.txt
コンパイル
結果をファイルにリダイレクト
偏微分方程式の数値解法 1
gnuplot を用いた描画方法の例
% gnuplot
gnuplot 起動
GNUPLOT
Version 4.0 patchlevel 0
last modified Thu Apr 15 14:44:22 CEST 2004
System: Linux 2.6.18-5-686
Copyright (C) 1986 - 1993, 1998, 2004
Thomas Williams, Colin Kelley and many others
（省略）
Terminal type set to 'x11'
gnuplot>
gnuplot> plot “result.txt” index 0 u 2:3 w l
gnuplot> replot “result.txt” index 1 u 2:3 w l
0 番目の時間の 2, 3 カラムのデータを描画
1 番目の時間の 2, 3 カラムのデータを重ね描き
偏微分方程式の数値解法 1
数値解の例
数値解は次のようになるだろう.
この場合は解析解も同じだろう.
偏微分方程式の数値解法 2
移流方程式
• ここでは例として以下に示す 1 次元移流方程
式を離散化を扱う.


 c
t
x
(1)
ただし, ここでは c は定数で, 0  x  2 とし, 以下の境界条件（周
期境界条件）を与える.
 (0, t )   (2 , t )
偏微分方程式の数値解法 2
離散化例 1
• (1)式は, 有限差分法（のとある方法）を用いると以下のよう
に離散化される.
 nj   ( x j , tn ), x j  jx, tn  nt
 nj 1  nj
t

n 1
j
 c
 nj1  nj1
2x
t
  c
 nj1  nj1
2x
n
j


ただし, この方法は数値的に安定ではありません.
（安定性については, 必要な人は各自調べてください.）
偏微分方程式の数値解法 2
プログラミング練習
プログラムはどうなるでしょう?
偏微分方程式の数値解法 2
数値解の例
例えば,
△x = 2 pi / 100, △t = 0.1
c = 0.1
の場合の数値解は以下のようになるだろう.
偏微分方程式の数値解法 2
離散化例 2
• 有限差分法による離散化には様々な（精度の）方法
がある.
– 先の例では
n
n
  j 1  j 1

x
2 x
としたが,
n
n
  j  j 1

x
x
n
n
n
n
  j  2  8 j 1  8 j 1  j  2

x
12x
としてもよい.
差分の諸々
導出 1
•
ここでは以下の差分式を導出する.
–  をテイラー展開すると
n
n
  j 1  j 1

x
2 x



 x 3  O x 4
x



 x 3  O x 4
x
1   2
  
 x  x    x   
 x    2
2  x
 x  x

1   3
2
 x    3
6  x
x
1   2
  
 x  x    x   
 x    2
2  x
 x  x

1   3
2
 x    3
6  x
x
となり, この 2 式より,
3
    x  x   x  x  1   
  3

 
2x
3  x
 x  x


 


 

 

 x 2  O x 4
x
となる.
2
つまり, この差分式は x  の誤差が生じるので, 2 次精度となる.

差分の諸々
導出 2
• 同様に式変形することで,
n
n
  j  j 1

x
x
n
n
n
n
  j  2  8 j 1  8 j 1  j  2

x
12x
の差分式はそれぞれ, 1 次精度, 4 次精度で
あることがわかる.
• 時間についての差分式も同様.
差分の諸々
プログラム実習
• 興味のある人は, 精度の異なる差分式を導
出し, プログラムを作って見てください.
まとめ
• Fortran 入門（の入門の入門…）
• 偏微分方程式の数値解法入門（の入門の入
門… ）

Download Report