10. 常微分方程式 (2)

誤差の推定

求めた数値解の誤差を推定する方法について説明する。計算結果が得られればそれで終わりというのではなく、その結果の誤差が推定できていれば、特異点などのために異常な解が得られていないかなどのチェックができ計算結果の利用価値は上がる。また、現在主流となっているきざみ幅の自動調節機能付きルーチンの実装には誤差推定が必須である。

1段法の誤差

初期値問題の数値解の誤差の原因には次の2つがある。

  • 離散化誤差 (打ち切り誤差ともいう)
    演算がすべて無限精度で行われたとしたときに、解法に起因して生じる誤差

    • 局所的離散化誤差
      それまでに得られている値が正確であるとして、次の1ステップで生じる誤差
    • 大域的離散化誤差
      初期値から計算される解と真の解との差
  • 丸め誤差
    数値を浮動小数で表し有限桁数で演算することにより生ずる誤差

いま、次の初期値問題を考える。

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

tnにおける計算値をynとし、un(tn) = yn とする。すなわち、計算値ynから決定される解曲線をun(t)とする。yn+1は yn, tn, f より計算される推定値である。例えば、オイラー法を使うならば、きざみ幅を h として、yn+1 = yn + hf(yn, tn) となる。

そうすると、局所的離散化誤差 dn は次のように表すことができる。

    dn = yn+1 - un(tn+1)

局所的離散化誤差は p次の解法において O(hp+1) である。すなわち、次式を満たす数 C が存在する。

    |dn| ≤ Chp+1

大域的離散化誤差 en は次のように表される。

    en = yn - y(tn)

以上の関係を下に図示する。これは、解曲線が少しずれたときに時間が進むにつれて離れていく関数の例である。すなわち、en > Σdi (i = 0 ~ n-1) となり、ステップが進むにつれて誤差が拡大していく。

 

ode_fig-3

下に他の例を示す。これは、解曲線が少しずれたときに時間が進むにつれて接近していく関数の例である。すなわち、en < Σdi (i = 0 ~ n-1) となり、ステップが進むにつれて誤差は縮小していく。

 

ode_fig-4

2つの例をあげたが、関数がどちらの傾向を示すかは ∂f/∂y の符号により決まり、正であれば前者、負であれば後者になる。もし、f が y に依存しない場合、すなわち、∂f/∂y = 0 ならば en = Σdi となる。また、∂f/∂y の符号が区間内で変化するような複雑な関数もある。∂f/∂y = 0 に近ければ、大域的離散化誤差はおおむね O(hp) といえる。

ここまでは丸め誤差を無視して考えてきたが、有限桁で計算する限り必ず丸め誤差 εn が発生する。例えば、オイラー法ならば、yn+1 = yn + hf(yn, tn) + εn となる。

ここで、丸め誤差 εn のノルムはある数 ε 以下であるとする。

    ||εn|| ≤ ε

いま、最終点を tN とすると、ステップ数 N = (tN – t0) / h である。最終的な大域的離散化誤差を eN ≈ Σdi (i = 0 ~ N-1) と近似し、丸め誤差を Nε と近似すると、最終的な全誤差は次のように近似することができる。

    全誤差 ≈ NChp+1 + Nε = (tN - t0)(Chp + ε/h)

実験 4

次の初期値問題をオイラー法により解く。

    y' = y, y(0) = 1

オイラー法の場合 p = 1 であるから、全誤差は次のように近似できる。

    全誤差 ≈ (tN - t0)(Ch + ε/h)

h が大きいときは第1項(離散化誤差)が優勢で、h が小さいときは第2項(丸め誤差)が優勢になると考えられる。すなわち、きざみ幅を小さくしていくと誤差は小さくなっていくが、ある値((ε/C)-2)から下では丸め誤差のためにかえって誤差が大きくなると予想される。

きざみ幅を変化させて、t = 1 における誤差をプロットした結果を下に示す。横軸はきざみ幅、縦軸は誤差である。なお、丸め誤差の影響を見やすくするために計算は単精度で行った。

