ルンゲ=クッタ法

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動先: 案内検索

ルンゲ=クッタ法: Runge–Kutta method)とは、数値解析において微分方程式の初期値問題に対して近似解を与える一連の方法である。この技法は1900年頃に数学者カール・ルンゲマルティン・ヴィルヘルム・クッタ英語: Martin Wilhelm Kutta)によって発展された。

古典的ルンゲ=クッタ法[編集]

一連のルンゲ=クッタ公式の中で最も広く知られているのが、古典的ルンゲ=クッタ法 (RK4、もしくは単に狭義の ルンゲ=クッタ法: the (classical) Runge–Kutta method) などと呼ばれる4次の公式である。

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

但し、y(t) が近似的に求めたい未知関数であり、その t における勾配は f (t, y) によって t 及び y(t) の関数として与えられている。時刻 t0 における初期値は y0 で与えられている。

今、時刻 tn における値 yn = y(tn) が既知のとき、十分に小さなステップ幅 h に対して yn+1, tn+1 を以下の式で与えると、yn+1y(tn+1) の 4次精度の近似になっている。このステップを逐次的に繰り返すことによって、初期値 y0 から任意の時刻 tn における近似値 yn が求められる。

ここで、

である。次の値 (yn+1) は、現在の値 (yn) に増分を加えたものであり、増分は勾配の推定値に間隔 h を乗じたものになっている。勾配の推定値は、k1, ..., k4 の4つの勾配の重み付け平均で求める。k1, ..., k4 のそれぞれの勾配は、特定の (t, y) に対する f によって与えられ、以下のように解釈できる。

  • k1 は区間の最初 tn における勾配である (オイラー法で用いる勾配に一致する)。
  • k2 は区間の中央 tn + h/2 における勾配の近似値である (中点法で用いる勾配)。計算に用いる中央の y の値は、初期位置の勾配 k1 を用いてオイラー法で推定する。
  • k3 は区間の中央における勾配のもう一つの近似値である。中央の yk2 の値から推定して用いる。
  • k4 は区間の最後 tn + h における勾配の近似値であり、k3 の値から推定された最後の点の y の値を用いる。

重み付き平均では、中央の勾配に対して大きな重みを用いる。シンプソン則を用いた平均と同等の形になる。

RK4は4次の方法である。厳密解とRK4のテイラー展開が4次の項まで一致し、1ステップの推定誤差は O(h5) のオーダーになる。目的の時刻の y を求めるのに必要なステップ数は O(1/h) になるので、全体の推定誤差は O(h4) になる。

ルンゲ゠クッタ法[編集]

前述の RK4 の一般化として、以下の形式を持つ s 段のルンゲ゠クッタ法を構成することができる。整数 s をそのルンゲ゠クッタ法の段数 (stage) という。

但し、

(文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する。)

具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度持つように適切に選ばれた係数 aij (ルンゲ゠クッタ行列)、bi (重み)、ci (節点) で指定される (但し、1 ≤ i, js)。特に ij に対して aij=0 を満たす方法が広く用いられ、総称して 陽的ルンゲ゠クッタ法 (ERK、: explicit Runge–Kutta methods) と呼ぶ。そうでないものを 陰的ルンゲ゠クッタ法 (IRK、: implicit Runge–Kutta methods) と呼ぶ。

近似値 ynyn+1 から計算するときに発生する誤差の大きさが O(hp+1) のとき、そのルンゲ゠クッタ公式は p 次精度を持つといい、p次数 (または 位数) と呼ぶ。p 次のルンゲ゠クッタ公式は、誤差の大きさの条件に誤差の表式を代入し、係数の条件を求めることによって得られる。例えば、2 段の陽的方法が2次精度を持つための係数に対する条件は、b1 + b2 = 1, b2c2 = 1/2, かつ b2a21 = 1/2 である。この条件を満たす範囲内で様々な方法を考え得る。

これらの係数を分かりやすく表記する方法として、以下のような形式のブッチャー配列 (Butcher tableau) が知られている:

実際には、それぞれのルンゲ゠クッタ法について各要素に具体的な値を入れて用いる。陽的な方法ではルンゲ゠クッタ行列の上三角成分は常に 0 になるので表記を省略する。例えば RK4 はブッチャー配列を用いて以下のように表現される。

ルンゲ゠クッタ法が一貫しているためには、次の条件が満たされている必要がある。

陽的ルンゲ゠クッタ法[編集]

陽的ルンゲ゠クッタ法では aij=0 (ij) が満たされるので、ブッチャー配列は以下の形になる。

