1. この記事について

1-1. 対象読者

物理シミュレーションに興味がある方
「物理シミュレーションを、オイラー法で実装すると、誤差が大きくなる」と言われてピンとくる方
ルンゲクッタ法を詳しく知らない方
python3 の利用者

1-2. この記事で分かること

ルンゲクッタ法はなぜ必要か
ルンゲクッタ法の成り立ち(軽く)
4次ルンゲクッタ法を利用した物理シミュレーションの例

次のことは、この記事では行わない。

ルンゲクッタ法の詳しい導出(気が向いたらやるかも)
安定性解析、エネルギー保存性の検証などの評価
ライブラリの紹介 (コメントいただきました。単にルンゲクッタ法を使うだけであれば、scipy.integrate.solve_ivp(関数, tの区間, 初期値, method="RK45") とするだけで簡単に利用できるようです。関数 は t と np.array([x,y,v_x,v_y])を受け取り、np.array([v_x,v_y,a_x,a_y])を返す関数だと思います。詳しい使い方は公式ドキュメントを読んでください。)

あくまで、物理シミュレーションの初心者が、オイラー法を卒業してルンゲクッタ法にたどり着くまでの一助としてこの記事を捉えてほしい。

1-3. 今回考える物理シミュレーションについて

今回考える物理シミュレーションでは、1つの点 がする 任意の運動 について再現することを考える。
任意の運動を再現するために、 加速度をコールバック関数 とし、ユーザが点の位置と速度をもとに直後のタイムステップにおける加速度を自由に決められるようにする。
$$\overrightarrow{a(t)}=\overrightarrow{{\rm callback}\left(t,\overrightarrow{r(t)},\overrightarrow{v(t)}\right)}
$$
($\overrightarrow{a(t)}$は時刻$t$における加速度、$\overrightarrow{v(t)}$は速度、$\overrightarrow{r(t)}$は位置)

python コードは次のとおり。(「基本コード」と呼ぶことにする)

基本コード

class シミュレータ型:
    def __init__(self, 初期位置, 初速度, 時刻幅, callback):
        (self.__r_x, self.__r_y) = 初期位置
        (self.__v_x, self.__v_y) = 初速度
        (self.__a_x, self.__a_y) = (0., 0.) # 初期の加速度
        self.__時刻幅 = 時刻幅
        self.__callback = callback
        self.__ステップ = -1
    @property
    def 位置(self): return (self.__r_x,  self.__r_y)
    @位置.setter
    def 位置(self, value):
        (self.__r_x, self.__r_y) = value
    @property
    def 速度(self): return (self.__v_x,  self.__v_y)
    @速度.setter
    def 速度(self, value):
        (self.__v_x,  self.__v_y) = value
    @property
    def 加速度(self): return (self.__a_x,  self.__a_y)
    @加速度.setter
    def 加速度(self, value):
        (self.__a_x,  self.__a_y) = value
    @property
    def 時刻(self): return self.__ステップ * self.__時刻幅

    def 進める(self):
        # 加速度
        self.加速度 = self.__callback(
            self.時刻, self.位置, self.速度
        )

        # 新しい速度と位置
        self.__速度と位置を更新する()

        # 時刻(ステップ)を更新
        self.__ステップ += 1

    def __速度と位置を更新する(self):
        [r_x, r_y, v_x, v_y] = [None]*4 # 未実装
        return ((r_x, r_y), (v_x, v_y))

2. オイラー法の復習

オイラー法とは、「時刻幅 $\Delta t$ の間、常に一定の加速度、速度が与えられている」とみなして、速度や位置を更新する方法だった。

$$\begin{array}{}
\overrightarrow{v(t)} = \overrightarrow{v(t-\Delta t)} + \overrightarrow{a(t-\Delta t)}\Delta t\\
\overrightarrow{r(t)} = \overrightarrow{r(t-\Delta t)} + \overrightarrow{v(t-\Delta t)}\Delta t\\
\end{array}$$

python で実装するなら、基本コードの __速度と位置を更新する メソッドを次の内容で書き換えればよい。

オイラー法

    def __速度と位置を更新する(self):
        x,y = 0,1
        r_before = self.位置
        v_before = self.速度
        a_before = self.加速度
        delta_t = self.__時刻幅

        r_x = r_before[x] + v_before[x] * delta_t
        r_y = r_before[y] + v_before[y] * delta_t
        v_x = v_before[x] + a_before[x] * delta_t
        v_y = v_before[y] + a_before[y] * delta_t

        self.位置 = (r_x, r_y)
        self.速度 = (v_x, v_y)

基本コードの __速度と位置を更新する メソッドを上記の内容で書き換えたものを euler.py として保存する。

