微分方程式を解くための4次のルンゲ＝クッタ法

目次

はじめに

2階微分方程式を解くために、[1]で紹介された4次のルンゲ＝クッタ法に基づく数値法と、解析的手法を使用した計算機が比較および教育目的で紹介されています。
最初に1階の微分方程式に対する4次のルンゲ＝クッタ法を紹介し、次に同じ方法を2階線形微分方程式に適用します。最後に、ルンゲ＝クッタ法で使用するステップサイズ \( h \) の影響を調べるインタラクティブなチュートリアルも用意されています。

1階微分方程式に対する4次のルンゲ＝クッタ法

4次のルンゲ＝クッタ法は、以下の形式の常微分方程式（ODE）を解くために広く使用される数値法です： \[ \dfrac{dy}{dt} = f(y, t) \] 初期条件は \( y(0) = y_0 \) です。 4次のルンゲ＝クッタ法は、将来の時刻 \( t + h \) での解 \( y(t) \) を、時刻 \( t \) での情報に基づいて近似します。この方法は、区間 \([t, t+h]\) 内の異なる点での傾きを4つの中間ステップを使って計算し、それらを組み合わせて加重平均を求めます。その傾きを使って \( y(t+h) \) を推定します。 以下は4次のルンゲ＝クッタ法の手順です：
ステップ1: \[ k_1 = h \cdot f(y_n, t_n) \] ここで、 \( y_n \) と \( t_n \) は \( y \) と \( t \) の現在の値です。
ステップ2: \[ k_2 = h \cdot f\left(y_n + \dfrac{k_1}{2}, t_n + \dfrac{h}{2}\right) \] ステップ3: \[ k_3 = h \cdot f\left(y_n + \dfrac{k_2}{2}, t_n + \dfrac{h}{2}\right) \] ステップ4: \[ k_4 = h \cdot f(y_n + k_3, t_n + h) \] ステップ5:
\( y \) と \( t \) の更新 上記の傾きを組み合わせて、区間 \([t, t+h]\) の終点での \( y \) の値を近似します。これらの傾きの加重平均は以下の式で計算されます： \[ \dfrac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \] この平均傾きを \( h \) で掛けて、 \( y_n \) に加えることで、\( t_{n+1} \) での \( y \) の値を更新します。 \[ y_{n+1} = y_n + h \left( \dfrac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \right) \] \[ t_{n+1} = t_n + h \]

2階微分方程式に対するルンゲ＝クッタ法

次に、以下の形式の2階微分方程式を考えます： \[ \dfrac{d^2y}{dt^2} = f(y, \dfrac{dy}{dt}, t) \] 初期条件： \[ y(0) = y_0 \] \[ \dfrac{dy}{dt}(0) = \dfrac{dy}{dt}_0 \] 新しい変数を導入して、これを2つの1階微分方程式の系に変換します。
次のようにします： \[ \dfrac{dy}{dt} = v \] したがって、2階の微分方程式は以下のように書けます： \[ \dfrac{dy}{dt} = v \] \[ \dfrac{dv}{dt} = \dfrac{d^2y}{dt^2} = f(y, v, t) \] これで、4次のルンゲ＝クッタ法を使用して、この1階微分方程式の系を解くことができます。
ステップ1:
\( k_{1y} \) と \( k_{1v} \) を計算 \[ k_{1y} = h \cdot v \] \[ k_{1v} = h \cdot f(y_n, v_n, t_n) \]
ステップ2:
\( k_{2y} \) と \( k_{2v} \) を計算 \[ k_{2y} = h \cdot (v_n + \dfrac{k_{1v}}{2}) \] \[ k_{2v} = h \cdot f\left(y_n + \dfrac{k_{1y}}{2}, v_n + \dfrac{k_{1v}}{2}, t_n + \dfrac{h}{2}\right) \]
ステップ3:
\( k_{3y} \) と \( k_{3v} \) を計算 \[ k_{3y} = h \cdot (v_n + \dfrac{k_{2v}}{2}) \] \[ k_{3v} = h \cdot f\left(y_n + \dfrac{k_{2y}}{2}, v_n + \dfrac{k_{2v}}{2}, t_n + \dfrac{h}{2}\right) \]
ステップ4:
\( k_{4y} \) と \( k_{4v} \) を計算 \[ k_{4y} = h \cdot (v_n + k_{3v}) \] \[ k_{4v} = h \cdot f(y_n + k_{3y}, v_n + k_{3v}, t_n + h) \]
ステップ5:
\( y \) と \( v \) を更新 \[ y_{n+1} = y_n + \dfrac{1}{6}(k_{1y} + 2k_{2y} + 2k_{3y} + k_{4y}) \] \[ v_{n+1} = v_n + \dfrac{1}{6}(k_{1v} + 2k_{2v} + 2k_{3v} + k_{4v}) \] \[ t_{n+1} = t_n + h \]

