減衰( a < 0 )または増大( a > 0 )現象を表わす微分方程式
d Z / d t = a Zの数値解を、差分法を、前方差分のほか、
グラフにすること、誤差の数値を出してみることなどの方法は、 前方差分の場合の記事を参照。
前方差分の例文のプログラムが、 前方差分の記事で紹介してあり、 詳しい説明は Fortran、 Pascal、 C、それぞれ別の記事に書いてある。
差分法を変える場合、このプログラムの大きな骨組み(時間ステップ による繰り返し)は変えず、ループの中の演算を変えればよい。 ただし、そのために変数を追加する必要があるかもしれない。 陽公式の場合は、すなおにプログラムの代入文に翻訳できると思う。 陰公式の場合、一般の g(z) について代入文に翻訳することはできないが、 今の例題( g(z) = a z )に限っては、代数的な式変形で、 未知数である時刻(t1 + Δt)でのzに関する値を左辺にまとめて、 「z(t1 + Δt) = ...」の形にすることができるので、陽公式と 同様に代入文にすることができる。
もう少し具体的に述べる。 前方差分の計算の主要部分は (プログラム forward.fでいえば)
z = zold + deltat * a * zold (1)
である。ここで、zold には時刻 t - Δt (前の説明では t1) での値が
はいっており、z として時刻 t (前の説明では t1 + Δt) での値が計算される。
後方差分の場合は、
z = zold + deltat * a * z (2)
と書きたいところであるが、Fortran(など手順だて型の言語)の代入文でこう
書くと、右辺の z は、計算結果(左辺の z)ではなく、過去に z という変数に
はいっていた値が使われてしまう。等式としての(2)をみたすような z を求める
ためには、z が左辺だけにくるように変形してやればよい。つまり、
z - deltat * a * z = zold (これはFortranの文ではない)
(1 - deltat * a) * z = zold ( " )
z = zold / (1 - deltat * a) (3)
のように変形してやれば、代入文として実行可能な形にできた。
ただし、この変形は、d Z / d t = a Z という問題だから簡単にできたので
あって、一般の d Z / d t = g(Z) では困難である。
そこで、後方差分の近似としてつぎのような方式を考えることがある。
zforw = zold + deltat * a * zold (4a)
z = zold + deltat * a * zforw (4b)
つまり、(2)式で両辺に出てくる z の値のうち、右辺のものだけを、前方差分で
求めた時刻 t での z の値で代用するのである。言いかえると、前方差分の結果
を後方差分の考え方で修正しているとも言えるので、「前方予測・後方修正型」
の差分方式と呼ぶことができる。(気象シミュレーションで、前方差分を使うと
振動の振幅が非現実的に大きくなってしまうとき、よく使われている。)
台形公式は、(2)と同様に代入文として正しくない表現をすれば
z = zold + deltat * a * (z + zold) / 2.0 (5)
である。代入文への変形は各自考えてほしい。
中央差分は、Δtの2次まで正しい精度がありながら、前方差分とほぼ同じ 計算量ですむので、よく使われている。
中央差分の場合は、時刻 t での値 z を求めるのに、時刻 t - Δt での値 (zmidとする)と、時刻 t - 2Δtでの値(zoldとする)の両方を必要とする。 差分の主要部は
z = zold + 2 * deltat * a * zmid (6)
である(2がかかっていることに注意)。
また、時間発展の計算を始めるにあたって、時刻 t_0 および t_0 + Δtの 2つの時点での値が必要となる。 t = Δt での値を求める第1ステップは、中央差分の公式どおりに計算する ことができない。その計算方法としては、次のようなものが考えられる。
中央差分の例文を、 center.p, center.f, center.c として用意する。 これでは、最初のステップに前方差分を使っている。 実際これを実行してみると、結果はめちゃくちゃ(正しい微分方程式の解と 全く対応しない)となる。これは、奇数番めと偶数番めのステップの値が だんだん離れていってしまうためである。 言いかえると、周期2Δtの振動の振幅が大きくなる。
このように、中央差分を使うのには注意が必要である。 特に、減衰問題(結果が指数関数になる問題) では、中央差分は避けるべきで、誤差が大きくなっても、 前方差分や後方差分のほうがよいとされている。(精度を要求するときは、 台形公式のたぐいや、予測・修正法を使う。) また、振動問題でも、奇数ばんめと偶数ばんめの情報が分離するのを避ける ため、中央差分に前方差分をいくらか混合することがある。
対策の一つとして、中央差分と前方差分をある比率で混ぜることがある。 例文を、 mixed.p, mixed.f, mixed.c として用意する。 混ぜる比率は仮に 9:1 とした。これで、あばれることは防げているが、 まだ多少、周期2Δtの計算上の振動が残っている。 うまく混ぜると、前方差分だけよりも誤差が小さくなるので、実際、気象の 数値予報にはこのような混合方式が使われることが多い。