🫁

odeintを使ってLorenz方程式を解いてみる

2025/03/09に公開

Python

はじめに

今回は、scipyのodeintというものを用いて、Lorenz方程式を解いてみます。odeintの使い方について詳しくまとまっているものがあまりなかったので、各パラメタについてもまとめていこうと思います。

使用するライブラリなど

今回使用するライブラリなどは以下の通りです。

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

odeintの使用方法

odeintの使用法はとても簡単で、以下のような引数があります。

sol = odeint(func, y0, t, args Dfun, col_deriv, full_output, printmessg, tfirst)

この中でよく設定して使うものは、func, y0, t, argsです。それ以降の引数はいじることが少なく、デフォルトのまま使用することが多いようです。

各パラメタの役割は以下の通りです。

func: 微分方程式 $\dfrac{dy}{dt} = f(y, t)$ において、右辺を計算する関数です。引数は通常 $(y, t)$ の順に定義することが多いですが、もし、 $(t, y)$ の順序で定義した際は、最後の引数tfirstについて、tfirst=Trueとする必要があります。
y0: 微分方程式の初期条件を指定するパラメタです。スカラーまたはベクトルで与えることができます。
t: 解を求める時間点を指定するものです。この配列は単調増加または単調減少である必要があります。
arg: func に渡す追加の引数を指定します。
Dfun: 関数 func のヤコビ行列（微分係数）を計算する関数です。これを指定することにより、高速化が期待できます。引数順序が (t, y) の場合は、tfirst=True を指定します。
col_deriv: ヤコビ行列が列方向に微分されている場合に True を指定します（高速化）。デフォルトは False（行方向に微分）になっています。
printmessg: Trueを指定すると、計算結果や収束についてのメッセージを出力してくれます。
tfirst: 先ほども説明したように、関数 func や Dfun の引数順序が (t, y) の場合に True を指定します。

これらは、以下のサイトを参考にしました。
scipy.integrate.odeint

Lorenz方程式を解いてみる

さて、説明が終わったので、実際に使ってみましょう。せっかくなので、全てのパラメタを使って実行してましょう。

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

# ローレンツ方程式の定義
def lorenz(y, t, sigma, beta, rho):
    x, y, z = y  # 状態変数を分解して表式を簡潔にする
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

# ヤコビ行列（ローレンツ方程式の偏微分）を定義
def jacobian(y, t, sigma, beta, rho):
    x, y, z = y
    return [
        [-sigma, sigma, 0],
        [rho - z, -1, -x],
        [y, x, -beta]
    ]

# 初期条件と時間配列
y0 = [1.0, 1.0, 1.0]  # 初期値
t = np.linspace(0, 50, 10000)  # 時間配列

# パラメータ
sigma = 10.0
beta = 8.0 / 3.0
rho = 28.0

# odeintの設定と実行
solution = odeint(
    func=lorenz,
    y0=y0,
    t=t,
    args=(sigma, beta, rho),  # 関数への追加引数
    Dfun=jacobian,           # ヤコビ行列を指定
    col_deriv=True,          # ヤコビ行列が列方向微分であることを指定
    printmessg=True         # 計算情報を取得するためにTrueを指定
)

# 結果と計算情報を分離
trajectory = solution[0]  # 解軌道
info_dict = solution[1]   # 計算情報

# 計算情報を出力（オプション）
print("計算情報:", info_dict)

# 結果をプロット（3次元プロット）
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(trajectory[:, 0], trajectory[:, 1], trajectory[:, 2])
ax.set_title("Lorenz Attractor")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()

これを実行すると、以下のような図が出力されます。

すごく綺麗に描けていますね

そして、printmessg=Trueにしたことによって出てくるメッセージは以下のようになります

計算情報: {'hu': array([0.00171406, 0.00171406, 0.00449938, ..., 0.00605863, 0.00605863,
       0.00605863]), 'tcur': array([5.93752319e-03, 1.10796990e-02, 1.72931351e-02, ...,
       4.99952095e+01, 4.99952095e+01, 5.00012681e+01]), 'tolsf': array([0., 0., 0., ..., 0., 0., 0.]), 'tsw': array([0., 0., 0., ..., 0., 0., 0.]), 'nst': array([  10,   13,   15, ..., 6691, 6691, 6692], dtype=int32), 'nfe': array([   23,    29,    33, ..., 13843, 13843, 13845], dtype=int32), 'nje': array([0, 0, 0, ..., 0, 0, 0], dtype=int32), 'nqu': array([4, 4, 5, ..., 6, 6, 6], dtype=int32), 'imxer': -1, 'lenrw': 68, 'leniw': 23, 'mused': array([1, 1, 1, ..., 1, 1, 1], dtype=int32), 'message': 'Integration successful.'}

基本的に、図を見たいことや、解を求めたいことがほとんどだと思うので、このメッセージについては、いじらずにデフォルトのFalseのままでもいいかもしれません

以上、odeintを紹介してみました。
何か使用の仕方についてのアドバイスやご指摘などありましたら、ぜひご連絡いただけると幸いです。

はじめに

使用するライブラリなど

odeintの使用方法

Lorenz方程式を解いてみる

Discussion