ルンゲ・クッタの方法（2次）で運動方程式を解きたい

#数学 #Qiita

https://qiita.com/hironarin/items/cc4e1c8f961155a8e511

# ルンゲクッタの方法（2次）

$$TeXmacro \def\DIS{\displaystyle} \def\fraction(#1/#2){\frac{#1}{#2}} \def\ff(#1/#2){\frac{#1}{#2}} \def\bm{\boldsymbol} \def\dd(#1/#2){\frac{d#1}{d#2}} \def\pd(#1/#2){\frac{\partial #1}{\partial #2}} \def\gl{{}^{\mathrm g}} \def\lo{{}^{\mathrm l}}$$

1回積分をして微分方程式を数値的に解く。

微分方程式， $dx/dt=f(t)$， $x⁡(t0)=x_0$を積分する方法を記述します。 $x$は $t$の関数で，結論は

$$x (t_0 + \varDelta t) = x_0 + f(t_0 + \ff(1/2) \varDelta t) \cdot \varDelta t$$ 

となります。ここで得られた $x(t_0 + \Delta t)$の値を $x_0$へと逐次当てはめることで，数値的な積分値が得られます。

1次までのテーラー展開と比較すると， $f(t_0) + f(t_0)\Delta t $ → $f(t_0)+f(t_0+\ff(1/2)\Delta t)\Delta t$と置き換わっています。

導出

関数 $x⁡(t)$ を $t_0$ まわりでテーラー展開します。

$$\begin{eqnarray} x(t_0 + \varDelta t) &=& x(t_0) + \frac{dx}{dt} \varDelta t + \ff(1/2) \frac{d^2x}{dt^2} \varDelta t^2 \cdots \frac{1}{n!}\frac{d^nx}{dt^n} \varDelta t^n + \cdots \end{eqnarray}$$

$\Delta t$の1次の係数に注目すると

$$\begin{eqnarray} x(t_0 + \varDelta t) &=& x(t_0) + f(t_0) \varDelta t + \ff(1/2) \frac{d^2x}{dt^2} \varDelta t^2 + \cdots + \frac{1}{n!}\frac{d^nx}{dt^n} \varDelta t^n + \cdots \end{eqnarray}$$
(1)

ここで，方程式右辺の関数 $f$ を用いて次の値を $\phi $ を考えます。

$$\phi(\varDelta t) := w_1 f(t_0) + w_2 f(t_0 + a\varDelta t)$$

ここで，$w_1$, $w_2$, $a$は定数です。次のステップ $x(t_0 +\Delta t)$ の値を $x(t_0) + \phi(\varDelta t) \varDelta t$ と近似します。つまり，この2つの値の差が小さくなるように定数 $w_1$, $w_2$, $a$ を定めていきます。

$\phi(\Delta t)$ を展開します。

$$\begin{eqnarray} \varDelta t \phi &=& \varDelta t \cdot \left( w_1 f(t_0) + w_2 f(t_0 + a\varDelta t) \right) \\ &=& \varDelta t \cdot \left( w_1 f(t_0) + w_2 f(t_0) + w_2 \fraction(d f(t_0)/dt) a\varDelta t + O(\varDelta t^2) \right) \\ &=& (w_1 + w_2)f(t_0) \varDelta t + aw_2 \fraction(d f(t_0)/dt) \varDelta t ^2 + O(\varDelta t ^3) \end{eqnarray}$$

となります。(1) の$\Delta t$の1次以降の項と比較すると，まず1次の項から$w_1+w_2=1$という条件が付くことが分かります。次に$\Delta t ^2$の係数を比較すると

$$aw_2 = \ff(1/2)$$

という関係を得ます。ここで$w_1=0$, $w_2=1$, $\DIS a=\ff(1/2)$とすることで目的の式を得ることができました。

関数fがx, tを変数に持つ場合

微分方程式に含まれている関数 $f$ が $t$, $x$ の2変数関数で，$\DIS \frac{dx}{dt} = f(t,x), x(t_0) = x_0$の場合に拡張して考えます。

物理的には，エネルギーが保存される系での微分方程式がこの形に相当します。1回だけの積分で事足りるので便利！まず，結論は

$$x (t_0 + \Delta t) = x_0 + f(t_0 + \ff(1/2) \varDelta t , x_0 + f(t_0,x_0) \ff(1/2) \varDelta t) \Delta t$$

となります。

導出

関数 $x(t)$ を $t_0$ まわりでテーラー展開します。

$$\begin{eqnarray} x(t_0 + \varDelta t) &=& x(t_0) + \frac{dx}{dt} \varDelta t + \ff(1/2) \frac{d^2x}{dt^2} \varDelta t^2 \cdots \frac{1}{n!}\frac{d^nx}{dt^n} \varDelta t^n + \cdots \end{eqnarray}$$

$\Delta t$の1次の係数は$f(t_0,x_0)$です。$\Delta t$の2次の係数に注目すると，

$$\frac{d^2x}{dt^2} = \ff({df}/{dt}) = \ff(\partial f/\partial t) + \ff(\partial f/\partial x) \ff(dx/dt) = \ff(\partial f/\partial t) + \ff(\partial f/\partial x) f$$

です。

$$\begin{eqnarray} x(t_0 + \varDelta t) &=& x(t_0) + f(t_0,x_0) \varDelta t + \ff(1/2) \left[ \ff(\partial f(t_0,x_0)/\partial t) + \ff(\partial f(t_0,x_0)/\partial x) f(t_0,x_0) \right] \varDelta t^2 + \cdots \end{eqnarray}$$

となります。

ここで，次の値 $\phi $ を考えます。

$$\phi(\varDelta t) := w_1 f(t_0) ,x_0+ w_2 f(t_0 + a\varDelta t,x_0 + b\varDelta t)$$

ここで，$w_1$, $w_2$, $a$, $b$は定数です。$\phi(\Delta t)$ を展開します。

$$\begin{eqnarray} \varDelta t \cdot \phi &=& \varDelta t \cdot \left( w_1 f(t_0,x_0) + w_2 f(t_0 + a\varDelta t,x_0+b\Delta t) \right) \\ &=& \varDelta t \cdot \left( w_1 f(t_0) + w_2 f(t_0) + w_2 \fraction(\partial f(t_0,x_0)/\partial t) a\varDelta t + w_2 \fraction(\partial f(t_0,x_0)/\partial x) b\varDelta t + O(\varDelta t^2) \right) \\ &=& (w_1 + w_2)f(t_0) \varDelta t + \left[ aw_2 \fraction(\partial f(t_0,x_0)/\partial t) + bw_2 \fraction(\partial f(t_0,x_0)/\partial x) \right] \varDelta t ^2 + O(\varDelta t ^3) \end{eqnarray}$$

$\Delta t$の1次の係数を比較すると，$w_1 + w_2 = 1$となり，2次の係数を比較すると

$$aw_2 = \ff(1/2), \quad bw_2 = \ff(1/2) f(t_0,x_0)$$

と分かります。$w_1=0$, $w_2=0$とすることで結論の式が得られます。

#### プログラム上の手順

数値的な手順では，一旦$f(t_0,x_0)$の値を計算し，それを再度$f$の値を計算するのに使用します。

$f$を定義

$t:=t_0$ %tの初期値 ふつうは0
$t:=x_0$ %xの初期値

$a:=f(t,x)$
$a:=f(t+0.5\Delta t,x+ 0.5a\Delta t)$

$t :=t+\Delta t$ %次のステップの時刻
$x := x + a \Delta t$

運動方程式を解く

運動方程式を数値的に解くことを記述します。運動方程式は通常，2階の微分方程式だから，2度積分をしなければなりません。まずは1次元の場合で考えていきます，運動方程式

$$\DIS \dd(^2x/t^2) = f(t,x,v)$$

から考えることにしよう。$f$は力に相当する項ですが，位置$x$や速度$dx/dt =v$の関数として取り扱うことにします。例えば，調和振動子は$f(x) = -kx$，空気抵抗やコリオリ力を考慮する場合は速度の関数になります。

積分の方法は，まず微分方程式を

$$\begin{eqnarray} \fraction(dv/dt) &=& f(t,x,v) \cdots (1) \\ \fraction(dx/dt) &=& v \cdots (2) \end{eqnarray}$$

と2つの方程式にして取り扱います。方程式(2)を積分するために，(1)も同時に積分するといった具合です。積分すると，

$$v (t_0 + \Delta t) = v_0 + f(t_0 + \fraction(1/2)\Delta t , x_0 + \fraction(1/2) v_0 \Delta t, v_0 + \fraction(1/2) f(t_0,x_0,v_0) \Delta t) \cdot \Delta t$$

$$x (t_0 + \Delta t) = x_0 + \left( v_0 + \fraction(1/2) f(t_0,x_0,v_0) \Delta t \right) \Delta t$$

となります。実際にプログラムを組むときは，次のようなより一般の微分方程式

$$\begin{eqnarray} \fraction(dv/dt) &=& F(t,x,v) \\ \fraction(dx/dt) &=& G(t,x,v) \end{eqnarray}$$

の$G(t,x,v)=v$という関数の場合とした方が便利です。また，これはハミルトンの正準方程式でも，そのまま使えます。積分をすると

$$v (t_0 + \Delta t) = v_0 + F(t_0 + \fraction(1/2)\Delta t , x_0 + \fraction(1/2) G_0 \Delta t, v_0 + \fraction(1/2) F_0 \Delta t) \cdot \Delta t$$

$$x (t_0 + \Delta t) = x_0 + G(t_0 + \fraction(1/2)\Delta t , x_0 + \fraction(1/2) G_0 \Delta t, v_0 + \fraction(1/2) F_0 \Delta t) \cdot \Delta t$$

ハミルトンの運動方程式（正準方程式）を解く

ニュートンの運動方程式の代わりにハミルトンの正準方程式を解いて運動を求めることができます。一座標と運動量の関数ハミルトニアン$H$を用いて，

$$\dot x = \pd(H/p) , \quad \dot p = -\pd(H/x)$$

と表される。力を表す関数の代わりに，2つの関数

$$\pd(H/p) ,\quad \pd(H/x)$$

が運動を決めます。自由度が2になる代わりに，1階の連立微分方程式を解くことになります。

TeXでは \def  で定義できる変数と， \tikzmath{}  内で使用する変数の2種類があります。通常の計算は \tikzmath{}  で行いますが， \loop  を使用すると， \tikzmath{}  内の変数は値を忘れてしまいます。そこで \def  も合わせて使用します。 \def  で定義する変数を$\gl a$

`\loop`
`\tikzmath`
$\DIS \lo t := \gl t$
$\DIS \lo x := \gl x$

計算をする。
次のステップの$t$, $x$を得る。

`\edef`
$\DIS \gl t := \lo t$
$\DIS \gl x := \lo x$