0


この時、各勾配 k1, ..., ks を以下のように逐次的に求めることができる。


陽的ルンゲ゠クッタ法が目的の精度 p になる係数を持つためには、十分に大きな段数 s が必要になる。少なくとも次数以上の段数が必要 (sp) であり、更にp = 5 次以上の場合には段数は次数よりも大きく取らなければならない (s > p) ことが知られている。各次数 p に対する最低段数 s の一般式は未解決問題である。具体的な値が幾つか知られている:


1段1次の公式[編集]

最も単純なルンゲ゠クッタ法が 1段1次 の方法であり、一意に定まる。(前進) オイラー法と等価になる。

以下のブッチャー配列で表現される。

0
1


2段2次の公式[編集]

2段2次の方法は1パラメータの自由度 α を持ち、以下のブッチャー配列で与えられる。

0


α = 1/2 の場合が中点法 (改良オイラー法)

に対応し、以下のブッチャー配列で与えられる。

0
1/2 1/2
0 1


α = 1 の場合はホイン法 (修正オイラー法)として知られ、ブッチャー配列は以下の通りである。

0
1 1
1/2 1/2


α = 2/3 の場合はRalston 法と呼ばれ、ブッチャー配列は以下の通りである。

0
2/3 2/3
1/4 3/4


4段4次の公式[編集]

古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる:

0
1/2 1/2
1/2 0 1/2
1 0 0 1
1/6 1/3 1/3 1/6


修正版として クッタの3/8公式 (: Kutta's 8/3-rule) が知られている。

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


対象の変数 y が多成分の巨大なベクトルの場合、計算機上で記憶するデータ量が最小で済むように構成された ルンゲ゠クッタ゠ギル法 (: Runge–Kutta–Gill method) も用いられる。

0
1/2 1/2
1/2 1/2 - α- α-
1 0 1 - α+ α+
1/6 α-/3 α+/3 1/6

但し、α± = 1 ± 1/√2 である。以下の手順で計算することにより、記録領域を3ベクトル (y, k, q) に収めることができる (RK4 や 3/8 公式では少なくとも4ベクトル必要になる)。但し、初期状態として yyn が保持されているとし、この手順により y に次の値 yn+1 が設定される。

  1. kh f (x, y),
  2. (y, q) ← (y + k/2, k),
  3. kh f (x + h/2, y),
  4. (y, q) ← (y + α- (k - q), q + α- (2 k - 3 q)),
  5. kh f (x + h/2, y),
  6. (y, q) ← (y + α+ (k - q), q + α+ (2 k - 3 q)),
  7. kh f (x + h, y),
  8. (y, q) ← (y + (1/6)(k - 2 q), 0).

誤差と収束性[編集]

関数階連続偏微分可能であるならば、次RK法から得られた近似解と厳密解の局所誤差にはステップ幅をパラメータとした誤差限界が存在する。

しかし、この誤差限界はが高次になると意味を成さなくなる。そのため、ステップ幅の選定には様々な補外法が提案されている。

計算例[編集]

Ralston法による例を示す。Ralston 法は以下のブッチャー配列で与えられる:

0
2/3 2/3
1/4 3/4

対応する方程式は以下のようになる。

初期値問題

をステップ幅 h = 0.025 のステップ 4 回で解く。

次の過程で計算が進められる。下線を付した値が数値解に対応する。

時刻 t0 = 1
y0 = 1 (初期条件).
時刻 t1 = 1.025
y0 = 1,
k1 = f (t0, y0) = 2.557407725,
k2 = f (t0 + (2/3) h, y0 + (2/3) h k1) = 2.7138981184,
y1 = y0 + h ((1/4) k1 + (3/4) k2) = 1.066869388.
時刻 t2 = 1.05
y1 = 1.066869388,
k1 = f (t1, y1) = 2.813524695,
k2 = f (t1 + (2/3) h, y1 + (2/3) h k1),
y2 = y1 + h ((1/4) k1 + (3/4) k2) = 1.141332181.
時刻 t3 = 1.075
y2 = 1.141332181,
k1 = f (t2, y2) = 3.183536647,
k2 = f (t2 + (2/3) h, y2 + (2/3) h k1),
y3 = y2 + h ((1/4) k1 + (3/4) k2) = 1.227417567.
時刻 t4 = 1.1
y3 = 1.227417567,
k1 = f (t3, y3) = 3.796866512,
k2 = f (t3 + (2/3) h, y3 + (2/3) h k1),
y4 = y3 + h ((1/4) k1 + (3/4) k2) = 1.335079087.