古典ルンゲ・クッタ法によるニュートンの運動方程式の数値計算

局所計算精度が４次である古典ルンゲ・クッタ法を用いて、２階常微分方程式であるニュートンの運動方程式を解いて、所定の精度が得られていることを確認します。質量m、バネ定数kの調和振動子（単振動）に対するニュートン運動方程式から得られる位置に対する方程式は次のとおりです。

$\frac{d^2x(t)}{dt^2} = - \frac{k}{m}\, x(t)$

ルンゲ・クッタ法は１階の微分方程式で成り立つ公式なので、上記の２階微分方程式を２つの１階微分方程式に分解して、連立微分方程式の形にします。

$\frac{dx(t)}{dt} = v(t)$
$\frac{dv(t)}{dt} = - \frac{k}{m}\, x(t)$

連立微分方程式と言っても、実質的には独立しているので簡単にルンゲ・クッタ法に適用することができます。次のグラフは $0\leq t \leq 1$ の範囲で、上記調和振動子の時間発展の計算結果（位置と速度）と、計算誤差（解析解との差）の分割数依存性です。

調和振動子（単振動）の時間発展

分割数Nの場合、時間間隔を dt=1/N として計算しています。分割数に対して、計算誤差は４次となることが確認できます。