微分方程式の数値解法

scipyのライブラリの利用

常微分方程式

簡単な例 (調和振動子)

調和振動子は2階微分方程式 $$ \frac{d^2 x}{dt^2} = - k x $$ で与えられる。速度$v$を定義すると、 $$ \begin{eqnarray} \frac{dx}{dt} & = & v\\ \frac{dv}{dt} & = & -kx \end{eqnarray} $$ のように２変数１階の微分方程式であらわされる。

一般に、複数の変数の高階の微分方程式でも

$$ \begin{eqnarray} \frac{d x_1}{dt} & = & f_1(x_1,x_2,\cdots, x_N) \\ \frac{d x_2}{dt} & = & f_2(x_1,x_2,\cdots, x_N) \\ & \vdots & \\ \frac{d x_N}{dt} & = & f_N(x_1,x_2,\cdots, x_N) \end{eqnarray} $$

のように１階の連立微分方程式で表すことができる。（右辺にあらわに$t$が入っている場合もある。）

pythonのscipyにはこのような微分方程式の数値解を求めるモジュールodeintが用意されている。

これを用いるには、あらかじめ微分方程式の右辺（上記の$f_1$から$f_N$に相当する関数：$\{x_1,\dots,x_N\}$や係数を引数とし、$N$個の関数の値をを返すもの）を定義しておき、odeintの引数に与える。

調和振動子の方程式に適用した例が次である。

In [17]:
import numpy as np
from scipy.integrate import odeint

def harmonic_eq(x,t, k):
    return x[1], - k*x[0]

dt = 0.01
n_steps = 5000
t = np.arange(0, n_steps*dt, dt)
x = np.zeros((2))
x[0] = 2.0    # initial position
x[1] = 0.0    # initial velocity
k = 0.5       # parameter in diff. eq.
y = odeint(harmonic_eq, (x[0], x[1]), t, args=(k,)) # The definition of the differential equation (1st argument) must be defined in the above.
  • 第２引数に初期値をタプル(tuple)あるいはリストで与える
  • 第３引数に変数$t$のリストを与える。上記では、刻み幅をdtとしてdt刻みでn_step個のリストを作成している
  • 第４引数に関数定義内で使うパラメータをタプルで与える (1個しかない場合は、カッコ内の末尾にカンマが必要）

以下は、結果$(t,y)$を図示する関数である。

In [19]:
import matplotlib.pyplot as plt

def plot_func(t,x):
    fig = plt.figure(figsize=(12,5))
    axes = fig.add_subplot(1,2,1)
    axes.plot(t, x[:,0], 'r', label="$x(t)$")
    axes.plot(t, x[:,1], 'b', label="$v(t)$")
    axes.set_xlabel("$t$")
    
    plt.legend(loc='upper left')
    #plot a trajectory
    axes = fig.add_subplot(1,2,2)
    axes.plot(x[:,0], x[:,1],'#ff8800')
    axes.set_xlabel("$x(t)$")
    axes.set_ylabel("$v(t)$")
    plt.show()

plot_func(t,y)

3変数以上の微分方程式と特徴的な力学系の例

https://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/simulationMethods/NumericalMethods2.html

に記した。

In [ ]: