差分法で偏微分方程式を解いてみる

(常)微分方程式を差分法で解く記事を投稿しました。

運動方程式とか、波動方程式とか、マクスウェル方程式など、身近な物理現象を記述してくれる式は全て微分方程式です。解が解析的に求まらなくても、数値的に微分方程式を満たす解を求める事が出来る場合もあります。今回の記事では、微分の差分化を使って簡単な運動方程式を解きます。

その続きとして、数学的にはとても難しい偏微分方程式を数値的に解いてみましょう。何も考えずに差分法で解くと、痛い目を見ます。痛い目をみて、何故そうなったか考えようという記事です。
元ネタは以下の書籍です。

差分の種類と近似精度
熱伝導方程式の数値解
Von Neumannの安定性解析
まとめ

差分の種類と近似精度

差分の種類と近似の精度について簡単にまとめておきます。
差分には以下の3種類があります。

前進差分 $ \frac{dx}{dt} \sim \frac{f(x+\Delta t) -f(x)}{\Delta t} $
後退差分 $ \frac{dx}{dt} \sim \frac{f(x) -f(x- \Delta t )}{\Delta t} $
中心差分 $ \frac{dx}{dt} \sim \frac{f(x +\Delta t ) -f(x- \Delta t )}{2 \Delta t} $

関数のテイラー展開を考えると、近似できる理由と、その精度が分かります。
点$a $の周りでの$f(x) $のテイラー展開をn次の項で打ち切ったものは以下の式で表されます。
$$\begin{eqnarray}
f(x) =\sum _{k=1}^n \frac{1}{k!} \frac{df^k}{d^k x}(a) (x-a)^k + O( (x-a)^{n+1} )
\end{eqnarray}$$
ただし、$0! =1 , \frac{df^0}{d^0 x}(a) = f(a) $です。
また、$ O(g(x)) $はランダウの記号です。
ランダウの記号 $ O(g(x)) $ の意味は、大雑把には$ g(x) $で割れば定数になるという事で、定数$C $が存在して以下の式が成り立つ事を言います。
$$\begin{eqnarray}
\left | \frac{f(x)}{g(x)} \right| \leq C
\end{eqnarray}$$
さて、1次の項で打ち切った$ f(x+ \Delta t), f(x-\Delta t )$ のテイラー展開を考えましょう。
$$\begin{eqnarray}
f(x+\Delta t) &=& f(x) +\frac{df}{dx}\Delta t + \frac{1}{2!} \frac{d^2 f }{dt^2}\Delta t ^2 + O(\Delta t^3 )\\
f(x-\Delta t) &=& f(x) – \frac{df}{dx}\Delta t + \frac{1}{2!} \frac{d^2 f }{dt^2}\Delta t ^2 + O(\Delta t^3 ) \\
f(x+ \Delta t)- f(x -\Delta t )&=& 2f(x) +O(\Delta t^3 )
\end{eqnarray}$$
以上より、前進差分(後退差分)、中心差分では以下の式が成り立ちます。
$$\begin{eqnarray}
\frac{dx}{dt} &=& \frac{f(x+\Delta t) -f(x) }{\Delta t} + O(\Delta t ) \\
\frac{dx}{dt} &=& \frac{f(x-\Delta t) – f(x)}{\Delta t} + O(\Delta t ) \\
\frac{dx}{dt} &=& \frac{f(x+ \Delta t)- f(x -\Delta t )}{2 \Delta t} + O(\Delta t^2 )
\end{eqnarray}$$
前進差分や後退差分は、ずれとして$ \Delta t$程度の量を含んでいるという解釈が出来ます。
一方で、中心差分のずれは $ \Delta t^2$程度です。$ \Delta t $ が1より小さい数だと思うと $ \Delta t ^2 < \Delta t $なので、近似の精度は中心差分の方が良い事が分かります。

熱伝導方程式の数値解