ode_test-4

 

図のように、この場合 10-4 ~ 10-5 付近を境に丸め誤差が支配的になることがわかる。

実験 5

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

    yn+1 = yn + k

例えばオイラー法では、k = hf(yn, tn) である。この計算を計算機の浮動小数データ(ワード長は有限)で行うと下図のようになる。

ode_fig-5

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

この方法を実験4のオイラー法に適用して計算してみると次の図のようになった。

ode_test-5

効果は抜群で、単精度の限界(おおよそ10-7)近くの精度で解を求めることができるようになった。

ただし、現代の計算環境においては倍精度で計算を行うのが普通(Excelもそうである)なので、この問題はあまり気にならない。倍精度計算において丸め誤差の影響が深刻になるほどきざみ幅を小さくすることはほとんどないからだと思われる。

リチャードソンの補外

p次の解法の与えられた点 tn での大域的離散化誤差は O(hp) となるので、次のように表すことができる。

    yn - y(tn) = chp + O(hq), q > p

きざみ幅 h で計算した yn を yn(h) と書くことにして、きざみ幅 2h でも計算すると次のようになる。

    yn(h) = y(tn) + chp + O(hq)
    yn(2h) = y(tn) + c(2h)p + O(hq)

これらの差をとると、

    yn(2h) - yn(h) = (2p - 1)chp + O(hq)

となり、yn(h) の誤差は O(hq) を無視すると次のように近似できる。

    yn(h)の誤差 ≈ chp = (yn(2h) - yn(h))/(2p - 1)

これは、ロンバーグ積分で使われるリチャードソンの補外 (Richardson extrapolation) と同じである。

実験 6

実験 1 の計算結果を使い、この式により誤差を推定してみると次のようになった。

ode_tab-1

 

(Err Est) の欄が推定誤差である。実際の誤差 (Err) と比較的よく一致しているのがわかる。

推定誤差が求められたら、きざみ幅 h を使った計算結果に対してそれを補正してやればより良い結果を得ることができるはずである。オイラー法、ホイン法、および、4次のRK法について実際にやってみると次のようになった。

ode_test-6

 

横軸はきざみ幅、縦軸は誤差である。それぞれの1本下の線 (凡例で+Rexと記した線) が補正を加えた計算結果である。かなり効果があることがわかる。

埋め込み型RK法

誤差の推定を行うのにリチャードソンの補外を使うためにはきざみ幅を変えて2回計算を行わなければならず、計算量が大きくなるのが欠点である。しかし、RK法自身に誤差推定機能を付加した公式を作ることができれば、1回の計算で済ますことができ有用と思われる。

RK法は次のように表された。ここで、この公式は p次であるとする。

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

同じ aij と ci を使い (すなわち、同じ ki を使い)、bi を b*i に変えた時の近似値 y*n+1 がp*次 (例えば、p* = p – 1、p* = p + 1 など) になっているように定数を決めることができれば、ki を計算し直すことなく次式により p*次の場合の近似値を求めることができる。

    y*n+1 = yn + hΣb*iki (Σ は i = 1, 2, ... , s)

このとき、次数が低い方の局所離散化誤差を次のように求めることができる。

    yn+1 - y*n+1 = hΣ(bi - b*i)ki (Σ は i = 1, 2, ... , s)

このような公式は実際に作ることができて、埋め込み型RK法 (embedded RK method) とよばれる。主な公式を以下に示す。

ルンゲ・クッタ・マーソン法 (Runge-Kutta-Merson method)

埋め込み型RK法で最初に提案されたのはマーソンの公式 (ルンゲ・クッタ・マーソン法) とされる。

    0   |
    1/3 | 1/3
    1/3 | 1/6 1/6
    1/2 | 1/8  0  3/8
    1   | 1/2  0 -3/2   2
    ---------------------------
    y   | 1/6  0   0   2/3 1/6 
    y*  | 1/2  0 -3/2   2   0 
    

