SciPyで常微分方程式の数値解を得る
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
Author And Source
この問題について(SciPyで常微分方程式の数値解を得る), 我々は、より多くの情報をここで見つけました https://qiita.com/sw1227/items/1ece810f93d1bd352b5e著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .