概要
2次元拡散方程式
\begin{align}
\frac{\partial h}{\partial t} = \alpha \left\{\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial z^2} \right\}h
\end{align}
の安定性条件を導出します。つまり、
\begin{align}
\Delta t \leq \frac{1}{2\alpha}\frac{\Delta x^2 \Delta z^2}{\Delta x^2 + \Delta z^2}
\end{align}
を導出します（ただし、$\alpha>0$）。3次元も同様に計算できます。
ノイマンの安定性解析
空間微分を差分化すると
\begin{align}
\frac{\partial^2}{\partial x^2}h&\approx \frac{h_{i+1,j}^n-2h_{i,j}^n+h_{i-1,j}^n}{\Delta x^2}\\
\frac{\partial^2}{\partial z^2}h&\approx \frac{h_{i,j+1}^n-2h_{i,j}^n+h_{i,j-1}^n}{\Delta z^2}\\
\end{align}
と書ける。$i,j$ は格子点であり、 $n$ は時間ステップ、$\Delta x , \Delta z$は $x,z$ 方向の格子間隔である。
差分化された変数 $h_{i,j}^n$ を
\begin{align}
h_{i,j}^n =&\sum_{k_x,k_z} A^n \exp\left\{\sqrt{-1}k_x(x_i) \right\}\exp\left\{\sqrt{-1}k_z(z_j) \right\}\\
=&\sum_{k_x,k_z} A^n \exp\left\{\sqrt{-1}k_x(i\Delta x) \right\}\exp\left\{\sqrt{-1}k_z(j\Delta z) \right\}\\
\end{align}
とする（フーリエ変換）。全ての波数 $k_x,k_z$ において、 振幅 $A^n$ が次のステップの振幅 $A^{n+1}$ より小さいとき、つまり
\begin{align}
|A| \equiv \left|\frac{A^{n+1} }{A^n} \right| \leq 1
\end{align}
のとき、計算が発散せず安定する。この方法は、ノイマンの安定性解析と呼ばれる。
$h_{i,j}^n$ を具体的に代入してみると、
\begin{align}
h_{i+1,j}^n-2h_{i,j}^n+h_{i-1,j}^n &= {A^n} e^{\sqrt{-1}k_x((i+1)\Delta x) +\sqrt{-1}k_z(j\Delta z) } \\
&- 2{A^n} e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }
+{A^n} e^{\sqrt{-1}k_x((i-1)\Delta x) +\sqrt{-1}k_z(j\Delta z) }
\\
&={A^n} \left(e^{\sqrt{-1}k_x(\Delta x)}-2+e^{-\sqrt{-1}k_x(\Delta x)} \right)
e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }\\
&=(2\sqrt{-1})^2{A^n} \frac{\left(e^{\sqrt{-1}k_x(\Delta x)/2}-e^{-\sqrt{-1}k_x(\Delta x)/2} \right)^2}{(2\sqrt{-1})^2}
e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }\\
&=-4{A^n} \sin^2\left(\frac{k_x\Delta x}{2} \right)
e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }\\
h_{i,j+1}^n-2h_{i,j}^n+h_{i,j-1}^n
&=-4{A^n} \sin^2\left(\frac{k_z\Delta z}{2} \right)
e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }\\
\end{align}
となる。$\sum_{k_x,k_z} $ は、後に任意の波数 $k_x,k_z$ において成立することを示すため省いた。
陽解法
陽解法の場合、時間微分の差分化は、$\sum_{k_x,k_z} $を省いて
\begin{align}
\frac{\partial h}{\partial t} &\approx \frac{h_{i,j}^{n+1}-h_{i,j}^n}{\Delta t}\\
&=\frac{A^{n+1}-A^n}{\Delta t} e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }
\end{align}
となる。よって、拡散方程式から
\begin{align}
\frac{A^{n+1}-A^n}{\Delta t} = -4\alpha{A^n}
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\}
\end{align}
が得られる。式を変形すると、
\begin{align}
A = \frac{A^{n+1}}{A^n} = 1-4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\}
\end{align}
となる。まず $A \leq 1 $ を考えると、
\begin{align}
&1-4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\}\leq 1\\
\therefore & -4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\} \leq 0
\end{align}
$\alpha>0 ,\Delta t>0 $ であり、他は二乗の値なので、$A \leq 1 $ を満たす。
次に $-1 \leq A $ を考えると、
\begin{align}
&-1\leq 1-4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\}\\
\therefore \ \ & 4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\} \leq 2
\end{align}
さらに、任意の波数 $k_x,k_z$ において $\sin^2(k_x\Delta x/2 )\leq 1,\sin^2(k_z\Delta z/2 )\leq 1 $ が成立するので、
\begin{align}
& 4\alpha\Delta t
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\} \leq
4\alpha\Delta t
\left\{\frac{1}{\Delta x^2}+\frac{1}{\Delta z^2} \right\} \leq 2 \\
\therefore \ \ &\Delta t \leq \frac{1}{2\alpha}\frac{\Delta x^2 \Delta z^2}{\Delta x^2 + \Delta z^2}
\end{align}
$\Delta t$ の安定性条件が求まる。
よって、$|A|\leq 1 $ を満たすには、
\begin{align}
\Delta t \leq \frac{1}{2\alpha}\frac{\Delta x^2 \Delta z^2}{\Delta x^2 + \Delta z^2}
\end{align}
を満たすように、$\Delta t$ を設定しなければならない。
陰解法
陰解法の場合、時間微分の差分化は、
\begin{align}
\frac{\partial h}{\partial t} &\approx \frac{h_{i,j}^{n}-h_{i,j}^{n-1}}{\Delta t}\\
&=\frac{A^{n}-A^{n-1}}{\Delta t} e^{\sqrt{-1}k_x(i\Delta x) +\sqrt{-1}k_z(j\Delta z) }
\end{align}
となる。よって、拡散方程式から
\begin{align}
\frac{A^n-A^{n-1}}{\Delta t} = -4\alpha{A^n}
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\}
\end{align}
が得られる。式を変形すると、
\begin{align}
A = \frac{A^{n+1}}{A^n} = \left(1+4\alpha
\left\{\frac{\sin^2\left(k_x\Delta x/2 \right)}{\Delta x^2}+\frac{\sin^2\left(k_z\Delta z/2 \right)}{\Delta z^2} \right\} \right)^{-1}
\end{align}
となる。
$\alpha>0 ,\Delta t>0 $ であり、他は二乗の値なので、明らかに、$0<A \leq 1 $ である。陰解法の場合は、$|A|\leq 1 $ を満たすので、陽解法のように $\Delta t$ を設定する必要がない。
最後に
陰解法は、陽解法のように $\Delta t$ を設定する必要がない。しかし、連立方程式を解く必要がある。実装は、逐次加速緩和法（SOR法）などが使わる。
拡散方程式を差分化し、まとめると、
\begin{align}
&\frac{h_{i,j}^{n}-h_{i,j}^{n-1}}{\Delta t}= \alpha \left\{\frac{h_{i+1,j}^n-2h_{i,j}^n+h_{i-1,j}^n}{\Delta x}
+ \frac{h_{i,j+1}^n-2h_{i,j}^n+h_{i,j-1}^n}{\Delta z} \right\}\\
\therefore & \left\{1+\frac{2\alpha\Delta t}{\Delta x}+\frac{2\alpha\Delta t}{\Delta z} \right\}h_{i,j}^{n}
= h_{i,j}^{n-1}+\frac{\alpha\Delta t}{\Delta x}h_{i+1,j}^n
+\frac{\alpha\Delta t}{\Delta x}h_{i-1,j}^n
+\frac{\alpha\Delta t}{\Delta z}h_{i,j+1}^n
+\frac{\alpha\Delta t}{\Delta z}h_{i,j-1}^n
\end{align}
\begin{align}
\therefore \ \ & h_{i,j}^{n} = \left(1+\frac{2\alpha\Delta t}{\Delta x}+\frac{2\alpha\Delta t}{\Delta z} \right)^{-1}
\left\{h_{i,j}^{n-1}+\frac{\alpha\Delta t}{\Delta x}(h_{i+1,j}^n
+h_{i-1,j}^n)
+\frac{\alpha\Delta t}{\Delta z}(h_{i,j+1}^n
+h_{i,j-1}^n) \right\}
\end{align}
となる。
SOR法は、繰り返しステップを $m$ とすると
\begin{align}
\tilde{h}_{i,j}^{(m+1)n} &= \left(1+\frac{2\alpha\Delta t}{\Delta x}+\frac{2\alpha\Delta t}{\Delta z} \right)^{-1}
\left\{h_{i,j}^{n-1}+\frac{\alpha\Delta t}{\Delta x}(h_{i+1,j}^{(l)n}
+h_{i-1,j}^{(l)n})
+\frac{\alpha\Delta t}{\Delta z}(h_{i,j+1}^{(l)n}
+h_{i,j-1}^{(l)n}) \right\}\\
h_{i,j}^{(m+1)n} &= h_{i,j}^{(m)n} +\omega \left(\tilde{h}_{i,j}^{(m+1)n} -h_{i,j}^{(m)n} \right) \ \ , \ \ (0<\omega<2)
\end{align}
と更新する。 $l$ は、更新の順番による。$h_{i,j}$ より先に更新された $h$ 場合は $l=m+1$ 、他の場合は $l=m$ である。