計算式は次のとおり。

    yn+1 = yn + h(k1 + 4k4 + k5)/6  (4次)
    y*n+1 = yn + h(k1 - 3k3 + 4k4)/2  (誤差評価用 4次)
    k1 = f(yn, tn)
    k2 = f(yn + hk1/3, tn + h/3)
    k3 = f(yn + h(k1/6 + k2/6), tn + h/3)
    k4 = f(yn + h(k1/8 + 3k3/8), tn + h/2)
    k5 = f(yn + h(k1/2 - 3k3/2 + 2k4), tn + h)

この公式による推定誤差は、f(y, t) が y と t に関して線形ならば 1/5 (yn+1 – y*n+1) となる。

1ステップあたり 5回の関数評価が必要である。すなわち、誤差評価を行うために、通常の 4次の RK法に比べ 1ステップあたり 1回余計に関数評価を行うという犠牲を払っている。

5(4)次 ルンゲ・クッタ・フェールベルグ法 (Runge-Kutta-Fehlberg method)

埋め込み型RK法の中でも有名な公式である。

    0     |
    1/4   | 1/4
    3/8   | 3/32       9/32
    12/13 | 1932/2197 -7200/2197 7296/2197
    1     | 439/216    -8        3680/513   -845/4104
    1/2   | -8/27       2       -3544/2565  1859/4104 -11/40
    ----------------------------------------------------------------
    y     | 16/135      0       6656/12825 28561/56430 -9/50  2/55 
    y*    | 25/216      0       1408/2565   2197/4104  -1/5    0 
    

この公式では、近似値 y の計算を 5次で行い、誤差評価用の値 y* の計算を 4次で行う。これを、5(4)次の公式と書く。1ステップあたり 6回の関数評価が必要である。

5(4)次 ドルマン・プリンス法 (Dormand-Prince method)

    0    |
    1/5  |     1/5
    3/10 |     3/40        9/40
    4/5  |    44/45      -56/15      32/9
    8/9  | 19372/6561 -25360/2187 64448/6561  -212/729
    1    |  9017/3168   -355/33   46732/5247    49/176  -5103/18656
    1    |    35/384        0       500/1113   125/192  -2187/6784    11/84
    --------------------------------------------------------------------------------
    y    |    35/384        0       500/1113   125/192  -2187/6784    11/84     0
    y*   |  5179/57600      0      7571/16695  393/640 -92097/339200 187/2100  1/40
    

5(4)次の公式で、1ステップあたり 7回の関数評価が必要である。この公式の特徴は、7番目(k7)の係数 a7j が近似値 y を計算するための係数 bi に等しいことである。これは FSAL (First Same As Last) とよばれるアイデアで、k7 が次の k1 = f(yn+1, tn+1) に等しくなるよう巧妙に作られている。

6(5)次 ルンゲ・クッタ・ヴァーナー法 (Runge-Kutta-Verner method)

    0    |
    1/6  |     1/6
    4/15 |     4/75     16/75
    2/3  |     5/6      -8/3       5/2
    5/6  |  -165/64     55/6    -425/64      85/96
    1    |    12/5       -8     4015/612    -11/36     88/255
    1/15 | -8263/15000 124/75   -643/680    -81/250  2484/10625 0
    1    |  3501/1720 -300/43 297275/52632 -319/2322 24068/84065  0  3850/26703
    -----------------------------------------------------------------------------------
    y    |     3/40       0     875/2244     23/72     264/1955   0  125/11592  43/616
    y*   |    13/160      0    2375/5984      5/16      12/85    3/44   0      0
    

6(5)次の公式で、1ステップあたり 8回の関数評価が必要である。長い間広く使われているサブルーチンDVERKに採用されている公式である。

多段法の誤差

多段法における局所誤差は次のように求められる。

解法 次数 局所誤差
アダムス-バシュフォース法 k γk hk+1 y(k+1)
アダムス-ムルトン法 k+1 γ*k+1 hk+2 y(k+2)
後退微分公式 (BDF法) k -1/(k+1) hk+1 y(k+1)

