Pythonで学ぶ拡散方程式

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使って、拡散方程式と呼ばれる偏微分方程式を調べたいと思います。偏微分方程式の数値的解法の導入、本記事で使っている記法の詳細については以下の記事で解説しています：

pianofisica.hatenablog.com

あわせて読んでみてください。

拡散方程式
熱核
熱核を用いた解
差分方程式
Pythonによる実装
解析解との比較

拡散方程式

上に引用した過去記事では、輸送方程式

$\qquad \displaystyle{\frac{\partial f}{\partial t}=-c\frac{\partial f}{\partial x}}$

を考え、その解析解と数値解を比較しつつ、偏微分方程式の数値計算の手法（差分法）について解説しました。本記事では、問題をもう１ステップ高度にして、拡散方程式（熱方程式、熱伝導方程式とも）と呼ばれる偏微分方程式

$\qquad \displaystyle{\frac{\partial f}{\partial t}=D\frac{\partial^2 f}{\partial x^2}\qquad(\ast)}$

を考えてみます。簡単のために、ここでは $D$ は定数としましょう。本質的な違いは、右辺の１階微分が２階微分に置き換わったことです。

熱核

方程式 $(\ast)$ を境界条件

$\qquad \displaystyle{f(0,x)=\varphi(x)}$

のもとで解くために、（唐突ではありますが）熱核（heat kernel） $K(t,x,y)$ を導入します：

$\qquad \displaystyle{K(t,x,y)=\frac{1}{\sqrt{4\pi D t}}\,e^{-\frac{(x-y)^2}{4Dt}} }$

性質を少し調べてみましょう。まず方程式 $(\ast)$ の微分演算子を作用させてみます：

import sympy as sp
x, y = sp.symbols('x, y', real = True)
t, D = sp.symbols('t, D', positive = True)

def K(t,x,y):
    return 1/sp.sqrt(4*sp.pi*D*t) * sp.exp(-(x-y)**2/(4*D*t)) 

Eq = sp.diff(K(t,x,y), (t, 1)) - D * sp.diff(K(t,x,y), (x, 2))
sp.simplify(Eq)

より

$\qquad \displaystyle{\left(\frac{\partial f}{\partial t}-D\frac{\partial^2 f}{\partial x^2}\right)K(t,x,y)=0\qquad(1)}$

であって、熱核は拡散方程式の"解"であることがわかります。次に、 $x$ を固定して、 $y$ に関して実軸全体にわたって積分してみます：

sp.integrate( K(t,x,y), (y, -sp.oo, sp.oo) )

よって

$\qquad \displaystyle{\int_{-\infty}^{\infty}K(t,x,y)dy=1 }$

です。また、 $y=0$ として、いくつかの $t$ の値に対する関数の概形を図示してみます（ $D=1$ とします）：

import matplotlib.pyplot as plt
import numpy as np

T = [1e0, 1e-1, 1e-2, 1e-3, 1e-4]

for t in T:
    X = np.linspace( -5, 5, int(1e5) )
    Y = 1/( 4 * np.pi * t ) * np.exp( -X**2 / (4*t) )
    plt.plot(X, Y)
    plt.show()

グラフは左から順に $t= 10^0,\,10^{-1},\,10^{-2},\,10^{-3},\,10^{-4}$ です。以上から、熱核 $K(t,x,y)$ は $t$ を小さく取るほど $x=0(=y)$ の点に局在する関数であり、 $t\to0$ の極限でDiracのデルタ関数を与えることがわかります：

$\qquad \displaystyle{\lim_{t\to0}K(t,x,y)=\delta(x-y) \qquad(2)}$

熱核を用いた解

さて、もともと考えていた問題だった、方程式 $(\ast)$ の境界条件 $f(0,x)=\varphi(x)$ のもとでの解ですが、それは熱核を用いて

$\qquad \displaystyle{f(t,x)= \int_{-\infty}^\infty K(t,x,y)\varphi(y)dy \qquad (\dagger)}$

で与えられます。これを確かめます。まず方程式 $(\ast)$ に代入してみます：

$\qquad \displaystyle{\begin{aligned}\left(\frac{\partial f}{\partial t}-D\frac{\partial^2 f}{\partial x^2}\right)f(t,x)&= \int_{-\infty}^\infty \left(\frac{\partial f}{\partial t}-D\frac{\partial^2 f}{\partial x^2}\right)K(t,x,y)\varphi(y)dy\\ &=0 \end{aligned}}$

ここで上で見た性質 $(1)$ を使いました。よって確かに方程式の解になっています。次に境界条件をチェックしてみます：

$\qquad \displaystyle{\begin{aligned}f(0,x)&= \int_{-\infty}^\infty \lim_{t\to0}K(t,x,y)\varphi(y)dy\\ &= \int_{-\infty}^\infty \delta(x-y)\varphi(y)dy\\ &=\varphi(x)\end{aligned}}$

ここで性質 $(2)$ およびデルタ関数の性質を使っています。よって境界条件も満たされています。

差分方程式

差分化の手法については過去記事で扱っているので、その導入については省略します。ポイントとなるのは２階偏導関数の差分化です。２つのTaylor展開

$\qquad \displaystyle{\begin{aligned} g(x+h)&=g(x)+hg'(x)+\frac{h^2}{2!}g''(x)+\frac{h^3}{3!}g''(x)+\mathcal{O}(h^4)\\ g(x-h)&=g(x)-hg'(x)+\frac{h^2}{2!}g''(x)-\frac{h^3}{3!}g''(x)+\mathcal{O}(h^4) \end{aligned}}$

の両辺を足し合わせると

$\qquad \displaystyle{g(x+h)+g(x-h)=2g(x)+h^2g''(x)+\mathcal{O}(h^4) }$