今、オイラー法で、次のような条件のシミュレーションを考えてみよう。

初速度: $\begin{pmatrix}10\\0\end{pmatrix}$
加速度: 進行方向に垂直な方向に大きさ$10$
時刻幅: 1, 0.1, 0.01, 0.001 をそれぞれ試す
終了時刻: 300

理論的には、進行方向に垂直な方向に一定の大きさの加速度を加えているため、等速円運動を行うはずである。
その半径は速さ^2 / 加速度の大きさ = $100/10 = 10$ となるはずである。
シミュレーションは、次のようなコード simulation.py で実行できる。(同ディレクトリに euler.py の他、回帰円を求める記事で紹介した fit_circle.py を保存しておく)

simulation.py

from euler import シミュレータ型
from fit_circle import fit_circle
import math
import matplotlib.pyplot as plt
def get_a(t, r, v):
    # 進行方向に垂直な方向に10だけ加速し続ける
    x,y = 0,1
    
    進行方向 = math.atan2(v[y], v[x])
    加速度の方向 = 進行方向 + math.pi / 2
    加速度の大きさ = 10

    a_x = 加速度の大きさ * math.cos(加速度の方向)
    a_y = 加速度の大きさ * math.sin(加速度の方向)
    return (a_x, a_y)

# シミュレーションの設定    
シミュレーション = シミュレータ型(
    初期位置 = ( 0, 0),
    初速度   = (10, 0),
    時刻幅   = 1, # あるいは 0.1, 0.01, 0.001
    callback = get_a
)
終了時刻 = 300

# シミュレーションの実行
ts, xs, ys = [], [], []
x,y = 0,1
while シミュレーション.時刻 <= 終了時刻:
    ts += [シミュレーション.時刻]
    xs += [シミュレーション.位置[x]]
    ys += [シミュレーション.位置[y]]
    シミュレーション.進める()


# 結果のグラフ
# 1. 用意
幅 = max([abs(x) for x in xs] + [abs(y) for y in ys])*2.5
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111)
plt.xlim(-幅/2, 幅/2)
plt.ylim(-幅/2, 幅/2)
plt.xlabel("X")
plt.ylabel("Y")
ax.set_aspect("equal", adjustable="box")

# 2. 軌跡
ax.plot(xs, ys, lw = 1, color = "blue") # lw は線の幅。適当にいじる

# 3. 時刻
for t, x, y in zip(ts, xs, ys):
    if t % 5 != 0: # 時刻表示の間隔。適当にいじる
        continue

    ax.text(
        x, y, str(int(t)), fontsize = 8, # fontsize は適当にいじる
        ha = "center", va = "center", color = "black"
    )

# 4. 回帰円 (t=295,296,297,298,299,300)
idxs = [i for i in range(len(ts)) if ts[i] in [295, 296, 297, 298, 299, 300]]
points = [(xs[i], ys[i]) for i in idxs]
(xe, ye, r) = fit_circle(points)
print(f"回帰円半径: {r}")
xs = []
ys = []
for theta1000 in range(int(2*math.pi*1000)):
    theta = theta1000 / 1000
    xs += [xe + r * math.cos(theta)]
    ys += [ye + r * math.sin(theta)]
ax.plot(xs, ys, lw = 1, color = "#888888", linestyle=":", zorder=-1)
    # lw は適当にいじる

# 5. 表示
ax.grid(True)
plt.show()

これを実行すると、次のようなグラフを得る。

シミュレーションでは、理論(半径10の等速円運動)どおりにならず、渦巻き運動が観測された。
渦巻き運動では、回帰円半径が時間と共に大きくなる。
時刻幅を小さくすればするほど、時刻295,296,297,298,299,300の点による回帰円半径が理論値10に近づいたことも確認できた。

時刻幅を小さくすれば確かに誤差も小さくなるのではあるが、同じ時間をシミュレートするのに必要なステップ数がその分増えてしまうという問題もある。

そこで、より大きな時刻幅でも誤差を小さくできるような方法があれば、それを採用したい。
そのような背景から様々な代替手法が提案されるが、その一つに「ルンゲクッタ法」があるというわけである。

3. ルンゲクッタ法の理論

位置と速度の両方についていきなり同時に議論すると混乱するので、3-1節から 3-6節まででは位置の更新式についてのみ議論する。

3-1. オイラー法の誤差の原因

