ガウス・ザイデル法

解は反復的に得られる $${\displaystyle \mathbf {L} \mathbf {x} ^{(k+1)}=\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)},}$$ 行列${\displaystyle \mathbf {A} }$下三角成分に分解される${\displaystyle \mathbf {L} }$、そして厳密に上三角関数の成分${\displaystyle \mathbf {U} }$そのため${\displaystyle \mathbf {A} =\mathbf {L} +\mathbf {U} }$^{[ 4 ]}より具体的には、${\displaystyle \mathbf {A} }$の中へ${\displaystyle \mathbf {L} }$そして${\displaystyle \mathbf {U} }$次のように与えられます。

$${\displaystyle \mathbf {A} =\underbrace {\begin{bmatrix}a_{11}&0&\cdots &0\\a_{21}&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}} _{\textstyle \mathbf {L} }+\underbrace {\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}} _{\textstyle \mathbf {U} }.}$$

しかし、三角形の形状を利用することで${\displaystyle \mathbf {L} }$要素${\displaystyle \mathbf {x} ^{(k+1)}}$各行について順次計算できます${\displaystyle i}$前方代入法を用いる：^{[ 5 ]}$${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.}$$

この式では、反復ごとに 2 つの総和が使用されていますが、これは 1 つの総和として表現できます。${\displaystyle \sum _{j\neq i}a_{ij}x_{j}}$これは、最近計算された反復処理を使用します。${\displaystyle x_{j}}$手順は一般的に、反復によって行われた変更が、十分に小さな残差などの許容範囲を下回るまで続けられます。

ガウス・ザイデル法の要素ごとの公式は、（反復）ヤコビ法の公式と関連しているが、重要な違いがある。

ガウス・ザイデル法では、${\displaystyle \mathbf {x} ^{(k+1)}}$要素を使用する${\displaystyle \mathbf {x} ^{(k+1)}}$すでに計算されている要素のみ${\displaystyle \mathbf {x} ^{(k)}}$計算されていない${\displaystyle (k+1)}$これは、ヤコビ法とは異なり、要素が計算されるにつれて上書きできるため、必要なストレージベクトルは1つだけで済むことを意味し、非常に大規模な問題に対して有利となる場合があります。

しかし、ヤコビ法とは異なり、各要素の計算は一般的に並列処理が非常に困難であり、クリティカルパスが非常に長くなる可能性があるため、疎行列の場合に最も適しています。さらに、各反復における値は、元の方程式の次数に依存します。

ガウス・ザイデル法は、次の過緩和法と同じです。${\displaystyle \omega =1}$。

線形システムとして以下のように表される${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$次のように与えられます。 $${\displaystyle \mathbf {A} ={\begin{bmatrix}16&3\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\end{bmatrix}}.}$$

方程式を使用する $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})}$$ 形式で $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c} }$$ どこ：

${\displaystyle \mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .}$

分解する${\displaystyle \mathbf {A} }$下三角成分の合計に${\displaystyle \mathbf {L} }$そして厳密な上三角成分${\displaystyle U}$: $${\displaystyle \mathbf {L} ={\begin{bmatrix}16&0\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\end{bmatrix}}.}$$

逆${\displaystyle \mathbf {L} }$は： $${\displaystyle \mathbf {L} ^{-1}={\begin{bmatrix}16&0\\7&-11\end{bmatrix}}^{-1}={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\\\end{bmatrix}}.}$$

次に以下を探してください： $${\displaystyle {\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\end{bmatrix}}={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1194\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}11\\13\end{bmatrix}}={\begin{bmatrix}0.6875\\-0.7439\end{bmatrix}}.\end{aligned}}}$$

と${\displaystyle \mathbf {T} }$そして${\displaystyle \mathbf {c} }$ベクトル${\displaystyle \mathbf {x} }$反復的に取得できる。

まず最初に、${\displaystyle \mathbf {x} ^{(0)}}$、 例えば$${\displaystyle \mathbf {x} ^{(0)}={\begin{bmatrix}1.0\\1.0\end{bmatrix}}.}$$推測値が最終解に近いほど、アルゴリズムに必要な反復回数は少なくなる。

