2 陰公式(implicit formula)

数値積分を実行する為に拡散方程式を差分化する必要がある。 差分化にはいろいろな方法があるが、 ここでは数値的な安定性を考えて、 陰公式(implicit formula) を用いることにする。 空間座標 $x$ について $J$ 等分する。 従って、隣り合う二点間の間隔は $\Delta x= 1/J$ で与えられる。 積分するときの時間間隔を $\Delta t$ として、 時刻 $t{}^n = t{}^0 + n\Delta t$ に於ける場所 $x_j = x_0 + j\Delta x$ での解を

$$U_j^n =u\left(x_j=x_0+j\Delta x,t{}^n=t{}^0 +n\Delta t\right)$$ (5)

と表すことにする。ここでは $x_0=0 t{}^0=0$ とする。(1)の差分化の一案として

$$\frac{U_j^{n+1} -U_j^n}{\Delta t} = p\frac{1}{\Delta x} \left( \frac{U_{j+1}^{n+1} -U_j^{n+1}}{\Del... ...{U_{j+1}^{n} -U_j^{n}}{\Delta x} -\frac{U_j^{n} -U_{j-1}^{n}}{\Delta x} \right)$$

$$= p\frac{U_{j+1}^{n+1} -2U_{j}^{n+1} +U_{j-1}^{n+1}}{(\Delta x)^2} +(1-p) \frac{U_{j+1}^{n} -2U_{j}^{n} +U_{j-1}^{n}}{(\Delta x)^2}$$ (6)

と書く。ここで $0\leq p\leq 1$ とする。 (6)は

$$\alpha =\frac{\Delta t}{(\Delta x)^2}$$ (7)

を用いて、

$$-\alpha p U_{j+1}^{n+1} +(1+2\alpha p)U_{j}^{n+1} -\alpha p U_{j-... ...V_j^n \equiv \alpha q U_{j+1}^{n} +(1-2\alpha q)U_{j}^{n} +\alpha q U_{j-1}^{n}$$ (8)

と変形できる。 ここで $q=1-p$ である。 これを連立一次方程式と考えて、大きな行列を使って表せば、 (8) の左辺は

$$\begin{pmatrix}1+2\alpha p & -\alpha p & & & & 0 -\alpha p & ... ... V_2^{n} \vdots \vdots V_{J-2}^{n} V_{J-1}^{n} \end{pmatrix}$$ (9)

となり、右辺は

$$\begin{pmatrix}V_1^{n} V_2^{n} \vdots \vdots V_{J... ... U_2^{n} \vdots \vdots U_{J-2}^{n} U_{J-1}^{n} \end{pmatrix}$$ (10)

と書くことができる。 $V_j^n$ は $U_j^n$ が与えられれば簡単に計算できる。 (9) と (10) は、 形式的に $A \vec{U}^{n+1} =\vec{V}^{n}$ と $\vec{V}^{n} = B \vec{U}^n$ と書くことができる。 $\vec{V}^n$ は $\vec{U}^n$ が与えられれば簡単に求められるので、 原理的には連立一次方程式 $A \vec{U}^{n+1} =\vec{V}^n$ を、 行列 $A$ の逆行列を用いて $\vec{U}^{n+1} = A^{-1} \vec{V}^n=A^{-1} B \vec{U}^n$ として解けばよいことになる。 これはつまり初期条件 $\vec{U}^n$ が与えられれば $\vec{U}^{n+1}$ が求められることを意味している。