オイラー法による位置$\overrightarrow{r(t)}$の誤差は、時間幅 $\Delta t$ の間、速度 $\frac{\rm d}{{\rm d}t}\overrightarrow{r(t)}$ が一定であるとみなしていることに起因するものである。
本当は$\overrightarrow{0}$ではないはずの加速度 $\left(\frac{\rm d}{{\rm d}t}\right)^2\overrightarrow{r(t)}$、加加速度$\left(\frac{\rm d}{{\rm d}t}\right)^3\overrightarrow{r(t)}$、加加加速度$\left(\frac{\rm d}{{\rm d}t}\right)^4\overrightarrow{r(t)}$、...がすべて$\overrightarrow{0}$とみなされていることで誤差が生まれている。

3-2. テイラーの定理

関数$f(x)$ が $ a \leq x \leq b $ の範囲で連続で、 $ a < x < b$ の範囲で $n$ 回微分可能であるとき、

$$f(b)= \sum_{k=0}^{n-1}\left(
\frac{
\left\{\left(\frac{\rm d}{{\rm d}x}\right)^k f\right\}(a)
}{k!}
(b-a)^k
\right)
+\frac{
\left\{\left(\frac{\rm d}{{\rm d}x}\right)^n f\right\}(c)
}{n!}
(b-a)^n
$$
を満たす $c$ が、$a < c < b$ の範囲で存在する。このことをテイラーの定理という。

但し、 $\left\{\left(\frac{\rm d}{{\rm d}x}\right)^k f\right\}(a)$は、 $\left(\frac{\rm d}{{\rm d}x}\right)^k f(x)$ に $x=a$ を代入したものを表す。

上記の定理について、$f$を$\overrightarrow{r}$、 $x$ を $t$ に置き換え、 $b=t, a=t-\Delta t, c=t-\delta$ を代入すると、

$$\overrightarrow{r(t)} = \sum_{k=0}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^k r\right\}(t-\Delta t)}
}{k!}
(\Delta t)^k
\right)
+\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}
}{n!}
(\Delta t)^n$$

となる。 ( $t-\Delta t < t-\delta < t$ なので、 $0 < \delta < \Delta t$ )

$n=2$ を代入すると、

$$\overrightarrow{r(t)}=\overrightarrow{r(t-\Delta t)} +
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right) r\right\}(t-\Delta t)} \Delta t +
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^2 r\right\}(t-\delta)}
}{2!}(\Delta t)^2
$$
つまり、

$$\overrightarrow{r(t)}=\overrightarrow{r(t-\Delta t)} +
\overrightarrow{v(t-\Delta t)} \Delta t +
\frac{
\overrightarrow{a(t-\delta)}
}{2!}(\Delta t)^2
$$

となる。

これをオイラー法による位置の更新式
$\overrightarrow{r(t)} = \overrightarrow{r(t-\Delta t)} + \overrightarrow{v(t-\Delta t)}\Delta t$

と比較すると、誤差は

$$\frac{
\overrightarrow{a(t-\delta)}
}{2!}(\Delta t)^2$$

となる。

3-3. ルンゲクッタ法は何をしてくれるのか

前節では、オイラー法により計算される位置の誤差が、$\frac{
\overrightarrow{a(t-\delta)}
}{2!}(\Delta t)^2$であることを、テイラーの定理によって説明した。

逆に言えば

$$\overrightarrow{r(t)} = \sum_{k=0}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^k r\right\}(t-\Delta t)}
}{k!}
(\Delta t)^k
\right)
+\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}
}{n!}
(\Delta t)^n$$

なのだから、位置の更新式を

$$\overrightarrow{r(t)} = \sum_{k=0}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^k r\right\}(t-\Delta t)}
}{k!}
(\Delta t)^k
\right)
$$

…①に出来れば、誤差は $\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}
}{n!} (\Delta t)^n$ に出来る。
この誤差は(基本的に) $n$ が大きくなればなるほど、小さくなる。

説明

$$\frac{
\left|\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{n+1} r\right\}(t-\delta)}\right|
}{
\left|\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}\right|
}
< \frac{n+1}{\Delta t}$$

が成り立つとき、 $n$ が大きくなればなるほど、誤差 $\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}
}{n!} (\Delta t)^n$ の大きさは小さくなる。

例えば $\left|\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^n r\right\}(t-\delta)}\right|=10^n, \Delta t = 0.1$ とすると、
$n=2$なら誤差の大きさ $\frac{10^2}{2!}0.1^2=0.5$
$n=3$なら誤差の大きさ $\frac{10^3}{3!}0.1^3=0.166...$

(説明終わり)

(なるべく大きな$n$で、)①の更新式を使うことで、誤差を小さくしようとするのが $n-1$ 次のルンゲクッタ法である。

3-4. ルンゲクッタ法の課題

加加速度を $\overrightarrow{a^{(2)}(t)}$, 加加加速度を $\overrightarrow{a^{(3)}(t)}$, ... 加×$i$速度を$\overrightarrow{a^{(i)}(t)}$ と書くことにする。

