数学とか語学とか楽しいよね

数学とか語学とか楽しいよね。ドイツ語とかフランス語とか数値計算とか勉強したことをまとめます。

【数値計算】Gauss-Seidel法(ガウス・ザイデル法)の説明 理論編

Gauss-Seidel法とは

連立一次方程式を解くための基本的な反復法です。反復法とは一回の計算で解を求めるのではなく、ある決まった手順を何度も繰り返しながら真の解へと近づいていく方法です。数値計算では有限桁までしか表現できないので、得られる解は誤差を含んだものになります。ちなみに、何故連立一次方程式を解くのかというと、有限要素法などで方程式を離散化すると最終的に連立一次方程式に帰着されるからです。

反復手順

まず解こうとする連立一次方程式を

 Ax=b

としましょう。ここで、  A n \times n の行列、 x b n 次元のベクトルです。 A b が既知で、 x が未知です。例えば、 3 \times 3 だと

 \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array} \right)

です。この行列  A

 A=D+L+U

と分解します。ここで、 D は対角行列(対角成分以外は成分が0の行列)、 L は下三角行列(対角線の下にだけ非零の成分が入っている行列)、 U は上三角行列(対角線の上にだけ非零の成分が入っている行列)です。いかめしく書いていますが、具体例を見ればすぐ分かると思います。例えば、

 A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\2 & 2 & 10 \\ \end{pmatrix}

としましょう。このとき、 A は、

 A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\ 2 & 2 & 10 \\ \end{pmatrix}=\begin{pmatrix} 10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 10 \\ \end{pmatrix}+\begin{pmatrix} 0 & 0 & 0 \\ 2 & 0 & 0 \\ 2 & 2 & 0 \\ \end{pmatrix}+\begin{pmatrix} 0 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{pmatrix}

のように分解できるよね、と言っているだけです。さて、この分解を最初の連立一次方程式に代入すると

 (D+L+U)x=b

となります。次に、 Ux を右辺に移項して

 (D+L)x=-Ux+b

です。最後に、この両辺に  D+L逆行列  (D+L)^{-1} を左からかけて

 x=-(D+L)^{-1}Ux+(D+L)^{-1}b

を得ます。これをもとにGauss-Seidel法の反復式を作ります。まず、 k ステップ目での  x の近似解を  x^k としましょう。次のステップ( k+1 ステップ目)での  x の近似解を先ほどの式を使って

 x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b

のように決めることにします。これがGauss-Seidel法の反復式です。適当な初期値  x^0 からスタートした数列  x^0, x^1,\cdots は、もし収束するなら  Ax=b を満たす  x へと収束します。なぜなら、数列  x^0, x^1,\cdots x^* に収束したとすると

 x^{*}=-(D+L)^{-1}Ux^*+(D+L)^{-1}b

が成り立ちます。後はこの式を逆に変形していけば

 Ax^*=b

となるので  x=x^* です。

実際の数値計算は有限桁なので、誤差が一定値以下になったら計算を切り上げます。誤差の決め方ですが、自分で決めた誤差の下限  \epsilon (小さい値)に対して

 \max_{1 \leq i \leq n} |x_i^{k+1} - x_i^{k}|\ < \epsilon

となったら反復を止めることにします。この式の意味は、 k+1 ステップ目の近似解と  k ステップ目の近似解の差を取ったとき、その差が最大になるものが、あらかじめ定めておいた誤差の下限を下回ったら反復を止めるということです。つまり、近似解の値があまり更新されなくなったら反復を止めよう、と言っているだけです。実は他にも止める基準はあるのですがそれはまた今度。

長くなってきたので一度切りますが、これだけではきっと分からないと思うので、具体的な行列に対する例は次回書きます。

注 Gauss-Seidel法の収束については、結構大変なので別の記事で触れます。

参考

英語で学ぶ数値解析

英語で学ぶ数値解析

楽しい反復法

楽しい反復法