ルンゲ＝クッタ法①～ルンゲ＝クッタ法を使うと物体の運動をコンピュータ上でシミュレーションすることができます

ルンゲ＝クッタ法と微分方程式

ルンゲ＝クッタ法は微分方程式の初期値問題に対して近似解を与える一連の方法です。

さて、物体の運動はニュートン方程式 $F=ma$ によって近似的に表されます。ここでFは力、mは質量、aは加速度であり、質量mの物体に力が働いた際にどのような加速度を受けるのか、反対に加速度があるときにどのような力が働いているのか、といった「力と加速度の関係」を表す方程式です。

そして、加速度aとは位置xの時間tによる２階微分（ $\frac{dx^2}{d^2t}$ ）ですが、このように、微分された関数による方程式を一般に微分方程式と呼んでいます。

ということは、ルンゲ＝クッタ法を使えばニュートン方程式の近似解を求めることができるのですから、物体の運動をコンピュータ上でシミュレーションすることができます。これは工作をする上で何かに役立ちそうです、というか本当のところを言うと、地震が起こったときに本棚がどのように動くのかをシミュレーションしたいがためにルンゲ＝クッタ法を調べたのでした。

微分方程式の解を数値的に計算する方法はたくさんありますが、ルンゲ＝クッタ法は計算公式がシンプルでプログラムに実装しやすく精度もよいため、一般的に使われているとのことです。

ルンゲ＝クッタ法の計算公式*1

ここでは微分方程式にyの関数やyの時間微分しか出てこない場合に、ルンゲ＝クッタ法の計算公式を書いてみます。要は「一変数の場合」ということですね～。

実際にシミュレーションするときには、縦×横の２次元空間や、縦×横×高さの３次元空間の中の物体の運動をシミュレーションすることが多そうなので、これだけでは実用に足らないかもしれません、多変数の場合は次回にします。

さて、初期条件と時刻 $t$ での位置 $y$ の時間変化が $f(t,y)$ で与えらている微分方程式を考えます。数式で書くと
$\frac{dy}{dt}=f(t,y),y(t_0)=y_0$
ですねー。

ルンゲ＝クッタ法の計算公式では時間の刻み幅を $h$ として、ある時刻の位置 $y_n$ から、次の時刻の位置 $y_{n+1}$ を近似的に計算するのですが、それは
$y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$
と計算されます。

ここで、 $k_1$ ～ $k_4$ は時刻 $t_n$ とそのときの位置 $y_n$ から $f$ によって計算される値であり、
$k_1=f(t_n,y_n)$
$k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)$
$k_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)$
$k_4=f(t_n+h,y_n+k_3)$
です。