10. 常微分方程式 (2) 誤差の推定, 補外法, ステップ幅の自動調節, 密出力, 2階常微分方程式, 遅延微分方程式

目次

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

誤差の推定補外法ステップ幅の自動調節密出力2階常微分方程式遅延微分方程式

誤差の推定

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

1段法の誤差

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

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

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

いま, 次の初期値問題を考える.

  
  dy/dt = f(t, y),  t = t0 において y = y0

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

そうすると, 局所的離散化誤差 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)

 

数値実験 (10)

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

  
  y' = y, y(0) = 1

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

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

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

ステップ幅を変化させて, t = 1 における相対誤差をプロットした結果を下に示す. 横軸はステップ幅(右にいくほど小さい), 縦軸は相対誤差である. なお, 丸め誤差の影響を見やすくするために計算は単精度で行った.

図中の青色のプロットのように, この場合 10-5 付近を境に丸め誤差が支配的になることがわかる.

ルンゲ・クッタ・ジル法の項目で丸め誤差を減らす方法を紹介したが, それをオイラー法に適用して計算した結果を赤色で表示した. 丸め誤差がなければhに比例して誤差が減っていくことがわかる.

リチャードソンの補外

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)と同じである.

数値実験 (11)

リチャードソンの補外
簡単な微分方程式の初期値問題の例をオイラー法, ホイン法, ルンゲ・クッタ法で計算する. 次に, リチャードソンの補外によりそれぞれの計算結果の誤差を推定してプロットする.

  
  y' = y, y(0) = 1

 

この例では全体として少なめに見積もられているが比較的近い値が得られている.

補外法

H = t – t0 とする. 正の整数列 n1 < n2 < n3 < … について hi = H/ni とすると, 整数列に対応するステップ幅 h1 > h2 > h3 > … が得られる.

p次の数値解法においてステップ幅hjを使ってnjステップ計算して求めたtにおける解をTj1とする. T11 ~ Tk1 (k <= j) の値を使って連立方程式を解けばhのk次補間多項式を作ることができる. そこで h = 0 として補外値を求めればよい近似を得ることができる.

実際の計算においては, 補間多項式そのものは不要で h = 0 における値のみあればよいので, 次のネヴェルのアルゴリズムとよばれる方法を使うことができる.

  
  Tjk = Tj(k-1) + (Tj(k-1) - T(j-1)(k-1))/(nj/nj-k+1 - 1)

 
この漸化式を使って求めた Tjk は解の良い近似となっている. この方法はp+k-1次の数値解法となっており, 補外法(Extrapolation method)とよばれる.

補外法の計算結果を表す次の表を補外表という.

  
  T11   |
  T21   |  T22
  T31   |  T32  T33
  T41   |  T42  T43  T44
  ...   |  ...  ...  ...  ...
 数値解 |  k-1個の補外値 →

 
T21についてはT22まで(k = 2), T31はT33まで(k = 3), …, Tj1はTjjまで(k = j)補外できる. すなわち, 最大のkを使用するならばjが増えるほど次数も大きくなっていく.

上において, k = 2, n1 = 1, n2 = 2 の場合はリチャードソンの補外に一致する.

整数列niのとり方であるが, 次のようにいくつかの方法が提案されている.

ロンバーク列: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, …
ブリアシュ列: 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, …
調和数列: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, …

数値実験 (12)

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

  
  y' = y, y(0) = 1

 
数値解法としてはオイラー法(p = 1次)を使った. 整数列としてはロンバーク列 1, 2, …, 128 を使った.

