メニュー

技術情報 Techinicalinfo

  1. ホーム
  2. 技術情報
  3. 【技術情報】有限要素法入門
  4. 6.4 有限要素法のプログラム(その6)

【技術情報】有限要素法入門

6.4 有限要素法のプログラム(その6)


この節では連立方程式の解法について述べます。有限要素法における連立方程式では係数行列の要素の多くがゼロとなっており、この性質を利用した解法が必要です。
また有限要素法ではこの係数行列が対称行列となっています。したがって係数行列の非対角項については半分だけ格納する配列を用意すればよいことになりメモリーの節約ができます。

ここで次の連立方程式を考えます。
\begin{equation}
A\boldsymbol{x}=\boldsymbol{b} \tag*{$(6.4-1)$}
\end{equation}
このときこの方程式を解くことは次の関数の最小値を求めることと同じです。
\begin{equation}
\begin{split}
F(\boldsymbol{x})&=\frac{1}{2}(\boldsymbol{x},A\boldsymbol{x})-(\boldsymbol{x},\boldsymbol{b}) \\
&=\sum_{ij}\frac{1}{2}x_iA_{ij}x_j-\sum_ix_ib_i
\end{split} \tag*{$(6.4-2)$}
\end{equation}
なぜなら、この微分をとると、
\begin{equation}
\begin{split}
\frac{\partial F(\boldsymbol{x})}{\partial x_i}&=\frac{1}{2}\sum_jA_{ij}x_j+\frac{1}{2}\sum_jx_jA_{ij}-b_i \\
&=\sum_jA_{ij}x_j-b_i
\end{split} \tag*{$(6.4-3)$}
\end{equation}
となりますが、最小値ではこれがゼロになるからです。
ここで $\boldsymbol{x}$ の適当な初期値 $\boldsymbol{x}^{(0)}$ から出発してこの関数 $F$ の最小値にたどり着くような探索を考えます。
まず探索方向 $\boldsymbol{p}$ を決めて、
\begin{equation}
\boldsymbol{x}^{(1)}=\boldsymbol{x}^{(0)}+\alpha^{(0)}\boldsymbol{p}^{(0)} \tag*{$(6.4-4)$}
\end{equation}
のときに関数 $F$ の値が最小になるように $\alpha^{(0)}$ を求めます。
\begin{equation}
\frac{\partial}{\partial\alpha^{(0)}}F(\boldsymbol{x}^{(0)}+\alpha^{(0)}\boldsymbol{p}^{(0)})=0 \notag
\end{equation}
この式の左辺を変形すると次のようになります。
\begin{equation}
\begin{split}
&\frac{\partial}{\partial\alpha^{(0)}}\bigl\{\sum_{ij}\frac{1}{2}(x_i^{(0)}+\alpha^{(0)}p_i^{(0)})A_{ij}(x_j^{(0)}+\alpha^{(0)}p_j^{(0)})
-\sum_i(x_i^{(0)}+\alpha^{(0)}p_i^{(0)})b_i\bigr\} \\
&=\sum_{ij}\frac{1}{2}p_i^{(0)}A_{ij}(x_j^{(0)}+\alpha^{(0)}p_j^{(0)})+\sum_{ij}\frac{1}{2}(x_i^{(0)}+\alpha^{(0)}p_i^{(0)})A_{ij}p_j^{(0)}
-\sum_ip_i^{(0)}b_i \\
&=(\boldsymbol{p}^{(0)},A\boldsymbol{x}^{(0)})+(\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})\alpha^{(0)}-(\boldsymbol{p}^{(0)},\boldsymbol{b})
\end{split} \notag
\end{equation}
これより $\alpha^{(0)}$ は次のようになります。
\begin{equation}
\alpha^{(0)}=\frac{(\boldsymbol{p}^{(0)},\boldsymbol{b}-A\boldsymbol{x}^{(0)})}{(\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})} \tag*{$(6.4-5)$}
\end{equation}
ここに、
\begin{equation}
\boldsymbol{r}^{(0)}\equiv\boldsymbol{b}-A\boldsymbol{x}^{(0)} \tag*{$(6.4-6)$}
\end{equation}
は残差ベクトルであり、(6.4-3)式より$F$の勾配ベクトルと逆向きのベクトルです。
これより(6.4-5)式は、
\begin{equation}
\alpha^{(0)}=\frac{(\boldsymbol{p}^{(0)},\boldsymbol{r}^{(0)})}{(\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})} \tag*{$(6.4-7)$}
\end{equation}
となります。さてこの探索方向のベクトルをどう決めるかですが、最初の探索方向としてこの残差ベクトルをとります。
\begin{equation}
\boldsymbol{p}^{(0)}=\boldsymbol{r}^{(0)} \tag*{$(6.4-8)$}
\end{equation}
これより、
\begin{equation}
\boldsymbol{x}^{(1)}=\boldsymbol{x}^{(0)}+\alpha^{(0)}\boldsymbol{p}^{(0)} \tag*{$(6.4-9)$}
\end{equation}
となります。これに対する残差ベクトルは次のように計算できます。
\begin{equation}
\begin{split}
\boldsymbol{r}^{(1)}&=\boldsymbol{b}-A\boldsymbol{x}^{(1)}=\boldsymbol{b}-A\boldsymbol{x}^{(0)}+A(\boldsymbol{x}^{(0)}-\boldsymbol{x}^{(1)}) \\
&=\boldsymbol{r}^{(0)}-\alpha^{(0)}A\boldsymbol{p}^{(0)}
\end{split} \tag*{$(6.4-10)$}
\end{equation}
次の探索方向のベクトルは前の方向と $A$ に関して共役方向となるように選びます。すなわち、
\begin{equation}
(\boldsymbol{p}^{(1)},A\boldsymbol{p}^{(0)})=0  \tag*{$(6.4-11)$}
\end{equation}
です。ここで、
\begin{equation}
\boldsymbol{p}^{(1)}=\boldsymbol{r}^{(1)}+\beta^{(0)}\boldsymbol{p}^{(0)} \tag*{$(6.4-12)$}
\end{equation}
とおけば上の式に代入して次のようになります。
\begin{equation}
(\boldsymbol{r}^{(1)}+\beta^{(0)}\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})=(\boldsymbol{r}^{(1)},A\boldsymbol{p}^{(0)})+\beta^{(0)}(\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})=0 \notag
\end{equation}
これより、
\begin{equation}
\beta^{(0)}=-\frac{(\boldsymbol{r}^{(1)},A\boldsymbol{p}^{(0)})}{(\boldsymbol{p}^{(0)},A\boldsymbol{p}^{(0)})} \tag*{$(6.4-13)$}
\end{equation}
となります。
これより次の $\alpha^{(1)}$ は次のようになります。
\begin{equation}
\alpha^{(1)}=\frac{(\boldsymbol{p}^{(1)},\boldsymbol{r}^{(1)})}{(\boldsymbol{p}^{(1)},A\boldsymbol{p}^{(1)})} \notag
\end{equation}
これより次の解が求まります。
\begin{equation}
\boldsymbol{x}^{(2)}=\boldsymbol{x}^{(1)}+\alpha^{(1)}\boldsymbol{p}^{(1)} \notag
\end{equation}
次の残差ベクトルは、
\begin{equation}
\boldsymbol{r}^{(2)}=\boldsymbol{r}^{(1)}-\alpha^{(1)}A\boldsymbol{p}^{(1)} \notag
\end{equation}
となり、
\begin{equation}
\beta^{(1)}=-\frac{(\boldsymbol{r}^{(2)},A\boldsymbol{p}^{(1)})}{(\boldsymbol{p}^{(1)},A\boldsymbol{p}^{(1)})} \notag
\end{equation}
です。これから次の探索方向ベクトルは次のように求まります。
\begin{equation}
\boldsymbol{p}^{(2)}=\boldsymbol{r}^{(2)}+\beta^{(1)}\boldsymbol{p}^{(1)} \notag
\end{equation}
後はこれの繰り返しなのでまとめると次のようになります。