例えば、4次のAdams法では、予測子(右肩にP)と修正子(右肩にC)のそれぞれの局所誤差は具体的には次のようになる。

    y(tn+1) - yPn+1 = γ4h5y(5) + O(h6) = 251/720 h5y(5) + O(h6)
    y(tn+1) - yCn+1 = γ*4h5y(5) + O(h6) = -19/720 h5y(5) + O(h6)

このように多段法では予測子と修正子で局所誤差が異なるので、それを使い誤差の推定をすることができる。

両辺の引き算をして高次の項を無視すると次の近似式が得られる。

    yCn+1 - yPn+1 ≈ 3/8 h5y(5)

これを上の修正子の式に代入して高次の項を無視すると次のような誤差評価式が得られる。

    y(tn+1) - yCn+1 ≈ -19/27(yCn+1 - yPn+1)

きざみ幅の自動調節

一定のきざみ幅で計算するのではなく、目標誤差を満たすように最適なきざみ幅に自動調節しながら計算するプログラムが現在の主流である。これは、目標誤差に対して計算量を最小化することと、計算の安定性を損なわないようにきざみ幅を維持することに効果がある。

上に説明した誤差推定機能を持つ公式を使えばそのようなプログラムを実装することができる。

きざみ幅の自動調節機能の実装

おおざっぱにいえば、きざみ幅の初期値と局所離散化誤差の許容限界を与え計算を開始し、ステップごとに推定誤差をチェックし、許容内であれば次ステップに進み計算を続ける(ステップ成功)。推定誤差が大きすぎる場合、きざみ幅を小さくしてやり直す(ステップ失敗)。また、推定誤差が小さすぎるとステップ数が不要に増えてしまうので、次ステップからきざみ幅を大きくするなどの調整を行う。このような処理を繰り返し、目標精度を満たしつつ最小のステップ数で計算を完了することを目指す。ただし、実際に実用に供するプログラムを実装するのは簡単ではなく多くの細かい工夫が必要である。

多段法の場合には、きざみ幅を変更すると以前の値の間隔が新しいきざみ幅と合わなくなり、公式がそのままでは使えなくなる。そこで、不等間隔のステップ点に対する公式に変更した可変ステップ幅アダムス法や可変ステップ幅BDF法などを使う必要がある。

実験 7

ルンゲ・クッタ・マーソン法(RKM)の公式を使い、簡単なきざみ幅自動調節機能付きプログラムを作成して実験してみる。

計算アルゴリズムは次のようにした。

  • きざみ幅の初期値と局所離散化誤差の許容限界を与え計算を開始する
  • RKM公式で1ステップ計算後、推定誤差が許容限界より大きければきざみ幅を1/2にしてやり直す
  • 推定誤差が許容限界より小さければ1ステップ進める。ただし、推定誤差が許容限界の1/32より小さいときには次ステップからきざみ幅を2倍にする

実験1と同じ問題を使って、作成したRKMルーチンの動作を確かめてみた。比較のために、4次のRK法およびXLPackの各ルーチン (Dderkf: フェールベルグ、Dopri5: ドルマン・プリンス、Dverk: ヴァーナー) による計算結果も示す。

RK法以外ではきざみ幅が変わるので、横軸は関数評価回数 (f(y, t) の呼び出し回数) とした。RK法ではきざみ幅を変えていき、そのときの誤差と関数評価回数をプロットした。その他のきざみ幅自動調節機能付きの各ルーチンでは、許容限界を 10-4, 10-5, … と変えていき、得られた結果の誤差とその時に実際に必要とした関数評価回数をプロットした。RKMではきざみ幅の初期値を 0.2 ((最終 t – 初期 t)の1/5) とした。XLPackの各ルーチンではきざみ幅の初期値は自動設定されるのでユーザーが与える必要はない。

ode_test-7

 

