差分法の全般的な説明の中で、 「前方差分はΔtの1次まで正しい、台形公式はΔtの2次まで正しい、...」 などと述べた。これは、差分法の式を、(無限回微分可能などいくつかの 条件のもとで)正確なテイラー(Taylor)展開の式と比較して、理論的に 導いたものである。
では、実際に差分計算をしてみると、Δtを変えるに従って誤差は どのように変わるだろうか。
前に行なった理論的な誤差の見積りでは、 無限に続くテイラー展開をある項数で 打ち切ることによる「打ち切り誤差」だけを考慮している。 実際の計算では、これに、 実数を有限桁数の浮動小数点数で近似することからくる「まるめ誤差」も 複合してくるので、 結果は必ずしもきれいに対応しないかもしれない。
問題としては、前と同様、指数関数型の減衰または増大現象を表わす微分方程式
d Z / d t = a Zをとりあげる。
前に述べた前方差分の例文プログラムforward.{p,f,c} (Fortran、 Pascal、 Cにそれぞれ説明をつけた)は、 ステップ数(nt)とステップの長さ(deltat)のひと組の値について計算し、 毎ステップ、結果を書き出していた。
今回の例文 deltat.p, deltat.f, deltat.c では、 計算の対象となる時間の長さ(tfinal = nt * deltat)を 一定にして、それをいくつのステップに分けるか(nt)を変えながら 何度も時間発展型計算を行ない、 最終時刻(t = tfinal)での解析解からの誤差を、nt, deltatとともに 書き出すようにする。 ntの値はキーボードから入れることにし、 ntにゼロまたは負の値を入れたら計算の終わりの印ということにする。 つまり、論理的に次のような構造になる。
deltat.{p,f,c}の例文では、差分の方式は前方差分を採用している。 内側のループの内容は、forward.{p,f,c}の例文と同じである。 ただし、ループの中での結果の書き出し(writeln, WRITE, printf)は 省略されている。
これを実行してみよう。 (オプション-lmはC言語の場合だけ必要。 ) 結果はあとでグラフにするため、画面に出さずにファイルに入れることにする。 結果を入れるファイル名は例のとおりでなくてもよいが、減衰問題 (指数関数型の問題)の前方差分であることがわかるようにしておく。 プログラム実行中に入れるntの値は、もちろんこのとおりでなくてもよい。
(deltat.c を用意する) % cc deltat.c -lm % a.out > deltat.gensui.zenpoo 10 20 30 50 100 200 300 500 1000 -1出力結果のファイルを見て、 数値が3列にならんでいることを確認してほしい。 それぞれ、nt, deltat, 誤差(z - ztrue)である。
% less deltat.gensui.zenpoo q ( ... lessを終わる)それをグラフにしてみよう。
% gnuplot gnuplot> plot "deltat.gensui.zenpoo" using 2:3これで、横軸にΔt、縦軸に誤差をとったグラフになる。 この例の場合、グラフを見ると、 Δtが0.1以下の部分で、誤差がΔtに対して 直線的に変化しているようである。この特徴は横軸の範囲を設定して グラフをかきなおしてみるとはっきりする。
gnuplot> set xrange [0:0.1] gnuplot> replotさらに、縦軸の0がグラフの中に表示されるように、次のようにしたほうが 見やすいだろう。
gnuplot> set yrange [:0] gnuplot> replot
このようにして見ると、この例題の場合、微分を前方差分で近似した ことによる誤差は、確かにΔtの1乗に比例しているようである。 ただし、Δtをもっと小さくしていくと、 (計算の時間がかかるという問題もあるが)、 結果を得るまでに積み重なる演算の回数が多くなるので、 まるめ誤差が無視できなくなり、今得られた直線の延長よりも 大きな誤差が出てくると予想される。 ただし、ntをどのくらい大きくするとまるめ誤差の効果が見えてくるかは まだ確認していない。
同じ微分方程式を、前方差分以外の差分方式で近似した場合の誤差が、 Δtを変えるに従ってどのように変わるか、 また、それが理論的に見た打ち切り誤差のΔt依存性と 対応しているか、ためしてみてほしい。 forward.{p,f,c}の例文の場合と同様に、deltat.{p,f,c}の例文でも、 時間ステップのループの中の、zの値を変えている代入文1つを書きかえれば よい。ただし、差分方式によっては代入文や宣言文を追加する必要があるかも しれない。
なお、今回の例文では、t = tfinalでの最終結果だけを書き出すように している。これは、差分による誤差が、時間発展とともに徐々に単調に (正または負の方向に絶対値が)増加する場合には適した方法なのだが、 誤差が時間とともに振動する場合には、たまたま見かけ上小さくなった ところを拾ってしまうおそれもある。 その場合は、複数のステップでの誤差を書き出しておき、その絶対値の 最大値を計算するなどのくふうが必要になるだろう。 差分方式の誤差の出かたを理解する上では、最終結果だけを見るのでは なく、forward.{p,f,c}の例文のように時間発展の経過を書き出し、 数値やグラフの形で見てみることも重要である。