Pythonによる計算物理（４）

常微分方程式の数値計算例題-ロジスティック方程式

ここでは、常微分方程式の数値計算の例題として、ロジスティック方程式を扱ってみる。

ロジスティック方程式

ロジスティック方程式は、生物の個体数の増加を表すモデルであり、以下のような形の方程式として表される。

ここで、 $y$ は生物の個体数を最大値で規格化した無次元数で、 $0 \le y \le 1$ の領域で与えられる。

解析解

ロジスティック方程式の解析解は以下の形で表される。

ここで、 $y_0=y(0)$ は初期条件であり、時刻 $t=0$ における個体数である。

数値計算

数値計算の手法として、オイラー法、中点法、ルンゲ・クッタ法の３手法を試し、解析解との比較をしてみた。

コード例（自作関数使用）

まず、Scipyなどの常微分方程式のソルバーを使わず、自分でオイラー法、中点法、ルンゲ・クッタ法を関数として自作したコード例を示す。

# このプログラムではロジスティック方程式をPythonで数値計算を行う
# Scipyなどのライブラリを使わないバージョン
# Coded by AKH1116  2024.2.3

import numpy as np
import matplotlib.pyplot as plt
import os

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

# オイラー法による時間発展
def euler(f, y, dt):
    return y + f(y) * dt

# 中点法による時間発展
def euler_midpoint(f, y, dt):
    k1 = f(y)
    k2 = f(y+k1 * 0.5 * dt)
    return y + k2 * dt

# ルンゲ・クッタ法による時間発展
def runge_kutta(f, y, dt):
    k1 = f(y)
    k2 = f(y + k1 * dt/2)
    k3 = f(y + k2 * dt/2)
    k4 = f(y + k3 * dt)
    return y + (k1 + 2*k2 + 2*k3 + k4) * dt / 6

# 常微分方程式を解く関数
# f: 関数f(y)
# y0: 初期値
# tmax: 終了時間
# nt: 時間の分割数
def solve_ode_euler(f, y0, tmax, nt):
    t = np.linspace(0, tmax, nt)
    dt = t[1] - t[0]
    y = [y0,]   # 結果を格納するリスト
    for _ in range(nt-1):   # 時間ループ
        y.append(euler(f, y[-1], dt))    # 時刻を1つ進めて結果をリストyに追加
    return t, np.array(y)                # 結果をNumpy配列として返す

def solve_ode_euler_midpoint(f, y0, tmax, nt):
    t = np.linspace(0, tmax, nt)
    dt = t[1] - t[0]
    y = [y0,]   # 結果を格納するリスト
    for _ in range(nt-1):   # 時間ループ
        y.append(euler_midpoint(f, y[-1], dt))    # 時刻を1つ進めて結果をリストyに追加
    return t, np.array(y)                         # 結果をNumpy配列として返す

def solve_ode_euler_rk(f, y0, tmax, nt):
    t = np.linspace(0, tmax, nt)
    dt = t[1] - t[0]
    y = [y0,]   # 結果を格納するリスト
    for _ in range(nt-1):   # 時間ループ
        y.append(runge_kutta(f, y[-1], dt))    # 時刻を1つ進めて結果をリストyに追加
    return t, np.array(y)                         # 結果をNumpy配列として返す

def exact_sol(y0, tmax, nt):
    t = np.linspace(0, tmax, nt)
    dt = t[1] - t[0]
    return t, y0 / (y0 + (1 - y0) * np.exp(-t))


def main():
    # 画像を保存するディレクトリ、ファイル名を指定
    dirname = "image/"
    filename = dirname + 'logistic.png'

    # ロジスティック方程式を解く
    t, y1 = solve_ode_euler(f_logistic, y0=1e-3, tmax=20, nt=101)
    t, y2 = solve_ode_euler_midpoint(f_logistic, y0=1e-3, tmax=20, nt=101)
    t, y3 = solve_ode_euler_rk(f_logistic, y0=1e-3, tmax=20, nt=101)
    t, y4 = exact_sol(y0=1e-3, tmax=20, nt=101)

    # グラフを描画
    fig, ax = plt.subplots()
    ax.plot(t, y1, 'x-', label='Euler')
    ax.plot(t, y2, 'o-', label='Euler Midpoint')
    ax.plot(t, y3, '+-', label='Runge-Kuuta')
    ax.plot(t, y4, ':', label='Excat')
    ax.set_xlabel('t')
    ax.set_ylabel('y')
    ax.legend()
    plt.show()
    fig.savefig(filename)
    print("Run!")

if __name__ == '__main__':
    main()

数値計算結果

以下に数値計算結果を示す。オイラー法よりも、中点法、ルンゲ・クッタ法といった高次精度の手法では解析解とほぼ同じ結果が得られている。

日々の雑記帳

日々感じたこと、学んだことの備忘録