前置き

　今回は常微分方程式の解法で有名な2次のルンゲ-クッタ法を実行するプログラムをpythonで作製しましたので紹介します。(n番煎じのような気はしますが、そこはご愛敬ということで……。)

プログラムの前に、標準形について

　プログラムを作製する前に、常微分方程式を計算するプログラムに必要な「標準形」について、ニュートンの運動方程式を例にあげて説明していきます。普通、ニュートンの運動方程式は以下のように表現されます。


\frac{d^2x}{d t^2} = \frac{1}{m} F (t,x,\frac{dx}{dt})

ここで$m$は物体の質量、$F$は物体に働く力で時間$t$・位置$x$・速度$dx/dt$に依存しています。
　さて、こちらを計算しやすいように標準形に変換していきましょう。まずは位置$x$を従属変数$y^{(0)}$、速度$dx/dt$を$y^{(1)}$とします。このようにおくと、上記の運動方程式は以下の2本の連立1次方程式に書き改めることができます。

\begin{eqnarray*}
\frac{dy^{(0)}}{dt} &=& y^{(1)} \\
\frac{dy^{(1)}}{dt} &=& \frac{1}{m} F (t,y^{(0)},y^{(1)}) \\
\end{eqnarray*}

　ここでさらに、$f^{(0)} = y^{(1)}$、$f^{(1)} = F/m$とおけば、以下のような非常にシンプルな形に表すことができます。


\frac{dy^{(i)}}{dt} = f^{(i)}(t, y^{(j)}) (i, j = 0,1)

　このように、常微分方程式を連立一次常微分方程式の形で表す方法を標準形と呼びます。初めにも言いましたが、プログラムの中で使用していますので確認してください。

　より一般的な常微分方程式の標準形は以下のように表されます。


\frac{dy^{(0)}}{dt} = f^{(0)}(t,y^{(i)}) \\
\frac{dy^{(1)}}{dt} = f^{(1)}(t,y^{(i)}) \\
\vdots \\
\frac{dy^{(N-1)}}{dt} = f^{(N-1)}(t,y^{(i)})

ここで右辺の$f^{(j)}$は、微分係数$dy^{(i)}/dt$は含みません。また、これらは$N$次元ベクトル${\bf y}$と${\bf f}$を用いて、以下のようにも書けます。

\frac{{\bf y}(t)}{dt} = {\bf f}(t,{\it y}), \\ 

  {\bf y} = \left(
    \begin{array}{c}
      y^{(0)}(t) \\
      y^{(1)}(t) \\
      \vdots \\
      y^{(N-1)}(t)
    \end{array}
  \right),

  {\bf f} = \left(
    \begin{array}{c}
      f^{(0)}(t,{\bf y}) \\
      f^{(1)}(t,{\bf y}) \\
      \vdots \\
      f^{(N-1)}(t,{\bf y})
    \end{array}
  \right),

プログラムの前に、ルンゲ-クッタ法について

　次に、常微分方程式を解く方法で有名な2次のルンゲ-クッタ法を説明します。(より精度の高い4次のルンゲ-クッタ法もありますが今回は割愛します。)
　ルンゲ-クッタ法は以下の微分方程式の積分から求めることができます。

\frac{dy}{dt}=f(t,y) \Rightarrow y = \int f(t,y) dt

ここで、時刻$t_n$の$y_n$から次の時刻$t_{n+1}$の$y_{n+1}$を求めるとすれば、上式は、

y_{n+1} = y_n + \int^{t_{n+1}}_{t_n} f(t,y) dt

と書けます。ここで、$f(t,y)$を積分区間の中点$t_{n+1/2}$のまわりでテイラー展開し、第2項まで残すと、以下となります。

f(t,y) \simeq f(t_{n+1/2}, y_{n+1/2}) + (t-t_{n+1/2})\frac{df}{dt}(t_{n+1/2}, y_{n+1/2}) + O(h^2)

ここで、$h$は時間の間隔$(t_{n+1}-t_n = h)$です。積分区間$t_n \le t \le t_{n+1}$では$(t-t_{n+1/2})$は正と負の値が等価に含まれているので定積分の値が消え、上式の$(t-t_{n+1/2})$の1次の項の積分は$0$となります。したがって、$f(t,y)$の積分、および、$y_{n+1}$を求める式は、

\begin{eqnarray*}
\int^{t_{n+1}}_{t_n} f(t,y) dt &\simeq& f(t_{n+1/2}, y_{n+1/2})h + O(h^3) \\
y_{n+1} &\simeq& y_n + f(t_{n+1/2}, y_{n+1/2})h + O(h^3)
\end{eqnarray*}

