10. 常微分方程式 (1) 概要, 1段法, 多段法, 使われなくなった解法

目次

常微分方程式 (1): 概要, 1段法, 多段法, 使われなくなった解法
常微分方程式 (2): 誤差の推定, 補外法, ステップ幅の自動調節, 密出力, 2階常微分方程式, 遅延微分方程式
常微分方程式 (3): 解法の安定性, スティフな方程式, 微分代数方程式
常微分方程式 (4): 実用ルーチンのベンチマーク (実用ルーチンの選択, スティフではない問題, スティフな問題, 微分代数方程式)

概要1段法多段法使われなくなった解法

常微分方程式の初期値問題について

初期値問題

1階の常微分方程式は次の形に書くことができる.

  
  y' = f(t, y(t))

 
ここで, y(t) は未知の関数である. 常微分方程式の初期値問題(initial value problem)は, 初期点 t0 における関数値 y(t0) が与えられたとき, 最終点 tf における値 y(tf)を求める問題である. 初期点における関数値を与える次式は初期条件(initial condition)とよばれる.

  
  y(t0) = y0

 

1階の連立方程式

1階の常微分方程式の初期値問題は, スカラー y をベクトル y に変えてそのまま連立方程式に拡張することができる.

  
  y' = f(t, y)
  y(t0) = y0

 
例えば, 2変数(ベクトルの要素数が2)の場合について書き下すと次のようになる.

  
  y1' = f1(t, y1, y2)
  y2' = f2(t, y1, y2)

 
この場合, 初期条件も2つ必要になる.

  
  y1(t0) = y1,0
  y2(t0) = y2,0

 

高階の常微分方程式

