プログラム forward.f に注釈をつけながら説明する。 Pascal、 Cの例文も同じ考え方で作ってある。
ここでは、z (時刻 t+Δtでの変数Zの値)と、 zold (時刻 t での変数Zの値)を 別々の(Fortranなどの言語の立場での)変数にしてみた。 前方差分だけに限れば、この2つを同じ変数にしてしまっても問題ない (なぜか考えてみよう)。しかし、後方差分その他への変更のしやすさを考えて、 区別しておくことにする。
PROGRAM decay
C
C Simple decay
C dZ/dt = aZ
C by forward difference
C written by K. Masuda, 1994-05-23
＊まず変数の宣言。
IMPLICIT NONE
これは、変数はすべて宣言する、という指定である。つまり、宣言しない変数を
使うとコンパイル段階でエラーになるようにしておく。
REAL*8 a, deltat, t, tzero,
& z, zold, zzero, ztrue
実数型変数はいずれも倍精度で宣言した。
INTEGER kt, nt
＊以下はFortran文法上は変数への代入文だが、実質的には定数を与えている。
. 方程式の係数 a
a = -1.0d0
. 時間発展の開始時刻
tzero = 0.0d0
. Zの初期値 (時刻 tzero での値)
zzero = 1.0d0
. 時間ステップ (解析解では結果の出力間隔)
deltat = 0.1d0
. 終了までの時間ステップ数
nt = 100
＊対象とする問題や計算の条件を示す定数の値を書き出しておく。
WRITE文のカッコの中
. 1つめの要素はファイル参照番号(星印なら標準出力)。
. 2つめは、FORMAT(=書式)文の番号(星印なら変数の型ごとに標準の書式を仮定)。
gnuplotでコメント扱いするため、各行の頭に「#」を入れている。
FORMAT文の要素 Aは文字型、Iは整数型、Fは実数型変数用の書式。
WRITE (*,600) 'Problem dZ/dt = aZ'
WRITE (*,600) 'Forward difference'
WRITE (*,610) 'a',a
WRITE (*,610) 'tzero',tzero
WRITE (*,610) 'zzero',zzero
WRITE (*,610) 'deltat',deltat
WRITE (*,620) 'nt',nt
WRITE (*,610) 'tfinal',tzero+nt*deltat
600 FORMAT('#',A)
610 FORMAT('#',A6,1X,F20.10)
620 FORMAT('#',A6,1X,I20)
＊シミュレーションの最初の時刻(tzero)での変数の値(初期値)の設定と出力。
t = tzero
z = zzero
ztrue = zzero
. 各列の見出しも、「#」つきで書き出しておく。
WRITE (*,700) 'time', 'simulated', 'true'
700 FORMAT('#',3A20)
WRITE (*,710) t, z, ztrue
710 FORMAT(' ',3F20.10)
. 時刻 Δt での計算の準備として、z の値を zold に入れる
zold = z
＊時間発展の計算。ステップ数 nt だけ、ループ内容を繰り返す。
DO..ENDDOは、Fortran77標準ではないが、90標準に含まれ、教室でも使える。
DO kt = 1, nt
. これから計算する時刻
t = tzero + kt * deltat
. これが前方差分の1ステップの計算
z = zold + deltat * a * zold
. この問題は真の解がわかっているので計算しておく
ztrue = zzero * exp( a * (t - tzero) )
. 値を書き出す
WRITE (*,710) t, z, ztrue
. 時刻 t (次のステップでは時刻 t - Δt になる)の z の値を zold に入れる
zold = z
END DO
＊終わり。
STOP
END