例えば、 $n=5$ として、 4次のルンゲクッタ法を考えると、①の更新式は
$$\overrightarrow{r(t)}=
\overrightarrow{r(t-\Delta t)} +
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right) r\right\}(t-\Delta t)} \Delta t
+
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^2 r\right\}(t-\Delta t)}
}{2!}(\Delta t)^2
+
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^3 r\right\}(t-\Delta t)}
}{3!}(\Delta t)^3
+
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^4 r\right\}(t-\Delta t)}
}{4!}(\Delta t)^4
$$

つまり

$$\overrightarrow{r(t)}=
\overrightarrow{r(t-\Delta t)} +
\overrightarrow{v(t-\Delta t)} \Delta t
+
\frac{
\overrightarrow{a(t-\Delta t)}
}{2!}(\Delta t)^2
+
\frac{
\overrightarrow{a^{(2)}(t-\Delta t)}
}{3!}(\Delta t)^3
+
\frac{
\overrightarrow{a^{(3)}(t-\Delta t)}
}{4!}(\Delta t)^4
$$

…②となり、誤差は$$
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^5 r\right\}(t-\delta)}
}{5!}(\Delta t)^5$$
つまり
$$\frac{
\overrightarrow{a^{(4)}(t-\delta)}
}{5!}(\Delta t)^5$$
となる。

しかし、シミュレーションでは、速度までしか考えていない※ため、加速度、加加速度、加加加速度を含む②の式は、そのままの形では利用できない。

※に対して「加速度は考えているだろ」と突っ込みたい方へ

鋭いご指摘です。ですが、位置の更新式に加速度を利用してはいけない理由があるのです。
今は「位置の更新式」について議論していますが、実際には「速度の更新式」についても議論する必要があります。
速度の更新式は、位置の更新式に対して、
位置→速度
速度→加速度
と置き換えただけのものにしたいです。
位置の更新式に「加速度」を含めてしまうと、速度の更新式では「加加速度」に置き換わることになります。
速度の更新式で加加速度を利用できないので、位置の更新式でも加速度を利用できない、ということになります。

(※に対して「加速度は考えているだろ」と突っ込みたい方へ終わり)

3-5. ルンゲクッタ法の工夫

そこで、ルンゲクッタ法では、 加加速度, 加加加速度, ... を速度の式で置き換える ということをする。
具体的には、次のようにする。

$$\begin{array}{l}
\sum_{k=0}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^k r\right\}(t-\Delta t)}
}{k!}
(\Delta t)^k
\right)
=\overrightarrow{r(t-\Delta t)}+\sum_{i=0}^{n-2}\left(b_i \overrightarrow{v(t-(1- c_i) \Delta t)}\right)\Delta t\\
ただし、c_0=0
\end{array}
$$

上記を成り立たせるパラメータ $b_i, c_i$が何故か※1あるので、それを何とかして※2見つける。

※1, ※2の解説は筆者の実力不足のため、行わない。

左辺の加速度、加加速度、加加加速度... が、右辺では速度で説明されていることに注目されたい。

3-6. ルンゲクッタ法の真の強み

さて、上記の式もルンゲクッタ法の式といって間違いではないのであるが、これではルンゲクッタ法の真の強みを発揮できていない。

上記の式では $\frac{\rm d}{{\rm d}t}\overrightarrow{r(t)}=\overrightarrow{v(t)}$ という前提で成り立っている。
しかし実際には、$\frac{\rm d}{{\rm d}t}\overrightarrow{r(t)}=\overrightarrow{v\left(t, \overrightarrow{r(t)}\right)}$ のように、速度は時刻だけでなく、位置にも依存するとして考える必要がある。

ルンゲクッタ法の真の強みは、速度が時刻と位置の両方に依存する場合にも対応できることである。

この場合、テイラーの定理による式は

$$
\overrightarrow{r(t)}=\overrightarrow{r(t-\Delta t)} +
\sum_{l=1}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{l-1} v\right\}\left(t-\Delta t,\overrightarrow{r(t-\Delta t)}\right)}
}{l!}
(\Delta t)^l
\right)
+\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{n-1} v\right\}\left(t-\delta,\overrightarrow{v(t - \delta)}\right)}
}{n!}
(\Delta t)^n
$$

となる。 (3-2節の式に対し、$\overrightarrow{\left\{
\left(
\frac{\rm d}{{\rm d}t}
\right)^l r
\right\}(t-\Delta t)}
$を、 $l=0$なら$\overrightarrow{r(t-\Delta t)}$に、それ以外なら$\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{l-1} v\right\}\left(t-\Delta t,\overrightarrow{r\left(t-\Delta t\right)}\right)}$に置き換えたものである。)

