$ \DeclareMathOperator{\arccosh}{arccosh} \DeclareMathOperator{\arcsinh}{arcsinh} \DeclareMathOperator{\rank}{rank} \DeclareMathOperator{\rot}{rot} \DeclareMathOperator{\grad}{grad} \DeclareMathOperator{\diver}{div} $

Runge-Kutta法による運動方程式の数値解析

運動方程式は時間に関する二階の常微分方程式なので、その解は一般的には求めることができません。
従って、コンピュータで数値計算を行い、解析的に運動方程式の数値解を求めることが重要になります。
しかし数値計算には必ず誤差の問題がつきまとうので、誤差をなるべく少なくするように計算をしないととんでもないことになってしまいます [1]
近似計算のうちでも誤差がかなり少なく、計算量もそこそこなためよく使われる計算方法にRunge-Kutta(ルンゲ-クッタ)というのがあります。 5次のRunge-Kutta法とか6次のRunge-Kutta法とかもあるらしいのですが、ここでは一般的に用いられる4次のRunge-Kutta法について解説していきます。
なお、Java でライブラリを作ってみたのでよかったら使ってみてください。

[1]最下部ののアプレット参照

1一変数ver.

一階の常微分方程式に対する Runge-Kutta 法から出発し、最終的に運動方程式に適用できる形にすることを目指します。
微分方程式、
\begin{equation} \frac{dx}{dt} = f(x, t) \end{equation}
(1)
と、初期値、
\begin{equation} x(t_0) = x_0 \end{equation}
(2)
が与えられている時、$h$ 秒後の $x$ の値を以下で近似します。
\begin{equation} x_{n+1} = x_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \end{equation}
(3)
ただし、
\begin{eqnarray*} \begin{cases} k_1 & = f\left(x_n, t_n\right) \\ k_2 & = f\left(x_n + \frac{h}{2}k_1, t_n + \frac{h}{2}\right) \\ k_3 & = f\left(x_n + \frac{h}{2}k_2, t_n + \frac{h}{2}\right) \\ k_4 & = f\left(x_n + hk_3, t_n + h\right) \end{cases} \end{eqnarray*}
(4)
とします。
これが一般に4次のRunge-Kutta法といわれる近似法です。 ちょっと複雑な感じがしますが、この方法だと誤差が $O(h^4)$ になるとか。

2多変数ver.

次に上の Runge-Kutta 法を多変数に拡張してみます。 …といっても太く(ベクトル表記に)するだけですけどw
変数が、
\begin{equation} \bm{x} = \left(\begin{array}{c} x^1 \\ \vdots \\ x^m \end{array}\right) \end{equation}
(5)
$m$ 個あり、微分方程式が、
\begin{equation} \frac{d\bm{x}}{dt} = \bm{f}(\bm{x}, t) = \left(\begin{array}{c} f^1(\bm{x}, t) \\ \vdots \\ f^m(\bm{x}, t) \end{array}\right) \end{equation}
(6)
が与えられている時、$h$ 秒後の $\bm{x}$ の値を以下で近似します。
\begin{equation} \bm{x}_{n+1} = \bm{x}_n + \frac{h}{6} (\bm{k}_1 + 2\bm{k}_2 + 2\bm{k}_3 + \bm{k}_4) \end{equation}
(7)
ただし、
\begin{eqnarray*} \begin{cases} \bm{k}_1 & = \bm{f}\left(\bm{x}_n, t_n\right) \\ \bm{k}_2 & = \bm{f}\left(\bm{x}_n + \frac{h}{2}\bm{k}_1, t_n + \frac{h}{2}\right) \\ \bm{k}_3 & = \bm{f}\left(\bm{x}_n + \frac{h}{2}\bm{k}_2, t_n + \frac{h}{2}\right) \\ \bm{k}_4 & = \bm{f}\left(\bm{x}_n + h\bm{k}_3, t_n + h\right) \end{cases} \end{eqnarray*}
(8)
とします。
以上です。 手抜きではありません。

3二階の常微分方程式への拡張

さて、いよいよ本命です。 上の多変数バージョンの Runge-Kutta 法を二階の常微分方程式へ対応できるように拡張しましょう。 なお、初めから変数は二個以上あるものだと思ってやっていきます。
…といってもこちらも簡単で、
\begin{equation} \frac{d^2\bm{x}}{dt^2} = \bm{f}\left(\bm{x}, \frac{d\bm{x}}{dt}, t\right) \end{equation}
(9)
という運動方程式を、
\begin{eqnarray*} \frac{d\bm{x}}{dt} & = & \bm{v}(\bm{x}, \bm{v}, t) = \bm{v} \\ \frac{d\bm{v}}{dt} & = & \bm{f}(\bm{x}, \bm{v}, t) \end{eqnarray*}
(10)
の二つに分けるだけです。 あとは、
\begin{equation} \bm{X} = \left(\begin{array}{c} x^1 \\ \vdots \\ x^m \\ v^1 \\ \vdots \\ v^m \end{array}\right) \end{equation}
(11)
として、この $\bm{X}$ を上の多変数バージョンの $\bm{x}$ だと思って計算すれば OK です。
簡単でしょ? …といいたいところですが、やってることは単純だけど変数が大量発生していて、いざプログラムを組もうとしてもなかなか面倒だったりします。
とりあえず Java でライブラリ的なものを作っておいたので、よかったら使ってみてください。

