Pythonによる計算物理（６）

前回、ロジスティック方程式について、数値解法を自作の関数を作成して解いた。

akh1116.hatenablog.com

今回は、scipyライブラリのintegrateモジュールを用いてコーディングをし、解いてみる。 自作で関数を作成するより、モジュールに渡す引数に注意すれば、簡潔かつバグも生じにくいと感じる。 以下、プログラム例である。

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

# dy/dt=f(t,y)の定義
def f_logistic(t, y):
    return y * (1.0 - y)    # ロジスティック方程式

# 解析解
def log_anlytical(y0, t):
    return y0 / (y0 + (1.0 - y0) * np.exp(-t))


# グラフ作成
def plot(t_sparse, y_sparse, t_dense, y_dense, y_exact, filename):
    fig, ax = plt.subplots()
    ax.plot(t_sparse, y_sparse, 'o', zorder=2, color='r', markerfacecolor='None', label="selected time points")
    ax.plot(t_dense, y_dense, '.', zorder=1, color='b', label="dense_output")
    ax.plot(t_dense, y_exact, '-', zorder=1.5, color='c', label="exact")
    ax.axhline(y=0, color='k', linestyle='dashed', zorder=0)   # 横線
    ax.axhline(y=1, color='k', linestyle='dashed', zorder=0)   # 横線
    ax.set_xlabel(r'$t$')  # xラベル(LaTeX表記)
    ax.set_ylabel(r'$y$')  # yラベル(LaTeX表記)
    ax.legend()            # 凡例
    plt.show()
    fig.savefig(filename)

# メイン処理
def main():
    y0 = np.array([1e-3,], dtype=float)      # 初期値(1成分のNumpy配列)
    t_start = 0     # 初期時刻
    t_end = 20.0    # 最終時刻

    # Scipy integrateモジュールで微分方程式を解く
    sol = solve_ivp(f_logistic, (t_start, t_end), y0, dense_output=True)
    print(sol.message)    # ソルバーのメッセージを表示
    print("sol.t.shape=", sol.t.shape)    # (n_points,)
    print("sol.y.shape=", sol.y.shape)    # (1, n_points)

    # グラフ描画用のデータ(dense_output)
    nt = 101     # 時刻の解像度を100分割
    t = np.linspace(t_start, t_end, nt)    # 時刻の等間隔メッシュを生成
    y = sol.sol(t)     # メッシュ点上の時刻tにおけるy(t)の値を取得

    print("t.shape=", t.shape)
    print("y.shape=", y.shape)

    # グラフ作成
    dirname = "image/"
    filename = dirname + 'logistic_scipy.png'
    plot(sol.t, sol.y[0], t, y[0], log_anlytical(y0, t), filename)



if __name__ == '__main__':
    main()

以下、計算結果のプロット図である。
dense_outputを有効にし、任意の時刻での結果を出力させている。 実際は、モジュール内で時間刻みを調整してるようで、より少ないステップ数で効率良く計算してようである。