1次元移流拡散方程式をPythonでインタラクティブに解く

🐍 Fortran経験者のためのPython入門

～1次元移流拡散方程式をインタラクティブに解いてみよう～

この記事はChatGPT4oにより生成されました

---

🎯 この記事の目的

• Pythonで 1次元移流拡散方程式を数値的に解く
• Fortranの知識を活かして、Pythonらしい書き方（OOP + NumPy）を理解する
• ipywidgets を使って、パラメータを自由に変えられるUI付きのシミュレーションを実装する！

---

👤 対象読者

• Fortranで数値計算はやっていたが、Pythonはまだ初心者
• @dataclass や ipywidgets に興味がある
• Pythonでも構造化して、きれいに数値コードを書いてみたい！

---

🧪 解く方程式：1次元移流拡散方程式

$$
\frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x} = D \frac{\partial^2 u}{\partial x^2}
$$

• $u(x, t)$：スカラー場（濃度や温度）
• $v$：移流速度（定数）
• $D$：拡散係数（定数）

数値スキームには FTCS（Forward Time Centered Space）法 を使います。

---

🏗️ コード全体構成

1. Grid1D：格子点・初期条件の管理
2. SolverConfig：物理パラメータの管理
3. AdvectionDiffusionSolver：1ステップの計算
4. plot_result：結果の描画
5. run_simulation：全体の統括処理
6. interact()：スライダーUIの追加（ipywidgets）

---

📦 フルコード

from dataclasses import dataclass
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact

# ===============================
# 1. 空間格子と初期条件の管理クラス
# ===============================
@dataclass
class Grid1D:
    nx: int
    length: float
    dx: float = None
    x: np.ndarray = None
    u: np.ndarray = None

    def __post_init__(self):
        self.dx = self.length / (self.nx - 1)
        self.x = np.linspace(0, self.length, self.nx)
        self.u = np.exp(- ((self.x - 0.5 * self.length)**2) / (2 * (0.1 * self.length)**2))

# ===============================
# 2. 物理パラメータを管理する設定クラス
# ===============================
@dataclass(frozen=True)
class SolverConfig:
    v: float
    D: float
    dt: float

# ===============================
# 3. 移流拡散の数値スキーム（1ステップ）
# ===============================
@dataclass
class AdvectionDiffusionSolver:
    grid: Grid1D
    config: SolverConfig

    def __call__(self):
        u = self.grid.u
        dx = self.grid.dx
        dt = self.config.dt
        v = self.config.v
        D = self.config.D

        advection = -v / (2 * dx) * (u[2:] - u[:-2])
        diffusion = D / dx**2 * (u[2:] - 2 * u[1:-1] + u[:-2])
        u_new = u.copy()
        u_new[1:-1] = u[1:-1] + dt * (advection + diffusion)

        self.grid.u = u_new

# ===============================
# 4. 結果のプロット
# ===============================
def plot_result(grid: Grid1D, title: str):
    plt.plot(grid.x, grid.u, label=title)
    plt.xlim(0, grid.length)
    plt.ylim(0, 1.1)
    plt.xlabel("x")
    plt.ylabel("u")
    plt.title(title)
    plt.grid(True)
    plt.legend(loc="lower left")
    plt.show()

# ===============================
# 5. シミュレーション実行関数
# ===============================
def run_simulation(nx=101, length=2.0, v=1.0, D=0.01, dt=0.001, steps=200):
    grid = Grid1D(nx=nx, length=length)
    config = SolverConfig(v=v, D=D, dt=dt)
    solver = AdvectionDiffusionSolver(grid=grid, config=config)

    for _ in range(steps):
        solver()

    plot_result(grid=grid,
                title=f"v={v}, D={D}, dt={dt}, length={length}, steps={steps}")

# ===============================
# 6. ウィジェットUIでパラメータ変更
# ===============================
interact(
    run_simulation,
    v=(0.0, 2.0, 0.1),
    D=(0.0, 0.1, 0.005),
    dt=(0.0001, 0.01, 0.0005),
    steps=(10, 300, 10),
    length=(0.5, 5.0, 0.1)
)