ルンゲ＝クッタ法 - trinoの走り書き

微分方程式を数値的に解くための方法です。実用的には４次のルンゲ＝クッタ法を用いるのが普通らしいです。まあ、詳しく書いてあるサイトはたくさんありますがなんとなく取り上げてみました(^^;

１階微分方程式
$\frac{dx(x)}{dx} = f(x,y(x))$
について解く場合、
$\Delta y_1 = \Delta x f(x_n,y_n)$
$\Delta y_2 = \Delta x f(x_n + \frac{\Delta x}{2},y_n + \frac{\Delta y_1}{2})$
$\Delta y_3 = \Delta x f(x_n + \frac{\Delta x}{2},y_n + \frac{\Delta y_2}{2})$
$\Delta y_4 = \Delta x f(x_n + \Delta,y_n + \Delta y_3)$
$y_{n+1} = y_n + \frac{1}{6}(\Delta y_1 + 2 \Delta y_2 + 2 \Delta y_3 + \Delta y_4 )$
という差分方程式になる。これでΔxを小さく取れば小さく取るほど精度は上がっていく。導出に関しては・・・ぐぐる様にでも聞いてください（ぉ

具体的にどういった場面で使うかというと、やはり一番は力学的シュミレーションであろう。しかし、力学で扱う基本方程式、いわゆるニュートンの運動方程式は２階の微分方程式であり、また、一般にはベクトルである。
$m\vec{a}=\vec{F}$
これは微分方程式の形にすると
$\frac{d^2\vec{r}(t)}{dt^2}=\vec{f}(t,\vec{r},\frac{d\vec{r}}{dt})$
この場合はどうしたらよいか？まず、ベクトルへの拡張は簡単である。ただ単にynやΔyをベクトルにすればよいだけの話である。次に肝心な２階にする方法。これは２元連立にして１階に落とせばよい。
$\vec{v}=\frac{d\vec{r}}{dt}$
と置いて、
$m\frac{d\vec{v}(t)}{dt}=\vec{f}(t,\vec{r},\vec{v})$
とすればよい。こうすれば、ルンゲ＝クッタ法は
$\Delta \vec{r}_1 = \Delta t \vec{v}_n$
$\Delta \vec{v}_1 = \Delta t \vec{f}(t_n,\vec{r}_n,\vec{v}_n)$
$\Delta \vec{x}_2 = \Delta t ( \vec{v}_n + \frac{\Delta\vec{v}_1}{2} )$
$\Delta \vec{v}_2 = \Delta t \vec{f}(t_n + \frac{\Delta t}{2},\vec{r}_n + \frac{\Delta\vec{r}_1}{2},\vec{v}_n + \frac{\Delta\vec{v}_1}{2})$
$\Delta \vec{x}_3 = \Delta t ( \vec{v}_n + \frac{\Delta\vec{v}_2}{2} )$
$\Delta \vec{v}_3 = \Delta t \vec{f}(t_n + \frac{\Delta t}{2},\vec{r}_n + \frac{\Delta\vec{r}_2}{2},\vec{v}_n + \frac{\Delta\vec{v}_2}{2})$
$\Delta \vec{x}_4 = \Delta t ( \vec{v}_n + \Delta\vec{v}_3 )$
$\Delta \vec{v}_4 = \Delta t \vec{f}(t_n + \Delta t,\vec{r}_n + \Delta\vec{r}_3,\vec{v}_n + \Delta\vec{v}_3)$
$\vec{r}_{n+1} = \vec{r}_n + \frac{1}{6}(\Delta \vec{r}_1 + 2 \Delta \vec{r}_2 + 2 \Delta \vec{r}_3 + \Delta \vec{r}_4 )$
$\vec{v}_{n+1} = \vec{v}_n + \frac{1}{6}(\Delta \vec{v}_1 + 2 \Delta \vec{v}_2 + 2 \Delta \vec{v}_3 + \Delta y_4 )$
となる。

さてはて、今まではいわゆる１粒子系での話。多粒子系においてはどうなるか？今度は外場だけでなく他の粒子からの位置や速度によって加速度が変わってくるしそいつらも運動している。こういう場合はどうしたらいいか？
$\frac{d^2\vec{r^{(n)}}(t)}{dt^2}=\vec{f^{(n)}}(t,\vec{r},\frac{d\vec{r}}{dt}) + \vec{g^{(n)}}(t,\vec{r^{(1)}},\vec{r^{(2)}},\cdots,\vec{v^{(1)}}\vec{v^{(2)}},\cdots )$
と、後半の他の粒子から来る項をどうするかという問題である。

・・・まあ実はその場合も簡単で、全粒子の位置や速度を一つのベクトルに表してしまえばOKである。どういうことかというと、３次元空間上で粒子がｎ個あったとしたら、自由度は全部で３×ｎであり、それを一つの空間に押し込める。これにより、ｎ個の粒子は３ｎ次元空間上の１個の粒子として扱える。いわゆる物理で言うところの位相空間である(正確にはちょい違うけど)。そうすると次元こそ違えど式の形は全く同じになり、先ほどの後半の項ｇは消える(ていうかｆに繰り込まれる)。そうすると計算方法も全く同じで迷うこともないだろう。

まあしかし、ゲームとかで使う分にはオイラー法で大抵は十分ですけど(^^;
$\Delta \vec{r}_{n+1} = \vec{r}_n + \Delta t \vec{v}_n$
$\Delta \vec{v}_{n+1} = \vec{v}_n + \Delta t \vec{f}(t_n,\vec{r}_n,\vec{v}_n)$