きざみ幅自動調節の各ルーチンはおおむね公式の次数なりの性能を示している。すなわち、自動調節がうまく機能しているようである。

作成したRKMルーチンも問題なく機能している。RKMの詳しい動作を下表に示す。

ode_tab-2

 

Tol は与えられた許容限界である。Error は結果の誤差、# of F call は関数評価回数、h (final) はきざみ幅の最終値、# of h dec. はきざみ幅を1/2にした回数、# of h inc. はきざみ幅を2倍にした回数である。

Tol=1e-4 ではきざみ幅を1回だけ2倍にしており、初期値 0.2では小さすぎたことがわかる。Tol=1e-7以下では1~6回きざみ幅を1/2にしており、初期値 0.2では大きすぎたことがわかる。きざみ幅を1/2にした回数はステップの失敗回数に等しい、すなわち無駄な関数評価を行ったステップ数を示すが、簡単な制御方法のわりに悪くない結果といえる。

これは簡単な関数の例であるが、途中で部分的にきざみ幅を小さくしなければ精度が得られないような難しい関数の場合に、きざみ幅の自動調節を行うルーチンの威力が発揮されるはずである。

安定性

初期値問題を解いている間に数値解が異常な値になったり発散したりすることがある。その原因は、微分方程式自体の問題である場合と、数値解法の問題である場合がある。

微分方程式自体の安定性

離散化誤差の説明図に示されるように、微分方程式の解曲線が時間が進むにつれて離れていく場合、初期値のわずかな差が大きく拡大する現象が生じることがある。これは不安定な微分方程式の例である。

逆に、微分方程式の解曲線が時間が進むにつれて近づいていく場合、むしろ誤差が縮小していくので安定な微分方程式の例である。

繰り返すと、関数がどちらの傾向を示すかは微分方程式の ∂f/∂y の符号により決まる。

数値解法の安定性

安定な微分方程式であっても、使用する数値解法によって不安定現象を示すことがある。これは数値解法の安定性の問題である。

きざみ幅の制約

いま、微分方程式の導関数 ∂f/∂y が一定値 λ で、λ < 0 とする(これは安定な微分方程式の初期値問題である)。また、各ステップで加えられる誤差はステップによらず一定であるとする。その場合、数値解法の安定性は次の線形テスト問題を調べることでわかる。

    dy/dt = λy (t > 0, λ(定数) < 0), y(0) = y0

この問題の解は、

    y(t) = y0exp(λt)

となるから、t → ∞ のとき y(t) → 0 となる。

(1) オイラー法

上のテスト問題をオイラー法で解くと、

    yn = y0(1 + hλ)n

となるから、n → ∞ のとき yn → 0 とならなければならない。そのためには、

    |1 + hλ| < 1

でなければならない。これを安定条件という。

λを複素数として z = hλ とおき、安定条件 |R(z)| < 1 (ただし、R(z) = 1 + z) を複素平面上に図示すると次のようになる。 ode_euler色付き部分 { z ∈ C; |R(z)| < 1 } を絶対安定領域という。λが実数であれば -2 < z = hλ < 0 である。区間 [-2, 0] を絶対安定区間という。このとき、きざみ幅は 0 < h < -2/λ としなければならない。

(2) 後退オイラー法

後退オイラー法ではR(z)は次のようになる。

    R(z) = 1/(1 - z)

絶対安定領域を図示すると次のようになる。

ode_beuler

このように、左半平面が絶対安定領域に含まれるときにA安定であるといい、きざみ幅 h に対する制限がなくなる。

(3) ルンゲ・クッタ法

s段p次のRK法では R(z) は次のようになる。

    R(z) = 1 + z + z2/2! + ... + zp/p! + γp+1zp+1 + ... + γszs

ここで、γp+1, … はRK法のパラメータにより決まる。

p = s = 2, 3, 4 について絶対安定領域を図示すると次のようになる。

ode_rk

