heat_wave_equation

熱伝導方程式と波動方程式を方程式を分数階微分でつなぎます。

\[ \frac{\partial^{1 + \alpha} u}{\partial t^{1 + \alpha}} = c^2 \frac{\partial^2 u}{\partial x^2} \]

ここでは、この式を熱伝導-波動方程式 (Heat-Wave Equation) と呼びます。拡散-波動方程式 (Diffusion-Wave Equation) と呼ばれることもあります。

Implcit FDM (Finite Difference Method) とGrünwald–Letnikovの分数階微分を組み合わせてシミュレーションします。

\[ d^{\alpha} f(x) = \lim_{h \to 0} \frac{1}{h^{\alpha}} \sum_{m=0}^{\infty} (-1)^{m} \frac{\Gamma(\alpha + 1)}{\Gamma(m + 1)\Gamma(\alpha - m + 1)} f(x - mh) \]

熱伝導-波動方程式を計算できる形にするために、左辺をGrünwald–Letnikovの分数階微分、右辺を有限差分として展開します。今回は \(u(\_, t)\) が今から計算する値で \(u(\_, t - dt)\) までの値が既に得られているものとしています。

\[ \begin{aligned} u(x, t) +& \sum_{m=1}^{\infty} (-1)^{m} \binom{1 + \alpha}{m} u(x, t - m\,dt)\\ =& \frac{c^2 dt^{1 + \alpha}}{dx^2} \left( u(x - dx, t) -2u(x, t) + u(x + dx, t) \right)\\ \end{aligned} \]

\[ \begin{aligned} C_1 u(x - dx, t) + C_2 u(x, t) + C_1 u(x + dx, t) =& \sum_{m=1}^{\infty} (-1)^{m} \binom{1 + \alpha}{m} u(x, t - m\,dt)\\ C_1 = \frac{c^2 dt^{1 + \alpha}}{dx^2},\quad C_2 = -(1 + 2C_1) \end{aligned} \]

分数階微分の演算をまとめるために関数 \(\hat{d}(\alpha, k, \mathbf{u}^{t})\) を定義します。

\[ \hat{d}(\alpha, k, \mathbf{u}^{t}) = \sum_{m=k}^{\infty} (-1)^{m} \binom{\alpha}{m} \mathbf{u}^{t - m\,dt} \]

\[ \mathbf{A} = \begin{bmatrix} C_2 & l & 0 & 0 & \cdots & 0 & 0\\ C_1 & C_2 & C_1 & 0 & \cdots & 0 & 0\\ 0 & C_1 & C_2 & C_1 & \cdots & 0 & 0\\ & \vdots & & & & \vdots &\\ 0 & 0 & 0 & 0 & \cdots & r & C_2 \end{bmatrix} ,\quad \mathbf{u}^t = \begin{bmatrix} u(x_0, t)\\ u(x_1, t)\\ u(x_2, t)\\ \vdots\\ u(x_{n-1}, t) \end{bmatrix} \]

\[ \begin{aligned} l =& C_1 \left (1 + L \right )\\ r =& C_1 \left (1 + R \right ) \end{aligned} \]

\(L\) と \(R\) の値は境界条件によって決まります。固定端のときは \(L = R = 0\) 、自由端のときは \(L = R = 1\) となります。

立てた連立方程式を左辺の \(\mathbf{u}^{t}\) について解くことでシミュレーションを1ステップ進めることができます。

コードを実行すると1次元の熱-波のシミュレーションを行い、結果を heat_wave1d.png に画像として書き出します。

連立方程式のソルバにnumpy.linalg.solve、画像の書き出しにImageioを使っています。

import numpy
import imageio
from scipy.special import binom