よって、位置の更新式は

$$
\overrightarrow{r(t)}=\overrightarrow{r(t-\Delta t)} +
\sum_{l=1}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{l-1} v\right\}\left(
t-\Delta t,\overrightarrow{r(t-\Delta t)}
\right)}
}{l!}
(\Delta t)^l
\right)
$$

となる。加速度, 加加速度, ... を速度の式で置き換える方法は、次のとおり。

$$\begin{array}{ll}
\sum_{l=1}^{n-1}\left(
\frac{
\overrightarrow{\left\{\left(\frac{\rm d}{{\rm d}t}\right)^{l-1} v\right\}\left(
t-\Delta t,\overrightarrow{r(t-\Delta t)}
\right)}
}{l!}
(\Delta t)^l
\right)=
\sum_{i=0}^{n-2}b_i\overrightarrow{k_i}\\
\overrightarrow{k_i} = \overrightarrow{v\left(
t-(1-c_i)\Delta t,
\overrightarrow{r(t-\Delta t)} + c_i\overrightarrow{k_{i-1}}\Delta t
\right)}\\
ただし、 c_0=0であり、 \overrightarrow{k_{-1}} は有限の大きさの成分を持つ任意のベクトル
\end{array}$$

上記を満たす $b_i, c_i$ を見つけて代入すればよい。

改めて、$n-1$次ルンゲクッタ法における位置の更新式は次のとおりになる。

$$\begin{array}{ll}
\overrightarrow{r(t)}=\overrightarrow{r(t-\Delta t)} +
\sum_{i=0}^{n-2}b_i\overrightarrow{k_i}\\
\overrightarrow{k_i} = \overrightarrow{v\left(
t-(1-c_i)\Delta t,
\overrightarrow{r(t-\Delta t)} + c_i\overrightarrow{k_{i-1}}\Delta t
\right)}\\
ただし、 c_0=0であり、 \overrightarrow{k_{-1}} は有限の大きさの成分を持つ任意の2次元ベクトル
\end{array}$$
…③

3-7. 位置と速度を同時に更新する

③の式で、ベクトルを行列に置き換えてやることで、位置と速度を同時に更新できる。
具体的には、 $$
\begin{array}{l}
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix} =
\begin{pmatrix}
時刻tでのx座標 & 時刻tでの速度のx成分\\
時刻tでのy座標 & 時刻tでの速度のy成分
\end{pmatrix}\\
\frac{\rm d}{{\rm d}t}\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}={\rm V\_and\_a}\left(
t,\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
\right)
\end{array}$$

として、③の式を次のように置き換えれやればよい。

ベクトル $\overrightarrow{r(t)}$ を行列 $
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
$ に置き換える
ベクトル $\overrightarrow{k_i}$ を行列 $K_i$ に置き換える
ベクトル $\overrightarrow{v\left(t, \overrightarrow{r(t)}\right)}$ を行列$
{\rm V\_and\_a}\left(
t,\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
\right)
$ に置き換える

すると、次のようになる。

$$\begin{array}{ll}
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
=
\begin{pmatrix}
\overrightarrow{r(t-\Delta t)} &
\overrightarrow{v\left(t-\Delta t,\overrightarrow{r(t - \Delta t)}\right)}
\end{pmatrix} +
\sum_{i=0}^{n-2}b_iK_i\\
K_i = {\rm V\_and\_a}\left(
t-(1-c_i)\Delta t,
\begin{pmatrix}
\overrightarrow{r(t-\Delta t)} &
\overrightarrow{v\left(t-\Delta t,\overrightarrow{r(t - \Delta t)}\right)}
\end{pmatrix} + c_iK_{i-1}\Delta t
\right)\\
ただし、 c_0=0であり、 K_{-1} は有限の大きさの成分を持つ任意の2次正方行列
\end{array}$$

…④

さて、④の式はだいぶ見た目がごちゃごちゃしているので、${\rm R\_and\_v}(t)=
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}$ と置いて整理すると、

$$\begin{array}{ll}
{\rm R\_and\_v}(t)={\rm R\_and\_v}(t-\Delta t) +
\sum_{i=0}^{n-2}b_iK_i\\
K_i = {\rm V\_and\_a}\left(
t-(1-c_i)\Delta t,
{\rm R\_and\_v}(t-\Delta t) + c_iK_{i-1}\Delta t
\right)\\
ただし、 c_0=0であり、 K_{-1} は有限の大きさの成分を持つ任意の2次正方行列
\end{array}$$

…⑤とできる。