次数が上がると安定領域も広くなることがわかる。絶対安定区間はそれぞれ、[-2, 0]、[-2.513, 0]、[-2.785, 0] である。ルンゲ・クッタ法の陽解法は A安定ではない。

(4) 多段法 (予測子・修正子法)

多段法も陽解法は A安定ではない。

多段法の陰解法では A安定な解法は2次以下である (A安定な多段法で最も高次なものは台形則である)。

(5) 台形則 (陰解法、アダムス・モルトン法の k=1)

台形則ではR(z)は次のようになる。

    R(z) = (1 + z/2)/(1 - z/2)

この絶対安定領域は左半平面全部となり、台形則はA安定である。

きざみ幅の制約 (連立方程式の場合)

連立方程式の場合、線形テスト問題で y をベクトルに拡張して扱う。

    dy/dt = Jy, y(0) = y0

J は y の元数と同じ次元数の正方行列となる。安定な微分方程式の初期値問題であるための条件は、J のすべての固有値 λi の実部が負であることである。

    Re(λi) < 0 (すべての i について)

また、z = hλi として、複素平面上で z が絶対安定領域に含まれていればその h について安定である。

実験 8

次の初期値問題をオイラー法により解く。

    y' = -10y, y(0) = 1

この方程式の解は、y(t) = exp(-10t) である。きざみ幅を 0.09, 0.19, 0.2, 0.21 と変えて計算してみた結果、次のようになった。

ode_test-8

 

h = 0.2 を境に、これよりきざみ幅を大きくすると発散した。0.2 は上の(1)で説明した絶対安定区間の上限にあたり、これを超えると不安定になることが確かめられた。なお、0.1 < h < 0.2 では振動しながら収束するが、0.1 < h では振動しない。

t = 0 の近くの値が急激に変化する部分を避けて、値がかなり落ち着いた t = 0.6 から計算を始めてみる。すなわち、初期条件を y(0.6) = exp(6) と変えて計算してみると、結果は次のようになった。

 

ode_test-8_2

同じような傾向を示し、発散は止まらないことがわかる。

A安定でない解法を使う限り、きざみ幅の値には上限がある。

不安定現象

ミルン法や中点法で生じる別の種類の不安定性について実験してみる。

実験 9

次の初期値問題をミルン法により解く。

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

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

ode_test-9

 

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

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

ode_test-9_2

 

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

この種の不安定性に対しては、適用する問題に対して発散していないことを確かめながら使うか、あきらめて解法を変える必要がある。なお、1段法ではこの種の不安定性は生じない。

スティフな(硬い)方程式

スティフな方程式の例

スティフな方程式は化学反応や電気回路などの問題で現れる。非常に急な変化をする解を持ち、通常の解法では解きにくいものをいう。

次の初期値問題はスティフな方程式の例である。

    dy/dt = ( -2    1    ) y + ( -cos(t)             )
            ( 1998 -1999 )     ( 1999cos(t) - sin(t) )

    y(0) = ( 1 )
           ( 2 )

上の方程式の y の係数行列の固有値は実数で、-1 および -2000 である。従って、上の安定性の議論から、オイラー法ならばきざみ幅 h < 0.001、4次のRK法ならば h < 0.00139 としなければ解くことができないはずである。実際に試してみると、h をこれらの制約条件よりも大きくすると発散して解くことができないことが確かめられた。

スティフな方程式は、このように安定性の問題により通常の解法ではきざみ幅に制約があるために、計算に非常に時間がかかったり、条件によっては振動や発散をしたりするものである。上の例のように、一見急な変化をするようには見えないものであっても、一般解に指数項を含んでいたりするとスティフな方程式となる。

実験 10

スティフな方程式であっても A安定な解法ならば解けるはずである。そこで、上の例題を後退オイラー法で解いてみる。

ode_test-10

上図は h = 0.2 として計算した y2 をプロットしたものである。このように後退オイラー法ならば問題なく解を求めることができた。

実験 11

XLPackに収録されている通常のきざみ幅自動調節機能付きのルーチンで同じ例題を解いてみる。要求精度は RTol = 0.01 とした。

