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 “おまじない” 変数宣言 （繰り返し）計算 ( 1 + 2 + 3 + … + 10 ) 画面に出力 繰り返し -------------------------------------- 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 write( 6, * ) i, num( i ) 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 題材とする方程式 は有限差分法（のとある方法）を用いると次のように離散化される. 初期条件 … … … … … 時間 … …

偏微分方程式の数値解法 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 psi( i ) = 5.0 + 2.0 * sin( x( i ) ) write( 6, * ) 0, x( i ), psi( i ) 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 psi( i ) = 5.0 + 2.0 * sin( x( i ) ) write( 6, * ) 0, x( i ), psi( i ) write( 6, *) 演算はどうすれば良いでしょう? end program partdiff1 do t = 1, 10 do i = 1, 100 psi_n( i ) = psi( i ) - 2 * x(i) * dt end do psi( i ) = psi_n( i ) write( 6, * ) t, x( i ), psi( i ) write( 6, * ) 計算 変数の入れ替え 結果の出力

偏微分方程式の数値解法 1 実行例 プログラム名を main.f90 とする. % ifort main.f90 % a.out > result.txt コンパイル 結果をファイルにリダイレクト

偏微分方程式の数値解法 1 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 gnuplot 起動 0 番目の時間の 2, 3 カラムのデータを描画 1 番目の時間の 2, 3 カラムのデータを重ね描き

偏微分方程式の数値解法 1 数値解の例 数値解は次のようになるだろう. この場合は解析解も同じだろう.

偏微分方程式の数値解法 2 移流方程式 ここでは例として以下に示す 1 次元移流方程式を離散化を扱う. (1) ただし, ここでは c は定数で, 0  x  2 とし, 以下の境界条件（周期境界条件）を与える.

偏微分方程式の数値解法 2 離散化例 1 (1)式は, 有限差分法（のとある方法）を用いると以下のように離散化される. ただし, この方法は数値的に安定ではありません. （安定性については, 必要な人は各自調べてください.）

偏微分方程式の数値解法 2 プログラミング練習 プログラムはどうなるでしょう?

偏微分方程式の数値解法 2 数値解の例 例えば, △x = 2 pi / 100, △t = 0.1 c = 0.1 の場合の数値解は以下のようになるだろう.

偏微分方程式の数値解法 2 離散化例 2 有限差分法による離散化には様々な（精度の）方法がある. 先の例では としたが, としてもよい.

差分の諸々 導出 1 ここでは以下の差分式を導出する. をテイラー展開すると となり, この 2 式より, となる. つまり, この差分式は の誤差が生じるので, 2 次精度となる.

差分の諸々 導出 2 同様に式変形することで, の差分式はそれぞれ, 1 次精度, 4 次精度であることがわかる. 時間についての差分式も同様.

差分の諸々 プログラム実習 興味のある人は, 精度の異なる差分式を導出し, プログラムを作って見てください.

まとめ Fortran 入門（の入門の入門…） 偏微分方程式の数値解法入門（の入門の入門… ）