拡散方程式 - 数理生物学演習

拡散方程式を離散化し，数値的に解く．1次元から始め，2次元へ拡張する．

import matplotlib.pyplot as plt
import numpy as np

01-01. １次元の拡散方程式¶

拡散方程式の空間2階微分を有限差分法 (finite difference method)で近似すると，各点の次の値を近傍の値から計算できる．時間方向はオイラー法，空間方向は中心差分で離散化する．

固定境界条件（左端を $u_\text{left}$ ，右端を 0 に固定）でシミュレーションしてみよう．

D = 1.0  # 拡散係数
dt = 0.05  # 時間方向の刻み幅
dh = 0.5  # 空間方向の刻み幅
u_left = 1.0  # 左端の固定値 u(0, t)
x_end = 10.0  # x の終端
t_end = 50.0  # t の終端

x = np.arange(0, x_end, dh)
n_steps = int(t_end / dt)

# u の初期化（左端のみ u_left，他は 0）
u = np.zeros(len(x))
u[0] = u_left

t_list = [0.0]  # 時刻を記録するリスト
u_list = [u.copy()]  # u を記録するリスト

for n in range(1, n_steps + 1):
    u_new = u.copy()
    # 内部の各点を中心差分で更新
    for i in range(1, len(u) - 1):
        u_new[i] = u[i] + D * (u[i - 1] + u[i + 1] - 2 * u[i]) / dh**2 * dt
    # 固定境界条件
    u_new[0] = u_left
    u_new[-1] = 0.0

    u = u_new
    t_list.append(n * dt)
    u_list.append(u.copy())

一定間隔で時刻を選び，濃度分布を重ねて描く．

plt.figure(dpi=100)
plt.ylim(-0.05, 1.05)
for i in range(0, len(u_list), 100):
    plt.plot(x, u_list[i], ".-", label=f"t = {t_list[i]:.1f}")
plt.xlabel("x")
plt.ylabel("u")
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
plt.show()

拡散性の物質が，左端から時間とともに右へ広がっていく．

01-02. ２次元への拡張と周期境界¶

2次元では上下左右の4近傍を使う．ここでは周期境界条件 (periodic boundary condition)を採用する．

まずは全セルを1つずつ更新する二重ループ版で実装しよう．隣接インデックスを剰余（%）で求めると，端が自動的に反対側へ折り返り，周期境界条件になる．

def diff_2d_loop(u, D, dh, dt):
    """2次元拡散方程式の1ステップ更新（二重ループ版，周期境界条件）."""
    n_y, n_x = u.shape
    u_new = np.zeros_like(u)
    for i in range(n_y):
        for j in range(n_x):
            # 剰余（%）で隣接インデックスを求めると，端は反対側へ折り返る
            lap = (
                u[(i - 1) % n_y, j]
                + u[(i + 1) % n_y, j]
                + u[i, (j - 1) % n_x]
                + u[i, (j + 1) % n_x]
                - 4 * u[i, j]
            ) / dh**2
            u_new[i, j] = u[i, j] + D * lap * dt
    return u_new

中央の1区画だけに初期濃度を置いて拡散させる．二重ループ版は遅いので，ここでは時間を短めにとる．

D = 1.0  # 拡散係数
dt = 0.05  # 時間方向の刻み幅
dh = 1.0  # 空間方向の刻み幅
u_peak = 100.0  # 初期濃度

x = np.arange(-25, 25, dh)
y = np.arange(-25, 25, dh)
xmesh, ymesh = np.meshgrid(x, y)

t_end = 10.0  # 二重ループ版は遅いので短めにする
n_steps = int(t_end / dt)

# 中央の1区画だけに初期濃度を置く
u = np.zeros((len(y), len(x)))
u[len(y) // 2, len(x) // 2] = u_peak

for n in range(n_steps):
    u = diff_2d_loop(u, D, dh, dt)

fig, ax = plt.subplots(dpi=120)
ax.set_aspect("equal")
pcm = ax.pcolormesh(xmesh, ymesh, u, vmin=0, vmax=1)
fig.colorbar(pcm, ax=ax, shrink=0.8)
plt.show()

中央から同心円状に広がる様子が得られる．

01-03. ラプラシアンによるベクトル化¶

拡散方程式の空間2階微分はラプラシアン (Laplacian)の差分で表される．同じ計算を NumPy のベクトル化で高速化する． 4近傍の和は，配列全体を上下左右にずらして足し合わせれば一度に計算できる． np.roll は配列をずらし，はみ出た要素を反対側へ回り込ませる関数で，これがそのまま周期境界条件になる．

def laplacian(f, dh):
    """周期境界条件でのラプラシアン（ベクトル化）.

    np.roll で配列を上下左右にずらして4近傍を取得する．
    はみ出た要素は反対側へ回り込むため，周期境界条件になる．
    """
    return (
        np.roll(f, 1, axis=0)
        + np.roll(f, -1, axis=0)
        + np.roll(f, 1, axis=1)
        + np.roll(f, -1, axis=1)
        - 4 * f
    ) / dh**2


def diff_2d(u, D, dh, dt):
    """2次元拡散方程式の1ステップ更新（ベクトル化，周期境界条件）."""
    return u + D * laplacian(u, dh) * dt

ベクトル化により高速になったので，全ステップの状態を記録しながら最後まで計算する．

D = 1.0  # 拡散係数
dt = 0.05  # 時間方向の刻み幅
dh = 1.0  # 空間方向の刻み幅
u_peak = 100.0  # 初期濃度

x = np.arange(-25, 25, dh)
y = np.arange(-25, 25, dh)
xmesh, ymesh = np.meshgrid(x, y)

t_end = 50.0  # t の終端
n_steps = int(t_end / dt)

# 中央の1区画だけに初期濃度を置く
u = np.zeros((len(y), len(x)))
u[len(y) // 2, len(x) // 2] = u_peak

t_list = [0.0]  # 時刻を記録するリスト
u_list = [u.copy()]  # u を記録するリスト
for n in range(1, n_steps + 1):
    u = diff_2d(u, D, dh, dt)
    t_list.append(n * dt)
    u_list.append(u.copy())

01-04. 可視化¶

記録した状態の1つを取り出して可視化する． pcolormesh で濃度を色で表す．

i = 200  # t = 10.0 の時点
fig, ax = plt.subplots(dpi=150)
ax.set_aspect("equal")
pcm = ax.pcolormesh(xmesh, ymesh, u_list[i], vmin=0, vmax=1)
fig.colorbar(pcm, ax=ax, shrink=0.8)
ax.set_title(f"t = {t_list[i]:.1f}")
plt.show()

同じ状態を3次元の曲面として描くと，濃度の山が広がっていく様子がよくわかる．

fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection="3d")
surf = ax.plot_surface(xmesh, ymesh, u_list[i], cmap="viridis", vmin=0, vmax=1)
fig.colorbar(surf, shrink=0.8)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u")
plt.show()

Solution to Exercise 1

拡散係数を大きくし，初期濃度を2箇所に置いてみる．

D = 2.0  # 拡散を速くする

# 初期濃度を2箇所に置く
u = np.zeros((len(y), len(x)))
u[15, 15] = u_peak
u[35, 35] = u_peak

for n in range(200):
    u = diff_2d(u, D, dh, dt)

fig, ax = plt.subplots(dpi=120)
ax.set_aspect("equal")
pcm = ax.pcolormesh(xmesh, ymesh, u, vmin=0, vmax=1)
fig.colorbar(pcm, ax=ax, shrink=0.8)
plt.show()