Dderkf (RKF45) と Ddeabm (アダムス法) では、スティフ性を検出(エラーコード=-4)して途中で停止した。

Dverk (ヴァーナー法) は計算を継続し、非常に細かいきざみ幅を使用して時間をかけて解を求めた。このときの y2 をプロットした図を下に示す。

ode_test-11

 

おおよそ 0.002 のきざみ幅が使われているため、点でプロットしているが連続した太い線に見える。t = 10 まで計算するのに、試行のためのものを含め関数評価回数は約44000回であった。スティフ性のために、きざみ幅を大きくしようとすると不安定性により誤差が大きく評価され、自動調節機能により細かいきざみ幅が使われてしまうものと思われる。逆に、不用意に大きなきざみ幅を使い発散などしてしまうことがないようにできているとも言える。

スティフな方程式の解法

スティフな方程式を解くためには、きざみ幅を小さくして力ずくで解くか、きざみ幅に制約がない解法、すなわち、A安定な解法を使う必要がある。

A安定な解法としては後退オイラー法や陰的台形則があるが、実用的にはもっと高次の公式が必要であり種々の解法が研究されてきた。以下、代表的な解法として、後退微分公式(BDF)および陰的ルンゲ・クッタ法について説明する。

後退微分公式 (BDF)

アダムス・モルトン法では (ti, fi) (i = n – k +1, …, n+1) を補間データとして f を補間したが、y そのものを補間する方法がある。すなわち、(ti, yi) (i = n – k +1, …, n+1) を補間データとして補間多項式 q(t) を求め、それが少なくても格子点の1か所(i=n+1-r の点)で微分方程式を満たすようにする。

    q'(tn+1-r) = f(yn+1-r, tn+1-r)

ここで、r=1 とすると陽解法になり、k=1 と k=2 の場合はそれぞれオイラー法と中点法になる。しかし、k >= 3 の場合には不安定になり微分方程式の解法として使うことはできない。

r=0 とすると陰解法になるが k≤6 で安定な解法が得られる。これは後退微分公式 (Backward differentiation formula (BDF) または ギヤ(Gear)法) とよばれ、k 次の公式になる。

    Σ (1/j)▽jyn+1 = hfn+1  (Σは j = 1, ... , k)

ただし、

0fn = fn,  ▽j+1fn = ▽jfn - ▽jfn-1  (後退差分)

k=1~6 について具体的な式を示すと次のようになる。

    k=1 : yn+1 - yn = hfn+1  (後退オイラー法に一致)
    k=2 : (3/2)yn+1 - 2yn + (1/2)yn-1 = hfn+1
    k=3 : (11/6)yn+1 - 3yn + (3/2)yn-1 - (1/3)yn-2 = hfn+1
    k=4 : (25/12)yn+1 - 4yn + 3yn-1 - (4/3)yn-2 + (1/4)yn-3 = hfn+1
    k=5 : (137/60)yn+1 - 5yn + 5yn-1 - (10/3)yn-2 + (5/4)yn-3 - (1/5)yn-4 = hfn+1
    k=6 : (147/60)yn+1 - 6yn + (15/2)yn-1 - (20/3)yn-2 + (15/4)yn-3 - (6/5)yn-4 + (1/6)yn-5 = hfn+1

BDF法の絶対安定領域を k=1~4 について下図に示す。

ode_bdf

図の曲線の外側が安定領域である。k=1 (後退オイラー法) と k=2 は A安定であるが、それ以外では少し左半平面にはみだしている部分があり A安定ではない。k が大きくなるほど安定領域は狭くなり、k≥7 では原点でも不安定になる。

スティフな方程式を解く際に、実際には多くの場合に完全に A安定である必要はない。そこで、条件を緩めた代わりに高次な解法としてBDF法がよく使われている。4~6次のBDF法は、A安定の条件を緩めたA(α)安定であり、また、硬安定である。

