5.1. 偏微分方程式の解法

5.1.1. 概観

この章では偏微分方程式（Partial Differential Equation; PDE）の数値解法を解説します。話を具体的にするために、1次元の拡散方程式を考えます。

(1)\frac{\partial u(x,t)}{\partial t} = D \frac{\partial^2 u(x,t)}{\partial x^2}

D は拡散係数です。 この方程式の解の性質を考えるために、u(x,t) が x=x_0 にピークを持っているとします。このとき、ピーク付近では (1) 式の右辺が負になるので、u(x_0,t) は時間と共に減少します。つまり、u(x,t) の凹凸がなくなる方向に時間発展します。u(x,t) が完全にフラットになると、(1) 式の右辺がゼロになり、それ以上時間変化をしなくなります。

PDEを数値的に解くには、空間座標 x を離散化して空間微分を差分で置き換えます。中心差分を適用すると (1) 式は

\frac{\partial u(x_i,t)}{\partial t} = \frac{D}{2 (\Delta x)^2} \left[ u(x_{i+1},t) - 2u(x_i,t) + u(x_{i-1},t) \right]

と表せます。 いま、[\bm{u}(t)]_i \equiv u(x_i, t) で定義されるベクトルを用いると、この方程式は

(2)\frac{d \bm{u}(t)}{dt} = f(\bm{u}(t))

の形の常微分方程式になっていることがわかります。こうすることで、時間微分に対して 常微分方程式の解法 で解説した方法がそのまま適用できます。このアプローチををFTCS法 (Forward-Time Central Space) と呼びます。

5.1.2. 不安定性

常微分方程式の解法 で、アルゴリズム特有の不安定性について解説しましたが、偏微分方程式ではこの不安定性が現れやすくなります。 これは「誤差の爆発」と言われたりします。時刻の刻み \Delta t を十分細かく取らないと、たとえプログラムが正しくても、正しい結果が得られません。一見するとプログラムのバグのように見えるので、これを知らないと、存在しないバグを探し続けることになってしまいます。

\Delta t の目安は、空間微分の階数の最大値を n として（(1) 式の拡散方程式は n=2）

(3)\frac{\Delta t}{\Delta x^n} \lesssim 1

です。分母の \Delta x^n は、空間微分の差分表示に起因します。これが分母にあるため、常微分方程式の場合よりもシビアです。 空間微分の精度を上げようと \Delta x を小さくしたら、それに応じて \Delta t も小さくしなければなりません。例えば、拡散方程式 (1) の場合、\Delta x を1/10にしたら \Delta t を1/100にしないと計算が不安定になる恐れがあります。

結果がおかしいと感じたら、\Delta t を小さしたりして、結果が変わるかを確認してください。おかしな結果の原因がバグか不安定性かを切り分けることが重要です。