時間発展問題の数値計算

差分法

○目次

注:

○ 対象として考える問題

微分方程式で表わされるような問題を、差分で近似することを考える。 ここでは、指数関数型減衰(あるいは増大)問題や、あとで扱う単振動問題を含む、 やや一般化された対象を考える。 次のように、量 Z の時間微分(時間あたりの変化)が、同じ時刻の Z 自体の 関数であるとしよう。 

dZ / dt = g(Z)                                          (1)
ここで g(Z) は、単に、「 Z の関数である」という意味である。 f をここの Z の意味で使っている文献が多いので、文字を変えた。 特に、 g(Z) = - a Z とすれば、減衰問題になる。 以下、 g(Z) のことを、 Z の時間微分であることを強調して、Zdot とも書く。 これは、Zの上に点を打った記号を想定している。

時刻 t1 での Z の値がわかっていたときに、 時刻 t1 + Δt での Z の値を求めることにする。 (以下の式では、時刻t1でのZの真の値から出発するかのように説明している。 実際の時間発展の計算では、真の値ではなく、初期時刻から差分で計算して きた結果得られた時刻t1でのZの値を使って1ステップを進めるのである。)

○Taylor (テイラー) 展開

もし、時刻 t1 での Z の t による1階微分 (d Z / d t) (t = t1) = Zdot(t1) だけではなく、 2階微分 (d2 Z / d t2) (t = t1) 、 3階微分 (d3 Z / d t3) (t = t1) ... もわかっていたとすれば、時刻 t1 + Δt での Z の値は 

Z(t1 + Δt) = 
Z(t1) +         Δt   (d   Z / d t) (t = t1)
       + 1/(2!) Δt2  (d2 Z / d t2) (t = t1)
       + 1/(3!) Δt3  (d3 Z / d t3) (t = t1) ...           (2)
という形で求めることができる。(2)のような形をTaylor展開という。 正確には、右辺の項が無限個続いた極限が左辺と等しくなるのだが、 Δt が ある程度小さければ、特殊な場合を除いて、右辺の項の絶対値は急速に小さくなる ので、求めたい有効数字を考慮して、適当な項数で打ち切って計算すればよい。

○前方差分 (forward difference; 前進差分、Euler (オイラー) 法)

しかし、ここでは1階微分しか計算できない場合を考えている。 そこで、式(2)を Δt の1次までで打ち切った 

ZF(t1 + Δt) = Z(t1) + Δt Zdot(t1)              (3)
のような形で求められる ZF を Z の近似として考える。 式(3)は、 
{ZF(t1 + Δt) - Z(t1)}
------------------------- = Zdot(t1)                (4)
           Δt
と書いても等式としては同じである。ただし、(3)では、求めたいものが左辺に あり、FortranやPascalなどのプログラム言語の代入文に近い形をしているので、 以下おもにこの形で議論する。

Δt → 0 の極限では、 ZF は Z に一致する。 つまり、差分式(4)は、Δt → 0 の極限では、微分式(1)に一致するのである。 しかし実際には、時間刻み Δt を無限に小さくはできないので、 有限の値のときの精度も問題にしなければならない。 差分式の精度は、Taylor展開(2)と比較して見積もることができる。 式(3)の場合、 ZF と Z の差は、 

  1          d2 Z
  - (Δt)2  ----- (t = t1) + ((Δt) の3次以上の項 ) 
  2          d t2
である。つまり、前方差分式は、Δt の1次まで微分式と一致する。

次の図では、横軸に t 、縦軸に Z をとって前方差分の概念を示した。

曲線 Z(t) を、 t = t1 での接線で置きかえて延長し、 t = t1 + Δt での Z の近似値 ZFを求めたことになる。

また、次の図では、 横軸に t 、縦軸に Zdot をとって示した。

左(または上):真の解 (true solution):::::: 右(または下):前方差分 (forward difference)

:::

この図では、曲線 Zdot(t) と、横軸、 t = t1、 t = t1 + Δt に囲まれた 面積が、 Z(t1 + Δt) - Z(t1) に相当する。 前方差分は、曲線の下側の面積を、Zdot = Zdot(t1) という水平の直線の 下側の長方形の面積で近似することに相当する。

t1 + Δt での Z の近似値が求まれば、今度はその時刻から Δt 後、つまり t1 + 2 Δt での値を同様な方法で求めることができる。 このようにして、初期時刻 t0 での Zの値が与えられれば、順に 時間ステップを重ねて計算をしていくことができる。 

    Z ---> Z ---> Z ---> Z ---> Z
     /    \ /    \ /    \ /    \ /
    +------+------+------+------+---------> t
    t0   t0+1dt t0+2dt t0+3dt t0+4dt
