概要

差分法で離散化した際に現れる連立一次方程式

\boldsymbol{A}\boldsymbol{f} = \boldsymbol{g}

をSOR法で解くときに，SOR法の漸化式は

\varDelta \boldsymbol{f} = \frac{\omega}{\boldsymbol{D_A}}\left(\boldsymbol{g}-\boldsymbol{A}\hat{\boldsymbol{f}}\right)

と記述できる．$\varDelta \boldsymbol{f}$はある反復における修正量，$\omega$は加速係数，$\boldsymbol{D_A}$は$\boldsymbol{A}$の対角成分である．$\hat{\boldsymbol{f}}$は，反復過程における更新前後の値が混在した未知数である．

差分法の場合，支配方程式を離散化した式から容易に漸化式を導出できるが，導出した漸化式はSOR法以外では利用しない．この記述では，支配方程式を離散化した式をそのまま用いて漸化式を構築できる．利点は覚えやすいということだけであり，計算量の点で有利であるとか，桁落ちに強いとかはない．

はじめに

差分法による非圧縮性流れの数値計算を学習する際，差分法による微分の近似を学ぶと，そこから自然な形でJacobi法系統の方法（Jacobi法，Gauss-Seidel法，SOR法）による連立方程式の求解法を学ぶことができる．その求解法は，係数行列を作る必要がなく，漸化式を実装すればよいため非常に簡単である．一方で，漸化式は解きたい方程式に固有の形になり，また非圧縮性流れの計算を行う過程でその式を流用することもない．

そのため，漸化式は覚えにくく，実装の度に導出し直す事になる．漸化式をよくよく眺めてみると，非常に覚えやすい形に整理できることがわかった．備忘録としてSOR法を例にその考察を示す．

SOR法

SOR法（Successive Over-Relaxation method, 逐次過緩和法）は，連立一次方程式を解く手法の一つである．定常反復法の一種であり，Gauss-Seidel法に加速係数（緩和係数）を導入した方法として知られている．

差分法におけるSOR法の実装例

漸化式の導出

1次元Poisson方程式

\begin{equation}
    \nabla^2 f = g \qquad (1)
\end{equation}

を用いて，差分法による実装の例を示す．

$0\le x \le L_x$の範囲に，等間隔（$\varDelta x$）に$N_x$個の点を配置し，その点上で$f, g$の値が定義されるとする．$N_x$個の各点に$1$から順番に番号を付け，番号$i$の点で定義される値を，下付添字$i$を用いて表すことにする．

Poisson方程式を2次精度中心差分で離散化すると，次式が得られる．

\begin{equation}
    \frac{f_{i+1}-2f_i+f_{i-1}}{\varDelta x^2} = g_i \qquad (2)
\end{equation}

このとき，既知の$f$から$g$を計算するのは至極簡単であるが，既知の$g$から$f$を求めるには，次のような連立方程式を解く必要がある．

\begin{equation}
\begin{pmatrix}
1&0&0&0&0\\
\frac{1}{\varDelta x^2}&\frac{-2}{\varDelta x^2}&\frac{1}{\varDelta x^2}&0&0\\
0&\frac{1}{\varDelta x^2}&\frac{-2}{\varDelta x^2}&\frac{1}{\varDelta x^2}&0\\
0&0&\frac{1}{\varDelta x^2}&\frac{-2}{\varDelta x^2}&\frac{1}{\varDelta x^2}\\
0&0&0&0&1\\
\end{pmatrix}
\begin{pmatrix}
f_1\\
f_2\\
f_3\\
f_4\\
f_5
\end{pmatrix}
=
\begin{pmatrix}
g_1\\
g_2\\
g_3\\
g_4\\
g_5
\end{pmatrix} \qquad (3)
\end{equation}

上式において，$g_1$と$g_5$は$f_1$と$f_5$の境界条件である．

反復法では，上式のような連立方程式を作らず，近似値$\tilde{f}$を計算してその値が真の値から遠ければ何らかの方法で$\tilde{f}$を更新する．更新された$\tilde{f}$が真の値に十分近くなるまで計算する手順を繰り返す．

Gauss-Seidel法では，離散化されたPoisson方程式$(2)$を変形し，$f_i$を

