【python】1次元拡散方程式の数値解法(陽解法)

はじめに

拡散方程式とは拡散が生じている物質あるいは物理量の密度のゆらぎを記述する偏微分方程式のことです。

ja.wikipedia.org

ちょっと何言ってるかわからないかもしれませんが、物体を熱が伝わる熱伝導の時間変化もこの方程式で表すことができます。

今回は例として1次元の熱伝導をシミュレーションしてみましょう。

はじめに
拡散方程式
- 安定化条件
実装
- 例題設定
- 実装
おわりに

拡散方程式

一般的な拡散方程式は以下の形で表されます。

$\dfrac{\partial\phi}{\partial t}=\nabla \cdot (D\nabla \phi)$

ここで $\phi$ は拡散物質の密度、 $D$ は拡散係数です。
1次元の熱伝導方程式は上の式から以下のように導出できます。

$\dfrac{\partial u}{\partial t}=\alpha \dfrac{\partial^2 u}{\partial x^2}$

ここで $u(x, t)$ は温度、 $x$ は位置、 $\alpha$ は熱拡散率を表します。

数値解析が行えるようにこれを差分化します。時間を $\Delta t$ 、空間を $\Delta x$ の長さの微小区間に等分して、上の式をテイラー展開してごにょごにょすると以下の式が得られます。

$\dfrac{u_i^{n+1}}{\Delta t}=\alpha \dfrac{u_{i+1}^n-2u_i^n+u_{i-1}^{n}}{(\Delta x)^2}$

ここで $u_i^n$ は $x = i \Delta x, t = n \Delta t$ での温度を表します。
式変形すると以下のようになります。

$u_i^{n+1}=u_i^n + \dfrac{\alpha \Delta t}{(\Delta x)^2}(u_{i+1}^n-2u_i^n+u_{i-1}^{n})$

この時左辺は時刻 $n+1$ での値、右辺は時刻 $n$ での値となるため現在時刻 $n$ での値を使うことで時刻 $n+1$ での値が求まるということを示しています。
こうして現在時刻の結果から未来の状態を予測する方法を陽解法といい、今回は陽解法をつかって解いてみます。

上の式より、境界条件と $\Delta t, \Delta x$ の値を決めてあげれば簡単にシミュレーションができそうです。

安定化条件

陽解法で解く場合は時間刻み $\Delta t$ を大きくしていくと計算が早く終わりますが、大きくしすぎると誤差の影響により解が発散して解けなくなります。
上の式の安定性の条件は以下であることが知られています。

$\dfrac{\alpha \Delta t}{(\Delta x)^2} < \dfrac{1}{2}$

空間刻みを小さく設定した場合はそれに伴って時間刻みも小さくしないといけなくなり、計算時間が非常に長くなるということです。

実装

例題設定

例として以下のような問題を考えます。

長さ $l=1$ の棒の片側を常に1度であっためて、逆側は断熱されていて熱の移動がないような状況を考えてみます。棒の初期温度は室温で0度としましょう。 $\alpha$ は棒の材質等で決まりますが、今回は簡単のため適当に1としておきます。

まとめると以下のようになります。

$\begin{cases} \dfrac{\partial u}{\partial t}=\alpha \dfrac{\partial^2 u}{\partial x^2} & \text{if $0 < x < l, t > 0$}\\ u = 1 & \text{if $x = 0, t \ge 0$}\\ \dfrac{\partial u}{\partial x} = 0 & \text{if $x = l, t \ge 0$}\\ u = 0 & \text{if $0 < x < l, t = 0$} \end{cases}$

ここで境界条件についてですが、 $x=1$ ではどの時刻においても同じ値になるように設定しており、そういう境界条件をディリクレ条件といい、 $x=l$ では温度の勾配がどの時刻においてもゼロになるように設定しており、こちらはノイマン条件と呼びます。シミュレーション対象に応じて使い分けます。

実装

それでは実際に実装していきましょう。

クラスをいくつか作りました。それぞれの役割は以下です。

Config: シミュレーションの設定値を管理する
Geom: シミュレーション対象の定義
BoundaryCondition: 境界条件の定義(オーバーライドしてカスタマイズできるように)
InitialCondition: 初期値の定義(オーバーライドしてカスタマイズできるように)
ExplicitSolver: 拡散方程式のソルバー

拡張性が高くなるように心がけました。

import dataclasses

import matplotlib.pyplot as plt
import numpy as np


@dataclasses.dataclass
class Config:
    """解析設定"""

    dx: float
    dt: float
    steps: int
    alpha: float = 1.0


class Geom:
    """ジオメトリ定義"""

    def __init__(self, length: float, dx: float):
        self.length = length
        self.dx = dx
        self.x = np.arange(0, length + dx, dx)
        self.N = len(self.x)


class BoundaryCondition:
    """境界条件の抽象クラス"""

    def apply(self, u: np.ndarray) -> None:
        raise NotImplementedError


class DirichletNeumannBC(BoundaryCondition):
    """左側ディリクレ条件、右側ノイマン条件の境界条件設定クラス"""

    def __init__(self, left: float = 0.0):
        self.left = left

    def apply(self, u: np.ndarray) -> None:
        u[0] = self.left
        u[-1] = u[-2]


class InitialCondition:
    """初期値設定の抽象クラス"""

    def apply(self, x: np.ndarray) -> np.ndarray:
        raise NotImplementedError


class UniformIC(InitialCondition):
    """すべて初期値を一定値に設定する"""

    def __init__(self, temp=0):
        self.temp = temp

    def apply(self, x: np.ndarray) -> np.ndarray:
        return x * 0 + self.temp


class ExplicitSolver:
    def __init__(
        self, config: Config, geom: Geom, bc: BoundaryCondition, ic: InitialCondition
    ):
        self.cfg = config
        self.geom = geom
        self.bc = bc
        self.u = ic.apply(geom.x)
        self.bc.apply(self.u)

        self.r = config.alpha * config.dt / config.dx**2
        if self.r > 0.5:
            print("Warning: 安定性に問題があります(r > 0.5)")

    def step(self):
        u = self.u
        r = self.r

        u_new = u.copy()
        u_new[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2])

        # 境界条件を適用
        self.bc.apply(u_new)
        self.u = u_new

    def run(self):
        for _ in range(self.cfg.steps):
            self.step()
        return self.u


def plot(x: np.ndarray, u: np.ndarray, title="Solution"):
    fig, ax = plt.subplots(figsize=(6, 4))
    ax.plot(x, u)
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_ylim(-0.2, 1.2)
    ax.set_title(title)
    plt.show()


cfg = Config(dx=0.01, dt=0.00005, steps=5000, alpha=1.0)
geom = Geom(length=1.0, dx=cfg.dx)
bc = DirichletNeumannBC(left=1.0)
ic = UniformIC(temp=0)
solver = ExplicitSolver(cfg, geom, bc, ic)
plot(geom.x, solver.u, title="Initial State")
u = solver.run()
plot(geom.x, u, title="Solution Steps: {}".format(cfg.steps))