次に計算します。 $${\displaystyle {\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}1.0\\1.0\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8494\\-0.6413\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8494\\-0.6413\\\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(4)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(5)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(6)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(7)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\end{aligned}}}$$

予想通り、アルゴリズムは解に収束する。 $${\displaystyle \mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} \approx {\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}}$$。

実際、行列Aは厳密には対角優位行列ではあるが、正定値行列ではない。

別の線形システムとして示されている${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$次のように与えられます。

$${\displaystyle \mathbf {A} ={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\\\end{bmatrix}}.}$$

方程式を使用する $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})}$$ 形式で $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c} }$$ どこ：

${\displaystyle \mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .}$

分解する${\displaystyle \mathbf {A} }$下三角成分の合計に${\displaystyle \mathbf {L} }$そして厳密な上三角成分${\displaystyle \mathbf {U} }$: $${\displaystyle \mathbf {L} ={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}.}$$

逆${\displaystyle \mathbf {L} }$は： $${\displaystyle \mathbf {L} ^{-1}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}^{-1}={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}.}$$

次に以下を探してください： $${\displaystyle {\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}={\begin{bmatrix}0.000&-1.500\\0.000&1.071\\\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}.\end{aligned}}}$$

と${\displaystyle \mathbf {T} }$そして${\displaystyle \mathbf {c} }$ベクトル${\displaystyle \mathbf {x} }$反復的に取得できる。

まず最初に、私たちは選択しなければなりません${\displaystyle \mathbf {x} ^{(0)}}$、 例えば$${\displaystyle \mathbf {x} ^{(0)}={\begin{bmatrix}1.1\\2.3\end{bmatrix}}}$$

次に計算します。 $${\displaystyle {\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}1.1\\2.3\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}4.911\\-1.651\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&=\cdots .\end{aligned}}}$$

収束テストでは、アルゴリズムが発散することがわかります。実際、行列は${\displaystyle \mathbf {A} }$対角優位でも正定値でもない。すると、厳密解への収束が起こる。 $${\displaystyle \mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} ={\begin{bmatrix}-38\\29\end{bmatrix}}}$$ 保証されるものではなく、この場合、発生しません。

与えられたと仮定する${\displaystyle n}$方程式と出発点${\displaystyle \mathbf {x} _{0}}$ガウス・ザイデル反復法のどのステップでも、最初の式を解いて${\displaystyle x_{1}}$に関しては${\displaystyle x_{2},\dots ,x_{n}}$;次に2番目の方程式を解いて${\displaystyle x_{2}}$に関しては${\displaystyle x_{1}}$発見されたばかりで、残りの${\displaystyle x_{3},\dots ,x_{n}}$; そして、${\displaystyle x_{n}}$そして、収束が達成されるまで反復計算を繰り返すか、解の発散が事前に定義されたレベルを超えて発散し始めた場合は処理を中断する。

例を挙げて考えてみましょう。 $${\displaystyle {\begin{array}{rrrrl}10x_{1}&-x_{2}&+2x_{3}&&=6,\\-x_{1}&+11x_{2}&-x_{3}&+3x_{4}&=25,\\2x_{1}&-x_{2}&+10x_{3}&-x_{4}&=-11,\\&3x_{2}&-x_{3}&+8x_{4}&=15.\end{array}}}$$

解決する${\displaystyle x_{1},x_{2},x_{3}}$そして${\displaystyle x_{4}}$与える： $${\displaystyle {\begin{aligned}x_{1}&=x_{2}/10-x_{3}/5+3/5,\\x_{2}&=x_{1}/11+x_{3}/11-3x_{4}/11+25/11,\\x_{3}&=-x_{1}/5+x_{2}/10+x_{4}/10-11/10,\\x_{4}&=-3x_{2}/8+x_{3}/8+15/8.\end{aligned}}}$$