差分法で解いてみる方程式は、1次元の熱伝導方程式です。
導出方法は熱力学の教科書を読むと大体載っていますが、インターネット上でも解説してくれてる人は多いです。例えばMITが無料で公開している講座の資料はわかりやすいです。
https://ocw.mit.edu/courses/mathematics/18-303-linear-partial-differential-equations-fall-2006/lecture-notes/heateqni.pdf
導出は別の資料を見てもらうことにして、扱いたい熱伝導方程式を書きます。
$$\begin{eqnarray}
\frac{\partial u}{ \partial t} = k \frac{\partial ^2 u}{\partial x^2}
\end{eqnarray}$$
早速差分化してみましょう。$t $については前進差分化を取り、$x$については中心差分を取ります。¹
$$\begin{eqnarray}
\frac{u(x,t+\Delta t) – u(x,t)}{\Delta t} = k \frac{u(x+\Delta x, t) -2u(x,t) +u(x-\Delta x, t)}{\Delta x^2}
\end{eqnarray}$$
$ c= k\Delta t/\Delta x^2 $²と置くことで、以下のようにまとまります。
$$\begin{eqnarray}
u(x,t+\Delta t) =(1-2c) u(x,t) +cu(x+\Delta x, t) +cu(x-\Delta x,t)
\end{eqnarray}$$
初期値と境界条件³を決めて、解いてみましょう。
$u(x,0) =0 $, $u(0,t) =u(10,t) =1$を初期条件と境界条件とします。⁴
$ k=1.0 $ として、プログラムを書いて解かせてみましょう。
$u $ は2変数関数ですが、$\Delta t$毎に$ u(x,t_i) $を計算して, $10/\Delta x $個の関数値を得ます。計算値を配列に格納し、計算を繰り返します。


import numpy as np
def calc_u(x_stride , t_stride, cur):
    """
    x_stride = Δx
    t_stride = Δt
    cur = u
    """
    n=len(cur)
    new_u = np.zeros(n)
    c= 1*t_stride/ (x_stride **2 ) 
    new_u[0]=1
    new_u[n-1] = 1
    for i in range(1, n-1):
        new_u[i] = (1-2*c)*cur[i] + c*(cur[i-1] +cur[i+1] )
    return new_u 

def diffusion(x_stride , t_stride):
    steps = round(2/t_stride ) #取りあえず2秒間分計算
    n = round(10 / x_stride )
    cur_u = np.zeros(n)
    cur_u[0] = 1 #境界条件
    cur_u[n-1] =1 #境界条件
    for i in range(0, steps):
        cur_u = calc_u(x_stride, t_stride, cur_u)
    return cur_u

calc_u が差分化した偏微分微分方程式を解く関数です。diffusion で$ u(x,t) $のいれものを作り、steps 回計算しています。この関数に適当な $ ( \Delta x , \Delta t ) $ を渡すと、長さ$10 $の筒の両端に同じ量の熱を加えたときの2秒後の熱分布(温度分布) $ u(x,2) $ を計算してくれます。
色々な$ ( \Delta x , \Delta t ) $で計算してみましょう。
初めに、$\Delta t $を固定して、$ \Delta x $を動かして計算してみます。


u_1 = diffusion(1,0.1)
u_2 = diffusion(0.9,0.1)
u_3 = diffusion(0.5, 0.1)
u_4 = diffusion(0.4, 0.1)

左側の図はなんとなく良い感じの図になっています。両端から同じ量の熱をかけたら一番冷えているのは真ん中で、温度はx軸に対して対称になっているはずなので。
しかし、右側の図は恐ろしい事になっています。物理的な条件はおかしくないので、数値計算が不安定になっているという結論になります。
次に、$ \Delta x $を固定して$ \Delta t $を変化させて計算させてみます。


v_1 = diffusion(1, 0.4)
v_2 = diffusion(1, 0.5)
v_3 = diffusion(1, 0.6)

またしても解が不安定になる条件が出てきました。どうやら(白々しいですが)数値計算の安定性には$ c$の値が関係しているようです。
なんとなく、クーラン数$c$ が0.5 以下では安定しているように見えます。少し考えてみましょう。