class HeatWave1D():
    def __init__(self, length, c, dx, dt, alpha, attenuation):
        """
        1次元の熱と波のシミュレータ。

        d^(1 + alpha) u / dt^(1 + alpha) = c^2 d^2 u / dx^2

        :param length: 波を表す1次元配列の長さ。
        :param c: 波の速度[m/s]。
        :param dx: 配列の要素間の距離[m]。
        :param dt: シミュレーションの1ステップで進む時間[s]。
        :param alpha: 分数階微分の階数。 [0, 1] の範囲。
        :param attenuation: 厳密でない波の減衰係数。 [0, 1] の範囲。
        """
        # u(x, t) -> self.wave[t][x]
        self.length = length
        self.fracDiffDepth = 256
        self.wave = numpy.zeros((self.fracDiffDepth, self.length))

        self.field = numpy.zeros(self.length)

        self.setParameters(c, dx, dt, alpha, attenuation)
        self.pick(0, 0)

    def value(self):
        """
        描画用に最新の波を返す。
        """
        return self.wave[0]

    def initBoundary(self):
        """
        境界条件を指定。 0 で固定端。 1 で自由端。
        """
        self.L = 0
        self.R = 0

    def setParameters(self, c, dx, dt, alpha, attenuation):
        """
        パラメータを更新する。
        シミュレーションの途中でパラメータが変更された時に呼び出す。
        """
        self.c = c
        self.dx = dx
        self.dt = dt
        self.alpha = alpha
        self.attenuation = attenuation

        self.refreshConstants()
        self.initBoundary()
        self.initMatrix()

    def refreshConstants(self):
        """
        シミュレーションで用いる定数を設定。
        """
        self.C1 = (self.c / self.dx)**2 * self.dt**(1 + self.alpha)
        self.C2 = -(1 + 2 * self.C1)

        self.fracCoefficients = [(-1)**m * binom(1 + self.alpha, m)
                                 for m in range(self.fracDiffDepth)]

    def initMatrix(self):
        """
        Implicit finite difference method で解く必要のある方程式の設定。
        ax = b の a。
        """
        mat = numpy.zeros((self.length, self.length))

        mat[0][0] = self.C2
        mat[0][1] = self.C1 * (1 + self.L)

        last = self.length - 1
        for i in range(1, last):
            mat[i][i - 1] = self.C1
            mat[i][i] = self.C2
            mat[i][i + 1] = self.C1

        mat[last][last - 1] = self.C1 * (1 + self.R)
        mat[last][last] = self.C2

        self.a = mat

    def step(self):
        """
        シミュレーションの1ステップ。
        """
        self.wave = numpy.roll(self.wave, 1, axis=0)
        if self.pickY != 0:
            self.wave[1][self.pickX] = self.pickY
        self.wave[0] = self.attenuation * numpy.linalg.solve(
            self.a,
            self.fractionalDifference(),
        )

    def fractionalDifference(self):
        """
        分数階微分の計算。方程式 ax = b の bを返す。
        """
        self.field.fill(0)
        for m in range(1, len(self.wave)):
            self.field += self.fracCoefficients[m] * self.wave[m]
        return self.field

    def reset(self):
        """
        波を 0 で埋めて初期状態に戻す。
        """
        self.wave.fill(0)

    def pick(self, x, y):
        """
        x, y で指定した位置の波をつまむ。

        :param x: 波をつまむ場所。[0, 1] の範囲。
        :param y: 波をつまむ高さ。任意の実数。
        """
        self.pickX = int((self.length - 1) * numpy.clip(x, 0.0, 1.0))
        self.pickY = y


if __name__ == "__main__":
    length = 512
    wave1d = HeatWave1D(length, 64, 0.1, 0.01, 0.5, 1)

    result = []
    wave1d.pick(0.5, 1)
    for t in range(0, length):
        wave1d.step()
        result.append(wave1d.value())

    imageio.imwrite("heat_wave1d.png", numpy.array(result))

実行結果の画像です。横はあるステップでの波の状態、縦は時間で上から下に向かって進んでいます。

\(\alpha\) の値を変えて熱伝導-波動をシミュレーションした動画です。

\[ \frac{\partial^{1 + \alpha} u}{\partial t^{1 + \alpha}} = c^2 \frac{\partial^2 u}{\partial x^2} \]

熱伝導-波動方程式の両辺の整数階微分を有限差分で変形します。

\[ \begin{aligned} D^{\alpha} \left( \frac{u(x, t) - u(x, t - dt)}{dt} \right) =& c^2 \left( \frac{u(x - dx, t) -2u(x, t) + u(x + dx, t)}{dx^2} \right)\\ \end{aligned} \]

\[ \begin{aligned} u(x, t) +& \sum_{m=1}^{\infty} (-1)^{m} \binom{\alpha}{m} u(x, t - m\,dt) - \sum_{m=0}^{\infty} (-1)^{m} \binom{\alpha}{m} u(x, t - (m + 1)\,dt)\\ =& \frac{c^2 dt^{1 + \alpha}}{dx^2} \left( u(x - dx, t) -2u(x, t) + u(x + dx, t) \right)\\ \end{aligned} \]

\[ \begin{aligned} & C_1 u(x - dx, t) + C_2 u(x, t) + C_1 u(x + dx, t)\\ &= \sum_{m=1}^{\infty} (-1)^{m} \left( \binom{\alpha}{m} + \binom{\alpha}{m - 1} \right) u(x, t - m\,dt)\\ & C_1 = \frac{c^2 dt^{1 + \alpha}}{dx^2},\quad C_2 = -(1 + 2A) \end{aligned} \]

ここで二項係数の加算についての性質 (Pascal’s rule) を使って二項係数を一つにまとめます。

「連立方程式を立てる」で変形した式と同じ式が出てきます。

\[ \begin{aligned} C_1 u(x - dx, t) + C_2 u(x, t) + C_1 u(x + dx, t) =& \sum_{m=1}^{\infty} (-1)^{m} \binom{1 + \alpha}{m} u(x, t - m\,dt) \end{aligned} \]

自由端にすると発散します。固定端でも \(0 < \alpha < 1\) かつ \(C_1\) が小さいときに発散することがあります。この発散がExplicit FDMの発散と関係がないことを確認するためにImplicit FDMで実装しました。

分数階微分で1次元の熱伝導-波動方程式

連立方程式を立てる

実装

デモ

整数階と小数階で分けて変形

問題点