プログラム forward.cに 注釈をつけながら説明する。 Fortran、 Pascalの例文も同じ考え方で作ってある。
ここでは、z (時刻 t+Δtでの変数Zの値)と、 zold (時刻 t での変数Zの値)を 別々の(Cなどの言語の立場での)変数にしてみた。 前方差分だけに限れば、この2つを同じ変数にしてしまっても問題ない (なぜか考えてみよう)。しかし、後方差分その他への変更のしやすさを考えて、 区別しておくことにする。
/*
Simple decay
dZ/dt = aZ
by forward difference
written by K. Masuda, 1994-05-23, 1996-05-29, 1996-06-03
*/
＊コンパイラとともに用意されたヘッダファイルのとりこみ
(printfなどを使う準備)
#include
(expなどを使う準備)
#include
int main(){
＊変数の宣言。
変数はすべて宣言する。
実数型変数は、倍精度を使うことにする。
double a, deltat, t, tzero,
z, zold, zzero, ztrue;
int kt, nt;
＊以下実行部。
＊以下はC文法上は変数への代入文だが、実質的には定数を与えている。
(注意: C言語では実数型定数は特に指定しなくても倍精度となる。)
. 方程式の係数 a
a = -1.0;
. 時間発展の開始時刻
tzero = 0.0;
. Zの初期値 (時刻 tzero での値)
zzero = 1.0;
. 時間ステップ (解析解では結果の出力間隔)
deltat = 0.1;
. 終了までの時間ステップ数
nt = 100;
＊対象とする問題や計算の条件を示す定数の値を書き出しておく。
gnuplotでコメント扱いするため、各行の頭に「#」を入れている。
printf関数のカッコの中の1つめの文字列の中に含まれた「%」で始まる
文字列は、書式指定である。%20.10lfの20は数値を表わす全体の桁数、
10は小数部の桁数を示す。
printf("#Problem dZ/dt = aZ\n");
printf("#Forward difference\n");
printf("# a %20.10lf\n",a);
printf("# tzero %20.10lf\n",tzero);
printf("# zzero %20.10lf\n",zzero);
printf("#deltat %20.10lf\n",deltat);
printf("# nt %20d\n",nt);
printf("#tfinal %20.10lf\n",tzero+nt*deltat);
＊シミュレーションの最初の時刻(tzero)での変数の値(初期値)の設定と出力。
t = tzero;
z = zzero;
ztrue = zzero;
. 各列の見出しも、「#」つきで書き出しておく。
printf("# time simulated true\n");
printf("%20.10lf%20.10lf%20.10lf\n", t, z, ztrue);
. 時刻 Δt での計算の準備として、z の値を zold に入れる
zold = z;
＊時間発展の計算。ステップ数 nt だけ、ループ内容を繰り返す。
for(kt = 1; kt <= nt; kt++)
{
. これから計算する時刻
t = tzero + kt * deltat;
. これが前方差分の1ステップの計算
z = zold + deltat * a * zold;
. この問題は真の解がわかっているので計算しておく
ztrue = zzero * exp( a * (t - tzero) );
. 値を書き出す
printf("%20.10lf%20.10lf%20.10lf\n", t, z, ztrue);
. 時刻 t (次のステップでは時刻 t - Δt になる)の z の値を zold に入れる
zold = z;
}
＊終わり。
}