2階の常微分方程式は1階の連立常微分方程式の形にして扱うことができる. 例えば, 次の方程式を考える.

  
  y'' = f(t, y, y')

 
これは, y1 = y, y2 = y’ とおけば, 次のように1階の連立常微分方程式になる.

  
  y1' = y2
  y2' = f(t, y1, y2)

 
この場合, 2つの初期条件が必要になる.

  
  y1(t0) = y(t0) = y0
  y2(t0) = y'(t0) = y'0

 
同様にしてさらに高階の常微分方程式を解くこともできる.

常微分方程式の初期値問題の数値解法

数値解法の分類

変数 t を離散化して, t0, t1, … , tn, … と表す. これらの点の間隔をステップ幅あるいはステップ幅(step size)といい, hn で表す.

  
  hn = tn+1 - tn

 
ステップ幅は必ずしも一定でなくてよいが, n によらず一定の場合には単に h と表すことがある. 未知関数 y(t) の tn における値 y(tn) の近似値を yn と表すことにする.

差分法(微分を差分で近似する方法)により数値解法を導いてみる. 代表的な差分の求め方には, 次のようなものがある.

  • 前進差分: (yn+1 – yn) / h
  • 後退差分: (yn – yn-1) / h
  • 中心差分: (yn+1 – yn-1) / 2h

まず, 前進差分を使って微分を近似する.

  
  y'n = f(tn, yn) ≒ (yn+1 - yn) / h

 
これより次の近似式が得られる.

  
  yn+1 = yn + hf(tn, yn)

 
この式を使い, 初期値 y0 から始めて, y1, y2, … と次々に各点における近似解を求めることができることがわかる. これをオイラー法 (Euler method) という.

後退差分を使った場合には, 同様にして次式により近似解を求めることができる.

  
  yn+1 = yn + hf(tn+1, yn+1)

 
これを後退オイラー法 (backward Euler method) という.

後退オイラー法では右辺の f() の中に yn+1 が含まれており, yn+1 を求めるために方程式を解かなければならない(一般的には非線形方程式になる). このような解法を陰解法 (陰的公式, implicit formula) という. これに対して, オイラー法の場合には方程式を解くことなく直接 yn+1 を計算することができた. そのような解法を陽解法 (陽的公式, explicit formula) という.

次に, 中心差分を使った場合には, 次式により近似解を求めることができる.

  
  yn+1 = yn-1 + 2hf(tn, yn)

 
これは中点則とよばれる.

中点則では, yn+1 の計算には yn の他に過去の値 yn-1 も必要で, 2つの値を使うので2段法とよばれ, 多段法 (多段階法, multistep method) とよばれる解法のひとつになる. これに対して, オイラー法と後退オイラー法は yn だけしか使わないので, 1段法 (1段階法, one-step method) とよばれる解法に分類される.

まとめると, 常微分方程式の初期値問題の解法を大きく分類すると1段法と多段法があり, それぞれに陽解法と陰解法がある. ここでは, 1段法の陽解法, および, 多段法について説明する. 1段法の陰解法についてはスティフな方程式の項で説明する.

数値解法の次数

関数 y(t) (真の解) を tn においてテイラー展開すると次のようになる.

  
  y(tn+1) = y(tn) + hy'(tn) + (1/2)h2y''(tn) + (1/6)h3y'''(tn) + ... + (1/p)hpy(p)(tn) + ...
          = y(tn) + hf(tn, y(tn)) + (1/2)h2f'(tn, y(tn))+ (1/6)h3f''(tn, y(tn)) + ... + (1/p)hpf(p-1)(tn, y(tn)) + ...

 
数値解法による近似解 yn+1 の右辺が上のテイラー展開 y(tn+1) のp次の項まで一致するとき, すなわち次式が成り立つときにその数値解法はp次の方法であるという.

  
  yn+1 - y(tn+1) = O(hp+1)

 
オイラー法はテイラー展開の1次(h)の項までで打ち切ったものに等しいから明らかに1次の方法である. また, 後退オイラー法は1次, 中点則は2次の方法であることを確かめることができる.

数値実験 (1)

数値解法の次数
最も簡単な微分方程式の初期値問題の例である.

  
  y' = y, y(0) = 1

 
解析解は y = et である. 上の3つの公式を用いてステップ幅を 0.1, 0.05, 0.01, … というように小さくしていきながら順に計算する. t = 1 において得られた近似解の真値 y = e からの誤差と必要としたステップ数の値を対数プロットすると下図のようになった. 横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数(= 関数評価(f()の計算)回数)である.
 

オイラー法と後退オイラー法では, 1,000ステップ(h = 1.0e-3)でおよそ3桁半の精度, 1,000,000ステップ(h = 1.0e-6)でおよそ6桁半の精度となり, 1次の方法であることが確かめられた. 一方, 中点則では, 1,000ステップ(h = 1.0e-3)でおよそ6桁半の精度, 1,000,000ステップ(h = 1.0e-6)でおよそ12桁半の精度があり, 2次の方法となっていることがわかる. 同じ計算量であれば中点則はオイラー法の倍の桁数を正しく求めることができる.

1段法

テイラー法

関数 y(t) の tn におけるテイラー展開を再度次に示す.

  
  y(tn+1) = y(tn) + hf(tn, y(tn)) + (1/2)h2f'(tn, y(tn))+ (1/6)h3f''(tn, y(tn)) + ... + (1/p)hpf(p-1)(tn, y(tn)) + ...

 
オイラー法はテイラー展開を1次の項までで打ち切って近似したものに相当した. そこで, より高次の項まで使った近似式を使えばオイラー法より精度のよい, すなわち, より次数の高い解法を得ることができると考えられる. これをテイラー法という.

上式でf'(t, y(t))とf”(t, y(t))は具体的には次のように計算する.

  
  f'(t, y(t)) = ft(t, y(t)) + fy(t, y(t))・f(t, y(t))
  f''(t, y(t)) = ftt(t, y(t)) + 2fty(t, y(t))・f(t, y(t)) + fyy(t, y(t))・f2(t, y(t)) + ft(t, y(t))・fy(t, y(t))
    + fy2(t, y(t))・f(t, y(t))

 
2次の項を計算するのでさえ結構複雑なことがわかる. そのため, f()の高階の微分が簡単に計算できる関数でない限りこの方法を実際に適用するのは難しい.

数値実験 (2)

テイラー法による計算例
数値解法の次数の実験と同じ最も簡単な微分方程式の初期値問題の例を使って1次(オイラー法), 2次, 3次のテイラー法による計算を実際に行って比較してみる.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. それぞれ1次~3次の傾きを示している. ただし, 1ステップの計算には上のf'(t, y(t))やf”(t, y(t))の計算が含まれるので, 1ステップあたりの計算量は, 2次のテイラー法ではf()の計算が3回相当, 3次のテイラー法ではf()の計算が6回相当ということになる. そのため計算量を指標とすると一般的にはオイラー法との差はやや縮まる. ただし, 上の例の方程式では微分係数が0か1と簡単なため1ステップあたりの計算量の差は少ない.

2次のルンゲ・クッタ法

ルンゲ・クッタ法 (Runge-Kutta method) は, f() の微分を計算することなく高次のテイラー展開に相当する精度を得る1段法である. いくつかのtとyについてf(t, y)を計算しそれらの重み付き平均を作ってテイラー展開のp次の項まで一致させることによりp次の公式を作る.

2次のルンゲ・クッタ法の一般形は次式のとおりである. a, b1, b2, c はパラメータである.

  
  k1 = f(tn, yn)
  k2 = f(tn + ch, yn + hak1)
  yn+1 = yn + h(b1k1 + b2k2)

 
k2を (tn, yn)のまわりでテイラー展開すると次のようになる.

  
  k2 = fn + chft,n + ak1fy,n + ...

 
これを上の式に代入して, 次式が得られる.

  
  yn+1 = yn + (b1 + b2)hfn + b2h2(cft,n + afy,nfn + ... )

 
この式とテイラー展開の式を比較して, hとh2の係数が等しいとすると次式が得られる.

  
  b1 + b2 = 1
  b2c = 1/2
  b2a = 1/2

 
未知数が4個で式が3本なので自由度が1ある. よく使われるものとしては, c = 1 としたものをホイン法(Heun method)とよぶ. また, c = 1/2 としたものを修正オイラー法とよぶことがある (2次のルンゲ・クッタ法とはあまりよばれないが, これらの総称ととらえてよい).

ホイン法の計算式は次のとおり.

  
  k1 = f(tn, yn)
  k2 = f(tn + h, yn + hk1)
  yn+1 = yn + h(k1 + k2)/2

 
修正オイラー法の計算式は次のとおり.

  
  k1 = f(tn, yn)
  k2 = f(tn + h/2, yn + hk1/2)
  yn+1 = yn + hk2

 

数値実験 (3)

ホイン法と修正オイラー法の計算例
再び最も簡単な微分方程式の初期値問題の例でホイン法と修正オイラー法の精度を確認してみる.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. ホイン法と修正オイラー法の計算結果はほぼ同じである. 2次の傾きを確認でき, 最初の方で示した中点則と似たグラフになった. ただし, 計算式からわかるようにホイン法と修正オイラー法では1ステップあたり2回の関数評価(f()の計算)を行っているのに対して, 中点則では1ステップあたり1回の関数評価しか行わなくてよい.

ルンゲ・クッタ法の一般形

上と同様に3次, 4次・・・と高次のルンゲ・クッタ法の公式を作ることができ, 多くの公式が提案されている.

ルンゲ・クッタ法は一般的に次のように表される.

  
  yn+1 = yn + hΣbiki (Σ は i = 1, 2, ... , s)
  ki = f(tn + cih, yn + hΣaijkj) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)

 