後は、 ${\rm V\_and\_a}(t, {\rm R\_and\_v}(t))$ と $b_i, c_i$ が具体的にわかれば、ルンゲクッタ法による位置と速度の更新が実現する。
${\rm V\_and\_a}(t, {\rm R\_and\_v}(t))$ については 3-8節、
$b_i,c_i$ の具体例については 3-9節でそれぞれ説明する。

3-8. 位置と速度を同時に微分する

$$
\begin{array}{lll}
{\rm V\_and\_a}(t, {\rm R\_and\_v}(t))
&=& \frac{\rm d}{{\rm d}t}{\rm R\_and\_v}(t)\\
&=& \frac{\rm d}{{\rm d}t}
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
\end{array}$$
なわけであるが、
$$\frac{\rm d}{{\rm d}t}\overrightarrow{r(t)}=\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}$$
であり、また
$$\frac{\rm d}{{\rm d}t}\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
=\overrightarrow{a\left(
t,\overrightarrow{r(t)},\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\right)}$$と書ける。
但し、1-3節で加速度をコールバック関数としたので
$$\frac{\rm d}{{\rm d}t}\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
=\overrightarrow{{\rm callback}\left(t,\overrightarrow{r(t)},\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}\right)}$$

よって、$$
\frac{\rm d}{{\rm d}t}
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}
=
\begin{pmatrix}
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)} &
\overrightarrow{
{\rm callback}\left(
t,\overrightarrow{r(t)},\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\right)
}
\end{pmatrix}
$$
なので
$${\rm V\_and\_a}\left(t, \begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix}\right)
=
\begin{pmatrix}
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)} &
\overrightarrow{
{\rm callback}\left(
t,\overrightarrow{r(t)},\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\right)
}
\end{pmatrix}
$$
である。
ややこしいので$
\begin{pmatrix}
\overrightarrow{r(t)} &
\overrightarrow{v\left(t,\overrightarrow{r(t)}\right)}
\end{pmatrix} =
\begin{pmatrix}
\overrightarrow{r} &
\overrightarrow{v}
\end{pmatrix}
$と簡略化して書くと

$${\rm V\_and\_a}\left(t, \begin{pmatrix}
\overrightarrow{r} &
\overrightarrow{v}
\end{pmatrix}\right)
=
\begin{pmatrix}
\overrightarrow{v} &
\overrightarrow{
{\rm callback}\left(
t,\overrightarrow{r},\overrightarrow{v}
\right)
}
\end{pmatrix}
$$
となる。

3-9. 4次ルンゲクッタ法 (RK4)

$n=5$とした4次のルンゲクッタ法については、次のパラメータがよく知られている。

$b_0=1/6, b_1=1/3, b_2=1/3, b_3=1/6$

$c_0=0, c_1=1/2, c_2=1/2, c_3=1$

このパラメータ使う場合、⑤に示した更新式は次のようになる。

$$\begin{array}{l}
{\rm R\_and\_v}(t)={\rm R\_and\_v}(t-\Delta t) + \left(
\frac{1}{6} K_0 +
\frac{1}{3} K_1 +
\frac{1}{3} K_2 +
\frac{1}{6} K_3
\right)\Delta t\\
K_0 = {\rm V\_and\_a}\left(
t-\Delta t,
{\rm R\_and\_v}(t-\Delta t)
\right)\\
K_1 = {\rm V\_and\_a}\left(
t-\frac{1}{2}\Delta t,
{\rm R\_and\_v}(t-\Delta t) + \frac{1}{2}K_0\Delta t
\right)\\
K_2 = {\rm V\_and\_a}\left(
t-\frac{1}{2}\Delta t,
{\rm R\_and\_v}(t-\Delta t) + \frac{1}{2}K_1\Delta t
\right)\\
K_3 = {\rm V\_and\_a}\left(
t,
{\rm R\_and\_v}(t-\Delta t) + K_2\Delta t
\right)\\
\end{array}
$$

ここで、 $
{\rm R\_and\_v}(t)=\begin{pmatrix}\overrightarrow{r} & \overrightarrow{v}\end{pmatrix}$、$
{\rm R\_and\_v}(t-\Delta t)=\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}\\
\end{pmatrix}$、$
K_i = \begin{pmatrix}\overrightarrow{k_{r,i}} & \overrightarrow{k_{v,i}}\end{pmatrix}$ のように書くと、

