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

1次元の熱伝導方程式です。

ut=c22ux2 \frac{\partial u}{\partial t} = c^2 \frac{\partial^2 u}{\partial x^2}

1次元の波動方程式です。

2ut2=c22ux2 \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}

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

1+αut1+α=c22ux2 \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の分数階微分を組み合わせてシミュレーションします。

Grünwald–Letnikovの分数階微分を再掲します。

dαf(x)=limh01hαm=0(1)mΓ(α+1)Γ(m+1)Γ(αm+1)f(xmh) 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) が今から計算する値で u(_,tdt)u(\_, t - dt) までの値が既に得られているものとしています。

u(x,t)+m=1(1)m(1+αm)u(x,tmdt)=c2dt1+αdx2(u(xdx,t)2u(x,t)+u(x+dx,t)) \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}

Implicit FDMの形に整理します。

C1u(xdx,t)+C2u(x,t)+C1u(x+dx,t)=m=1(1)m(1+αm)u(x,tmdt)C1=c2dt1+αdx2,C2=(1+2C1) \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}

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

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

連立方程式を立てます。

Aut=d^(1+α,1,ut) \mathbf{A} \mathbf{u}^{t} = \hat{d} (1 + \alpha, 1, \mathbf{u}^{t})

A\mathbf{A}ut\mathbf{u}^t です。

A=[C2l0000C1C2C10000C1C2C1000000rC2],ut=[u(x0,t)u(x1,t)u(x2,t)u(xn1,t)] \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}

A\mathbf{A} に含まれる llrr です。

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

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

立てた連立方程式を左辺の ut\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))

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

Image of a result of 1d heat-wave simulation.

デモ

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

デモのコードは次のリンクから読むこともできます。

デモのコードを見る (github.com)

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

熱伝導-波動方程式を再掲します。

1+αut1+α=c22ux2 \frac{\partial^{1 + \alpha} u}{\partial t^{1 + \alpha}} = c^2 \frac{\partial^2 u}{\partial x^2}

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

Dα(u(x,t)u(x,tdt)dt)=c2(u(xdx,t)2u(x,t)+u(x+dx,t)dx2) \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}

左辺の小数階の微分を変形します。

u(x,t)+m=1(1)m(αm)u(x,tmdt)m=0(1)m(αm)u(x,t(m+1)dt)=c2dt1+αdx2(u(xdx,t)2u(x,t)+u(x+dx,t)) \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}

Implicit FDMの形に整理します。

C1u(xdx,t)+C2u(x,t)+C1u(x+dx,t)=m=1(1)m((αm)+(αm1))u(x,tmdt)C1=c2dt1+αdx2,C2=(1+2A) \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) を使って二項係数を一つにまとめます。

(αm)+(αm1)=(α+1m) \binom{\alpha}{m} + \binom{\alpha}{m - 1} = \binom{\alpha + 1}{m}

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

C1u(xdx,t)+C2u(x,t)+C1u(x+dx,t)=m=1(1)m(1+αm)u(x,tmdt) \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<α<10 < \alpha < 1 かつ C1C_1 が小さいときに発散することがあります。この発散がExplicit FDMの発散と関係がないことを確認するためにImplicit FDMで実装しました。