横軸は相対誤差(右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットとした. 補外に必要な計算量は無視した. 一番左のラインがオイラー法による数値解(Tj1)である. 一番右の赤いラインは最高次数で求めた補外値(Tjj)である.

数値実験 (13)

補外法の計算例 (2)
同じ簡単な微分方程式の初期値問題の例で補外法のパラメータ(kの値, 整数列)を変えた場合の挙動を確認する.

  
  y' = y, y(0) = 1

 
数値解法としてはオイラー法(p = 1次)を使った. ロンバーク列で k = j と k = 4 の場合, ブリアシュ列で k = j, 調和数列で k = j の場合を比較する.

横軸は相対誤差(右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットとした. 補外に必要な計算量は無視した. ロンバーク列で k = 4 とした(補外表の横幅を4に制限した)場合には途中からこの(対数)プロットで直線となり次数が増えなくなったのがわかる. ブリアシュ列と調和数列はこの例では最初はよかったが途中で失速して誤差が減らなくなった.

GBSアルゴリズム

補間多項式が偶数のべき乗のみで表される場合, ネヴェルのアルゴリズムは次のようになり収束が速くなる.

  
  Tjk = Tj(k-1) + (Tj(k-1) - T(j-1)(k-1))/((nj/nj-k+1)^2 - 1)

 
GBS(Gragg-Bulirsch-Stoer)アルゴリズムはこれを利用したもので, 次の手順で計算を行う.

  
  y1 = y0 + hjf(t0, y0)                           [オイラー法]
  yi+1 = yi-1 + 2hjf(ti, yi) (i = 1, 2, ..., nj)  [中点則]

  Tj1 = (1/4)(ynj-1 + 2ynj + ynj+1)

 
このようにすると Tj1 は偶数のべき乗のみで表され, 上式を適用して効率よく補外値を求めることができる. なお, 中点則には不安定性の問題があることを「使われなくなった解法」の項で説明したが, この場合は問題は生じない.

整数列niは偶数になるようにする.

ロンバーク列: 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, …
ブリアシュ列: 2, 4, 6, 8, 12, 16, 24, 32, 48, 64, …
調和数列: 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, …

数値実験 (14)

GBSアルゴリズムによる補外法の計算例
最も簡単な微分方程式の初期値問題の例でGBSアルゴリズムによる補外法と通常の補外法の精度を比較する.

  
  y' = y, y(0) = 1

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットである. GBSアルゴリズムの方が明らかに速く誤差が減少している. なお, 図で線が切れているところは誤差が0になったためプロットされていない区間であり, 精度いっぱいで解が求められていることを示している.

ステップ幅の自動調節

埋め込み型ルンゲ・クッタ法

誤差の推定を行うのにリチャードソンの補外を使うためにはステップ幅を変えて2回計算を行わなければならず, 計算量が大きくなるのが欠点である. しかし, ルンゲ・クッタ法自身に誤差推定機能を付加した公式を作ることができれば, 1回の計算で誤差推定値も得ることができ有用である. 特に, あらかじめ要求精度を設定しておき, それに必要なステップ幅に自動調節するプログラムを作るのに必須である.

ルンゲ・クッタ法は次のように表された. ここで, この公式は p次であるとする.

  
  ki = f(yn + hΣaijkj, tn + cih) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
  yn+1 = yn + hΣbiki (Σ は 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)

 
このような公式は実際に作ることができて, 埋め込み型ルンゲ・クッタ法(embedded Runge-Kutta method) とよばれる. 主な公式を以下に示す.

ルンゲ・クッタ・マーソン法

埋め込み型ルンゲ・クッタ法で最初に提案されたのはマーソンの公式 (ルンゲ・クッタ・マーソン法(Runge-Kutta-Merson method)) とされる. 係数は次のとおりである.

  
  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 
    

 
通常のルンゲ・クッタ法のブッチャー配列に対してy*のための係数が1行追加されている.

計算式は次のとおりである.

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

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

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

5(4)次 ルンゲ・クッタ・フェールベルグ法

5(4)次 ルンゲ・クッタ・フェールベルグ法(Runge-Kutta-Fehlberg method)は埋め込み型ルンゲ・クッタ法の中でも有名な公式である. RKF45というプログラムが広く使われている.

  
  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)次 ドルマン・プリンス法

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)次 ルンゲ・クッタ・ヴァーナー法

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次のアダムス法では, 予測子(右肩に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法などを使う必要がある.

数値実験 (15)

ステップ幅自動調節機能付きプログラムの計算例
ルンゲ・クッタ・マーソン法を使い簡単なステップ幅自動調節機能付きプログラム(RKM)を作成して実験してみる.

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

  • ステップ幅の初期値と局所離散化誤差の許容限界を与え計算を開始する
  • マーソンの公式で1ステップ計算後, 推定誤差が許容限界より大きければステップ幅を1/2にしてやり直す(失敗ステップ)
  • 推定誤差が許容限界より小さければ1ステップ進める(成功ステップ). ただし, 推定誤差が許容限界の1/32より小さいときには次ステップからステップ幅を2倍にする

2種類の微分方程式の初期値問題例を使ってステップ幅の自動調節機能の動作を確認する. 比較のために, 4次のルンゲ・クッタ法(RK4)および5(4)次のルンゲ・クッタ・フェールベルグ法(RKF45)(XLPackのDerkfを使用)による計算結果も示す.

横軸は相対誤差(右に行くほど精度がよい. n = 2 の場合はそれぞれの値の平均), 縦軸は関数評価回数の対数プロットである. RKMとRKF45では要求精度(Tol)を少しづつ変えて計算を行い, 得られた結果の実際の誤差と関数評価回数をプロットした. RK4ではステップ幅を変えていき, そのときの誤差と関数評価回数をプロットした. なお, RKMではステップ幅の初期値を50*Tolとした.

(例1)

  
  y' = y, y(0) = 1

 
例1は変化が少ない方程式の例で, 解析解は y = et である. t = 1 における値 (= e) を求める.

(例2)

  
  y1' = 1 + y12y2 - 4y1
  y2' = 3y1 - y12y2
  y1(0) = 1.5, y2(0) = 3

 
例2は文献[3]のブラセレータとよばれる変化の激しい方程式の例で n = 2 である. t = 16 における値を求める.

RKMとRK4はどちらも4次であるが, 例1ではRK4の方が速く, 例2ではRKMの方が速い. これは, 変化が少ない関数においてはステップ幅の調節を行う必要がないが, RKMでは調節のために無駄なステップを費やしているためと思われる. 数値実験(16)で詳細を確認する.

数値実験 (16)

ステップ幅自動調節の動作の詳細
上の2つの例を計算したときのRKMの動作を詳細に見てみる.

(例1)

Tol = 1.0e-7 としたときの動作. 上段は計算値, 下段はステップ幅である. ステップ幅の初期値は 50*Tol とした. 丸は各ステップを示す.

関数が滑らかなので h = 0.082 になった後は変更されなかった. ステップ幅の初期値が小さすぎたようなので工夫が必要であろう. なお, 最後の1ステップでステップ幅が変化しているのは誤差の問題ではなく, 終点をちょうど1に合わせるために調整されたものである.

(例2)

Tol = 1.0e-4 としたときの動作. 上段は計算値, 下段はステップ幅である. ステップ幅の初期値は 50*Tol とした. 丸は各ステップを示す.

急な変化をするところではステップ幅を小さく, ゆるやかな変化をするところではステップ幅を大きく調節しているのがわかる. 成功ステップが97回で計算を終えているが, 失敗ステップが25回あった. すなわち, 2割ほど無駄ステップを踏んでいるがそれでもRK4よりは効率がよかったことになる.

この例でRK4の方が遅かったのは, 全体の精度を確保するためには急な変化をするところで必要なステップ幅に合わせる必要があり, ゆるやかな変化をするところで必要以上にステップ幅が小さくなり無駄な計算を行ったためと考えられる.

密出力

一定間隔あるいは決められた点で結果を出力したい場合には, ステップ幅をそれに合わせたり, ステップ幅自動調節機能付きプログラムであれば出力が必要な点ごとに積分を終了させたりしなければならないが, 特に出力点の間隔が狭い場合には効率が悪くなってしまう. そこでステップ幅は必要精度の面から最適なものを選び, ステップに一致しない出力点においては計算を行うことを避け補間値を使う方法が考えられる. これを密出力(dense output)などとよぶ.

多段法

アダムス・モルトン法を例にとる. アダムス・モルトン法で使用した説明図を再掲する.

tnの次のステップ tn+1 = tn + h における値を求める際には図のように補間多項式p*(t)を求め, これを区間[tn, tn+1]で積分することにより値を求めた.

  
  yn+1 = yn + ∫ p*(x) dx [tn, tn+1]

 
同様に, tn + θh (0 ≦ θ ≦ 1) における値 y(tn + θh) は同じ補間式を区間[tn, tn + θh]で積分すれば求めることができる. 新たに関数評価を行う必要がないし精度が落ちることもない.

数値実験 (17)

アダムス・モルトン法における密出力
微分方程式の初期値問題の例で4次アダムス・モルトン法で t = 0.4, 0.5, …, 1.0 におけるyの値を求める. t = 0.1, 0.2, 0.3 における値(出発値)は4次のルンゲ・クッタ法で求めた. t = 0.35, 0.45, …, 0.95 における値は上で説明した補間により求める.

  
  y' = y, y(0) = 1

 

青い点がアダムス・モルトン法による計算値で y = exp(t) のグラフとなる. 緑の点は補間により求めた点を示す. 赤は計算値および補間値の相対誤差を表し, 値は右側の対数目盛りによる. 期待通り誤差が増えることなく補間値が求められている.

ルンゲ・クッタ法

ルンゲ・クッタ法は1段法なので過去の値は使用しないためすぐに使える補間式はない. そこで, 負担の少ない, すなわちこのための追加の関数評価がないまたはごく少数回であるような, 密出力のための公式が研究されてきた. 例えば, tn + θh (0 ≦ θ ≦ 1) における値 y(tn + θh) を求める次の形の式である.

  
  y(tn + θh) = yn + hΣbi(θ)ki (Σ は i = 1, 2, ... , s*)
  ki = f(tn + cih, yn + hΣaijkj) (Σ は j = 1, 2, ... , i - 1) (i = 1, 2, ... , s*)

 
ルンゲ・クッタ法の公式と同じ形をしているが, 係数biは定数ではなくθに依存する. このようなルンゲ・クッタ法の自然な拡張を連続ルンゲ・クッタ法(Continuous Runge-Kutta (CRK) method)ともいう.

4次のルンゲ・クッタ法については, s* = s = 4 に対して密出力公式の次数 p* = 3 となる次の解が得られている.

  
  b1(θ) = θ - (3/2)θ2 + (2/3)θ3
  b2(θ) = b3(θ) = θ2 - (2/3)θ3
  b4(θ) = -(1/2)θ2 + (2/3)θ3

 
kiは4次ルンゲ・クッタ法の公式と同じ, すなわち, 新たな関数評価は必要としない. なお, 4次ルンゲ・クッタ法は関数評価を追加しなければ p* = 4 とはならない.

通常は p* = p – 1 であれば十分である. 導関数 y'(tn + θh) が正確な必要がある場合には p* = p でなければならない.

数値実験 (18)

ルンゲ・クッタ法における密出力
微分方程式の初期値問題の例で4次ルンゲ・クッタ法で t = 0.1, 0.2, …, 1.0 におけるyの値を求める. t = 0.05, 0.15, …, 0.95 における値は上で示した3次の密出力公式を使って求める.

  
  y' = y, y(0) = 1

 

青い点がルンゲ・クッタ法による計算値で y = exp(t) のグラフとなる. 緑の点は補間により求めた点を示す. 赤は計算値および補間値の相対誤差を表し, 値は右側の対数目盛りによる. 補間値の誤差が少し大きい傾向が見られた.

次に同様の計算を Dopri5 (5(4)次ドルマン・プリンス法) により行う. ドルマン・プリンス法では5次の密出力公式も提案されているが, Dopri5では別の4次の密出力公式を採用している.

ここではあえて補間値の誤差が少し大きい傾向が見られる Tol = 1e-8 の例を示した. Tolをもっと小さくしていくと誤差の差はずっと小さくなった.

2階常微分方程式

実際に応用で現れる方程式は2階常微分方程式が多い.

  
  y'' = f(t, y, y')
  y(t0) = y0, y'(t0) = y'0 (初期条件)

 
これは「概要」の項で説明したように1階の連立常微分方程式に変換して解くことができる.

ここでは, 次のようにf()がy’に依存しない特殊な場合に効率のよい方法を紹介する.

  
  y'' = f(t, y)
  y(t0) = y0, y'(t0) = y'0 (初期条件)

 

ニュストレム法

ルンゲ・クッタ法に似ているが係数の決め方が異なるニュストレム法(Nystrom method)が知られており, 次のように計算を行う.

  
  k'i = f(tn + cih, yn + cihy'n + h2Σaijk'j) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
  yn+1 = yn + hy'n + h2Σbik'i (Σ は i = 1, 2, ... , s)
  y'n+1 = y'n + hΣbik'i (Σ は i = 1, 2, ... , s)

 
係数は次のとおりである.

4次ニュストレム法

  
  0   |
  1/2 | 1/8         aij
  2/3 |  0   1/2
  ---------------------
  bi  | 1/6  1/3  0 
  bi  | 1/6  4/6  1/6
    

 

5次ニュストレム法

  
  0   |
  1/5 |  1/50         aij
  2/3 | -1/27   7/27
  1   |  3/10  -2/35    9/35
  -------------------------------------
  bi  | 14/336 100/336 54/336  0 
  bi  | 14/336 125/336 162/336 35/336
    

 
ルンゲ・クッタ法に比べて段数が少ないことがわかる(4次: 4段→3段, 5次: 6段→4段).

補外法

1階の連立常微分方程式に変換してGBSアルゴリズムを適用する.

  
  y1 = y0 + hjy'0
  y'1 = y'0 + hjf(t0, y0)

  yi+1 = yi-1 + 2hjy'i (i = 1, 2, ..., nj)
  y'i+1 = y'i-1 + 2hjf(ti, yi) (i = 1, 2, ..., nj)

  Tj1 = (1/4)(ynj-1 + 2ynj + ynj+1)
  T'j1 = (1/4)(y'nj-1 + 2y'nj + y'nj+1)

 
この場合実は次のようになるので半分だけ(奇数添字のときだけ)関数評価すればよくなる.

  
  Tj1 = ynj
  T'j1 = (1/2)(y'nj-1 + y'nj+1)

 

数値実験 (19)

2階常微分方程式(y'' = f(t, y))の計算例
次の2階常微分方程式の初期値問題の例を4次, 5次ニュストレム法および補外法(GBS, ロンバーク列)により計算する.

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

 

横軸は相対誤差(右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 4次のルンゲ・クッタ法の計算結果も示した.

4次ニュストレム法は段数は少ないが4次のルンゲ・クッタ法と同程度の精度であることが確認できる. 関数評価回数が少ない分だけプロットは下にきている.

遅延微分方程式

遅延微分方程式(Delay differential equations (DDE))とは次のような微分方程式をいう.

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

 
τは遅延時間を表す正定数である. 複数の遅延時間が含まれるなどより一般的なものも考えられる.

これをルンゲ・クッタ法で解くためにはステップ幅を τ = kh (kは整数)となるように選んで, 計算したy(t – τ)の値をとっておくようにすればよい. しかし, ステップ幅を自動調節するルーチンの場合にはうまくいかない. これは密出力を行う方法と同じく, 多段法か密出力機能を持ったルンゲ・クッタ法ならば計算できる. XLPackに収録されているRetardはドルマン・プリンス法ルーチンDopri5を遅延微分方程式用に修正したものである. Dopri5では1ステップ前までの過去の値を求められたが, Retardではずっと前までの値を得ることができる.

遅延微分方程式の計算例は「実用ルーチンのベンチマーク」を参照せよ.


参考文献

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