より

$\qquad \displaystyle{g''(x)=\frac{g(x+h)-2g(x)+g(x-h)}{h^2}+\mathcal{O}(h^2) }$

です。よって２階微分の差分化として

$\qquad \displaystyle{g''(x_n)=\frac{g(x_{n+1})-2g(x_n)+g(x_{n-1})}{h^2} }$

とすることができます。以上から $(\ast)$ に対応する差分方程式として

$\qquad \displaystyle{\frac{f^{(i+1)}_n-f^{(i-1)}_n}{2h_t}-D\, \frac{f^{(i)}_{n+1}-2f^{(i)}_n+f^{(i)}_{n-1}}{h_x^2}=0 }$

を立式できます。ただし、時間微分に関する１階偏導関数については中心差分をとりました。式を整理すると

$\qquad \displaystyle{f^{(i+1)}_n=\frac{2Dh_t}{h_x^2}\, \left(f^{(i)}_{n+1}-2f^{(i)}_n+f^{(i)}_{n-1}\right)+f^{(i-1)}_n }$

となります。境界条件から $f^{(0)}_n=\varphi(x_n)$ は既知です。しかしこの漸化式を解くためには、 $f^{(1)}_n$ の情報も必要になります。そこで、 $f^{(1)}_n$ についてのみ、時間微分について前進差分を用いることにします：

$\qquad \displaystyle{f^{(1)}_n=\frac{Dh_t}{h_x^2}\, \left(f^{(0)}_{n+1}-2f^{(0)}_n+f^{(0)}_{n-1}\right)+f^{(0)}_n }$

Pythonによる実装

境界条件を

$\qquad \displaystyle{f(0,x)=\varphi(x)=e^{-x^2/2}}$

として、上記の定式化をPythonで実装してみましょう。まず入力値です。

Nt = 100
Nx = 200
D = 0.5

tmin = 0.0
tmax = 3.0
xmin = -15.0
xmax = 15.0

ht = ( tmax - tmin ) / Nt
hx = ( xmax - xmin ) / Nx

次に境界条件の関数の配位を表示してみます。

X = []
f0 = []

for n in range(Nx):
    X.append( xmin + n*hx )
    f0.append(np.exp(-X[n]**2/2))
    n += 1

import matplotlib.pyplot as plt

plt.xlim(-2.0, 2.0)  # x軸の表示範囲
plt.ylim(-.25, 1.25)  # y軸の表示範囲

plt.plot( X, f0 )
plt.show()

より次のグラフを得ます。

次に $f^{(1)}_n$ を求めます：

f1 = [(f0[0]+f0[1])/2]    # 差分方程式から求まらない点は前の値を引き継ぐことにしました

for n in range(1,Nx-1):
    f1.append( D*ht/(hx**2) * (f0[n+1] - 2*f0[n] + f0[n-1]) + f0[n] )
    n += 1
f1.append((f0[Nx-1]+f0[Nx-2])/2)    # 差分方程式から求まらない点は前の値を引き継ぐことにしました

$f^{(0)}_n$ と合わせて図示してみます。

plt.xlim(-2.0, 2.0)
plt.ylim(-.25, 1.25)

plt.plot( X, f0 )
plt.plot( X, f1 )
plt.show()

逐次的に差分方程式から $f^{(i)}_n$ を求めます：

F = [f0,f1]

for m in range(2,Nt-2):
    f = [(F[m-1][0]+F[m-1][1])/2]    # 差分方程式から求まらない点は前の値を引き継ぐことにしました
    for n in range(1, Nx-1):
        f.append( 2*D*ht/(hx**2) * (F[m-1][n+1] - 2*F[m-1][n] + F[m-1][n-1]) + F[m-2][n] )
        n += 1
    f.append((F[m-1][Nx-1]+F[m-1][Nx-2])/2)
    F.append(f)
    m += 1

関数の時間発展のタイムラプスを作成してみます。

for k in range(5):
    plt.xlim(-2.0, 2.0)
    plt.ylim(-.25, 1.25)
    plt.plot( X, F[4*k] )
    k += 1
plt.show()

から次のグラフを得ます：

初期状態のガウス関数（青線）から橙→緑→赤→紫と、時間の経過とともに山の高さが徐々に低くなって裾が広がっていく様子がわかります。これは温度の高い部分が、低い部分にじわじわと移っていく現象を表しています（熱伝導）。

解析解との比較

さて、解析解は $(\dagger)$ で与えられました。ゆえに、境界条件 $f(0,x)=\varphi(x)=e^{-x^2/2}$ のもとでは

x, y = sp.symbols('x, y', real = True)
t, D = sp.symbols('t, D', positive = True)

def K(t,x,y):
    return 1/sp.sqrt(4*sp.pi*D*t) * sp.exp(-(x-y)**2/(4*D*t)) 
    
def phi(y):
    return sp.exp(-y**2/2) 

sol = sp.integrate( K(t,x,y) * phi(y), (y, -sp.oo, sp.oo) )
sp.simplify(sol)

より

$\qquad \displaystyle{f(t,x)=\frac{1}{\sqrt{2Dt+1}}\exp\left(-\frac{x^2}{2(2Dt+1)}\right) }$

となります。これを数値計算の重ね合わせてグラフに表示し、見比べてみます。

D = 0.5
T = []
for k in range(5):
    T.append( tmin + 4*ht*k )

for n in range( len(T) ):
    X = np.linspace( xmin, xmax, Nx )
    Y = 1/np.sqrt( 2 * D * T[n] + 1 ) * np.exp( - X**2 / ( 2 * ( 2 * D * T[n] + 1 ) ) )
    plt.plot(X, Y)
    plt.plot( X, F[4*n] )
    plt.show()