まず $k$ 番目までの探索方向のベクトル $\boldsymbol{p}^{(k)}$ と残差ベクトル $\boldsymbol{r}^{(k)}$ が求まっている場合次のように求めて行きます。
\begin{equation}
\alpha^{(k)}=\frac{(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})} \tag*{$(6.4-14)$}
\end{equation}
\begin{equation}
\hspace{5mm} \boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha^{(k)}\boldsymbol{p}^{(k)} \tag*{$(6.4-15)$}
\end{equation}
\begin{equation}
\hspace{7mm} \boldsymbol{r}^{(k+1)}=\boldsymbol{r}^{(k)}-\alpha^{(k)}A\boldsymbol{p}^{(k)} \tag*{$(6.4-16)$}
\end{equation}
\begin{equation}
\hspace{4mm} \beta^{(k)}=-\frac{(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})} \tag*{$(6.4-17)$}
\end{equation}
\begin{equation}
\hspace{8mm} \boldsymbol{p}^{(k+1)}=\boldsymbol{r}^{(k+1)}+\beta^{(k)}\boldsymbol{p}^{(k)} \tag*{$(6.4-18)$}
\end{equation}

このように求めていった残差ベクトルはすべて直交しています。
\begin{equation}
(\boldsymbol{r}^{(i)},\boldsymbol{r}^{(j)})=0 \hspace{5mm} (i\ne j) \tag*{$(6.4-19)$}
\end{equation}
また探索方向ベクトルはすべて $A$ に関して共役です。
\begin{equation}
(\boldsymbol{p}^{(i)},A\boldsymbol{p}^{(j)})=0 \hspace{5mm} (i\ne j) \tag*{$(6.4-20)$}
\end{equation}
次にこれを示します。まず $\boldsymbol{p}^{(0)}=\boldsymbol{r}^{(0)}$ なので、
\begin{equation}
\begin{split}
(\boldsymbol{r}^{(1)},\boldsymbol{r}^{(0)})&=(\boldsymbol{r}^{(0)}-\alpha^{(0)}A\boldsymbol{p}^{(0)},\boldsymbol{r}^{(0)}) \\
&=(\boldsymbol{r}^{(0)},\boldsymbol{r}^{(0)})-\frac{(\boldsymbol{p}^{(0)},\boldsymbol{r}^{(0)})}{(\boldsymbol{p}^{(0)},A\boldsymbol{P}^{(0)})}(A\boldsymbol{p}^{(0)},\boldsymbol{r}^{(0)}) \\
&=0
\end{split} \notag
\end{equation}
となります。次に、
\begin{equation}
\begin{split}
(\boldsymbol{r}^{(2)},\boldsymbol{r}^{(0)})&=(\boldsymbol{r}^{(1)}-\alpha^{(1)}A\boldsymbol{p}^{(1)},\boldsymbol{r}^{(0)})=-\alpha^{(1)}(A\boldsymbol{p}^{(1)},\boldsymbol{p}^{(0)})=0 \\
(\boldsymbol{r}^{(2)},\boldsymbol{r}^{(1)})&=(\boldsymbol{r}^{(1)}-\alpha^{(1)}A\boldsymbol{p}^{(1)},\boldsymbol{r}^{(1)}) \\
&=(\boldsymbol{r}^{(1)},\boldsymbol{r}^{(1)})-\frac{(\boldsymbol{p}^{(1)},\boldsymbol{r}^{(1)})}{(\boldsymbol{p}^{(1)},A\boldsymbol{P}^{(1)})}(A\boldsymbol{p}^{(1)},\boldsymbol{p}^{(1)}-\beta^{(0)}\boldsymbol{p}^{(0)}) \\
&=(\boldsymbol{r}^{(1)},\boldsymbol{r}^{(1)})-(\boldsymbol{r}^{(1)}+\beta^{(0)}\boldsymbol{p}^{(0)},\boldsymbol{r}^{(1)}) \\
&=0
\end{split} \notag
\end{equation}
です。また $\beta^{(k)}$ は、
\begin{equation}
(\boldsymbol{p}^{(k+1)},A\boldsymbol{p}^{(k)})=0 \notag
\end{equation}
が成立するように決めたので、番号が一つ違いのベクトルについては(6.4-20)式は成り立っています。
\begin{equation}
\begin{split}
(\boldsymbol{p}^{(1)},A\boldsymbol{p}^{(0)})&=0 \\
(\boldsymbol{p}^{(2)},A\boldsymbol{p}^{(1)})&=0 \\
\end{split} \notag
\end{equation}
つぎに、
\begin{equation}
\begin{split}
(\boldsymbol{p}^{(2)},A\boldsymbol{p}^{(0)})&=(\boldsymbol{r}^{(2)}+\beta^{(1)}\boldsymbol{p}^{(1)},A\boldsymbol{p}^{(0)})=(\boldsymbol{r}^{(2)},A\boldsymbol{p}^{(0)}) \\
&=(\boldsymbol{r}^{(2)},\frac{1}{\alpha^{(0)}}(\boldsymbol{r}^{(0)}-\boldsymbol{r}^{(1)}))=0
\end{split} \notag
\end{equation}
となります。

