2次元拡散方程式の安定性条件のメモ

# 概要

　2次元拡散方程式

```math
\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}
```
の安定性条件を導出します。つまり、

```math
\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次元も同様に計算できます。


# ノイマンの安定性解析

　空間微分を差分化すると

```math
\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$ を

```math
\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}$ より小さいとき、つまり

```math
\begin{align}
|A| \equiv \left|\frac{A^{n+1} }{A^n} \right| \leq 1
\end{align}
```
のとき、計算が発散せず安定する。この方法は、ノイマンの安定性解析と呼ばれる。

　$h_{i,j}^n$ を具体的に代入してみると、

```math
\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} $を省いて

```math
\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}
```
となる。よって、拡散方程式から

```math
\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}
```
が得られる。式を変形すると、

```math
\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 $ を考えると、

```math
\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 $ を考えると、

```math
\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 $ が成立するので、

```math
\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 $ を満たすには、

```math
\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$ を設定しなければならない。



# 陰解法


　陰解法の場合、時間微分の差分化は、

```math
\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}
```
となる。よって、拡散方程式から

```math
\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}
```
が得られる。式を変形すると、

```math
\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法）などが使わる。


　拡散方程式を差分化し、まとめると、

```math
\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}
```

```math
\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$ とすると
```math
\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$ である。