今回は常微分方程式の解法に広く使われているルンゲ・クッタ法について実習を行う．
前回の実習で行ったEular（オイラー）法による数値積分は，刻み幅を非常に小さく設定しないと精度が悪い問題があった．
今日は，比較的大きな刻み幅（即ち，少ない計算回数）でも精度のよいルンゲ・クッタ法について学ぶ．
関数 $x=x(t)$ にたいして，微分方程式 $dx/dt = f(t, x)$ が定義されているとする．
$t$ の刻み幅を $h$ とし，
現在の時刻を $t_i$，次の時刻を $t_{i+1}= t_i+h$ とする
ルンゲクッタ法では，微小時間 $h$ 経過後の関数の値 $x_{i+1}$ を次のように算出する．
$x_{i+1}=x_i+k$， $k=(h/6) (k_0+2k_1+2k_2+k_3)$
ただし，係数 $k_0$ から $k_3$ はそれぞれ，
$k_0= f(t_i , x_i )$
$k_1= f(t_i+h/2, x_i+k_0 h/2)$
$k_2= f(t_i+h/2, x_i+k_1 h/2)$
$k_3= f(t_i+h , x_i+k_2 h )$
とする．
複雑な計算をしているように見えるが，$t$ や $x$ の異なる位置での微係数を計算し，重みを付けて次の時刻の関数値を計算しているだけである．
計算手順をまとめると，
となり，これをを繰り返すことによって任意の時刻まで値を順次求める（＝数値積分する）ことができる． 参考までに，以下にソースファイルを示す．適宜改変して使用して良い．
ルンゲ・クッタ法→rktest.cpp
ルンゲ・クッタ法以外にも，ルンゲ・クッタ・ギル法，ミルンの方法，アダムス法などがあるので，興味がある人は，各自で調べても良い．
機械工学で利用する微分方程式は，二階，もしくはそれ以上の高階の方程式が頻繁に登場する．
ルンゲ・クッタ法は基本的に一階の微分方程式を解くための方法だが，少し工夫すれば，二階微分方程式を解くことができる．
具体的にどうするかと言うと，一階微分方程式を連立して解く方法を利用する．
いま，$dx/dt = z$ とおいて，導関数を新たに変数とすれば，時間微分を「'」で表して，$x'= z$ となる．
これをもう一回微分すると， $x '' = z'$ が得られる．
新たに変数が発生したため，
\begin{align*}
\frac{dx}{dt} &= f(t, x, z)\\
\frac{d^2x}{dt^2} &= \frac{dz}{dt} = g(t, x, z)
\end{align*}
の連立方程式を立てて同時に解けば良い．
具体的には，
\begin{align}
t_{i+1}&=t_i+h \notag\\
x_{i+1}&=x_i+k, \ \ \ \ k = (h/6) (k_0+2k_1+2k_2+k_3)\notag\\
z_{i+1}&=z_i+l, \ \ \ \ l = (h/6) (l_0+2l_1+2l_2+l_3)
\label{eqn:runge1}
\end{align}
但し
\begin{align}
k_0 &= f (t_i, x_i, z_i)\notag\\
l_0 &= g (t_i, x_i, z_i)\notag\\
k_1 &= f (t_i + h/2, x_i + k_0 h/2, z_i + l_0 h/2)\notag\\
l_1 &= g (t_i + h/2, x_i + k_0 h/2, z_i + l_0 h/2)\notag\\
k_2 &= f (t_i + h/2, x_i + k_1 h/2, z_i + l_1 h/2)\notag\\
l_2 &= g (t_i + h/2, x_i + k_1 h/2, z_i + l_1 h/2)\notag\\
k_3 &= f (t_i + h , x_i + k_2 h , z_i + l_2 h )\notag\\
l_3 &= g (t_i + h , x_i + k_2 h , z_i + l_2 h )
\label{eqn:runge2}
\end{align}
\eqref{eqn:runge2}の計算結果を\eqref{eqn:runge1}に代入して，計算をすれば良い．
ルンゲ・クッタ法は Eular 法に比べて計算式が増えており，一見すると煩雑に見えるが，基本的な処理としてある点での関数値に勾配（＝微係数）を足して，次の関数値を計算することに変わりは無い．