ただし、 差分が微分の近似として完全でないため、正しい微分の解からはしだいに離れて いくことになる。

○後方差分 (backward difference、後退差分) と 台形公式 (trapezoidal formula)

平均変化率を区間の初めの瞬間変化率で代用するのは不公平とも考えられよう。 代わりに、区間の終わり t1 + Δt での瞬間変化率で代用すれば、 

ZB(t1 + Δt) = Z(t1) + Δt ZdotB(t1 + Δt)           (5)
のような形が得られる。 ただし、ZdotB(t1 + Δt)は、Zdot = g(Z)の Z に ZB(t1 + Δt)を 代入したものである。 (5)の形を後方差分という。 Taylor展開の式(2)で、t1、Δt とあるところをそれぞれ t1 + Δt 、-Δt と置きかえたものと比較すれば、 後方差分も Δt の1次まで正しい微分の差分近似であることがわかる。 次の図では、 t = t1 + Δt での接線と平行な(同じ傾きをもつ)直線で、 t1 から t1 + Δt までの Z の変化を近似している。

また、次の図では、 t = t1 + Δt での Zdot の高さの長方形で曲線の下の面積を 近似している。

また、区間の両端を公平に考慮するために、算術平均をとれば、 

ZT(t1 + Δt) = 

Z(t1) + (1/2)Δt { Zdot(t1) + ZdotT(t1 + Δt) }
                                                          (6)
のような形が得られる。

次の図を見ればわかるように、これは数値積分でいう台形公式である。

図でも長方形の場合よりも曲線に近づいていることがわかるが、Taylor展開との 比較をしてみると、台形公式は Δt の2次まで正しい微分の近似式になっている (確かめてみてほしい)。

○陽公式(explicit scheme) と陰公式 (implicit scheme)

今考えている状況は、時刻 t1 での Z 、 Zdot がわかっていて、 時刻 t1 + Δt での Z を求めるものである。

前方差分の式(3)は、右辺はすべて既知の量なのでそのままFortranなどの代入文 に置きかえることができる。 このような式を、陽公式(陽解法、explicit scheme)という。 (未知数が複数あっても、それぞれが既知の数で表わされていれば陽公式である)。

一方、後方差分(5)や、台形公式(6)では、右辺に時刻 t1 + Δt での Zdotが はいっているが、これは、関係式 Zdot = g(Z) で左辺と結びついた未知数で ある。(そのため、式の中では、下つき添え字 B 、T をつけておいた)。 このように、仮に左辺に1つの未知数を出しても、右辺に未知数が残ってしまう ような式、つまり複数の未知数を連立で求める必要がある式を、陰公式 (陰解法、implicit scheme)という。

陽公式、陰公式は、差分の問題に限ったものではない。 たとえば、 (x, y) 平面上の円周の式は、ふつう 

x2 + y2 = r2                                  (7)
と書かれるが、これは x から y 、あるいは y から x を求める形を していないので陰公式である。 
y1 = + sqrt(r2 - x2)
y2 = - sqrt(r2 - x2)                          (8)
は陽公式である。[ sqrt() は平方根の記号を表わす ]

○予測修正法

具体的な g(Z) の形にもよるが、陰公式を数値計算可能な形にするのは むずかしいことが多い。そこで、厳密に表現すると陰公式になってしまうもの を、近似して陽公式の形に持ちこむ必要が生じる。その代表的方法が、 各時間ステップごとに、まず陽公式で(t1 + Δt)での仮の値を求めておき、 陰公式の右辺で(t1 + Δt)での値を必要とするところではその仮の値で 代用することである。言いかえると、まず陽公式で予測し、それを陰公式を 使って修正することになっているので、「予測・修正法」という。

まず、「前方差分で予測した仮の値を、後方差分で修正する」 ことを考えよう。 後方差分の式(5)の右辺にある未知の量 ZdotB(t1 + Δt) を、前方差分で求められた ZF(t1 + Δt) を Zdot = g(Z) に代入して得られる ZdotF(t1 + Δt) で代用する。つまり、まず(3)を計算したあと、 

ZFB(t1 + Δt) = Z(t1) + Δt ZdotF(t1 + Δt) (9)
を計算するのである。(記号は「前方・後方」の略でFBとした。) 計算結果は、陰公式である本来の後方差分とやや数値は違うが、 誤差の傾向は似ており、Δt の1次まで正しいという意味でも同じである。