いま $k$ 番目までの残差ベクトルに対して(6.4-19)式が成り立っているとしたとき、$k$ より小さな $i$ 番目の残差ベクトルに対して、
\begin{equation}
\begin{split}
(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(i)})&=(\boldsymbol{r}^{(k)}-\alpha^{(k)}A\boldsymbol{p}^{(k)},\boldsymbol{r}^{(i)}) \\
&=-\alpha^{(k)}(A\boldsymbol{p}^{(k)},\boldsymbol{p}^{(i)}-\beta^{(i-1)}\boldsymbol{p}^{(i-1)}) \\
&= 0
\end{split} \notag
\end{equation}
となります。$i$ が $k$ のときは次のようになります。
\begin{equation}
\begin{split}
(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(k)})&=(\boldsymbol{r}^{(k)}-\alpha^{(k)}A\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)}) \\
&=(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})-\alpha^{(k)}(A\boldsymbol{p}^{(k)},\boldsymbol{p}^{(k)}-\beta^{(k-1)}\boldsymbol{p}^{(k-1)}) \\
&=(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})-\frac{(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})}(A\boldsymbol{p}^{(k)},\boldsymbol{p}^{(k)}) \\
&=(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})-(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)}) \\
&=(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})-(\boldsymbol{r}^{(k)}-\beta^{(k-1)}\boldsymbol{p}^{(k-1)},\boldsymbol{r}^{(k)}) \\
&=-\beta^{(k-1)}(\boldsymbol{p}^{(k-1)},\boldsymbol{r}^{(k)}) \\
&=-\beta^{(k-1)}(\boldsymbol{p}^{(k-1)},\boldsymbol{r}^{(k-1)}-\alpha^{(k-1)}A\boldsymbol{p}^{(k-1)}) \\
&=-\beta^{(k-1)}\bigl\{(\boldsymbol{p}^{(k-1)},\boldsymbol{r}^{(k-1)})-\alpha^{(k-1)}(\boldsymbol{p}^{(k-1)},A\boldsymbol{p}^{(k-1)})\bigr\} \\
&=0
\end{split} \tag*{$(6.4-21)$}
\end{equation}

