Pythonとscipy.integrate.odeintの使い方

Pythonのscipy.integrate.odeintは、常微分方程式の数値解を求めるための関数です。この記事では、その基本的な使い方について説明します。

常微分方程式の解法

まずは、単純な一階の常微分方程式を解いてみましょう。以下のような微分方程式を考えます。

$$\frac{dy}{dt} = -y$$

この微分方程式は、odeintを使って以下のように解くことができます。

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

# 1階常微分方程式（放射性崩壊）
def func_dydt(y, t):
    dydt = -y
    return dydt

# 2d可視化
def plot2d(t_list, y_list, t_label, y_label):
    plt.xlabel(t_label) #x軸の名前
    plt.ylabel(y_label) #y軸の名前
    plt.grid() #点線の目盛りを表示
    plt.plot(t_list, y_list)
    plt.show()

# メイン実行部
if (__name__ == '__main__'):
    # 常微分方程式（放射性崩壊）
    t_list = np.linspace(0.0, 10.0, 1000)
    y_init = 1.0 #初期値
    y_list = odeint(func_dydt, y_init, t_list)
    print(y_list)

    # 可視化
    plot2d(t_list, y_list[:, 0], "$t$", "$y(t)$")

このコードは、時間に対するyの変化を計算し、結果をプロットします。結果は、放射性物質が指数的に減少するグラフになります。

連立常微分方程式の解法

次に、連立常微分方程式の一種である、Lorenz方程式を解いてみましょう。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 1階連立常微分方程式（Lorenz方程式）
def func_lorenz(var, t, p, r, b):
    dxdt = -p*var[0] +p*var[1]
    dydt = -var[0]*var[2] +r*var[0] -var[1]
    dzdt = var[0]*var[1] -b*var[2]
    return [dxdt, dydt, dzdt]

# 3d可視化
def plot3d(t_list, var_list):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.set_xlabel("$x$") #x軸の名前
    ax.set_ylabel("$y$") #y軸の名前
    ax.set_zlabel("$z$") #z軸の名前
    ax.plot(var_list[:, 0], var_list[:, 1], var_list[:, 2])
    plt.show()

# メイン実行部
if (__name__ == '__main__'):
    # 1階連立常微分方程式（Lorenz方程式）
    t_list = np.linspace(0.0, 100.0, 10000)
    p = 10
    r = 28
    b = 8/3
    var_init = [0.1, 0.1, 0.1] #3次元座標上での初期値
    var_list = odeint(func_lorenz, var_init, t_list, args=(p, r, b))
    print(var_list)

    # 可視化
    plot3d(t_list, var_list)

このコードは、Lorenz方程式の解を計算し、結果を3Dプロットします。結果は、ローレンツ・アトラクタと呼ばれる複雑な形状になります。

以上が、Pythonのscipy.integrate.odeintの基本的な使い方です。この関数を使うことで、様々な常微分方程式の数値解を求めることができます。ただし、常微分方程式の数値計算は、初期値だけでなく刻み幅の影響も受けるため、注意が必要です。.