\begin{equation}
f_i = \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \left(\frac{f_{i+1}+f_{i-1}}{\varDelta x^2}\right)\right] \qquad (4)
\end{equation}

から計算する．繰り返しの階数を右肩に括弧でくくって表示すると，次のような手順で$f_2,f_3,f_4$を更新する．（式$(3)$では$f_1,f_5$は境界条件によって値が決まるため，更新の必要がない．）

\left\{
    \begin{align}
    f_2^{(1)} &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_2 - \left(\frac{f_{1}+f_{3}^{(0)}}{\varDelta x^2}\right)\right]\\
    f_3^{(1)} &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_3 - \left(\frac{f_{2}^{(1)}+f_{4}^{(0)}}{\varDelta x^2}\right)\right]\\
    f_4^{(1)} &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_4 - \left(\frac{f_{3}^{(1)}+f_{5}}{\varDelta x^2}\right)\right]\\
    \end{align}
\right.

$f_3$を更新する式をみると，更新された$f_2$と更新前の$f_4$の値が用いられていることがわかる．Gauss-Seidel法では，ある反復で更新された値をその反復の中で利用する．式$(2)$のように，ある点$i$に対する式として書くと，ある反復$(m)$において更新される$f_i$の値は

\begin{equation}
    f_i^{(m+1)} = \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}+f_{i+1}^{(m)}}{\varDelta x^2}\right)\right]
\end{equation}

と書ける．

SOR法は，係数を導入してGauss-Seidel法の各反復における修正量を大きくし，収束に要する繰返し回数を削減する．Gauss-Seidel法における各反復の更新は

\begin{align}
    d f_i &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}+f_{i+1}^{(m)}}{\varDelta x^2}\right)\right] - f_i^{(m)}\\
    f_i^{(m+1)} &= f_i^{(m)} + d f_i
\end{align}

と表すことができる．SOR法では，$d f$を計算する式はGauss-Seidel法と同じであるが，$f$の更新は加速係数$\omega$を導入して

\begin{align}
    d f_i &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}+f_{i+1}^{(m)}}{\varDelta x^2}\right)\right] - f_i^{(m)}\\
    f_i^{(m+1)} &= f_i^{(m)} + \omega d f_i
\end{align}

とする．加速係数もまとめて更新量$\varDelta f_i$として書けば

\varDelta f_i := \omega df_i

となる．

いつ反復を止めるか

反復法では，近似値が真の値に近くなるまで反復を繰り返す．近似解と真の値との近さの指標として，残差と変化量の二つが知られている．

残差

連立方程式

\boldsymbol{A}\boldsymbol{f} = \boldsymbol{g}

を

\boldsymbol{r} = \boldsymbol{g} - \boldsymbol{A}\boldsymbol{f}

と変形し，$\boldsymbol{f}$に$f^{(m+1)}$を代入して計算する．ここで，本記事の場合，$\boldsymbol{A}\boldsymbol{f}$は式$(2)$の左辺と一致する．$\boldsymbol{r}$を残差とよび，残差が予め決めた値よりも小さくなるまで計算を繰り返すことで，真の値に近い近似値を得る．反復を止めるかどうかの指標としては，残差の$l2$ノルム$|\boldsymbol{r}|$や，その相対値$|\boldsymbol{r}|/|\boldsymbol{g}|$が利用される．

変化量

変化量は，$(m)$回目の繰り返しで得られた$\boldsymbol{f}^{(m)}$と$(m+1)$回目の繰り返しで得られた$\boldsymbol{f}^{(m+1)}$を計算する．変化量が小さいと，$f$がこれ以上更新されない＝真の解に近いと判断する．

変化量が小さいから真の値に近いとする根拠は弱いが，反復を止めるかどうかの指標として相対変化量$|\boldsymbol{f}^{(m+1)}-\boldsymbol{f}^{(m)}|/|\boldsymbol{f}^{(m+1)}|$を用いる場合は，$|\boldsymbol{f}^{(m+1)}|$が真の解になっていれば相対誤差とみなせる．