次に $k$ 番目までの探索方向ベクトルに対して(6.4-20)式が成り立っているとしたとき、$k$ より小さな $i$ 番目の探索方向ベクトルに対して、
\begin{equation}
\begin{split}
(\boldsymbol{p}^{(k+1)},A\boldsymbol{p}^{(i)})&=(\boldsymbol{r}^{(k+1)}+\beta^{(k)}\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(i)}) \\
&=(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(i)}) \\
&=\bigl(\boldsymbol{r}^{(k+1)},\frac{1}{\alpha^{(i)}}(\boldsymbol{r}^{(i)}-\boldsymbol{r}^{(i+1)})\bigr) \\
&=0
\end{split} \notag
\end{equation}
となります。$i$ が $k$ のときは次のようになります。
\begin{equation}
\begin{split}
(\boldsymbol{p}^{(k+1)},A\boldsymbol{p}^{(k)})&=(\boldsymbol{r}^{(k+1)}+\beta^{(k)}\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)}) \\
&=(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(i)})-\frac{(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})}(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)}) \\
&=0
\end{split} \notag
\end{equation}

以上(6.4-19)、(6.4-20)式を示すことができました。したがって出発点 $\boldsymbol{x}^{(0)}$ からこの探索を続けてゆけば残差ベクトルはすべて直交しているので、この方程式の次数 $n$ 回繰り返すと次に直交することはできずゼロとならざるをえません。つまり方程式の解に到達することになります。
しかし、コンピュータによる計算では丸め誤差があるのでこの回数で解を得ることができない場合が多くあります。

これまでの計算結果を使うと(6.4-14)式の $\alpha^{(k)}$ と(6.4-17)式の $\beta^{(k)}$ は次のようにかきなおすことができます。
\begin{equation}
\alpha^{(k)}=\frac{(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})}
=\frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})} \tag*{$(6.4-22)$}
\end{equation}
\begin{equation}
\hspace{8mm} \beta^{(k)}=-\frac{(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})}
=\frac{(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(k+1)})}{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})} \tag*{$(6.4-23)$}
\end{equation}
(6.4-22)式は(6.4-21)式から、
\begin{equation}
(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})-(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)})=0 \notag
\end{equation}
であるからです。また(6.4-23)式は、
\begin{equation}
\begin{split}
-\frac{(\boldsymbol{r}^{(k+1)},A\boldsymbol{p}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})}
&=-\frac{\frac{1}{\alpha^{(k)}}(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(k)}-\boldsymbol{r}^{(k+1)})}{\frac{1}{\alpha^{(k)}}(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)}-\boldsymbol{r}^{(k+1)})} \\
&=\frac{(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(k+1)})}{(\boldsymbol{p}^{(k)},\boldsymbol{r}^{(k)}-\boldsymbol{r}^{(k+1)})} \\
&=\frac{(\boldsymbol{r}^{(k+1)},\boldsymbol{r}^{(k+1)})}{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}
\end{split} \notag
\end{equation}
となるからです。ここで(6.4-21)式の、
\begin{equation}
-\beta^{(k-1)}(\boldsymbol{p}^{(k-1)},\boldsymbol{r}^{(k)})=0 \notag
\end{equation}
を使いました。

このようにして連立方程式を解く方法を共役勾配法といいます。ところがこの方法では連立方程式の次数が大きくなってくると収束までに非常に多くの繰り返しが必要になる場合があり、有限要素法の解法には次に述べる不完全コレスキー分解による前処理を行った後この共役勾配法を使うという方法が使われています。
この方法は ICCG 法( Incomplete Cholesky Conjugate Gradient method )とよばれています。

連立方程式(6.4-1)式の係数行列を次のように修正コレスキー分解します。
\begin{equation}
A=LDL^T \tag*{$(6.4-24)$}
\end{equation}
このとき、(6.4-1)式を次のように変形します。ここに$D$は対角行列です。
\begin{equation}
(LD^{1/2})^{-1}A((LD^{1/2})^T)^{-1}(LD^{1/2})^T\boldsymbol{x}=(LD^{1/2})^{-1}\boldsymbol{b} \tag*{$(6.4-25)$}
\end{equation}
ここで、
\begin{equation}
\begin{split}
(LD^{1/2})^{-1}A((LD^{1/2})^T)^{-1}&=(LD^{1/2})^{-1}LDL^T((LD^{1/2})^T)^{-1} \\
&=D^{-1/2}L^{-1}LDL^T(D^{1/2}L^T)^{-1} \\
&=D^{-1/2}DD^{-1/2} \\
&=I
\end{split} \notag
\end{equation}
と単位行列 $I$ となるので修正コレスキー分解を行うことは連立方程式を解くのと同じでです。