ルンゲ＝クッタ法計算機

次に、以下の形式の2階微分方程式の解を近似するための計算機を紹介します。 \[ a \dfrac{d^2y}{dt^2} + b \dfrac{dy}{dt} + c y = 0 \] 4次のルンゲ＝クッタ法を使用して、数値的に解を求め、さらに解析的な方法で得られた解と比較します。
上記の微分方程式は解析的にも解かれ、比較が行われます。
Let \[ \dfrac{dy}{dt} = v \] したがって、以下の系に4次のルンゲ＝クッタ法を適用する必要があります。 \[ \dfrac{dy}{dt} = v \] \[ \dfrac{dv}{dt} = \dfrac{d^2y}{dt^2} = \dfrac{1}{a} (- b \; v - c \; y ) \] この計算機は4次のルンゲ＝クッタ法を使用して同じ微分方程式を解析的に解き、結果をグラフで比較します。

係数と初期条件を入力してください

係数 \(a, b \) と \( c \)、および初期条件 \( y(0) \) と \( dy/dt(0) \)、ステップサイズ \( h \) を入力し、数値解法と解析解のグラフを比較してください。
注：
1) グラフの上にマウスを置くと、グラフ上の点の読み取りができますので、解の比較が可能です。
グラフ上の点の読み取り
2) マウスを使って、グラフの凡例をクリックすることで、1つまたは両方のグラフを表示できます。
グラフの凡例
3) ルンゲ＝クッタ法の誤差はステップサイズ \( h \) に依存します。この計算機では \( h \) を変更できるため、数値法と解析法を比較し、4次ルンゲ＝クッタ法の近似誤差に対する \( h \) の影響を調べることができます。
4) 計算に使用するステップ数 \( n \) も変更できます。グラフは \( x = 0 \) から \( x = n \times h \) までの範囲を表示します。








解：

インタラクティブなチュートリアル：ルンゲ＝クッタ法による近似に対するステップサイズ \( h \) の影響

定数 \( b, c, d \) の値として以下を入力してください：
\( a = 1 \), \( b = 0.5 \), \( c = 5 \) は減衰する振幅を持つ振動関数を与えるはずです。
次に、 \( h = 0.01 \) から始め、 \( h = 1 \) まで増やしてみてください。どの \( h \) の値が解析的手法と比較してより良い近似を与えるかを確認してください。

参考文献およびリンク

1 - https://personal.math.ubc.ca/~cbm/aands/abramowitz_and_stegun.pdf - 数学関数ハンドブック - A. Abramowitz と I. Stegun - ページ 875
2 - University Calculus - Early Transcendental - Joel Hass, Maurice D. Weir, George B. Thomas, Jr., Christopher Heil - ISBN-13 : 978-0134995540
3 - Calculus - Gilbert Strang - MIT - ISBN-13 : 978-0961408824
4 - Calculus - Early Transcendental - James Stewart - ISBN-13: 978-0-495-01166-8