Von Neumannの安定性解析

全ての説明をするとそれだけで記事になるので、波動方程式とかを解いたことある人向けに説明します。⁵手持ちの本にもVon Neumannの安定性解析の議論が書いてました。⁶

熱伝導方程式は、変数分離法とか、初期値のフーリエ変換で解が具体的に求まりますが、取りあえず解の形を以下に仮定します。
$$\begin{eqnarray}
u(x,t) = g(t)\exp(i\theta (x))
\end{eqnarray}$$
この時に、$\exp(i\theta) $の部分が波の挙動を表し、$g (t) $が振幅を表します。$ g(t) $の増加率が大きい⁷ と解が振動してしまうイメージがわくと思います。
この解を差分方程式に代入します。
$$\begin{eqnarray}
u(x,t+\Delta t) &=&(1-2c) u(x,t) +cu(x+\Delta x, t) +cu(x-\Delta x,t) \\
g(t+\Delta t) &=& g(t) -2c g(t) (1- \cos(\theta) )
\end{eqnarray}$$
式変形ではオイラーの公式$ \exp(i\theta)= \cos \theta +i \sin \theta $を使います。この時、増加率$ g= g(t+\Delta t) /g(t) $の絶対値が1以下であるためには、以下の式を満たす⁸必要があります。
$$\begin{eqnarray}
0 \leq c (1-\cos \theta ) \leq 1
\end{eqnarray}$$
$1-\cos \theta = 2 \sin ^2 \frac{\theta}{2} $に注意すると、
$$\begin{eqnarray}
0 \leq c \sin ^2 \frac{\theta}{2} \leq 1 /2
\end{eqnarray}$$
が成り立たなくてはいけません。$ \sin ^2 $は1以下の正の数なので、結局
$$\begin{eqnarray}
0 \leq c \leq 1 /2
\end{eqnarray}$$
が成り立つ時、数値解が安定しているという事になります。このように、振幅の増幅率から解の安定性を見定める方法を Von Neumann の安定性解析と呼びます。
前に計算してみた$ c=0.5 $というのは本当にギリギリを責めていたわけでした。今回の場合$c= k\frac{\Delta t}{\Delta x^2} $です。$\Delta x $を細かくして計算したいなら、$ \Delta t $も上手く決めてやる必要があります。
例えば、以下のようにしなくてはいけません。

偏微分方程式を解く際に、 FTCSスキームを使うと、$c =k \frac{\Delta t}{\Delta x^2} $で解の安定性がコントロールされていることが分かりました。差分の取り方を変えると、そもそも不安定な解しか求まらないという場合もあるので注意が必要です。
例えば$ \frac{ \partial u} {\partial t} =k \frac{\partial u}{\partial x} $をFTCSスキームで解くと、不安定な解しか求まりません。興味があれば Von Neumann の安定性解析をしてみてください。
今回の記事はこれで終わりにします。この記事の続きとしては、簡単なFEM をpythonで実装するとか、流体の式を解かせてみるとかやってみようと思います。

まとめ

差分には種類があり、精度が違う
偏微分方程式も差分化出来る
適当に差分を取ると、解が不安定になる
解の不安定さはクーラン数でコントロールされている

FTCSスキーム( Forward in Time and Centered Difference in Space ) と呼ぶみたいです。
この定数をクーラン数と呼びます。
初期値は $ u(x,0) $の事です。境界条件は、方程式を$ [x_0 , x_1 ] $ で解こうとしている時は、$ u(x_0, t) , u(x_1 , t) $ の事です。初期値と言ってますが、xの関数ですし、境界条件は一般にはtの関数です。
長さ10の筒の両端から同じ熱量を流すようなイメージですね。
振動　波動とかでGoogle に聞いてみてください。
最近第3版が出たらしいです。
$ \| g(t+\Delta t)/g(t) \| > 1 $という事です。
式変形を端折りました。