ここで、不完全コレスキー分解、
\begin{equation}
A\sim LDL^T \tag*{$(6.4-26)$}
\end{equation}
を行うと、
\begin{equation}
B=(LD^{1/2})^{-1}A((LD^{1/2})^T)^{-1} \tag*{$(6.4-27)$}
\end{equation}
は、単位行列とはならないがある程度それに近くなることが期待できます。
\begin{equation}
\begin{split}
&\boldsymbol{y}=(LD^{1/2})^T\boldsymbol{x} \\
&\boldsymbol{c}=(LD^{1/2})^{-1}\boldsymbol{b}
\end{split} \tag*{$(6.4-28)$}
\end{equation}
とかくと、(6.4-25)式は次のようにかけます。
\begin{equation}
B\boldsymbol{y}=\boldsymbol{c} \tag*{$(6.4-29)$}
\end{equation}
この方程式に共役勾配法を適用します。
残差ベクトルを $\boldsymbol{s}$ とすれば、
\begin{equation}
\boldsymbol{s}=\boldsymbol{c}-B\boldsymbol{y}=(LD^{1/2})^{-1}(\boldsymbol{b}-A\boldsymbol{x})=(LD^{1/2})^{-1}\boldsymbol{r} \tag*{$(6.4-30)$}
\end{equation}
です。

いま $k$ 回までの探査で $\boldsymbol{y}^{(k)}$ が求まっている場合、次の解を次のようにかきます。
\begin{equation}
\boldsymbol{y}^{(k+1)}=\boldsymbol{y}^{(k)}+\alpha^{(k)}\boldsymbol{q}^{(k)} \tag*{$(6.4-31)$}
\end{equation}
ここで $\boldsymbol{q}^{(k)}$ は探索方向のベクトルです。この方向で関数値を最小にすることにより $\alpha^{(k)}$ を求めると(6.4-22)式を求めたのと同様にして次のようになります。
\begin{equation}
\alpha^{(k)}=\frac{(\boldsymbol{s}^{(k)},\boldsymbol{s}^{(k)})}{(\boldsymbol{q}^{(k)},B\boldsymbol{q}^{(k)})} \tag*{$(6.4-32)$}
\end{equation}
この式の右辺は(6.4-27)式、(6.4-30)式を使うと次のように変形できます。
\begin{equation}
\frac{\bigl((LD^{1/2})^{-1}\boldsymbol{r}^{(k)},(LD^{1/2})^{-1}\boldsymbol{r}^{(k)}\bigr)}{(\boldsymbol{q}^{(k)},(LD^{1/2})^{-1}A((LD^{1/2})^T)^{-1}\boldsymbol{q}^{(k)}\bigr)}
=\frac{\bigl(\boldsymbol{r}^{(k)},((LD^{1/2})^{-1})^{T}(LD^{1/2})^{-1}\boldsymbol{r}^{(k)}\bigr)}
{(((LD^{1/2})^{-1})^T\boldsymbol{q}^{(k)},A((LD^{1/2})^T)^{-1}\boldsymbol{q}^{(k)}\bigr)} \notag
\end{equation}
ここで、
\begin{equation}
((LD^{1/2})^{-1})^T(LD^{1/2})^{-1}=((LD^{1/2})^T)^{-1}(LD^{1/2})^{-1}
=(LD^{1/2}(LD^{1/2})^T)^{-1}=(LDL^T)^{-1} \notag
\end{equation}
と変形できます。また、
\begin{equation}
\boldsymbol{q}^{(k)}=(LD^{1/2})^T\boldsymbol{p}^{(k)} \tag*{$(6.4-33)$}
\end{equation}
とおけば次のようになります。
\begin{equation}
\alpha^{(k)}=\frac{\bigl(\boldsymbol{r}^{(k)},(LDL^T)^{-1}\boldsymbol{r}^{(k)}\bigr)}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})} \tag*{$(6.4-34)$}
\end{equation}
これより、
\begin{equation}
\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha^{(k)}\boldsymbol{p}^{(k)} \tag*{$(6.4-35)$}
\end{equation}
および、
\begin{equation}
\boldsymbol{r}^{(k+1)}=\boldsymbol{r}^{(k)}-\alpha^{(k)}A\boldsymbol{p}^{(k)} \tag*{$(6.4-36)$}
\end{equation}
となります。

