【科学技術計算講座4-8】共役勾配法で線形方程式系を解く

前回は有限体積法の計算のもとになる線形方程式を作成するPythonプログラムを作りました。

今回はその線形方程式を解くための解法について説明したいと思います。

線形方程式系

これまで見てきたように、セル $C$ とその隣接セル $F$ に対して、線形方程式を書くと次のような式になります。

$\begin{matrix} (8-1) & a_{C} T_{C} + \sum_{F} a_{F} T_{F} = b_{C} \end{matrix}$

実際にはこの式がセル数分あり連立方程式を形成しています。このような線形方程式の連立方程式を線形方程式系（system of linear algebraic equations）といいます。

これらの式をセル数分書くのは大変なので、次のように表します。

$\begin{matrix} (8-2) & A T = b \end{matrix}$

ここで、 $T$ は各セルの温度 $T_{C}$ をセル数分格納したベクトルを表します。 $b$ は同様に $b_{C}$ のベクトル量です。 $A$ は係数 $a_{C}$ 、 $a_{F}$ を表す行列で係数行列と呼びます。この式は次の図のような構造をしています。

行列 $A$ は ● がセル $C$ の係数 $a_{C}$ を表していて、行列の対角成分に位置しています。○ は隣接セルの係数 $a_{F}$ を表しています。この行列は対角成分を挟んで折り返すと係数が一致します。このような行列は対称行列（symmetric matrix）と呼ばれます。

また、●○以外の成分は全てゼロです。セルの数がいくら多くても、一つの行に高々５つの成分（自分自身と隣接セル４つ）しかなく、ほとんどがゼロになっています。このような行列を疎行列（sparse matrix）といいます。

線形方程式系の解法

(8-2)式の線形方程式系を解くと、求めたい温度の解が得られます。この解法には実に様々な方法があります。

差分法のときに紹介したヤコビ法やガウス=ザイデル法もその一つです。今回の問題もガウス=ザイデル法を使って解くことは可能です。 $T_{C} =$ の式に変形して反復計算するだけなので、プログラムは簡単に作れます。ただし、実際にガウス=ザイデル法のプログラムで計算してみると、収束までの反復回数が多く、なかなか収束しないことがわかります。

そこで今回は、ガウス=ザイデル法と同じ反復解法の一つである共役勾配法（conjugate gradient method、CG法）で解いてみたいと思います。

共役勾配法とは対称正定値行列（※）を解くための反復解法の一種で、収束が早く今回のような大規模な疎行列が係数行列となる線形方程式系の解法としてよく使われます。共役勾配法にはいくつかの種類がありますが、汎用の流体解析ソフトなどでもこれらの解法が使われています。当サイトのCATCFDzeroはBiCGSTABという共役勾配法の一種を使っています。

今回はベースとなるオリジナルの共役勾配法で解いてみたいと思います。

※　行列が正定値であるとは、対称行列 $A$ が $0$ でない任意のベクトル $x$ に対して $(x, A x) > 0$ を満たすことをいいます。。。ここではこれ以上深入りしないことにしましょう。

共役勾配法

共役勾配法のコンセプトは、次の式の最小値を求めるというものです。

$\begin{matrix} (8-3) & f (T) = \frac{1}{2} (T, A T) - (b, T) \end{matrix}$

ここで、 $(b, T)$ はベクトル $b$ と $T$ の内積を表します。