Gauss-Seidel法やSOR法の場合，計算過程で$\boldsymbol{f}^{(m)}$と$\boldsymbol{f}^{(m+1)}$を計算するため，変化量を用いる方が実装は簡単である．残差を用いる場合は，残差を計算するための処理（2階微分など）を新たに実装する必要がある（ようにみえる）．

漸化式の整理

式$(11)$をみると，角括弧外側の係数は，係数行列の対角項である．

\begin{align}
    d f_i = \frac{1}{D_i}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}+f_{i+1}^{(m)}}{\varDelta x^2}\right)\right] - f_i^{(m)}
\end{align}

また，$f_i^{(m)}$を丸括弧内に入れると

\begin{align}
    d f_i &= \frac{1}{D_i}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}+f_{i+1}^{(m)}}{\varDelta x^2}  + D_i f_i^{(m)}\right)\right]\\
    &= \frac{1}{D_i}\left[g_i - \left(\frac{f_{i-1}^{(m+1)}-2f_i^{(m)}+f_{i+1}^{(m)}}{\varDelta x^2}\right)\right] \qquad (5)
\end{align}

となり，式$(5)$の角括弧内は残差の定義$\boldsymbol{g} - \boldsymbol{A}\boldsymbol{f}$そのものであるように見える．ただし，未知数$\boldsymbol{f}$には更新された値と未更新の値が混在している．（区別のためにハットを付けて$\hat{\boldsymbol{f}}$と書くことにする）

差分法によって離散化した際に現れる連立方程式をSOR法で計算する場合，係数行列式を作る必要はないが，導出される漸化式が覚えにくいという問題があった．漸化式をよくよく見てみれば，解きたい式を離散化した式の形と，着目している離散点$i$における未知数$f_i$の係数のみで漸化式は簡単に作れることがわかる．$(m)$回目の値を使うか，$(m+1)$回目の値を使うかは，実装時に離散点を計算していく方向に応じて決めればよい．例えば，$i$を1から走査する場合は，

\hat{\boldsymbol{f}} =
\begin{pmatrix}
f_1^{(m+1)}\\
\vdots\\
f_{i-1}^{(m+1)}\\
f_i^{(m)}\\
f_{i+1}^{(m)}\\
\vdots\\
f_{N_x}^{(m)}\\
\end{pmatrix}

となる．

反復を止める指標の一つである相対変化量について見直してみると，相対変化量の分子は

\|\boldsymbol{f}^{(m+1)}-\boldsymbol{f}^{(m)}\|=\|\omega d\boldsymbol{f}\|=\left\|\frac{\omega}{\boldsymbol{D_A}}\left(\boldsymbol{g}-\boldsymbol{A}\boldsymbol{f}\right)\right\|

なので，係数$\frac{\omega}{\boldsymbol{D_A}}$がかかっていることや規格化に使う値（$\boldsymbol{g}$か$\boldsymbol{f}^{(m+1)}$か）が異なるものの，結局は残差を計算していることがわかる．そのため，反復を止める指標として相対変化量を用いたとしても，致命的な問題とはならない．

まとめ

SOR法の漸化式は，以下の手順で簡単に導出できる．
1 解きたい式を差分法で離散化する．

\begin{equation}
    \frac{f_{i+1}-2f_i+f_{i-1}}{\varDelta x^2} = g_i
\end{equation}

2 離散化した式を右辺に集め，添字$i$の未知数にかかる係数で式全体を除す．

\begin{equation}
    \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \frac{f_{i+1}-2f_i+f_{i-1}}{\varDelta x^2}\right]
\end{equation}

3 どの反復における値かを決定する．（配列を1から走査する場合は，$i-1$が$(m+1)$回目の反復の値となる）

\begin{equation}
    \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \frac{f_{i+1}^{(m)}-2f_i^{(m)}+f_{i-1}^{(m+1)}}{\varDelta x^2}\right]
\end{equation}

4 得られた式を，SOR法の反復における未知数の更新量$df_i$とする．

\begin{align}
    d f_i &= \frac{1}{\frac{-2}{\varDelta x^2}}\left[g_i - \frac{f_{i+1}^{(m)}-2f_i^{(m)}+f_{i-1}^{(m+1)}}{\varDelta x^2}\right]\\
    f_i^{(m+1)} &= f_i^{(m)} + \omega d f_i