$$\begin{pmatrix}\overrightarrow{r} & \overrightarrow{v}\end{pmatrix}
=\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}
\end{pmatrix} + \left(
\frac{1}{6}
\begin{pmatrix}\overrightarrow{k_{r,0}} & \overrightarrow{k_{v,0}}\end{pmatrix} +
\frac{1}{3}
\begin{pmatrix}\overrightarrow{k_{r,1}} & \overrightarrow{k_{v,1}}\end{pmatrix} +
\frac{1}{3}
\begin{pmatrix}\overrightarrow{k_{r,2}} & \overrightarrow{k_{v,2}}\end{pmatrix} +
\frac{1}{6}
\begin{pmatrix}\overrightarrow{k_{r,3}} & \overrightarrow{k_{v,3}}\end{pmatrix}
\right)\Delta t
$$
となる。ただし、

$$\begin{array}{lll}
\begin{pmatrix}\overrightarrow{k_{r,0}} & \overrightarrow{k_{v,0}}\end{pmatrix}
&=&
{\rm V\_and\_a}\left(t-\Delta t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}\\
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
\overrightarrow{v_{\rm before}} &
\overrightarrow{{\rm callback}\left(
t-\Delta t, \overrightarrow{r_{\rm before}}, \overrightarrow{v_{\rm before}}
\right)}
\end{pmatrix}
\end{array}
$$

$$\begin{array}{lll}
\begin{pmatrix}\overrightarrow{k_{r,1}} & \overrightarrow{k_{v,1}}\end{pmatrix}
&=&
{\rm V\_and\_a}\left(t-\frac{1}{2}\Delta t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}\\
\end{pmatrix} +
\frac{\Delta t}{2}
\begin{pmatrix}\overrightarrow{k_{r,0}} & \overrightarrow{k_{v,0}}\end{pmatrix}
\right)\\
&=&
{\rm V\_and\_a}\left(t-\frac{1}{2}\Delta t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{r,0}} &
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,0}}
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,0}} &
\overrightarrow{{\rm callback}\left(
t-\frac{1}{2}\Delta t,
\overrightarrow{r_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{r,0}},
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,0}}
\right)}
\end{pmatrix}
\end{array}
$$

$$\begin{array}{lll}
\begin{pmatrix}\overrightarrow{k_{r,2}} & \overrightarrow{k_{v,2}}\end{pmatrix}
&=&
{\rm V\_and\_a}\left(t-\frac{1}{2}\Delta t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}\\
\end{pmatrix} +
\frac{\Delta t}{2}
\begin{pmatrix}\overrightarrow{k_{r,1}} & \overrightarrow{k_{v,1}}\end{pmatrix}
\right)\\
&=&
{\rm V\_and\_a}\left(t-\frac{1}{2}\Delta t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{r,1}} &
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,1}}
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,1}} &
\overrightarrow{{\rm callback}\left(
t-\frac{1}{2}\Delta t,
\overrightarrow{r_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{r,1}},
\overrightarrow{v_{\rm before}} + \frac{\Delta t}{2}\overrightarrow{k_{v,1}}
\right)}
\end{pmatrix}
\end{array}
$$