4Java による実装

Java でライブラリのようなものを作ってみました。
「シミュレーションをやりたいけど、Runge-Kutta 法を実装するのは面倒!」という人は使ってみるときっと幸せになれます。

4.1ソースコード

Java ソースコードです。 使い方は一つ下のセクションを参照してください。 Javadoc もあります。
なお、ライセンスは GNU LGPL v3 としておきます。

4.1.1多変数ver.

4.1.2運動方程式ver.

※上の RungeKutta.java, RungeKuttaFunc.java も必要です

4.2サンプルコード

サンプルコードです。適宜 Javadoc も参照してください。

4.2.1多変数ver.

多変数バージョンの Runge-Kutta 法のサンプルです。
class RKTest{
	public static void main(String args[]){
		// 初期値を配列で指定します。
		double x0[] = {1D};
		// 開始時刻です。
		double t0 = 0D;
		// RungeKutta インスタンスを生成します。
		RungeKutta rk = new RungeKutta(
			x0, t0,
			new RungeKuttaFunc(){
				public double[] calc(double x[], double t){
					// ここで dx/dt の値を計算し、返します。
					// このサンプルの場合、 x をそのまま返しているため、
					//     dx/dt = x
					// という微分方程式を表します。
					return x;
				}
			}
		);
		// 計算を行っていきます。
		for(;;){
			// 何秒後の状態を計算するのかを表します。
			double dt = 0.01D;
			// 計算した値を取得します。
			double xCalc = rk.get(0);
			// dx/dt = x の解は、
			//     x = e^t
			// なので、この式から理論計算による値を計算します。
			double xReal = Math.exp(rk.getTime());
			// 左から順に、現在時刻、理論値、計算値、(理論値-計算値) を標準出力に出力します。
			System.out.printf(
				"\rt=%.2f real=%.2f calc=%.2f diff=%.10f",
				rk.getTime(), xReal, xCalc, xReal - xCalc
			);
			// dt 秒後の状態に移行します。
			rk.calcNext(dt);
			// 速すぎて確認できないので 10ms ほど休憩 ZZzz
			try{Thread.sleep(10);}catch(Exception e){}
		}
	}
}

4.2.2運動方程式ver.

多変数バージョンの Runge-Kutta 法のサンプルです。 上のとほぼ同じですが…。
class RKETest{
	public static final double K = 1D;
	public static void main(String args[]){
		// 位置の初期値を配列で指定します。
		double x0[] = {0D};
		// 速度の初期値を配列で指定します。
		double v0[] = {1D};
		// 開始時刻です。
		double t0 = 0D;
		// RungeKuttaEOM インスタンスを生成します。
		RungeKuttaEOM rke = new RungeKuttaEOM(
			x0, v0, t0,
			new RungeKuttaEOMFunc(){
				public double[] calc(double x[], double v[], double t){
					// ここで d^2x/dt^2 の値を計算し、返します。
					// このサンプルの場合、
					//     d^2x/dt^2 = - K * x
					// という微分方程式を表します。
					// 単振動ですね。
					double xx[] = new double[1];
					xx[0] = - K * x[0];
					return xx;
				}
			}
		);
		// 計算を行っていきます。
		for(;;){
			// 何秒後の状態を計算するのかを表します。
			double dt = 0.01D;
			// 現在時刻です。
			double t = rke.getTime();
			// 計算した値を取得します。
			double xCalc = rke.getX(0);
			// d^2x/dt^2 = - K * x の解は、
			//     x = v0 / sqrt(K) * sin(sqrt(K) * t)
			// なので、この式から理論計算による値を計算します。
			double xReal = (v0[0] / Math.sqrt(K)) * Math.sin(Math.sqrt(K) * t);
			// 左から順に、現在時刻、位置の理論値、位置の計算値、(理論値-計算値) を標準出力に出力します。
			System.out.printf(
				"\rt=%.2f real=%.2f calc=%.2f diff=%.10f",
				t, xReal, xCalc, xReal - xCalc
			);
			// dt 秒後の状態に移行します。
			rke.calcNext(dt);
			// 速すぎて確認できないので 10ms ほど休憩 ZZzz
			try{Thread.sleep(10);}catch(Exception e){}
		}
	}
}

4.3デモ - 振り子

振り子は、振り幅が小さい場合には高校でも習うように、その運動は単振動に近似できますが、 振れ幅が大きい場合には近似が使えないため運動方程式を解くことができません[1]
そこで振り子の運動を Runge-Kutta 法を用いてシミュレートしてみます。
なお、同時に Euler 法という近似方法を使って計算したものも右に載せておきます。 こちらは近似誤差が大きいため、すぐに振れ幅が大きくなって、しまいにはぐるぐる回転しはじめてしまいます。
なお Euler 法とは、
\begin{equation} \bm{x}_{n+1} = \bm{x}_n + h \frac{d\bm{x}}{dt}(\bm{x}_n, t_n) \end{equation}
(12)
という、安直な近似法であり、誤差は $O(h^2)$ です[2]。 テイラー展開して $h^2$ 以降の項を消し去っただけですね。
※クリックすると始めからシミュレーションを行います

[1]よく知らないけど、実は楕円函数とかいうのを使えば解けるとか何とか。でも結局計算は無理…なんじゃないかなぁ
[2]Runge-Kutta法では $O(h^4)$ でしたね