\end{align}

1次元のPoisson方程式だけでなく，他の式でも同様に導出できる．過去に書いた記事（差分法によるPoisson方程式の離散化とJacobi法で用いる漸化式）で取り扱った2次元準調和方程式

\nabla\cdot\left(\frac{1}{\rho}\nabla p\right) = \frac{1}{\Delta t}\theta^*

をSOR法で解く漸化式を，同様の手順計算する．

まず，準調和方程式を2次精度中心差分で離散化する．

\frac{1}{\Delta x}\left(\frac{2}{\rho_{i+1,j}+\rho_{i,j}}\frac{p_{i+1,j}-p_{i,j}}{\Delta x}-\frac{2}{\rho_{i,j}+\rho_{i-1,j}}\frac{p_{i,j}-p_{i-1,j}}{\Delta x}\right)+\frac{1}{\Delta y}\left(\frac{2}{\rho_{i,j+1}+\rho_{i,j}}\frac{p_{i,j+1}-p_{i,j}}{\Delta y}-\frac{2}{\rho_{i,j}+\rho_{i,j-1}}\frac{p_{i,j}-p_{i,j-1}}{\Delta y}\right) = \frac{1}{\Delta t}\theta^*_{i,j}

離散化した式を右辺に集め，添字$i$の未知数にかかる係数で式全体を除す．

\frac{\frac{1}{\Delta t}\theta^*_{i,j} -\left\{ \frac{1}{\Delta x}\left(\frac{2}{\rho_{i+1,j}+\rho_{i,j}}\frac{p_{i+1,j}-p_{i,j}}{\Delta x}-\frac{2}{\rho_{i,j}+\rho_{i-1,j}}\frac{p_{i,j}-p_{i-1,j}}{\Delta x}\right)+\frac{1}{\Delta y}\left(\frac{2}{\rho_{i,j+1}+\rho_{i,j}}\frac{p_{i,j+1}-p_{i,j}}{\Delta y}-\frac{2}{\rho_{i,j}+\rho_{i,j-1}}\frac{p_{i,j}-p_{i,j-1}}{\Delta y}\right)\right\}}{-\frac{1}{\Delta x^2}\left(\frac{2}{\rho_{i+1,j}+\rho_{i,j}}+\frac{2}{\rho_{i,j}+\rho_{i-1,j}}\right)-\frac{1}{\Delta y^2}\left(\frac{2}{\rho_{i,j+1}+\rho_{i,j}}+\frac{2}{\rho_{i,j}+\rho_{i,j-1}}\right)}

配列を$i=1, j=1$から走査することとし，得られた式をSOR法の反復における未知数の更新量$dp_{i,j}$とする．

\begin{align}
dp_{i,j} &= \frac{\frac{1}{\Delta t}\theta^*_{i,j} -\left\{ \frac{1}{\Delta x}\left(\frac{2}{\rho_{i+1,j}+\rho_{i,j}}\frac{p_{i+1,j}^{(m)}-p_{i,j}^{(m)}}{\Delta x}-\frac{2}{\rho_{i,j}+\rho_{i-1,j}}\frac{p_{i,j}^{(m)}-p_{i-1,j}^{(m+1)}}{\Delta x}\right)+\frac{1}{\Delta y}\left(\frac{2}{\rho_{i,j+1}+\rho_{i,j}}\frac{p_{i,j+1}^{(m)}-p_{i,j}^{(m)}}{\Delta y}-\frac{2}{\rho_{i,j}+\rho_{i,j-1}}\frac{p_{i,j}^{(m)}-p_{i,j-1}^{(m+1)}}{\Delta y}\right)\right\}}{-\frac{1}{\Delta x^2}\left(\frac{2}{\rho_{i+1,j}+\rho_{i,j}}+\frac{2}{\rho_{i,j}+\rho_{i-1,j}}\right)-\frac{1}{\Delta y^2}\left(\frac{2}{\rho_{i,j+1}+\rho_{i,j}}+\frac{2}{\rho_{i,j}+\rho_{i,j-1}}\right)}\\
p_{i,j}^{(m+1)} &= p_{i,j}^{(m)} + \omega d p_{i,j}
\end{align}