なぜこの最小値が(8-2)式を解くことになるのかというと、 $f (T)$ の最小値のところでは $f (T)$ を $T$ で微分した値（勾配） $f^{'} (T)$ がゼロになることに由来します。すなわち、

$\begin{matrix} (8-4) & f^{'} (T) = A T - b = 0 \end{matrix}$

となります。これは(8-2)式を示しています。つまり、 $f (T)$ が最小になるときの $T$ が、求める解となるわけです。

共役勾配法は反復法の一種です。反復法とは、適当な初期値から始めて、次のステップの解を近似しながら繰り返し計算を行い、目的とする解を求める方法です。

今、反復ステップ $k$ の温度を $T_{k}$ とすると、次のステップの近似解を

$\begin{matrix} (8-5) & T_{k + 1} = T_{k} + α_{k} d_{k} \end{matrix}$

として求めることにします。ここでベクトル $d_{k}$ は次の解への方向ベクトルです。関数 $f (T)$ の最小値（極値）を求める方法の一つは、次の図の青色の経路をたどる方法がイメージしやすいです。青色の経路は関数の負の勾配 $- f^{'} (T)$ に沿って進んでいく方法です。（図は関数 $f (T)$ の等高線を表していて、その中央が最小値（ちょうどお椀の底）になっています。その最小点で真の解 $T$ をとります。）

ここで $r_{k}$ は残差ベクトルといいます。

$\begin{matrix} (8-6) & r_{k} = b - A T_{k} = - f^{'} (T_{k}) \end{matrix}$

この方法は、最急降下法（steepest descent method）と呼ばれています。(8-5)式で $d_{k} = r_{k}$ とした場合です。勾配 $r_{k}$ に沿って下っていき、その線上の最低点まできたら再びその地点での勾配方向に沿って下っていきます。ただし、この方法は図からもイメージできるように真の解 $T$ にたどり着くまでに反復回数が多そうです。

実際の真の解 $T$ の方向は図の赤いベクトルの方向であり、 $r_{k}$ の方向から少しずれています。これを次式で求めることにします。

$\begin{matrix} (8-7) & d_{k} = r_{k} + β_{k - 1} d_{k - 1} \end{matrix}$

これは前のステップの $d_{k - 1}$ を使って $r_{k}$ から少しずらすことを意味します。このとき、 $d_{k}$ と $A d_{k - 1}$ が直交するように $β$ を選んでやります。このようにすると、これまでの方向と重ならない方向で経路を決めることができて効率的です。このベクトル $d$ の方向を行列 $A$ に対して共役な方向といいます。

この考え方から $α$ と $β$ は次のように計算されます。

$\begin{matrix} (8-8) & α_{k} = \frac{(d_{k}, r_{k})}{(d_{k}, A d_{k})} \end{matrix}$

$\begin{matrix} (8-9) & β_{k} = \frac{(r_{k + 1}, r_{k + 1})}{(r_{k}, r_{k})} \end{matrix}$

共役勾配法のアルゴリズム

数値計算をするために共役勾配法のアルゴリズムを整理しておきます。

　適当な初期値 $T_{0}$ を決める。
　 $r_{0} = b - A T_{0}$ 、 $d_{0} = r_{0}$ とする。
　収束するまで以下の計算を繰り返す。
①　 $α_{k} = (d_{k}, r_{k}) / (d_{k}, A d_{k})$
②　 $T_{k + 1} = T_{k} + α_{k} d_{k}$
③　 $r_{k + 1} = r_{k} - α_{k} A d_{k}$
④　残差 $r_{k + 1}$ で収束判定。収束なら終了。
⑤　 $β_{k} = (r_{k + 1}, r_{k + 1}) / (r_{k}, r_{k})$
⑥　 $d_{k + 1} = r_{k + 1} + β_{k} d_{k}$

共役勾配法の導出過程はややこしいですが、アルゴリズムはとてもシンプルです。

残差と収束判定

反復法はどこかで収束を判定し、計算を打ち切ってやる必要があります。以前、ガウス=ザイデル法では、前のステップとの差をとってそれが小さくなった段階で収束と判定しました。

$\begin{matrix} (8-10) & | | T_{k + 1} - T_{k} | | \leq ε \end{matrix}$

$| | x | |$ はベクトル $x$ の絶対値（ノルム）を表します。これは誤差基準の収束判定です。

今回の共役勾配法では残差 $r$ を計算しています。残差とは、(8-6)式を見ればわかるように、解くべき線形方程式系(8-2)式の右辺と左辺の差を表しています。残差がゼロになると、(8-2)式が満たされ真の解が求まったことになります。

この残差がゼロになる（実際には小さな数値 $ε$ 以下になる）のを収束判定に使います。残差を使う収束判定の式もいろいろありますが、今回は右辺 $b$ の絶対値で割った相対残差という指標を使うことにします。

$\begin{matrix} (8-11) & \frac{| | r_{k + 1} | |}{| | b | |} \leq ε \end{matrix}$

(8-11)式を満たせば共役勾配法の反復を終了します。

まとめ

今回は解くべき線形方程式系を行列で表し、共役勾配法という解法について説明しました。共役勾配法は、数学的な話は難しいですが、プログラムにするためのアルゴリズムはとてもシンプルです。

次回は共役勾配法をPythonでプログラミングして、有限体積法のプログラムを完成させたいと思います。

←前回　Pythonで有限体積法プログラミング～線形方程式の作成

→次回　有限体積法プログラミング～２次元熱伝導問題

全体の目次