1. 概要

scipy.integrate.odeintを用いることで、常微分方程式の数値解を得る。

2. 問題設定

今回は以下のような微分方程式を題材とする。

2重井戸ポテンシャル $V(x) = - \frac{1}{2}x^2 + \frac{1}{4} x^4$ の中を動く質量1の粒子を考えると、運動方程式は以下の形を取ることが知られている。

\ddot{x} = x - x^3

位置 $x$ と速度 $\dot{x}$ の初期値が与えられたとき、解析的に解くことなしに粒子の挙動を知りたいとする。

3. 方法

3.1. 一階への変形

先述の運動方程式は二階の微分方程式だが、新たな変数 $v$ を $v = \dot{x}$ として定義してやることにより、以下のように一階の形に変形できる。この変形はscipy.integrate.odeintを利用する際に重要である。

\begin{eqnarray}
  \begin{cases}
    \dot{x} = v & \\
    \dot{v} = x - x^3 &
  \end{cases}
\end{eqnarray}

3.2. 実装

まずは必要なライブラリをimportする。ここではJupyter Notebookでのグラフ描画用も含む。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
%matplotlib inline

次に、先ほどの微分方程式をもとに、「従属変数と独立変数を受けて導関数を返す関数」を実装する。従属変数は配列であってよい。今回の例では $x, v$ の時間微分を考えているため、x, vの配列を第一引数・時間tを第二引数とするような関数を実装する。

def f(y, t):
    x, v = y
    return [
        v,
        x - x**3
    ]

ほとんど微分方程式そのままであることに注意。今回は第二引数tを使っていないが、強制外力を受ける振動子のように非自律的な方程式を扱う場合には用いることができる。

あとは、考えたい時間ステップと初期値を与えてscipy.integrate.odeintを利用するだけ。
今回は時間を0-10で200分割し、初期値は位置が2, 速度0とした。

t = np.linspace(0, 10, 200)
y0 = [2, 0]
result = odeint(f, y0, t)

odeintの返り値（今回はresult変数）には数値計算の結果が収められており、shapeは(200, 2)となっている。$x, v$ それぞれについて各step（時間を200分割していた）の値が入っている。

3.3. 可視化

result変数の中身をplotする。周期的になっていることがわかる。

fig, axes = plt.subplots(1, 2)
axes[0].plot(result[:, 0])
axes[0].set_xlabel("step")
axes[0].set_ylabel("x")
axes[1].plot(result[:, 1])
axes[1].set_xlabel("step")
axes[1].set_ylabel("v")

本筋とは関係ないが、ついでにstreamplotを行ってベクトル場を可視化してみよう。横軸・縦軸に $x, v$ をとり、エネルギー $E = \frac{1}{2} \dot{x}^2 - \frac{1}{2} x^2 + \frac{1}{4} x^4$ で色分けした。軌道はエネルギーの等高線になることが知られている。

w = 3
Y, X = np.mgrid[-w:w:100j, -w:w:100j]
U = Y
V = X - X**3

E = 0.5 * Y**2 - 0.5 * X**2 + 0.25 * X**4

plt.streamplot(X, Y, U, V, color=E)

4. 参考

ストロガッツ非線形ダイナミクスとカオス (p177)
scipy.integrate.odeint
Streamplot - Matplotlib

SciPyで常微分方程式の数値解を得る