ルンゲ＝クッタ法
ルンゲ=クッタ法の局所誤差を精確に計算することが難しいので、実践では誤差を一定の範囲にコントロールするのが望ましいである。そのために開発されたのが 埋め込み型ルンゲ=クッタ法(Embedded Runge-Kutta method)である。適応型ルンゲ=クッタ法 (adaptive Runge–Kutta method) とも呼ばれる。線型多段法にもミルンデバイスと呼ばれる相似の方法が存在する。
埋め込み方法は陽的ルンゲ=クッタ法二つを利用する(下の方法を上の方法に埋め込むように見えるので埋め込み型と呼ばれる)。ブッチャー配列の以下のように拡張する。
ここで、bi は p 次陽的方法に対応し、b *i は p-1 次陽的方法に対応する。二つの方法は係数 aij と ci を共用する。正しく係数を選ぶと、二つの方法をともに収束させることができる。そのとき、埋め込み方法の時刻 tn での相対(局所)誤差 en は次の公式で与えられる。
e n = y n − y n ∗ = h ∑ i = 1 s ( b i − b i ∗ ) k i = O ( h p ) {\displaystyle e_{n}=y_{n}-y_{n}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i}=O(h^{p})}
ルンゲ=クッタ法の収束性から、相対誤差も0に収束することがわかる。埋め込みルンゲ=クッタ法は、アルゴリズムを用いて一時刻ごと(自動的)に刻み幅 h を調整し、誤差をコントロールすることができる(故に適応型と呼ばれる)。
y ( t n + 1 ) = y n + ∫ y n y n + 1 f ( t , y ( t ) ) d t = y n + h ∫ 0 1 y ( t n + h τ , y ( t n + h τ ) ) d τ = y n + h ∑ i = 1 s b i f ( t n + c i h , y ( t n + c i h ) ) {\displaystyle y(t_{n+1})=y_{n}+\int _{y_{n}}^{y_{n+1}}f(t,y(t))dt=y_{n}+h\int _{0}^{1}y(t_{n}+h\tau ,y(t_{n}+h\tau ))d\tau =y_{n}+h\sum _{i=1}^{s}b_{i}f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))}
k i = f ( t n + c i h , y ( t n + c i h ) ) {\displaystyle k_{i}=f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))} とおく。ki あるいは y(tn + ci h) を適切に(線型結合として)近似することでルンゲ=クッタ法の公式となる。その上、係数をテイラー展開より正しく選択すると、方法の収束性も求積法の収束性より保証される。しかし、局所誤差のオーダーや上界は、方法によって大きく異なるので、方法別に計算しなければならない。
ここで、
k 1 = f ( t n , y n ) , k 2 = f ( t n + h 2 , y n + h 2 k 1 ) , k 3 = f ( t n + h 2 , y n + h 2 k 2 ) , k 4 = f ( t n + h , y n + h k 3 ) {\displaystyle {\begin{aligned}k_{1}&=f\left(t_{n},y_{n}\right),\\k_{2}&=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{1}\right),\\k_{3}&=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{2}\right),\\k_{4}&=f\left(t_{n}+h,y_{n}+hk_{3}\right)\end{aligned}}}
である。次の値 (yn+1) は、現在の値 (yn) に増分を加えたものであり、増分は勾配の推定値に間隔 h を乗じたものになっている。勾配の推定値は、k1, ..., k4 の4つの勾配の重み付け平均で求める。k1, ..., k4 のそれぞれの勾配は、特定の (t, y) に対する f によって与えられ、以下のように解釈できる。