SciPy - 数理生物学演習

SciPy は NumPy の上に構築された科学技術計算ライブラリで，微分方程式・最適化・統計・信号処理などのアルゴリズムを提供する．ここでは常微分方程式の数値解法 scipy.integrate.solve_ivp を主に扱う．格子・画像の処理をおこなう scipy.ndimage は別ページ（SciPy：ndimage）にまとめた．

import numpy as np

微分方程式を解く：solve_ivp¶

P1-04 ではオイラー法を自分で実装したが，SciPy には高精度な数値解法が用意されている． solve_ivp は，導関数を返す関数と初期値を渡すだけで常微分方程式 (ordinary differential equation)（ODE）を解く．

解きたいのは「ある時刻 $t$ と状態 $x$ における変化率 $dx/dt$ 」を与える関数である．例として，ロジスティック成長 $\frac{dx}{dt} = r x (1 - x/K)$ を解く．

from scipy.integrate import solve_ivp


def logistic(t, x, r, K):
    """ロジスティック成長の導関数 dx/dt を返す"""
    return r * x * (1 - x / K)


# solve_ivp(導関数, [開始時刻, 終了時刻], [初期値], ...)
t_eval = np.linspace(0, 20, 200)
sol = solve_ivp(logistic, [0, 20], [0.1], args=(0.5, 1.0), t_eval=t_eval)

返り値の sol.t に時刻，sol.y に各変数の値（行ごとに1変数）が入る．

import matplotlib.pyplot as plt

plt.figure()
plt.plot(sol.t, sol.y[0])
plt.xlabel("time")
plt.ylabel("population")
plt.title("Logistic growth (solve_ivp)")
plt.show()

args で導関数に追加のパラメータ（ここでは r，K）を渡せる．t_eval は解を記録する時刻を指定する．

連立方程式（ベクトル形式）¶

変数が複数あるときは，状態をベクトルとして受け取り，変化率のベクトルを返す．ロトカ-ヴォルテラモデル（P1-05）を解いてみる．

def lotka_volterra(t, state, a, b, c, d):
    prey, predator = state
    d_prey = a * prey - b * prey * predator
    d_predator = -c * predator + d * prey * predator
    return [d_prey, d_predator]


t_eval = np.linspace(0, 50, 1000)
sol = solve_ivp(
    lotka_volterra,
    [0, 50],
    [10, 5],
    args=(1.0, 0.1, 1.5, 0.075),
    t_eval=t_eval,
    rtol=1e-8,
    atol=1e-10,
)

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(sol.t, sol.y[0], label="prey")
ax[0].plot(sol.t, sol.y[1], label="predator")
ax[0].set_xlabel("time")
ax[0].legend()
ax[1].plot(sol.y[0], sol.y[1])  # 相図
ax[1].set_xlabel("prey")
ax[1].set_ylabel("predator")
plt.show()

ロトカ-ヴォルテラ系の軌道は本来，相平面上で閉じた周期軌道になる．既定の許容誤差のままだと数値誤差が蓄積して軌道がわずかにずれて閉じないため，rtol（相対誤差）と atol（絶対誤差）を小さく設定して精度を上げている．

さらに進んだ機能¶

本演習では詳しく扱わないが，SciPy には次のようなモジュールもある．必要に応じて公式ドキュメントを参照するとよい．

scipy.linalg：np.linalg を拡張した線形代数（行列分解など）
scipy.optimize：関数の最小化・方程式の求根・曲線あてはめ
scipy.stats：確率分布・統計検定