SIP法によるPoisson方程式の数値計算
注意
以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。
線形方程式の反復解法
次の線形方程式を考える。
\begin{align}
\label{eq:lieq}
\mathbf{A} \mathbf{u} = \mathbf{f}
\end{align}ここで、$\mathbf{A}$は係数行列、$\mathbf{u}$は変数ベクトル、$\mathbf{f}$は係数ベクトルである。もし$\mathbf{A}$が疎行列であるなら、上式の計算は反復解法によって行う必要がある。そこで、$\mathbf{A}$が次のように近似できるとする。
\begin{align}
\mathbf{M} = \mathbf{A} + \mathbf{N}
\end{align}ここで、$\mathbf{M}$は反復行列、$\mathbf{N}$は誤差行列であり、$\| \mathbf{M} \| \gg \| \mathbf{N} \|$とする。これを元の式に代入すると、反復計算の式を得る。
\begin{align}
\label{eq:ite}
\mathbf{M} \mathbf{u}^{(k+1)} = \mathbf{N} \mathbf{u}^{(k)} + \mathbf{f}
\end{align}ここで、$k$は反復計算のステップ数である。両辺から$\mathbf{M} \mathbf{u}^{(k)} $を引くと、
\begin{align}
\mathbf{M} \paren{ \mathbf{u}^{(k+1)} - \mathbf{u}^{(k)} } = \mathbf{f} - \paren{ \mathbf{M} - \mathbf{N} } \mathbf{u}^{(k)}
\end{align}ここで、$\mathbf{A} = \mathbf{M} - \mathbf{N}$より、修正量を$\mathbf{\delta}^{(k+1)} = \mathbf{u}^{(k+1)} - \mathbf{u}^{(k)}$、残差を$\mathbf{\rho}^{(k)} = \mathbf{f} - \mathbf{A} \mathbf{u}^{(k)}$とおけば、次式の修正方程式を得る。
\begin{align}
\label{eq:M}
\mathbf{M} \mathbf{\delta}^{(k+1)} = \mathbf{\rho}^{(k)}
\end{align}上式を解いた後、次式によって解を更新する。
\begin{align}
\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + \mathbf{\delta}^{(k+1)}
\end{align}反復法では式\eqref{eq:M}を解が収束するまで繰り返し計算する。そのため、$\mathbf{M}$が$\mathbf{A}$のよい近似であり、かつ式\eqref{eq:M}の求解が容易かつ効率的に行えることが望ましい。LU分解は一つの選択肢であろう。しかし、$\mathbf{A}$が疎行列である場合、LU分解によって得られる行列は元の行列の疎構造を継承せず、0だった部分に0でない値が現れる、いわゆる「fill-in」が生じる。これを防ぐため、行列の疎構造を保存しつつ近似的に分解を行う不完全LU分解（ILU分解）を適用することができる。しかしながら、標準的なILU分解は解の収束が比較的遅い。
SIP法
もしも式\eqref{eq:lieq}の係数行列$\mathbf{A}$が楕円形偏微分方程式を離散化して得られる疎行列であるなら、Stone(1968)によって提案されたILU分解法の一種であるStrongly Implicit Procedure（強陰的解法、以下、SIP法）を適用することができる。例として、次式のような直交座標系におけるPoisson方程式を考える。
\begin{align}
\pddif{u}{x} + \pddif{u}{y} = f
\end{align}ここで、$u$は求める未知関数、$f$は既知関数とする。上式を等間隔格子上で2次精度中心差分（5点スキーム）によって離散化すると、次の差分方程式を得る。
\begin{align}
\frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{\paren{\Delta x}^{2}} + \frac{u_{i,j+1} - 2 u_{i,j} + u_{i,j-1}}{\paren{\Delta y}^{2}} = f_{i,j}
\end{align}ここで、$\delta x$および$\delta y$は各方向の格子幅であり、添字$i$は$x$方向、$j$は$y$方向のインデックスであるとする。上式は次のようにまとめることができる。
\begin{align}
\label{eq:fd eq}
A_{W} u_{i-1,j} + A_{S} u_{i,j-1} + A_{P} u_{i,j} + A_{N} u_{i,j+1} + A_{E} u_{i+1,j} = f_{i,j}
\end{align}ここで、係数$A$の添え字は格子点上の現在位置（$P$）と東西南北（$EWSN$）の隣接点における量であることを表す（下図を参照）。
格子点上の変数を1次元配列に格納する順番は列優先（column-wise）とする。すなわち、2次元格子点上の$u_{i,j}$を列を固定して行方向に$u_{1,1}, u_{1,2}, u_{1,3}, \cdots, u_{2,1} , u_{2,2}, \cdots$と並べる（配列は左下隅から始める）。このような順番付けはプログラムの中の2重ループを組むときに$i$を外側、$j$を内側のループとして回すことに対応する。
こうして得られる係数行列$\mathbf{A}$は0でない要素が5本の対角線上に$WSPNE$の順に並んだ構造を持つ帯行列となる（行優先で変数を格納した場合は$SWPEN$の順番になる）。次式は係数行列のイメージである。
\begin{align}
\mathbf{A} =
\left(\begin{array}{ccccccc}
A_{P} & A_{N} & & & A_{E} & & \\
A_{S} & A_{P} & A_{N} & & & A_{E} & \\
& A_{S} & A_{P} & A_{N} & & & A_{E} \\
& & A_{S} & A_{P} & A_{N} & & \\
A_{W} & & & A_{S} & A_{P} & A_{N} & \\
& A_{W} & & & A_{S} & A_{P} & A_{N} \\
& & A_{W} & & & A_{S} & A_{P}
\end{array}\right)
\end{align}$\mathbf{A}$が帯行列であることから、$\mathbf{M}$を$\mathbf{A}$の対角行の内$WSP$成分と等しい値を持つ下三角行列$\mathbf{L}$と、$NE$成分が$\mathbf{A}$と等しく対角成分が1である（単位）上三角行列$\mathbf{U}$に分解することを考える。しかし、$\mathbf{L}$と$\mathbf{U}$の積を計算すると、$NW$と$SE$の部分に非零要素が現れる。
\begin{align}
\mathbf{L} &=
\left(\begin{array}{ccccccc}
L_{P} & & & & & & \\
L_{S} & L_{P} & & & & & \\
& L_{S} & L_{P} & & & & \\
& & L_{S} & L_{P} & & & \\
L_{W} & & & L_{S} & L_{P} & & \\
& L_{W} & & & L_{S} & L_{P} & \\
& & L_{W} & & & L_{S} & L_{P}
\end{array}\right)
\\ \nonumber
\\
\mathbf{U} &=
\left(\begin{array}{ccccccc}
1 & U_{N} & & & U_{E} & & \\
& 1 & U_{N} & & & U_{E} & \\
& & 1 & U_{N} & & & U_{E} \\
& & & 1 & U_{N} & & \\
& & & & 1 & U_{N} & \\
& & & & & 1 & U_{N} \\
& & & & & & 1
\end{array}\right)
\\ \nonumber
\\
\mathbf{LU} = \mathbf{M} &=
\left(\begin{array}{ccccccc}
M_{P} & M_{N} & & & M_{E} & & \\
M_{S} & M_{P} & M_{N} & & M_{NW} & M_{E} & \\
& M_{S} & M_{P} & M_{N} & & M_{NW} & M_{E} \\
& & M_{S} & M_{P} & M_{N} & & M_{NW} \\
M_{W} & M_{SE} & & M_{S} & M_{P} & M_{N} & \\
& M_{W} & M_{SE} & & M_{S} & M_{P} & M_{N} \\
& & M_{W} & M_{SE} & & M_{S} & M_{P}
\end{array}\right)
\end{align}この分解によって得られる方程式を式\eqref{eq:fd eq}のように表すと次式のようになる。
\begin{align}
\begin{split}
A_{W} u_{i-1,j} + N_{NW} u_{i-1,j+1} + A_{S} u_{i,j-1} + A_{P} u_{i,j} \\
+ A_{N} u_{i,j+1} + N_{SE}u_{i+1,j-1} + A_{E} u_{i+1,j} = f_{i,j}
\end{split}
\label{eq:mod fd eq}
\end{align}上式は式\eqref{eq:ite}に対応するが、これをそのまま計算しても収束が遅い。これは誤差項$\mathbf{N}$が最小化されていないためである。
Stoneは元々の楕円形方程式が空間的に滑らかに分布する解を持つことを利用して、誤差項を次のように近似すれば収束を加速できることを示した。
\begin{align}
u_{i-1,j+1} &\approx - u_{i,j} + u_{i,j+1} + u_{i-1,j} \\
u_{i+1,j-1} &\approx - u_{i,j} + u_{i+1,j} + u_{i,j-1}
\end{align}この近似は当然、一般の行列に対しては成り立たないため、SIP法を汎用の線形計算法として利用することはできないことに注意すべきである（別の近似式を用いることもできる）。上式を元の式\eqref{eq:mod fd eq}から差し引くと、
\begin{align}
\begin{split}
&A_{W} u_{i-1,j} + A_{S} u_{i,j-1} + A_{P} u_{i,j} + A_{N} u_{i,j+1} + A_{E} u_{i+1,j} \\
& + N_{NW}\bracket{ u_{i-1,j+1} - \alpha \paren{ - u_{i,j} + u_{i,j+1} + u_{i-1,j} } } \\
& + N_{SE}\bracket{ u_{i+1,j-1} - \alpha \paren{ - u_{i,j} + u_{i+1,j} + u_{i,j-1} } } = f_{i,j}
\end{split}
\end{align}あるいは
\begin{align}
\begin{split}
&\paren{ A_{W} - \alpha N_{NW}} u_{i-1,j} + \paren{ A_{S} - \alpha N_{SE} } u_{i,j-1} \\
&+ \paren{ A_{P} + \alpha N_{NW} + \alpha N_{SE} } u_{i,j} \\
&+ \paren{ A_{N} - \alpha N_{NW} } u_{i,j+1} + \paren{ A_{E} - \alpha N_{SE} } u_{i+1,j} \\
&+ N_{NW} u_{i-1,j+1} + N_{SE}u_{i+1,j-1} = f_{i,j}
\end{split}
\end{align}ここで、$\alpha$は収束を制御するパラメータであり、$0 < \alpha < 1$である。$\alpha$が$1$に近いほど収束を速めることができるが、具体的な最適値は問題による。この工夫によって反復誤差の影響が反復行列のすべての要素に反映されることにより、標準的なILU分解よりも収束性を改善することができる。
ILU分解と修正方程式の反復アルゴリズム
具体的な計算は以下のように行う。
1. 初期値$\mathbf{u}^{(0)}$から初期残差$\mathbf{\rho}^{(0)} = \mathbf{f} - \mathbf{A} \mathbf{u}^{(0)}$を計算する。$0 < \alpha < 1$の範囲でパラメータを設定する。
2. 反復計算の前に$\mathbf{L}$と$\mathbf{U}$を一度だけ計算する。これらは反復の度に更新する必要はない。
\begin{align}
L_{W}^{l} &= A_{W}^{l} / \paren{ 1 + \alpha U_{N}^{l - N_{j}} } \\
L_{S}^{l} &= A_{S}^{l} / \paren{ 1 + \alpha U_{E}^{l - 1} } \\
L_{P}^{l} &= A_{P}^{l} + \alpha \paren{ L_{W}^{l} U_{N}^{l - N_{j}} + L_{S}^{l} U_{E}^{l - 1} } - L_{W}^{l} U_{E}^{l - N_{j}} - L_{S}^{l} U_{N}^{l - 1} \\
U_{N}^{l} &= \paren{ A_{N}^{l} - \alpha L_{W}^{l} U_{N}^{l - N_{j}} } / L_{P}^{l} \\
U_{E}^{l} &= \paren{ A_{E}^{l} - \alpha L_{S}^{l} U_{E}^{l - 1} } / L_{P}^{l}
\end{align}ここで、$l$は行列ないしベクトルのインデックスである。ただし、上式の計算は上からこの順に実行しなければならない。また、境界付近では境界より外側の値は0として計算する。
3. 各反復において、残差と修正量の関係を表す修正方程式を解く。$\mathbf{M} = \mathbf{LU}$より、
\begin{align}
\mathbf{LU} \mathbf{\delta}^{(k+1)} = \mathbf{\rho}^{(k)}
\end{align}上式はLU分解によって解くことができる。まず、後退代入によって下式を解く。
\begin{align}
\mathbf{U} \mathbf{\delta}^{(k+1)} = \mathbf{L}^{-1} \mathbf{\rho}^{(k)} = \mathbf{R}^{(k)}
\end{align}$\mathbf{R}^{(k)}$は$l$の昇順にしたがって次のように計算できる。
\begin{align}
R^{l} = \paren{ \rho^{l} - L_{S}^{l} R^{l - 1} - L_{W}^{l} R^{l - N_{j}} } / L_{P}^{l}
\end{align}$\mathbf{R}^{(k)}$が求まれば、$\mathbf{\delta}^{(k+1)}$は前進代入によって$l$の降順で解くことができる。
\begin{align}
\delta^{l} = R^{l} - U_{N}^{l} \delta^{l+1} - U_{E}^{l} \delta^{l+N_{j}}
\end{align}上式を解いた後、次式によって解を更新する。
\begin{align}
\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + \mathbf{\delta}^{(k+1)}
\end{align}プログラムで計算するときには、行列要素はすべて境界の外側の位置まで含めて配列を確保する。たとえば、LU分解の要素の計算は次のように書ける。
int NN=(NX+2)*(NY+2); int Nj=NY+2; vector<double> AW(NN,0.0); vector<double> LW(NN,0.0); vector<double> UN(NN,0.0); for(int i=1;i<=NX;i++){ for(int j=1;j<=NY;j++){ int l = i * Nj + j; LW[l] = AW[l] / ( 1.0 + alpha * UN[l - Nj] ); } }
例題
テスト計算として次の極座標系におけるPoisson方程式を解いた。
\begin{align}
\pddif{u}{r} + \frac{1}{r} \pdif{u}{r} + \frac{1}{r^{2}} \pddif{u}{\theta} = - 4 \sin\paren{\pi r} \sin\paren{2 \theta}
\end{align}解析条件は以下のとおり。
- 解析領域：穴のある円形（$d_{0} \leq r \leq 1$)、$d_{0}=0.1$：穴の半径
- 境界条件：上下（内外周）で$u = 0$、周期境界条件$u(0)=u(2\pi)$
- 分割数は$64\times64$、$128\times128$の2種類
- 収束判定：初期残差$\mathbf{\rho}^{(0)}$で規格化した残差の$L^{1}$ノルム（各要素の絶対値和）で評価
\[ \displaystyle \frac{\| \mathbf{\rho}^{(k)} \|_{1}}{\mathbf{\rho}^{(0)}} < 10^{-4} \]
- 反復回数の上限は$10^{4}$回
- 比較としてGuss-SeidelS法でも計算
計算結果のカラーマップ図を以下に示す。
収束履歴は下図のようになった。残差の値は相対値ではなく各ステップ毎の減少量としてプロットした方がわかりやすかったように思うが、面倒なのでこれでよしとする。
収束に要した反復回数は下表のようになった。
| Grid | GS | SIP |
|---|---|---|
| $64 \times 64$ | 1895 | 58 |
| $128 \times 128$ | 7509 | 190 |
反復回数で見ると、GS法に比べてSIP法の方が10倍以上速い。とはいえ、計算負荷や使用メモリ量はSIP法の方が大きいため、状況によってはGS法（ないしそれを加速したSOR法）で十分という場合も考えられる。問題の性質や規模に応じた解法の選択が重要である。
おわりに
日本語で書かれたSIP法の解説がネットにはあまりないようなので備忘録として残すことにした。基本的には参考文献(1)ないし(2)の書き写しである。サンプルプログラムが手に入るURLなども載っているので、詳しくは(1)を参照されたい（10年以上前の本だがプログラムはまだ入手できる）。LU分解部分の具体的な計算は(2)にある。5点スキーム以外の問題の行列要素を求めるときに参考になるだろう（交差微分などを含む一般の楕円形方程式には9点スキームが必要）。また、NAG計算ライブラリにはSIP法が収録されているようなので、それを利用できる立場の人は自分でプログラムを作る必要はないかもしれない（もっというなら普通の人は共役勾配法系のライブラリを使うんじゃなかろうか……）。
参考文献
(1) J. H. Ferziger, M. Peric, 小林敏雄他訳，コンピュータによる流体力学，シュプリンガー・フェアラーク，2003年
(2) http://www.engr.uky.edu/~acfd/me690-lctr-nts.pdf
(3) http://www.nag.co.uk/numeric/fl/manual/pdf/D03/d03ebf.pdf
(4) Stone, H. L. "Iterative Solution of Implicit Approximations of Multidimensional Partial Differential Equations". SIAM Journal of Numerical Analysis 5 (3): 530–538., (1968)
(5) G. E. Schneider and M. Zedan, "A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems", Numer. Heat Transf. 4, 1–19, (1981).
