SciPy - 数理生物学演習

SciPy は NumPy の上に構築された科学技術計算ライブラリで，微分方程式・最適化・統計・信号処理などのアルゴリズムを提供する．ここでは本演習で特に有用な2つの機能を扱う．

scipy.integrate.solve_ivp：常微分方程式（ODE）の数値解法
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（絶対誤差）を小さく設定して精度を上げている．

格子上の場の処理：ndimage¶

scipy.ndimage は格子（配列）に対する各種の処理を提供する．拡散方程式や反応拡散モデル（P2-03）で必要なラプラシアン (Laplacian)は ndimage.laplace で一度に計算できる．自分で近傍を取り出して差分を書く代わりに，配列全体を1行で処理できる．

from scipy.ndimage import laplace

f = np.zeros((5, 5))
f[2, 2] = 1.0  # 中央に1つだけ値を置く

# 離散ラプラシアン（上下左右の和 - 4×中央）
# mode="wrap" は周期境界条件（端が反対側につながる）
print(laplace(f, mode="wrap"))

[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0.  1. -4.  1.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.]]

手書きの差分計算と比べて簡潔で，しかも高速である．拡散方程式の更新は次のように書ける．

# u の拡散による1ステップ更新（D: 拡散係数, dt: 時間刻み, dh: 空間刻み）
def diffuse(u, D, dt, dh):
    return u + D * laplace(u, mode="wrap") / dh**2 * dt


u = np.zeros((20, 20))
u[10, 10] = 100.0
for _ in range(50):
    u = diffuse(u, D=1.0, dt=0.1, dh=1.0)

plt.figure()
plt.pcolormesh(u)
plt.gca().set_aspect("equal")
plt.colorbar()
plt.title("diffusion via ndimage.laplace")
plt.show()

任意のカーネルとの畳み込み (convolution)は ndimage.convolve で行える（近傍の重み付き和を計算する）．

さらに進んだ機能¶

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

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