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

誤差の推定

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

1段法の誤差

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

• 離散化誤差 (打ち切り誤差ともいう)
  演算がすべて無限精度で行われたとしたときに, 解法に起因して生じる誤差
  • 局所的離散化誤差
    それまでに得られている値が正確であるとして, 次の1ステップで生じる誤差
  • 大域的離散化誤差
    初期値から計算される解と真の解との差
• 丸め誤差
  数値を浮動小数で表し有限桁数で演算することにより生ずる誤差

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

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

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

そうすると, 局所的離散化誤差 d_n は次のように表すことができる.

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

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

  
  |dn| ≦ Chp+1

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

  
  en = yn - y(tn)

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

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

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

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

ここで, 丸め誤差 ε_n のノルムはある数 ε 以下であるとする.

  
  ||εn|| ≦ ε

 
いま, 最終点を t_N とすると, ステップ数 N = (t_N – t₀) / h である. 最終的な大域的離散化誤差を e_N ≒ Σd_i (i = 0 ～ N-1) と近似し, 丸め誤差を Nε と近似すると, 最終的な全誤差は次のように近似することができる.

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

数値実験 (10)

  
  y' = y, y(0) = 1

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

リチャードソンの補外

p次の数値解法の与えられた点 t_n での大域的離散化誤差は O(h^p) となるので, 次のように表すことができる.

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

 
ステップ幅 h で計算した y_n を y_n^(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)

 
となり, y_n^(h) の誤差は O(h^q) を無視すると次のように近似できる.

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

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

数値実験 (11)

  
  y' = y, y(0) = 1

補外法

H = t – t₀ とする. 正の整数列 n₁ < n₂ < n₃ < … について h_i = H/n_i とすると, 整数列に対応するステップ幅 h₁ > h₂ > h₃ > … が得られる.

p次の数値解法においてステップ幅h_jを使ってn_jステップ計算して求めたtにおける解をT_j1とする. T₁₁ ～ T_k1 (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)

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

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

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

 
T₂₁についてはT₂₂まで(k = 2), T₃₁はT₃₃まで(k = 3), …, T_j1はT_jjまで(k = j)補外できる. すなわち, 最大のkを使用するならばjが増えるほど次数も大きくなっていく.

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

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

ロンバーク列: 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

数値実験 (13)

  
  y' = y, y(0) = 1

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)

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

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

ロンバーク列: 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)

  
  y' = y, y(0) = 1

ステップ幅の自動調節

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

誤差の推定を行うのにリチャードソンの補外を使うためにはステップ幅を変えて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)

 
同じ a_ij と c_i を使い (すなわち, 同じ k_i を使い), b_i を b^*_i に変えた時の近似値 y^*_n+1 がp^*次 (例えば, p^* = p – 1, p^* = p + 1 など) になっているように定数を決めることができれば, k_i を計算し直すことなく次式により 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 (y_n+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番目(k₇)の係数 a_7j が近似値 y を計算するための係数 b_i に等しいことである. これは FSAL (First Same As Last) とよばれるアイデアで, k₇ が次の k₁ = f(y_n+1, t_n+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に採用されている公式である.

多段法の誤差

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

例えば, 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)

  

  y' = y, y(0) = 1

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

数値実験 (16)

密出力

次の初期値問題の解を求める.

  
  y' = y, y(0) = 1

 
解析解は y = e^t である.

ステップ幅を自動調節するプログラム Dopri5(5(4)次 ドルマン・プリンス法) と Derkf(5(4)次 ルンゲ・クッタ・フェールベルグ法) を使って, 要求精度(RTol, ATol) = 1.0e-7 として, t = 0 ～ 1 で 0.05 ごとに 0.0 ～ 0.05, 0.05 ～ 0.1, …, 0.95 ～ 1.0 というように終点を変えながら繰り返し解を求めた. このときの誤差をプロットすると次のようになった. 青の実線と左側の目盛は関数形を表す. 点線と右側の目盛は各解法による誤差を表す.

Dopri5もDerkfも 0.05 幅の各区間では1ステップで解を求め, 誤差は 1.0e-10 付近と要求精度よりかなりよい精度で解を求めた. これは, 自動調節されるべき要求精度に見合ったステップ幅よりも 0.05 という幅が小さく, 終点に合わせるためにステップ幅が常に 0.05 に調節されたためである. すなわち, 要求精度に比べて公式の精度が良すぎるのである. 比較のために4次のルンゲ・クッタ法(RK4)でステップ幅を 0.05 として計算してみると, 誤差はちょうど 1.0e-7 付近になり, RK4の精度程度でよいことがわかる. 関数評価回数は, Dopri5が128回, Derkfが124回で, RK4は80回となった. Dopri5とDerkfはこの場合には高精度なゆえに無駄な計算をしていることになる.

そこで, ステップ幅は必要精度の面から最適なものを選び, ステップに一致しない出力点においては計算を行わずに補間値を使うことにより無駄な計算を避ける方法が考えられる. これを密出力(dense output)とよぶ.

多段法の場合

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

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

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

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

数値実験 (17)

ルンゲ・クッタ法の場合(連続ルンゲ・クッタ法)

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

 
ルンゲ・クッタ法の公式と同じ形をしているが, 係数b_iは定数ではなくθに依存する. このようなルンゲ・クッタ法の自然な拡張を連続ルンゲ・クッタ法(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

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

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

数値実験 (18)

ステップ幅の自動調節を行う場合

上の2つの例でみた補間を使う方法は, ステップ幅の自動調節を行うプログラムに適用すると効率を上げることができる. すなわち, ステップ幅は必要精度の面から最適なものを選び, ステップに一致しない出力点においては補間値を使うことにより無駄な関数計算を避けることができる.

多段法の例として, XLPackにも収録されている Deabm (アダムス・バシュフォース・ムルトン法) は, 補間値を使うことができる場合には補間値を使うように作られている.

ルンゲ・クッタ法の例として, Dopri5には4次の補間公式を使った密出力機能が組み込まれており, XLPackのDerkfには5次の補間公式を使った密出力機能が組み込まれている. Dopri5とDerkfの場合には密出力のための追加関数評価を必要としない公式を採用している. XLPackに収録されている他の多くのプログラムでも密出力機能が組み込まれている.

Dopri5とDerkfを使って要求精度(RTol, ATol) = 1.0e-7 で計算を行い, 0.02ごとに細かく結果の出力を行った例を以下に示す.

緑の丸が計算ステップで, 紫の三角がその時の誤差である. 要求精度を満たすにはステップ幅は0.1～0.15程度で十分であることがわかる. Dopri5の場合はこの例では精度 1.0e-8 程度を達成しており, やや安全側に出ている. Derkfではちょうど 1.0e-7 程度に調節されている. 必要としたステップ数は, Dopri5が9ステップ, Derkfが7ステップであった. 関数評価回数は, Dopri5が56回, Derkfが46回であった.

青の丸は補間値(密出力), 赤の三角はその誤差である. 補間値も十分な精度で求められている. 密出力においては補間値をどれだけ細かな間隔で求めても関数評価回数は増えない(補間計算程度の計算量しか増えない)ので, 精密なグラフを描くために細かく値を出力したい場合などに有効である. 逆に, 必要なステップ幅が出力点の間隔よりも小さくなる場合には密出力の効果は小さくなるので, 単純に繰り返し計算を行ってもよい.

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)

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

遅延微分方程式

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

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

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

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

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