と表すことができます。さらに、時刻$t_{n+1/2} = t_n + h/2$の時の$y_{n+1/2}$をオイラー法によって近似的に求めると、

y_{n+1/2} \simeq y_n + \frac{h}{2} \frac{dy}{dt} = y_n + \frac{h}{2}f(t_n,y_n)

となります。したがって、時刻$t_n$の$y_n$から$y_{n+1}$を求める式は以下のようにまとめることができます。

\begin{eqnarray*}
y_{n+1} &\simeq & y_n + k_2 \\ \\
k_1  & = & hf(t_n,y_n) \\
k_2 & = & hf(t_n + \frac{h}{2}, y_n+ \frac{k_1}{2})
\end{eqnarray*}

上記のような数式を用いて常微分方程式の解を求めるアルゴリズムを、2次のルンゲクッタ法と言います。

ようやくルンゲ-クッタ法のプログラム

さて、前置きが長くなりましたが、ルンゲクッタ法のプログラムは非常に簡単なものとなります。使用したパッケージはnumpyだけです。

import numpy as np

#時間間隔h、時刻tn、その時のynの値のリストlist_yn、標準形の右辺funcを与えることで
#時刻tn+1の値yn+1のリストを返すルンゲクッタ法のアルゴリズム
def rk2(h,tn,list_yn, func):
    k1 = h*func(tn,list_yn)
    k2 = h*func(tn+h/2, list_yn+k1/2)
    return list_yn + k2

#初期時刻t_ini、終了時刻t_last、分割点数n、初期条件list_y_ini、
#標準形の右辺funcを与えることで
#時刻のリストlist_t、ルンゲクッタ法の解のリストlist_ynを返す。
def runge_kutta2(t_ini, t_last, n, list_y_ini, func):
    #時間間隔hを求める
    h = (t_last-t_ini)/n

    #list_t, list_ynの入れ物を作製
    list_t = np.zeros(n+1)
    list_yn = np.zeros([n+1,list_y_ini.shape[0]])
    
    #初期位置の格納
    list_t[0] = t_ini
    list_yn[0,:] = list_y_ini

    #ルンゲクッタ法で随時計算し結果を格納
    for i in range(1,n+1):
        list_yn[i,:] = rk2(h, list_t[i-1] ,list_yn[i-1,:] ,func)
        list_t[i] = list_t[i-1] + h

    return list_t, list_yn

ルンゲ-クッタ法のプログラムのテスト

　例の如く、プログラムのテストをしてみます。調和振動子の運動方程式を使って確認していきましょう。調和振動子の運動方程式の標準形は、

\begin{eqnarray*}
\frac{dy^{(0)}}{dt} &=& y^{(1)} \\
\frac{dy^{(1)}}{dt} &=& - \frac{k}{m}y^{(0)} \\
\end{eqnarray*}

と書けるので、プログラムにおいては、

def harmonic(t, list_y):
    return np.array([list_y[1], -4*np.pi**2*list_y[0]])

のように定義します。list_yは0列に位置、1列に速度が格納された1行2列のリストです。また、$k/m = 4 \pi ^2$としていますが、これは調和振動子の厳密解が、

\begin{eqnarray*}
x(t) &=& A \sin (\omega_0 t) \\
v(t) &=& \omega_0 A \cos (\omega_0 t) \\
\omega_0 &=& \sqrt{\frac{k}{m}}
\end{eqnarray*}

であるため、$\omega_0 = 2 \pi$とし、プログラムの確認を用意にするためです。

　さて、このharmonicを使用してルンゲ-クッタ法のプログラムをテストしていきます。時刻は0~1、分割点数を100、初期位置を0、初速を1として計算しました。


num = 100
t_ini = 0
t_last = 1
list_y_ini = np.array([0,1])

list_t2H, list_yn2H = runge_kutta2(t_ini, t_last, num, list_y_ini, harmonic)

さて上記の計算結果は以下のような図になりました。青い実線が位置、赤い破線が速度の計算結果です。また、同じ色の丸はそれぞれの厳密解を表していますが、なかなかに一致していることが分かります。実際に計算を行うと1E-3程度まで一致していることが分かりました。そこそこの精度ですね。

#感想
　今回は、2次のルンゲ-クッタ法を実行するプログラムを作製しました。数式を以下にコンピュータで取り扱っていくかの工夫は勉強していて面白いですね。
　数式やプログラムに間違いなど、ありましたらお知らせいただけると幸いです。

参考書籍
R.H.Landau・他 (2018)『実践Pythonライブラリー計算物理学I -数値計算の基礎/HPC/フーリエウェーブレット解析』 (小柳義夫・他訳) 朝倉書店