原点を頂点として、左側に非有界であり、実軸とある角度αをもった扇型の領域が絶対安定領域に含まれる場合、A(α)安定であるという(A(90°)安定はA安定に等しい)。方程式の係数行列の固有値がすべてこの中に入っているならば、実質的に A安定と同じである。3次のBDF法はA(86°)安定、4次のBDF法はA(73°)安定、5次のBDF法はA(52°)安定、6次のBDF法はA(18°)安定である。

左半平面がほとんど絶対安定領域に含まれるが、原点付近を含まない虚軸近傍だけが絶対安定領域に含まれない場合を硬安定という。虚軸に近いところに係数行列の固有値ない限り、これも実質的に A安定と同じである。

実験 12

実験 10、11 と同じ例題をXLPackに収録されているBDF法のルーチン Ddebdf を使って解いてみる。要求精度を RTol = 0.01 としたときの y2 をプロットした図を下に示す。

ode_test-12

全く問題なく解いていることがわかる。t = 10 まで計算するのに、試行のためのものも含め関数評価回数はわずか71回であった。

陰的ルンゲ・クッタ法 (Implicit Runge-Kutta (IRK) method)

s段ルンゲ・クッタ法の一般形は次のようなものであった。

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

aij, bi および ci はパラメータで、bi および ci は次式を満たす。

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

aij=0 (i≤j) の場合、陽解法になる。それ以外の場合、陰的ルンゲ・クッタ(IRK)法という。IRK法は A安定なものを作れるが、1 ステップごとに連立非線形方程式を解かなければならないため計算の手間がかかる。

手間を軽減するために aij=0 (i<j) で aii≠0 (少なくても 1つの i について) としたものを半陰的ルンゲ・クッタ法という。対角陰的ルンゲ・クッタ(DIRK)法ともいう。さらに、すべての対角成分が同一の値である場合、単純対角陰的ルンゲ・クッタ(SDIRK)法という。

多くのIRK法の公式が提案されているが、安定性に優れたものとして、A安定な3段ラダウIIA法を下に示す。s段ラダウIIA法は 2s-1次になるため、下の公式は5次である。

    (4-√6)/10 |   (88-7√6)/360    (296-169√6)/1800  (-2+3√6)/225
    (4+√6)/10 | (296+169√6)/1800    (88+7√6)/360    (-2-3√6)/225
        1     |    (16-√6)/36        (16+√6)/36           1/9
    -------------------------------------------------------------
              |    (16-√6)/36        (16+√6)/36           1/9
    

XLPackに収録されている Radau5 ルーチンはこの公式を実装したものである。

実験 13

実験 10~12 と同じ例題をXLPackの Radau5 を使って解いてみる。要求精度を RTol = 0.01 としたときの y2 をプロットした図を下に示す。

ode_test-13

Ddebdf に比べて関数評価回数は多く、t = 10 まで計算するのに、試行のためのものも含め関数評価回数は1500回であった。ただし、求められた解の精度は10-8近くで、要求精度よりもかなり良い。要求精度を10-8にしても関数評価回数は1900回程度とあまり増えないので、このルーチンは高精度計算に適しているのかもしれない。

スティフな方程式を解くためには、スティフ用の解法が使われる。非スティフな問題にも同じ解法を使うことができればよいが、今のところ、効率の観点などから非スティフな問題には通常の解法を使い分けるのが普通である。


参考文献

[1] E. Hairer、他「常微分方程式の数値解法 I」(2007) シュプリンガー・ジャパン
[2] E. Hairer、他「常微分方程式の数値解法 II」(2008) シュプリンガー・ジャパン
[3] 三井斌友「岩波講座 応⽤数学13 [⽅法3] 微分⽅程式の数値解法 I」(1993) 岩波書店
[4] 森正武、「計算機のための数値計算法」科学技術出版社 (1978)
[5] 名取亮「数値解析とその応⽤」(1990) コロナ社
[6] 森⼝繁⼀「数値計算⼯学」(1989) 岩波書店

Top