次の探索ベクトルを、
\begin{equation}
\boldsymbol{q}^{(k+1)}=\boldsymbol{s}^{(k+1)}+\beta^{(k)}\boldsymbol{q}^{(k)} \tag*{$(6.4-37)$}
\end{equation}
とおき、
\begin{equation}
(\boldsymbol{q}^{(k+1)},B\boldsymbol{q}^{(k)})=0 \tag*{$(6.4-38)$}
\end{equation}
を要請すると、(6.4-23)式と同様に次のようになります。
\begin{equation}
\beta^{(k)}=\frac{(\boldsymbol{s}^{(k+1)},\boldsymbol{s}^{(k+1)})}{(\boldsymbol{s}^{(k)},\boldsymbol{s}^{(k)})}
=\frac{(\boldsymbol{r}^{(k+1)},(LDL^T)^{-1}\boldsymbol{r}^{(k+1)})}{(\boldsymbol{r}^{(k)},(LDL^T)^{-1}\boldsymbol{r}^{(k)})} \tag*{$(6.4-39)$}
\end{equation}
これより(6.4-37)式は次のようになります。
\begin{equation}
(LD^{1/2})^T\boldsymbol{p}^{(k+1)}=(LD^{1/2})^{-1}\boldsymbol{r}^{(k+1)}-\beta^{(k)}(LD^{1/2})^T\boldsymbol{p}^{(k)} \notag
\end{equation}
これより、
\begin{equation}
\boldsymbol{p}^{(k+1)}=(LDL^T)^{-1}\boldsymbol{r}^{(k+1)}-\beta^{(k)}\boldsymbol{p}^{(k)} \tag*{$(6.4-40)$}
\end{equation}
となります。

ここで具体的に修正コレスキー分解を行う方法について述べます。係数行列が対称な場合次のように分解できます。
\begin{equation}
A=LDL^T \tag*{$(6.4-41)$}
\end{equation}
ここに $L$ は対角項より右の要素がすべてゼロの左下三角行列で対角項はすべて $1$ です。また $D$ は対角行列です。
これらの要素を $l_{ij}$、$d_{ij}$ とかけば、
\begin{equation}
\begin{split}
&l_{ij}=0 \hspace{5mm} (i<j) \\
&l_{ii}=1 \hspace{5mm} (i=1,2,\cdots,n)
\end{split} \tag*{$(6.4-42)$}
\end{equation}
\begin{equation}
d_{ij}=d_i\delta_{ij} \tag*{$(6.4-43)$}
\end{equation}
です。
これより(6.4-41)式は次のようにかけます。
\begin{equation}
a_{ij}=\sum_k^nl_{ik}d_kl_{jk} \tag*{$(6.4-44)$}
\end{equation}
少し具体的に書き下すと、
\begin{equation}
\begin{split}
&a_{11}=l_{11}d_1l_{11}=d_1 \\
&a_{21}=l_{21}d_1l_{11}=l_{21}d_1 \\
&a_{22}=l_{21}d_1l_{21}+l_{22}d_2l_{22}=l_{21}d_1l_{21}+d_2 \\
&a_{31}=l_{31}d_1l_{11}=l_{31}d_1 \\
&a_{32}=l_{31}d_1l_{21}+l_{32}d_2l_{22}=l_{31}d_1l_{21}+l_{32}d_2 \\
&a_{33}=l_{31}d_1l_{31}+l_{32}d_2l_{32}+d_3 \\
\end{split} \tag*{$(6.4-45)$}
\end{equation}
ですが、変形すれば次のようにかけます。
\begin{equation}
\begin{split}
&l_{11}d_1=a_{11} \\
&l_{21}d_1=a_{21} \\
&l_{22}d_2=a_{22}-l_{21}d_1l_{21} \\
&l_{31}d_1=a_{31} \\
&l_{32}d_2=a_{32}-l_{31}d_1l_{21} \\
&l_{33}d_3=a_{33}-l_{31}d_1l_{31}-l_{32}d_2l_{32}
\end{split} \tag*{$(6.4-46)$}
\end{equation}
ここで、
\begin{equation}
w_{ij}=l_{ij}d_j \tag*{$(6.4-47)$}
\end{equation}
とおけば、これに関しては次のように求めていくことができます。
\begin{equation}
\begin{split}
&w_{11}=a_{11} \\
&w_{21}=a_{21} \\
&w_{22}=a_{22}-l_{21}d_1l_{21}=a_{22}-w_{21}w_{21}/w_{11} \\
&w_{31}=a_{31} \\
&w_{32}=a_{32}-l_{31}d_1l_{21}=a_{32}-w_{31}w_{21}/w{11} \\
&w_{33}=a_{33}-l_{31}d_1l_{31}-l_{32}d_2l_{32}=a_{33}-w_{31}w_{31}/w_{11}-w_{32}w_{32}/w_{22}
\end{split} \tag*{$(6.4-48)$}
\end{equation}
一般的に次のようにかけます。
\begin{equation}
w_{ij}=a_{ij}-\sum_{k=1}^{j-1}w_{i_k}w_{jk}/w_{kk} \tag*{$(6.4-49)$}
\end{equation}