同様に、「前方差分で予測した仮の値を、台形公式で修正する」 こともできる。 台形公式 式(6)の右辺にある未知の量 ZdotT(t1 + Δt) を、前方差分で求められた ZF(t1 + Δt) を Zdot=g(Z) に代入して得られる ZdotF(t1 + Δt) で代用することが考えられる。つまり、まず(3)を計算したあと、 

ZH(t1 + Δt) = Z(t1) 
  + (1/2) Δt { Zdot(t1) + ZdotF(t1 + Δt) }   (10)
を計算するのである。 この方式をHeun(人名)方式ともいうので、記号はHとした。 精度は、陰公式である本来の台形公式に比べて、定量的にはやや悪いが、 Δt の2次まで正しいという意味では同じである。

○Runge-Kutta法

予測・修正を繰りかえすことによって、 陽解法でありながら精度のよい公式が作られている。 有名なものとして、Runge-Kutta法がある。 Runge, Kuttaはそれぞれ人名で、この名前は一連の公式の総称だが、 そのうちでよく使われるのは次の4次のRunge-Kutta法 (Δtの4次まで正しい)である。 1ステップごとに、Zの1回微分の仮の値を4回計算する。 ここでは詳しい説明は省略する。 興味のある人は、数値計算の教科書の常微分方程式の章を参照してほしい。 

Z0 = Z( t1 )
G0 = g(Z0, t1)
Z1 = Z0 + (Δt/2) G0
G1 = g(Z1, t1 + (1/2) Δt)
Z2 = Z0 + (Δt/2) G1
G2 = g(Z2, t1 + (1/2) Δt)
Z3 = Z0 + Δt G2
G3 = g(Z3, t1 + Δt)
ZR( t1 + Δt ) = Z0 + (Δt/6) ( G0 + 2 G1 + 2 G2 + G3 )   (11)

○中央差分 (centered difference)

前方差分の不公平さを避けるため、台形公式のように両端を平均する方法のほかに、 区間のほうを t1 - Δt から t1 + Δt までに取りなおすことが考えられる。 

ZC(t1 + Δt) = Z(t1 - Δt) + 2 Δt Zdot(t1)         (12)
これは(次に述べる条件を満たすならば)陽公式であり、 Taylor展開(2)と比較すると Δt の2次まで正しい。

:::

しかし、これまでにあげた差分公式と違うのは、 t1 と t1 - Δt の 2つの時刻での Z の値を必要とすることである。 最初に2つの時刻 t0 と t0 + Δt での値がわかっていれば、 次の図のように、順次計算していくことができる。 

    Z ----------> Z ----------> Z
     /   dZ/dt   \ /   dZ/dt   \ /
    +------+------+------+------+---------> t
            \   dZ/dt   / \   dZ/dt
           Z ----------> Z ----------> 
    t0   t0+1dt t0+2dt t0+3dt t0+4dt
つまり、Δtの偶数倍の時刻のZを求めるには、その2Δt前(偶数倍)のZと、 1Δt前(奇数倍)のZdotを、Δtの奇数倍の時刻のZを求めるには、奇数倍の 時刻のZと偶数倍の時刻のZdotを使う。 (このような時間発展計算の方法をleap-frogということがある。直訳すれば カエル(蛙)とびであるが、日本でいう「うまとび」に相当し、 「うま」になる人ととぶ人の役割が交互に代わることをさしているらしい。)

時間発展問題では多くの場合、初期値としてある瞬間の値だけがわかっている ことが多いので、中央差分を使う場合は、 Δt だけ離れたもう1つの 瞬間の値を、前方差分などを使って推定しなければならない。 また、問題の性質によっては、中央差分で計算すると、偶数番めと奇数番め のステップの値が相互に離れていってしまい、微分方程式の近似解が 得られないことがある。 こういった問題点は、例題に関連して述べる予定である。

○参考文献

別記事参照。 特に、「常微分方程式」を扱っている、 高橋(1996)5章、川上(1989)6章、戸川(1991)第9章。

ただし、わたしの教材の用語は、気象学関係の習慣に従っているので これらの教科書にある概念を別の名前をつけているところがある。 (「修正オイラー法」とか「改良オイラー法」のようなまぎらわしい 表現を避けて、なるべく内容を直接表現する用語を選んだつもりである。)

気象学の人がよく使う文献も、別記事にあげておいた。

1996-06-16; 2000-04-18; 2001-02-13; 2002-07-11
増田 耕一 (MASUDA Kooiti)

[この文書の最初へ] [目次へ]