s は段数(stage)とよばれ, 1ステップ進むのに必要な f(t, y) の計算回数を表す(1段法や多段法の”段”とは異なる). aij, bi および ci はパラメータで, bi および ci は次式を満たす.

  
  ci = Σaij (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
  Σbi = 1 (Σ は i = 1, 2, ... , s)

 
ルンゲ・クッタ法のパラメータには次のような表記法がある (ブッチャー配列という).

  
  c1 | a11  a12 ... a1s
  c2 | a21  a22 ... a2s
  c3 | a31  a32 ... a3s
  .  | .  .
  cs | as1  as2 ... ass
  ----------------------------
     | b1   b2  ... bs  
    

 
aijの値によって3つの場合に分けられる.

(1) aij = 0 (i ≤ j) の場合 (右上三角部分と対角部分が0)
陽解法. (単に「ルンゲ・クッタ法」あるいは「陽的(explicit)ルンゲ・クッタ法」という)
c1 = 0 である.
ブッチャー配列ではaijの右上三角部分と対角部分を書かない(空白にする).

(2) aij = 0 (i < j) の場合 (右上三角部分が0)
半陰解法. (半陰的(semi-implicit)ルンゲ・クッタ法という)

(3) aij ≠ 0 (i < j) の場合 (右上三角部分に0でないものがある)
陰解法. (陰的(implicit)ルンゲ・クッタ法という)

ここでは(1)についてのみ説明する. (2), (3)については別途スティフな方程式の項で触れる.

上で示したホイン法と修正オイラー法をブッチャー配列で表すと次のようになる.

ホイン法

  
  0   |
  1   | 1
  ----------------
      | 1/2  1/2  
    

 

修正オイラー法

  
  0   |
  1/2 | 1/2
  -------------
      | 0   1  
    

 
(陽的)ルンゲ・クッタ法ではできるだけ高次の公式になるようにパラメータを決めるが, s段ならばs次になるわけではなく, s段の公式が達成しうる最大の次数pは次のようになることがわかっている.

段数 s 1 2 3 4 5 6 7 8 9 10 11
最大次数 p 1 2 3 4 4 5 6 6 7 7 8

4次のルンゲ・クッタ法

4段4次の公式を作ろうとすると, パラメータが12個で自由度が4なので様々な公式を作ることが可能である. その中で最良のものを選ぶための多くの研究が行われた. 中でもクッタが選んだ以下の2組の係数が有名である.

古典的ルンゲ・クッタ法

単に「ルンゲ・クッタ法」というとこの公式を指す.

  
  0   |
  1/2 | 1/2
  1/2 | 0   1/2
  1   | 0   0    1
  -----------------------
      | 1/6 1/3 1/3 1/6 
    

 
計算式は次のとおり.

  
  k1 = f(tn, yn)
  k2 = f(tn + h/2, yn + hk1/2)
  k3 = f(tn + h/2, yn + hk2/2)
  k4 = f(tn + h, yn + hk3)
  yn+1 = yn + h(k1 + 2k2 + 2k3 + k4)/6

 

クッタの3/8公式

  
  0   |
  1/3 |  1/3
  2/3 | -1/3  1
  1   |  1   -1   1
  -----------------------
      |  1/8 3/8 3/8 1/8 
    

 
計算式は次のとおり.

  
  k1 = f(tn, yn)
  k2 = f(tn + h/2, yn + hk1/3)
  k3 = f(tn + h/2, yn + h(-k1/3 + k2))
  k4 = f(, tn + h, yn + h(k1 - k2 + k3))
  yn+1 = yn + h(k1 + 3k2 + 3k3 + k4)/8

 
ルンゲ・クッタ法一族において最もよく使われているのは4次の古典的ルンゲ・クッタ法である. いくつかある4段4次の公式間での精度の差はほとんどないが, この公式は係数aijに0が多く, すべての係数が正であるなどの理由により他のものより使いやすい. また, 5次の公式は6段以上でないと作れないため, これ以上高次のものを求めると複雑になるためこれが普及したのかもしれない.

数値実験 (4)

4次のルンゲ・クッタ法の計算例
最も簡単な微分方程式の初期値問題の例で4次の古典的ルンゲ・クッタ法および3/8公式の精度を確認してみる.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 古典的ルンゲ・クッタ法および3/8公式の計算結果はほぼ同じである. 両方とも4次の傾きをしていることが確認できる. ただし, 1ステップあたり4回の関数評価(f()の計算)が必要である.

多段法

1段法は yn+1 の計算に現在の値 yn のみを使う解法であったが, 過去の値(それまでに計算済の値) yn-1, yn-2, … も使った方がよい結果が得られそうに思われる. このアイデアを実現したのが多段法である.

1階の常微分方程式の積分形を考える.

  
  y(tn+1) = y(tn) + ∫ f(t, y(t)) dt [tn, tn+1]

 
ここで, 被積分関数f()を補間多項式で近似して右辺の積分を求める方法を一般にアダムス法(Adams method)という.

アダムス・バシュフォース法

y(tn+1) と y(tn)は近似値 yn+1 と yn を使い, 右辺の積分は (ti, fi) (i = n – k +1, …, n) を通る補間多項式 p(t) で f(t, y(t)) を置き換えると, 次のように k段の陽解法の公式が得られ, アダムス・バシュフォース法 (Adams-Bashforth method) とよばれる. ただし, fi = f(ti, yi) を表し, yi (i = n – k + 1, …, n) が計算済なので fi はすぐに求めることができる. この公式は, k次の公式になる.

  
  yn+1 = yn + ∫ p(t) dt [tn, tn+1]
       = yn + hΣγjjfn  (Σは j = 0, ... , k-1)

 
ただし,

0fn = fn,  ▽j+1fn = ▽jfn - ▽jfn-1  (後退差分)
    γj = (-1)j ∫ C(-s, j) ds [0, 1]

 
である. ただし, C(n, k)は二項係数である. なお, γj は直接計算するよりも次の漸化式を用いて求めるとよい.

  
  γ0 = 1,  γj + (1/2)γj-1 + (1/3)γj-2 + ... + (1/(j+1))γ0 = 1

 
j = 0 ~ 7 のときの値は次のとおりである.

j 0 1 2 3 4 5 6 7
γj 1 1/2 5/12 3/8 251/720 95/288 19087/60480 5257/17280

この方法を図に示すと次のようになる. 区間 [tn, tn+1] の積分を計算するために, 区間 [tn-k+1, tn] の補間データを使用した補間関数 p(t) を使う. 図で, 黒丸は補間点, 色つきの部分は積分範囲を表す.
 
ode_fig-1

k = 1 ~ 5 について, fn, fn-1, fn-2, … について整理すると次式のようになる.

  
  k = 1 : yn+1 = yn + hfn
  k = 2 : yn+1 = yn + h/2(3fn - fn-1)
  k = 3 : yn+1 = yn + h/12(23fn - 16fn-1 + 5fn-2)
  k = 4 : yn+1 = yn + h/24(55fn - 59fn-1 + 37fn-2 - 9fn-3)
  k = 5 : yn+1 = yn + h/720(1901fn - 2774fn-1 + 2616fn-2 - 1274fn-3 + 251fn-4)

 
k = 1 の場合は1段法とみることができオイラー法に一致する.
なお, アダムス・バシュフォース法はk次である.

アダムス・モルトン法

アダムス・バシュフォース法では補間関数を使って外挿により積分を求めたが, 精度を上げるために tn+1を加えて区間 [tn-k+1, tn+1] の補間データを使う, すなわち, (ti, fi) (i = n – k + 1, …, n, n + 1) を通る補間多項式 p*(t) を使い, 外挿することを避けたのがアダムス・モルトン法 (Adams-Moulton method) である. こちらの方が精度がよさそうであるが, その代わり右辺に yn+1 が現れ陰解法になる. 式は次のとおりである. この公式は, k+1次の公式になる.

  
  yn+1 = yn + ∫ p*(x) dx [tn, tn+1]
          = yn + hΣγ*jjfn+1  (Σは j = 0, ... , k)

 
ただし,

0fn = fn,  ▽j+1fn = ▽jfn - ▽jfn-1  (後退差分)
  γ*j = (-1)j ∫ C(-s+1, j) ds [0, 1]

 
である. ただし, C(n, k)は二項係数である. なお, γ*j は次の漸化式を用いて計算するとよい.

  
  γ*0 = 1,  γ*j + (1/2)γ*j-1 + (1/3)γ*j-2 + ... + (1/(j+1))γ*0 = 0

 
j = 0 ~ 7 のときの値は次のとおりである.

j 0 1 2 3 4 5 6 7
γ*j 1 -1/2 -1/12 -1/24 -19/720 -3/160 -863/60480 -275/24192

この方法を図に示すと次のようになる. 下図で, 黒丸は補間点, 色つきの部分は積分範囲を表す. 積分範囲は補間区間に含まれており, 外挿されていない.

ode_fig-2

 
k = 0 ~ 5 について, fn+1, fn, fn-1, … について整理すると次式のようになる.

  
  k = 0 : yn+1 = yn + hfn+1
  k = 1 : yn+1 = yn + h/2(fn+1 + fn)
  k = 2 : yn+1 = yn + h/12(5fn+1 + 8fn - fn-1)
  k = 3 : yn+1 = yn + h/24(9fn+1 + 19fn - 5fn-1 + fn-2)
  k = 4 : yn+1 = yn + h/720(251fn+1 + 646fn - 264fn-1 + 106fn-2 - 19fn-3)
  k = 5 : yn+1 = yn + h/1440 (475fn+1 + 1427fn - 798fn-1 + 482fn-2 - 173fn-3 + 27fn-4)

 
k = 0 の場合は1段法とみることができ後退オイラー法に一致する. k = 1 の場合は台形則とよばれる公式に一致する.
なお, アダムス・モルトン法はk+1次である.

数値実験 (5)

アダムス・バシュフォース法とアダムス・モルトン法の計算例
最も簡単な微分方程式の初期値問題の例で k = 4 のアダムス・バシュフォース法(陽解法)とアダムス・モルトン法(陰解法)単体(注)で計算を行い精度を確認する. (注: 実務では後で説明するように組み合わせて使用し, 単体では用いない)

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 4次のルンゲ・クッタ法の計算結果も示した. 多段法では陰解法の方が精度がよいことが確認できる. 多段法のほうが同じ次数では少し精度が悪いように感じられるが, ステップ数で表示している縦軸を関数評価回数(計算量)で表示し直すとほとんど差がなくなる.

多段法では, はじめに出発値 y1, …, yk-1 の値を求めておく必要があるが, そのために1段解法などを補助的に使わなければならない. この精度が悪いと以降の計算値に影響するため, 使用する多段法の公式の精度に十分見合うだけの精度を持った方法を使用する必要がある. ここでは最初の4ステップ目まではルンゲ・クッタ法(4次)により計算を行った.

予測子・修正子法

多段法においては, 陽解法は計算が容易であるが満足のいく精度が得られない. 一方, 陰解法は精度はよいが一般的には非線形方程式を毎回解かなければならないため計算が難しい.

そこで広く使われている方法が予測子・修正子法(predictor-corrector method)である. まず予測子(predictor)とよばれる陽解法によりyn+1の第1近似を求める. 次にその値を使って修正子(corrector)とよばれる陰解法を1回~数回適用しこれを改良する方法である.

予測子と修正子の組み合わせにはバリエーションがあるが, 予測子にアダムス・バシュフォース法, 修正子にアダムス・モルトン法を使った予測子・修正子法をアダムス・バシュフォース・モルトン法(あるいは単にアダムス法)とよび, 広く使われている.

予測子・修正子法の計算手順は, (1) 予測子を適用して Yn+1を求める(P), (2) それを使って関数値 fn+1を計算する(E), (3) 修正子を適用して Yn+1を更新する(C), (4) 更新した Yn+1を使って fn+1を更新する(E) となる. これを PECE と表すことにすると, 最後の E を省略して PEC とする方法, 修正子の適用を 2回繰り返して PECECE とする方法など, バリエーションがある.

数値実験 (6)

予測子・修正子法(アダムス・バシュフォース・モルトン法)の計算例
最も簡単な微分方程式の初期値問題の例で k = 4 の予測子・修正子法(アダムス・バシュフォース・モルトン法)の精度を確認してみる.

  
  y' = y, y(0) = 1

 
ここではPECとPECEの2種類について試してみた. 出発値を求めるためには4次のルンゲ・クッタ法を使用した.
 

横軸は相対誤差(右に行くほど精度がよい), 縦軸は今回は関数評価回数(ステップ数ではない)の対数プロットとした. ステップあたりの関数評価回数がPECとPECEで異なるため, 総合的に評価するように縦軸のプロットを変えた.
そのPECとPECEであるが, この例は簡単な関数であるためかほとんど計算結果に差が出なかった. そのため, 2回目のEの計算のぶんだけPECEの関数評価回数が増えただけという結果になった. なお, 別途難しい関数で試してみるとPECEのほうがよくなることがあった.

使われなくなった解法

よく知られているが現在では使われなくなった解法についてここで触れておく.

ルンゲ・クッタ・ジル法

ルンゲ・クッタ・ジル法(Runge-Kutta-Gill method)は係数が次のように表される4次のルンゲ・クッタ法である.

  
  0   |
  1/2 |    1/2
  1/2 | (√2-1)/2 (2-√2)/2
  1   |     0     (-√2)/2  1+(√2)/2
  ---------------------------------------
      |    1/6   (2-√2)/6  (2+√2)/6  1/6 
    

 
Gillの提案は係数の違いだけを見てもよくわからないが, その特長は次の2点である.

(1) 普通の4次のルンゲ・クッタ法では計算のために4個の変数(連立方程式の場合には配列)が必要なのに対し, この方法では係数を工夫して3個の変数(配列)だけで計算できるようになっていてその手順も示されている. [メモリの節約]

(2) 丸め誤差を減らすための手順が含まれている(これは他のルンゲ・クッタ法にも適用可能). [精度の向上]

メモリが少なくおそらく単精度計算が主流であった初期の計算機におけるプログラムに適した方法だったと思われる. かつて広く使われていたが, メモリが豊富で倍精度計算があたりまえの現代では必要性が薄れ使われなくなった.

ここで, (2)について少し説明する.
常微分方程式の基本的な数値解法では各ステップにおいて関数値の更新差分(ここでは kと表す)を計算し, 次のようにして関数値を更新する.

  
  yn+1 = yn + k

 
この計算を計算機の浮動小数データ(ワード長は有限)で行うと下図のようになる.

ode_fig-5

ここで, k をワード長いっぱい正しく求めたとしても, この更新に使用されるのは色のついた一部だけで残りの部分は捨てられてしまう(丸め誤差という). そこで, この捨てられる部分を別の場所に累積しておき, 色のついた部分に頭が出るほどに貯まるたびにkに加えて関数値の更新に使用することにより丸め誤差を減らす方法が考案された. この方法がルンゲ・クッタ・ジル法に使われている.

数値実験 (7)

ルンゲ・クッタ・ジル法の計算例
最も簡単な微分方程式の初期値問題の例でルンゲ・クッタ・ジル法の精度を確認してみる.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸は関数評価回数を対数プロットとした. 比較のため普通のルンゲ・クッタ法もプロットしている.

相対誤差が1e-14付近までは差がないが, 普通のルンゲ・クッタ法では1e-15あたりで頭打ちになり, その先は刻み幅を小さくするとむしろ誤差が大きくなってしまう. これは, 64ビット倍精度では1回の演算の精度の限界は2e-16くらいなので, この場合は総合的には1e-15くらいが最高精度だと思われ, それ以上計算を繰り返すと丸め誤差が積み重なっていきかえって精度が落ちてしまうと考えられる.

一方で, ルンゲ・クッタ・ジル法では最後のほうの誤差は1.6e-16または0となった(図で線が切れているところは完全に0のためプロットされていない). つまり, 最終1ビットが違うかどうかのめいっぱいの精度で計算されたことを示している.

ミルン法

アダムス法と似ているが, 積分範囲が異なる次の積分形を考える.

  
  y(tn+1) = y(tn-1) + ∫ f(t, y(t)) dt [tn-1, tn+1]

 
アダムス法と同様に (ti, fi) (i = n – k + 1, …, n) を通る補間多項式p(t)で積分を近似すると, 次の陽解法の公式が得られる.

  
  yn+1 = yn-1 + hΣκjjfn  (Σは j = 0, ... , k-1)

 
ただし,

  
  κj = (-1)j ∫ C(-s, j) ds [-1, 1]

 
また, (ti, fi) (i = n – k +1, …, n+1) を通る補間多項式 p*(x) で積分を近似すると, 次の陰解法の公式が得られる.

  
  yn+1 = yn-1 + hΣκ*jjfn+1  (Σは j = 0, ... , k)

 
ただし,

  
  κ*j = (-1)j ∫ C(-s+1, j) ds [-1, 1]

 
予測子・修正子法の一種にミルン法(Milne method)とよばれる方法があり, 上の陰解法の式において k = 2 とした次の式がその修正子として使われている.

  
  yn+1 = yn-1 + h/3 (fn+1 + 4fn + fn-1)

 
ミルン法の予測子には, 積分範囲をさらに [tn-3, tn+1] に広げた積分形より導かれる次式が使われている.

  
  yn+1 = yn-3 + h/3 (8fn - 4fn-1 + 8fn-2)

 
ミルン法は4次の公式で, コンピュータが普及する以前の手計算の時代によく使われていた. しかし, コンピュータで使われるようになると安定性の面で問題があることがわかり, 現在では使われなくなった.

なお, 安定性の問題を解決するために, ミルン法の予測子と組み合わせて使う安定な修正子として次のハミング(Hamming) の修正子が提案されている.

  
  yn+1 = 9/8 yn - 1/8 yn-2 + h/8 (3fn+1 + 6fn - 3fn-1)

 

数値実験 (8)

ミルン法の計算例
最も簡単な微分方程式の初期値問題の例でミルン法の精度を確認してみる.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸は関数評価回数を対数プロットとした. 比較のためルンゲ・クッタ法とアダムス・バシュフォース・モルトン法もプロットしている. なお, 出発値を求めるためには4次のルンゲ・クッタ法を使用した.

ミルン法は公式が比較的簡単ではあるが図のように性能がよい解法である.

数値実験 (9)

ミルン法の不安定性
次の初期値問題をミルン法により解く.

  
  y' = -y + 1, y(0) = 2

 
この方程式の解は y(t) = exp(-t) + 1 である. ステップ幅を 0.001 および 0.0001 として計算した結果, 次のようになった(グラフ作成の都合上データ点を間引きして表示).

ode_test-9

ステップ幅が十分小さいにもかかわらず計算が進んでステップ数が多くなると振動を始め発散した. これは, 解法の差分方程式の一般解が真の解の近似(有効成分)の他に必要のない成分を含み, 条件によってそれが有効成分よりも大きくなったときに起きる現象である.

ミルン法ではこの不安定性の対策としてハミングの修正子があるので, それを使って同じ計算を行うと次のようになった(グラフ作成の都合上データ点を間引きして表示).

ode_test-9_2

 

ハミングの修正子を使うことにより不安定性が収まることが確認された.

数値実験(1)で中点則の性能がよいことに気づいたかもしれないが, 中点則にもこの種の不安定性があることがわかっている. なお, 1段法ではこの種の不安定性は生じない.


参考文献

[1] 名取亮「数値解析とその応用」(1990) コロナ社
[2] 三井斌友「常微分方程式の数値解法」(2003) 岩波書店
[3] E. Hairer, 他「常微分方程式の数値解法 I」(2007) シュプリンガー・ジャパン