修正コレスキー分解ができると(6.4-1)式は次のようにかけます。
\begin{equation}
LDL^T\boldsymbol{x}=\boldsymbol{b} \tag*{$(6.4-50)$}
\end{equation}
この式は次の2段階に分けて解くことができます。
\begin{equation}
\begin{split}
&L\boldsymbol{y}=\boldsymbol{b} \\
&DL^T\boldsymbol{x}=\boldsymbol{y}
\end{split} \tag*{$(6.4-51)$}
\end{equation}
まず最初の式は次のようにかけます。
\begin{equation}
\sum_{j=1}^{i-1}l_{ij}y_j=b_i \tag*{$(6.4-52)$}
\end{equation}
(6.4-47)式よりこの式は次のようにかけます。
\begin{equation}
\sum_{j=1}^{i-1}w_{ij}\bigl(\frac{y_j}{d_j}\bigr)=b_i \tag*{$(6.4-53)$}
\end{equation}
(6.4-49)式によって $w_{ij}$ が求まっているのでこの式を解いた解は、
\begin{equation}
z_i=y_i/d_i \tag*{$(6.4-54)$}
\end{equation}
となります。
2番目の式は、
\begin{equation}
\sum_{j=i}^nd_il_{ji}x_j=d_iz_i \notag
\end{equation}
両辺を $d_i$ で割って変形すれば次のようになります。
\begin{equation}
\sum_{j=i}^nl_{ji}x_j=\sum_{j=i}^n\frac{w_{ji}}{d_i}x_j=z_i \tag*{$(6.4-55)$}
\end{equation}

次にこの解法のためのプログラムを作成します。
ファイル solver.h をホルダー Include の中に作成し、ファイル solver.cpp をホルダー Source の中に作成します。
ファイル solver.h に次の ICCG 法に関するクラスを定義します。


#include <math.h>
#include <stdlib.h>
#include <fstream>
#include <iostream>
#include <time.h>

using namespace std;

//ICCG法による連立方程式の解法
class Iccg
{
protected:
    int      nfree;                     //方程式の次数
    int*     num_column;                //行の非ゼロ要素数
    int**    index_matrix;              //行の非ゼロ要素の列番号
    double** matrix;                    //係数行列の下三角行列の配列
    double*  vector;                    //方程式の右辺ベクトル

    //解ベクトル用
    double*  x;                         //解ベクトル

    //ICCG法で使うパラメータ
    int      iccg_max;
    double   iccg_eps;

    //ICCG法で使う作業用配列
    double** lmatrix;
    double   *r,*p,*ur,*ap;

    double*  lmatrix_inv;               //高速化のための配列

    ofstream* fp_chk;                   //チェックライト
    bool checkwriteflag;                //チェックライト用フラグ

public:
    Iccg(int nfree, int* num_column, int** index_matrix, double** matrix, double* vector,
         int iccg_max, double iccg_eps, ofstream* fp_chk)
    {
        this->nfree = nfree;
        this->num_column = num_column;
        this->index_matrix = index_matrix;
        this->matrix = matrix;
        this->vector = vector;
        this->iccg_max = iccg_max;
        this->iccg_eps = iccg_eps;
        this->fp_chk = fp_chk;
        checkwriteflag = 0;
        x  = new double[nfree];
        r  = new double[nfree];
        p  = new double[nfree];
        ur = new double[nfree];
        ap = new double[nfree];
        lmatrix = new double*[nfree];
        for (int n=0; n<nfree; n++)
            lmatrix[n] = new double[num_column[n]]; 
    }

    virtual ~Iccg()
    {
        if(nfree > 0)
        {
            delete[] x;
            delete[] r;
            delete[] p;
            delete[] ur;
            delete[] ap;
            for (int n=0; n<nfree; n++)
                delete[] lmatrix[n];
            delete[] lmatrix;
        }
    }
    double getX(int n)
    {
        return x[n];
    }
    void solver();
    void Ic_decomp ();
    void cg ();
    void Ic_solver(double* b, double* x);
    void productAP ();
    int getNfree() { return nfree; }
};

ファイルの先頭にはヘッダーのインクルード文を入れておいてください。

次に、ファイル solver.cpp に次のようなコードを入力します。


#include "../include/solver.h"

