プログラム forward.pに 注釈をつけながら説明する。 Fortran、 Cの例文も同じ考え方で作ってある。 ただし、このプログラムは 1996年度東京都立大学ワークステーション教室のDEC Pascalで動作を 確認したが、他の機種やコンパイラでは変更が必要かもしれない。
ここでは、z (時刻 t+Δtでの変数Zの値)と、 zold (時刻 t での変数Zの値)を 別々の(Pascalなどの言語の立場での)変数にしてみた。 前方差分だけに限れば、この2つを同じ変数にしてしまっても問題ない (なぜか考えてみよう)。しかし、後方差分その他への変更のしやすさを考えて、 区別しておくことにする。
PROGRAM decay(output);
{ Simple decay
dZ/dt = aZ
by forward difference
written by K. Masuda, 1994-05-31, revised 1995-04-26, 1996-06-03
}
＊まず変数の宣言。
Pascalでは、変数はすべて宣言する。
VAR
a, deltat, t, tzero, z, zold, zzero, ztrue: DOUBLE;
実数型変数はいずれも倍精度で宣言した。
kt, nt: INTEGER;
＊以下実行部。
BEGIN
＊以下はPascal文法上は変数への代入文だが、実質的には定数を与えている。
(注意: 定数のあとについている「d0」は、倍精度の定数であることを示す
DEC Pascalコンパイラでの約束である。Pascal言語一般に通用するものでは
ない。DEC Pascalでは、これをつけないと、単精度の定数になってしまう。)
. 方程式の係数 a
a := -1.0d0;
. 時間発展の開始時刻
tzero := 0.0d0;
. Zの初期値 (時刻 tzero での値)
zzero := 1.0d0;
. 時間ステップ (解析解では結果の出力間隔)
deltat := 0.1d0;
. 終了までの時間ステップ数
nt := 100;
＊対象とする問題や計算の条件を示す定数の値を書き出しておく。
gnuplotでコメント扱いするため、各行の頭に「#」を入れている。
writeln手続きのカッコの中は、いずれも書き出す変数あるいは定数である。
なお、数値の後にコロンでくぎってつけてある数字は、書式指定で、
1つめは数値を表わす全体の桁数、2つめは小数部の桁数を示す。
writeln( '#', 'Problem dZ/dt = aZ');
writeln( '#', 'Forward difference');
writeln( '#', 'a ',' ', a:20:10 );
writeln( '#', 'tzero ',' ', tzero:20:10 );
writeln( '#', 'zzero ',' ', zzero:20:10 );
writeln( '#', 'deltat',' ',deltat:20:10 );
writeln( '#', 'nt ',' ', nt:20 );
writeln( '#', 'tfinal',' ',(tzero+nt*deltat):20:10 );
＊シミュレーションの最初の時刻(tzero)での変数の値(初期値)の設定と出力。
t := tzero;
z := zzero;
ztrue := zzero;
. 各列の見出しも、「#」つきで書き出しておく。
writeln( '#','time ',
'simulated ',
'true ' );
writeln( ' ', t:20:10, z:20:10, ztrue:20:10 );
. 時刻 Δt での計算の準備として、z の値を zold に入れる
zold := z;
＊時間発展の計算。ステップ数 nt だけ、ループ内容を繰り返す。
FOR kt := 1 TO nt DO
BEGIN
. これから計算する時刻
t := tzero + kt * deltat;
. これが前方差分の1ステップの計算
z := zold + deltat * a * zold;
. この問題は真の解がわかっているので計算しておく
ztrue := zzero * exp( a * (t - tzero) );
. 値を書き出す
writeln( ' ', t:20:10, z:20:10, ztrue:20:10 );
. 時刻 t (次のステップでは時刻 t - Δt になる)の z の値を zold に入れる
zold := z;
END;
＊終わり。
END.