ルンゲ=クッタ法

出典: フリー百科事典『ウィキペディア(Wikipedia)』
ルンゲクッタ法から転送)
移動: 案内, 検索

ルンゲ=クッタ法(Runge-Kutta method)とは、数値解析において常微分方程式の近似解を求める一連の方法である。この技法は1900年頃に数学者 C. RungeM.W. Kutta によって発展された。

目次

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

一般に用いられているルンゲ=クッタ法は4次のルンゲ=クッタ法(RK4)と呼ばれるものである。

初期値問題を次のように設定する。

 y' = f(t, y), \quad y(t_0) = y_0

次に、RK4ではこの問題に対して次式を与える。

 y_{n+1} = y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4)

ここで、

 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)

である。すなわち、次の値(yn+1)は、現在の値(yn)に間隔(h)と推定された勾配の積を加えたものである。その勾配は次の4つの勾配の重み付け平均で求める。

  • k1は初期値における勾配である。
  • k2は区間の中央における勾配であり、勾配k1を用いてtn + h/2におけるyの値をオイラー法により決定したものである。
  • k3は区間の中央における勾配を再計算したものであり、k2の値から決められたyの値を用いる。
  • k4は区間の最後における勾配であり、k3の値から決められたyの値を用いる。

これら4つの平均を取るには、中央の勾配に対して大きな重み付けを与える。

\mbox{slope} = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}.

RK4は4次の方法である。すなわち、厳密解とRK4のテイラー展開が4次の項まで一致する(したがって全体の推定誤差はh5オーダーになる)。

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

陽的ルンゲ=クッタ法は前述したRK4の一般化である。次式で与えられる。

 y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i,

ここで、

 k_1 = f(t_n, y_n), \,
 k_2 = f(t_n+c_2h, y_n+a_{21}hk_1), \,
 k_3 = f(t_n+c_3h, y_n+a_{31}hk_1+a_{32}hk_2), \,
 \vdots
 k_s = f(t_n+c_sh, y_n+a_{s1}hk_1+a_{s2}hk_2+\cdots+a_{s,s-1}hk_{s-1}).

(注:これらの式は差分を含んでいるが、異なる書式による等価の定義である。)

個々の問題としては、整数s(段数)と、係数aij (for 1 ≤ j < is), bi (for i = 1, 2, ..., s) および ci (for i = 2, 3, ..., s)を与える必要がある。それらのデータはRunge-Kutta tableauとして、記憶を助ける配置として知られている。:

0
c2 a21
c3 a31 a32
 \vdots  \vdots  \ddots
cs as1 as2  \cdots as,s − 1
b1 b2  \cdots bs − 1 bs

ルンゲ=クッタ法は次のように一貫したものである、

\sum_{j=1}^{i-1} a_{ij} = c_i\ \mathrm{for}\ i=2, \ldots, s.

私達があるオーダーpを持つ方法を必要とする場合、切断誤差が O(hp+1)となる付随する条件がある。切断誤差自体の定義からこれらを得ることができる。例えば、2次の方法はb1 + b2 = 1, b2c2 = 1/2, and b2a21 = 1/2であれば2次のオーダーを持つ。

[編集]

RK4をこの枠組みに当てはめると、その tableauは:

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

ところで、最も単純なルンゲ=クッタ法は(前方)オイラー法であり、公式yn + 1 = yn + hf(tn,yn)で与えられる。これは1段の陽的ルンゲ=クッタ法から構成される。対応するtableauは:

0
1

2次の方法による例は中間値法により与えられる。

 y_{n+1} = y_n + hf\left(t_n+\frac{h}{2},y_n+\frac{h}{2}f(t_n, y_n)\right).

対応するtableauは:

0
1/2 1/2
0 1


[編集] 誤差と収束性

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

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

[編集] 用法

2段の陽的ルンゲ=クッタ法の例を示す。

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

次の初期値問題を解くものとする。

 y' = (\tan{y})+1,\quad y(1)=1,\ t\in [1, 1.1]

ステップの値をh=0.025.とする。

tableauは、方法の定義のもとで同等な対応をする方程式をもたらす。:

 u_1 = y_n \,
 u_2 = y_n + 2/3hf(t_n, u_1) \,
 y_{n+1} = y_n + h(1/4f(t_n,u_1)+3/4f(t_n+2/3h,u_2))\,
t0 = 1
y0 = 1
t1 = 1.025
u1 = y0 = 1 f(t0,u1) = 2.557407725 u2 = y0 + 2 / 3hf(t0,u1) = 1.042623462
y1 = y0 + h(1 / 4f(t0,u1) + 3 / 4f(t0 + 2 / 3h,u2) = 1.066869388
t2 = 1.05
u1 = y1 = 1.066869388 f(t1,u1) = 2.813524695 u2 = y1 + 2 / 3hf(t1,u1) = 1.113761467
y2 = y1 + h(1 / 4f(t1,u1) + 3 / 4f(t1 + 2 / 3h,u2) = 1.141332181
t3 = 1.075
u1 = y2 = 1.141332181 f(t2,u1) = 3.183536647 u2 = y2 + 2 / 3hf(t2,u1) = 1.194391125
y3 = y2 + h(1 / 4f(t2,u1) + 3 / 4f(t2 + 2 / 3h,u2) = 1.227417567
t4 = 1.1
u1 = y3 = 1.227417567 f(t3,u3) = 3.796866512 u2 = y3 + 2 / 3hf(t3,u1) = 1.290698676
y4 = y3 + h(1 / 4f(t3,u1) + 3 / 4f(t3 + 2 / 3h,u2) = 1.335079087

数値解は下線を付した値に対応する。yiの再計算を避けるためにf(ti,u1)を計算する。

個人用ツール
名前空間
変種
操作
案内
ヘルプ
ツールボックス
他の言語