Gauß-Seidel Iteration

Gegeben: Gleichungssystem mit $n$ Gleichungen und $n$ Unbekannten
Gesucht: Näherungslösung, Konvergenzrate



In [1]:

    
# input: number of iterations
n = 2
# input: equations in the form Ax = b
A = matrix([
    [5, -3],
    [4, -6]
])
b = vector([-2, 1])
# input: start vector
x0 = vector([1, 1])
#

Matrix Form der Gauß-Seidel Iteration:

$L$ ist eine strikt linke untere Dreiecksmatrix (Diagonalelemente sind 0)
$R$ ist eine strikt rechte obere Dreiecksmatrix (Diagonalelemente sind 0)
$D$ ist eine Diagonalmatrix



In [2]:

    
L = matrix([[A[i,j] if j < i else 0
             for j in range(A.dimensions()[1])] 
             for i in range(A.dimensions()[0])
            ])
R = matrix([[A[i,j] if j > i else 0
             for j in range(A.dimensions()[1])] 
             for i in range(A.dimensions()[0])
            ])
D = matrix([[A[i,j] if j == i else 0
             for j in range(A.dimensions()[1])] 
             for i in range(A.dimensions()[0])
             ])
show(LatexExpr("L = "), L, LatexExpr(r"\quad D = "), D, LatexExpr(r"\quad R = "), R)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}L = \left(\begin{array}{rr} 0 & 0 \\ 4 & 0 \end{array}\right) \quad D = \left(\begin{array}{rr} 5 & 0 \\ 0 & -6 \end{array}\right) \quad R = \left(\begin{array}{rr} 0 & -3 \\ 0 & 0 \end{array}\right)$

$$ T_{GS} = -(D + L)^{-1} \cdot R $$

Gauss-Seidel-Schritt: $$ x_{k+1} = T_{GS} \cdot x_k + (D + L)^{-1} b $$



In [3]:

    
# compute intermediate values
TGS = -(D + L)^-1 * R
DLb = (D+L)^-1 * b
show(TGS, DLb)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 0 & \frac{3}{5} \\ 0 & \frac{2}{5} \end{array}\right) \left(-\frac{2}{5},\,-\frac{13}{30}\right)$



In [4]:

    
# perform n steps of the iteration, saving the result
x = [x0]
show(LatexExpr("x_0 = "), x0)
for i in range(n):
    xk = TGS * x[-1] + DLb
    x.append(xk)
    show(LatexExpr("x_{} = ".format(i + 1)), xk)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}x_0 = \left(1,\,1\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}x_1 = \left(\frac{1}{5},\,-\frac{1}{30}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}x_2 = \left(-\frac{21}{50},\,-\frac{67}{150}\right)$



In [5]:

    
show(LatexExpr(r"x_n \approx "), vector(RDF, x[-1]))









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}x_n \approx \left(-0.42,\,-0.446666666667\right)$

Aufteilen von $A$ in $$A = S - T$$ wobei $S$ eine linke untere Dreiecksmatrix ist



In [6]:

    
S = matrix([[A[i,j] if j <= i else 0
             for j in range(A.dimensions()[1])] 
             for i in range(A.dimensions()[0])
            ])
T = matrix([[A[i,j] if j > i else 0
             for j in range(A.dimensions()[1])] 
             for i in range(A.dimensions()[0])
            ])
assert A == S + T
show(S, T)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 5 & 0 \\ 4 & -6 \end{array}\right) \left(\begin{array}{rr} 0 & -3 \\ 0 & 0 \end{array}\right)$



In [7]:

    
sit = S^-1 * T
show(sit)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 0 & -\frac{3}{5} \\ 0 & -\frac{2}{5} \end{array}\right)$

Spekatralradius $\rho(S^{-1}T) = \max |\lambda_i|$



In [8]:

    
spektral_radius = lambda M: max(map(abs, M.eigenvalues()))
show(spektral_radius(sit))









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{5}$

Konvergenzrate $r = - \log_{10} \rho(S^{-1} T)$



In [9]:

    
r = - log(spektral_radius(sit), 10)
show(LatexExpr("r = "), r, LatexExpr(r"\approx"), RDF(r))









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}r = -\frac{\log\left(\frac{2}{5}\right)}{\log\left(10\right)} \approx 0.397940008672$

Exakte Lösung zum Vergleich



In [10]:

    
show(A^-1 * b)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(-\frac{5}{6},\,-\frac{13}{18}\right)$