$$\begin{array}{lll}
\begin{pmatrix}\overrightarrow{k_{r,3}} & \overrightarrow{k_{v,3}}\end{pmatrix}
&=&
{\rm V\_and\_a}\left(t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} & \overrightarrow{v_{\rm before}}\\
\end{pmatrix} +
\Delta t
\begin{pmatrix}\overrightarrow{k_{r,2}} & \overrightarrow{k_{v,2}}\end{pmatrix}
\right)\\
&=&
{\rm V\_and\_a}\left(t,
\begin{pmatrix}
\overrightarrow{r_{\rm before}} + \Delta t\overrightarrow{k_{r,2}} &
\overrightarrow{v_{\rm before}} + \Delta t\overrightarrow{k_{v,2}}
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
\overrightarrow{v_{\rm before}} + \Delta t\overrightarrow{k_{v,2}} &
\overrightarrow{{\rm callback}\left(
t,
\overrightarrow{r_{\rm before}} + \Delta t\overrightarrow{k_{r,2}},
\overrightarrow{v_{\rm before}} + \Delta t\overrightarrow{k_{v,2}}
\right)}
\end{pmatrix}
\end{array}
$$

4. 4次ルンゲクッタ法の実践

基本コードの __速度と位置を更新する メソッドを次の内容で書き換えればよい。

4次ルンゲクッタ法

    def __速度と位置を更新する(self):
        import numpy as np
        
        r_before = np.array(self.位置)
        v_before = np.array(self.速度)
        delta_t = self.__時刻幅
        t = self.時刻 + delta_t 
        callback = lambda *args: np.array(self.__callback(*args))

        # ルンゲクッタ係数の計算
        k_r_0 = v_before
        k_v_0 = callback(t - delta_t, r_before, v_before)

        k_r_1 = v_before + delta_t / 2 * k_v_0
        k_v_1 = callback(t - delta_t / 2, r_before + delta_t / 2 * k_r_0, k_r_1)

        k_r_2 = v_before + delta_t / 2 * k_v_1
        k_v_2 = callback(t - delta_t / 2, r_before + delta_t / 2 * k_r_1, k_r_2)

        k_r_3 = v_before + delta_t * k_v_2
        k_v_3 = callback(t, r_before + delta_t * k_r_2, k_r_3)

        # 次の位置と速度の計算
        r = r_before + delta_t / 6 * (k_r_0 + 2 * k_r_1 + 2 * k_r_2 + k_r_3)
        v = v_before + delta_t / 6 * (k_v_0 + 2 * k_v_1 + 2 * k_v_2 + k_v_3)

        # 更新
        self.位置 = r
        self.速度 = v

基本コードの __速度と位置を更新する メソッドを上記の内容で書き換えたものを rk4.py として保存する。

2章(オイラー法)でもやったシミュレーションを、4次ルンゲクッタ法でももう一度やってみよう。

(再掲)

初速度: $\begin{pmatrix}10\\0\end{pmatrix}$

加速度: 進行方向に垂直な方向に大きさ$10$

時刻幅: 1, 0.1, 0.01, 0.001 をそれぞれ試す

終了時刻: 300

理論的には、進行方向に垂直な方向に一定の大きさの加速度を加えているため、等速円運動を行うはずである。
その半径は速さ^2 / 加速度の大きさ = $100/10 = 10$ となるはずである。

シミュレーションは、2章の simulation.py について、
from euler import シミュレータ型
を
from rk4 import シミュレータ型
に置き換えた、次のコードで実行できる

simulation.py(ほぼ再掲)

from rk4 import シミュレータ型 # ここだけ置き換えた
from fit_circle import fit_circle
import math
import matplotlib.pyplot as plt
def get_a(t, r, v):
    # 進行方向に垂直な方向に10だけ加速し続ける
    x,y = 0,1
    
    進行方向 = math.atan2(v[y], v[x])
    加速度の方向 = 進行方向 + math.pi / 2
    加速度の大きさ = 10

    a_x = 加速度の大きさ * math.cos(加速度の方向)
    a_y = 加速度の大きさ * math.sin(加速度の方向)
    return (a_x, a_y)

# シミュレーションの設定    
シミュレーション = シミュレータ型(
    初期位置 = ( 0, 0),
    初速度   = (10, 0),
    時刻幅   = 1, # あるいは 0.1, 0.01, 0.001
    callback = get_a
)
終了時刻 = 300

# シミュレーションの実行
ts, xs, ys = [], [], []
x,y = 0,1
while シミュレーション.時刻 <= 終了時刻:
    ts += [シミュレーション.時刻]
    xs += [シミュレーション.位置[x]]
    ys += [シミュレーション.位置[y]]
    シミュレーション.進める()


# 結果のグラフ
# 1. 用意
幅 = max([abs(x) for x in xs] + [abs(y) for y in ys])*2.5
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111)
plt.xlim(-幅/2, 幅/2)
plt.ylim(-幅/2, 幅/2)
plt.xlabel("X")
plt.ylabel("Y")
ax.set_aspect("equal", adjustable="box")

# 2. 軌跡
ax.plot(xs, ys, lw = 1, color = "blue") # lw は線の幅。適当にいじる

# 3. 時刻
for t, x, y in zip(ts, xs, ys):
    if t % 5 != 0: # 時刻表示の間隔。適当にいじる
        continue

    ax.text(
        x, y, str(int(t)), fontsize = 8, # fontsize は適当にいじる
        ha = "center", va = "center", color = "black"
    )

# 4. 回帰円 (t=295,296,297,298,299,300)
idxs = [i for i in range(len(ts)) if ts[i] in [295, 296, 297, 298, 299, 300]]
points = [(xs[i], ys[i]) for i in idxs]
(xe, ye, r) = fit_circle(points)
print(f"回帰円半径: {r}")
xs = []
ys = []
for theta1000 in range(int(2*math.pi*1000)):
    theta = theta1000 / 1000
    xs += [xe + r * math.cos(theta)]
    ys += [ye + r * math.sin(theta)]
ax.plot(xs, ys, lw = 1, color = "#888888", linestyle=":", zorder=-1)
    # lw は適当にいじる

# 5. 表示
ax.grid(True)
plt.show()

これを実行すると、次の図の「rk4」の結果が得られた。

オイラー法 (euler) と比較して、遥かに半径10の等速円運動に近い結果を得ることが出来ている。

他にも、 simulation.py の get_a関数をいじることで、様々なシミュレーションができる。