//ICCG法による連立方程式の解法
void Iccg::solver()
{
    double vectorEps = 1.0e-20;       //右辺ベクトルのゼロ判定

    for (int n=0; n<nfree; n++)
    {
        x[n]       = 0.0;
        r[n]       = 0.0;
        p[n]       = 0.0;
        ur[n]      = 0.0;
        ap[n]      = 0.0;
    }

    double vector2 = 0.0;
    for (int n=0; n<nfree; n++)
        vector2 = vector2 + vector[n]*vector[n];
    vector2 = sqrt(vector2);

    if(vector2 < vectorEps)
    {
        *fp_chk << "vector**2 = " << vector2 << "  in ICCG" << endl;
        return;
    }

    Ic_decomp ();                    //不完全コレスキー分解

    cg ();                           //CG法
}
void Iccg::Ic_decomp ()
{
    lmatrix[0][0] = matrix[0][0];
    for (int n=1; n<nfree; n++)
    {
        for (int i=0; i<num_column[n]-1; i++)
        {
            int ii = index_matrix[n][i] - 1;
            if(ii >= 0)
            {
                lmatrix[n][i] = matrix[n][i];
                int ij = 0;
                int nj = 0;
                int j1 = index_matrix[ii][ij] - 1;
                int j2 = index_matrix[n][nj] - 1;
                while (j1 < ii && j2 < ii)
                {
                    if (j1 < j2)
                    {
                        ij++;    j1 = index_matrix[ii][ij] - 1;
                    }
                    else if (j1 > j2)
                    {
                        nj++;    j2 = index_matrix[n][nj] - 1;
                    }
                    else
                    {
                        int jj = num_column[j1] - 1;
                        lmatrix[n][i] -= lmatrix[n][nj]*lmatrix[ii][ij]/lmatrix[j1][jj];
                        ij++;    j1 = index_matrix[ii][ij] - 1;
                        nj++;    j2 = index_matrix[n][nj] - 1;
                    }
                }
            }
        }
        lmatrix[n][num_column[n]-1] = matrix[n][num_column[n]-1];
        for (int i=0; i<num_column[n]-1; i++)
        {
            int ii = index_matrix[n][i] - 1;
            if(ii >= 0)
                lmatrix[n][num_column[n]-1] -= lmatrix[n][i]*lmatrix[n][i]/lmatrix[ii][num_column[ii]-1];
        }
    }
} 
void Iccg::cg ()
{
    productAP ();

    double vector2 = 0.0;
    for (int i=0; i<nfree; i++)
    {
        r[i] = vector[i] - ap[i];
        vector2 = vector2 + vector[i]*vector[i];
    }

    Ic_solver(r,ur);

    double rur0 = 0.0;
    for (int i=0; i<nfree; i++)
    {
        p[i] = ur[i];
        rur0 = rur0 + r[i]*ur[i];
    }

    bool convergence = false;
    double  ratMin = 0.0;                  //収束トレランス最小
    int     itMin  = 0;
    for (int it=1; it<=iccg_max; it++)
    {
        productAP ();

        double pap  = 0.0;
        for (int i=0; i<nfree; i++)
        {
            pap = pap + p[i]*ap[i];
        }

        double alpha = rur0/pap;

        for (int i=0; i<nfree; i++)
        {
            x[i] = x[i] + alpha*p[i];
            r[i] = r[i] - alpha*ap[i];
        }
        double rr = 0.0;
        for (int i=0; i<nfree; i++)
            rr += r[i]*r[i];

        double rat = sqrt(rr/vector2);

        if(checkwriteflag == 0)
        {
         // if(it%100 == 0)
                cout << "             ICCG IT. NO. " << it << "      RAT = " << rat << endl;
        }

        //収束した場合
        if (rat < iccg_eps)
        {
            if(checkwriteflag == 0)
            {
                cout << "             ICCG IT. NO. " << it << "      RAT = " << rat << endl;
            }
            convergence = true;
            break;
        }

        //収束しない場合
        Ic_solver(r,ur);

        double rur1 = 0.0;
        for (int i=0; i<nfree; i++)
            rur1 = rur1 + r[i]*ur[i];
 
        double beta = rur1/rur0;
        rur0 = rur1;
 
        for (int i=0; i<nfree; i++)
            p[i] = ur[i] + beta*p[i];
    }

    if(!convergence)
    {
         *fp_chk << "ICCGが発散しました!  IT. NO. " << itMin << "      RAT = " << ratMin << endl;
    }
}
void Iccg::Ic_solver(double* b, double* x)
{
    x[0] = b[0]/lmatrix[0][num_column[0]-1];
    for (int n=1; n<nfree; n++)
    {
        x[n] = b[n];
        for (int i=0; i<num_column[n]-1; i++)
        {
            int ii = index_matrix[n][i] - 1;
            x[n] -= lmatrix[n][i]*x[ii];
        }
        x[n] /= lmatrix[n][num_column[n]-1];
    }
    for (int n=nfree-1; n<0; n--)
    {
        for (int i=0; i<num_column[n]-1; i++)
        {
            int ii = index_matrix[n][i] - 1;
            x[ii] -= lmatrix[n][i]*x[n]/lmatrix[ii][num_column[ii]-1];
        }
    }
}
void Iccg::productAP ()
{        
    for (int n=0; n<nfree; n++)
    {
        ap[n] = 0.0;
        for (int i=0; i<num_column[n]; i++)
        {
            int ii = index_matrix[n][i] - 1;
            ap[n] += matrix[n][i]*p[ii];
        }
    }
    for (int n=1; n<nfree; n++)
    {
        for (int i=0; i<num_column[n]-1; i++)
        {
            int ii = index_matrix[n][i] - 1;
            ap[ii] += matrix[n][i]*p[n];
        }
    }
}