初期近似値が (0, 0, 0, 0)であると仮定すると、最初の近似解は次のように与えられます。$${\displaystyle {\begin{aligned}x_{1}&=3/5=0.6,\\x_{2}&=(3/5)/11+25/11=3/55+25/11=2.3272,\\x_{3}&=-(3/5)/5+(2.3272)/10-11/10=-3/25+0.23272-1.1=-0.9873,\\x_{4}&=-3(2.3272)/8+(-0.9873)/8+15/8=0.8789.\end{aligned}}}$$

得られた近似値を用いて、所望の精度に達するまで反復手順を繰り返す。以下は、4回の反復後の近似解である。

このシステムの正確な解は(1, 2, −1, 1)です。

以下の反復手順により、線形方程式系の解ベクトルが得られます。

import numpy as np反復回数制限= 1000# 行列を初期化しますA = np.array ( [ [ 10.0 , -1.0 , 2.0 , 0.0 ] , [ -1.0 , 11.0 , -1.0 , 3.0 ] , [ 2.0 , -1.0 , 10.0 , -1.0 ] , [ 0.0 , 3.0 , -1.0 , 8.0 ] , ] ) # 右辺ベクトルを初期化しますb = np.array ( [ 6.0 , 25.0 , -11.0 , 15.0 ] )print ( "連立方程式:" ) for i in range ( A.shape [ 0 ]): row = " + " .join ( f " { A [ i , j ] : 3g } *x { j + 1 } " for j in range ( A.shape [ 1 ] ) ) print ( f " [ { row } ] = [ { b [ i ] : 3g } ] " )x = np.zeros_like ( b ) for it_count in range ( 1 , ITERATION_LIMIT ) : x_new = np.zeros_like ( x ) print ( f " Iteration { it_count } : { x } " ) for i in range ( A.shape [ 0 ] ) : s1 = A [ i , : i ] @ x_new [: i ] s2 = A [ i , i + 1 : ] @ x [ i + 1 :] x_new [ i ] = ( b [ i ] - s1 - s2 ) / A [ i , i ] if np.allclose ( x , x_new , rtol = 1e - 8 ) : break x = x_newprint ( f "解答: { x } " ) error = A @ x - b print ( f "エラー: { error } " )

出力結果は以下のとおりです。

連立方程式: [ 10*x1 + -1*x2 + 2*x3 + 0*x4] = [ 6] [ -1*x1 + 11*x2 + -1*x3 + 3*x4] = [ 25] [ 2*x1 + -1*x2 + 10*x3 + -1*x4] = [-11] [ 0*x1 + 3*x2 + -1*x3 + 8*x4] = [ 15]反復 1: [ 0. 0. 0. 0.]反復 2: [ 0.6 2.32727273 -0.98727273 0.87886364]反復 3: [ 1.03018182 2.03693802 -1.0144562 0.98434122]反復4: [ 1.00658504 2.00355502 -1.00252738 0.99835095]反復 5: [ 1.00086098 2.00029825 -1.00030728 0.99984975]反復 6: [ 1.00009128 2.00002134 -1.00003115 0.9999881 ]反復 7: [ 1.00000836 2.00000117 -1.00000275 0.99999922]反復 8: [ 1.00000067 2.00000002 -1.00000021 0.99999996]反復 9: [ 1.00000004 1.99999999 -1.00000001 1. ]反復 10: [ 1. 2. -1. 1.]解: [ 1. 2. -1. 1.]エラー: [ 2.06480930e-08 -1.25551054e-08 3.61417563e-11 0.00000000e+00]

以下のコードは、数式を使用します。 $${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad {\begin{array}{l}i=1,2,\ldots ,n\\k=0,1,2,\ldots \end{array}}}$$

function x = gauss_seidel ( A, b, x, iters ) for i = 1 : iters for j = 1 : size ( A , 1 ) x ( j ) = ( b ( j ) - sum ( A ( j ,:) '.* x ) + A ( j , j ) * x ( j )) / A ( j , j );終わり終わり終わり