# Example: compare GS and SOR

Let's look at 

$$
\begin{equation}
\begin{pmatrix}
  4  &  3 & 0 \\
  3  &  4 & -1 \\
  0  & -1 & 4
\end{pmatrix}
%
\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} 
%
= \begin{pmatrix} 24 \\ 30 \\ -24 \end{pmatrix} \nonumber
\end{equation}
$$

with both GS and SOR. In  both cases we'll start with $\vec{x}^{(0)} = [1, 1, 1]^T$. 

The exact solution is $[3, 4, -5]^T$.

## Gauss Seidel

Recall:

for $i = 1, \dots, n$:
$$\qquad x^{(k+1)}_i = \frac{1}{a_{ii}}(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)}) \:.$$

Writing this out for our case gives

\begin{align}
x_1^{(k+1)} &= \frac{1}{4}(24 - 3 x_2^{(k)} - 0 x_3^{(k)}) 
= 6 - \frac{3}{4}x_2^{(k)} \nonumber \\
%
x_2^{(k+1)} &= \frac{1}{4}(30 - 3 x_1^{(k+1)} - (-1) x_3^{(k)}) 
= \frac{30}{4} - \frac{3}{4}x_1^{(k+1)} + \frac{1}{4}x_3^{(k)}\nonumber \\
%
x_3^{(k+1)} &= \frac{1}{4}(-24 - 0 x_1^{(k+1)} - (-1) x_2^{(k+1)}) 
= -6 + \frac{1}{4}x_2^{(k+1)}\nonumber
\end{align}


```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def GS(A, b, x_0, K, ref):

    # Define the size of our system
    n = len(x_0)
    
    # We've got an initial guess, x_0, and an iteration count to hit, K
    # To track progress through the iterations, we'll store all of our
    # estimates. You wouldn't typically do this, but we will for illustration.
    x = np.zeros((K+1, n))
    x[0] = x_0
    # We'd also like to track our error (since we know the real solution)
    error = np.zeros((K+1, n))
    error[0] = np.abs(ref - x[0])
    
    for k in range(0,K):

        for i in range(0, n):

            # Initialize the summation terms to zero on each iteration
            sum_1 = 0.0
            sum_2 = 0.0

            # the k+1 term from 0 to i-1 
            for j in range(i):

                sum_1 = sum_1 + A[i][j] * x[k+1][j]

            # the k term from i+1 to n-1
            for j in range(i+1, n):

                sum_2 = sum_2 + A[i][j] * x[k][j]

            x[k+1][i] = 1.0/A[i][i] * (b[i] - sum_1 - sum_2)
            
        # compute the error
        error[k+1] = np.abs(ref - x[k+1])

    return x, error
    

x_0 = np.ones(3)
A = [[ 4,  3, 0], [3,  4, -1] ,[0, -1, 4]]
b = [24, 30, -24]
ref = [3, 4, -5]

x_GS, error_GS = GS(A, b, x_0, 4, ref)
print x_GS, "\n"
print error_GS

```

    [[ 1.          1.          1.        ]
     [ 5.25        3.8125     -5.046875  ]
     [ 3.140625    3.8828125  -5.02929688]
     [ 3.08789062  3.92675781 -5.01831055]
     [ 3.05493164  3.95422363 -5.01144409]] 
    
    [[ 2.          3.          6.        ]
     [ 2.25        0.1875      0.046875  ]
     [ 0.140625    0.1171875   0.02929688]
     [ 0.08789062  0.07324219  0.01831055]
     [ 0.05493164  0.04577637  0.01144409]]


## Successive Over Relaxation

Recall:

for $i = 1, \dots, n$:
$$\qquad x^{(k+1)}_i = (1-\omega)x_i^{(k)} + \frac{\omega}{a_{ii}}(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)}) \:.$$

Let $\omega = 1.25 = 5/4$

Writing this out for our case gives

$$\begin{align}
x_1^{(k+1)} &= (1-\frac{5}{4})x_1^{(k)} + \frac{\frac{5}{4}}{4}(24 - 3 x_2^{(k)} - 0 x_3^{(k)}) 
= \frac{15}{2} - \frac{1}{4}x_1^{(k)} - \frac{15}{16}x_2^{(k)} \nonumber \\
%
x_2^{(k+1)} &= (1-\frac{5}{4})x_2^{(k)} + \frac{\frac{5}{4}}{4}(30 - 3 x_1^{(k+1)} - (-1) x_3^{(k)}) 
= \frac{150}{16} - \frac{15}{16}x_1^{(k+1)} - \frac{1}{4}x_2^{(k)} + \frac{5}{16}x_3^{(k)}\nonumber \\
%
x_3^{(k+1)} &= (1-\frac{5}{4})x_3^{(k)} + \frac{\frac{5}{4}}{4}(-24 - 0 x_1^{(k+1)} - (-1) x_2^{(k+1)}) 
= -\frac{15}{2} + \frac{5}{16}x_2^{(k+1)} - \frac{1}{4}x_3^{(k)}\nonumber
\end{align}$$


```
def SOR(A, b, x_0, w, K, ref):

    # Define the size of our system
    n = len(x_0)
    
    # We've got an initial guess, x_0, and an iteration count to hit, K
    # To track progress through the iterations, we'll store all of our
    # estimates. You wouldn't typically do this, but we will for illustration.
    x = np.zeros((K+1, n))
    x[0] = x_0
    # We'd also like to track our error (since we know the real solution)
    error = np.zeros((K+1, n))
    error[0] = np.abs(ref - x[0])
    
    for k in range(0,K):

        for i in range(0, n):

            # Initialize the summation terms to zero on each iteration
            sum_1 = 0.0
            sum_2 = 0.0

            # the k+1 term from 0 to i-1 
            for j in range(i):

                sum_1 = sum_1 + A[i][j] * x[k+1][j]

            # the k term from i+1 to n-1
            for j in range(i+1, n):

                sum_2 = sum_2 + A[i][j] * x[k][j]

            x[k+1][i] = (1 - w) * x[k][i] + w/A[i][i] * (b[i] - sum_1 - sum_2)
            
        # compute the error
        error[k+1] = np.abs(ref - x[k+1])

    return x, error

x_SOR, error_SOR = SOR(A, b, x_0, 1.25, 4, ref)
print x_SOR, "\n"
print error_SOR
```

    [[ 1.          1.          1.        ]
     [ 6.3125      3.51953125 -6.65014648]
     [ 2.62231445  3.95852661 -4.60042381]
     [ 3.13330269  4.01026464 -5.09668635]
     [ 2.95705123  4.00748383 -4.97348972]] 
    
    [[ 2.          3.          6.        ]
     [ 3.3125      0.48046875  1.65014648]
     [ 0.37768555  0.04147339  0.39957619]
     [ 0.13330269  0.01026464  0.09668635]
     [ 0.04294877  0.00748383  0.02651028]]


Let's compare with a plot or two


```
idx = [1, 2, 3]
plt.plot(idx, x_GS[4], 'b*')
plt.plot(idx, x_SOR[4], 'g+')
plt.plot(idx, ref, 'ko')
plt.axis([0, 4, -6, 6])
plt.legend(("GS", "SOR", "ref"), loc="lower left")
```


```
plt.axis([0.5, 4.5, 1e-2, 1e1])
for k in range(1,5):
    #print error_GS[k][0]
    plt.semilogy(k, error_GS[k][0], 'g*')    
plt.title("GS error at point 1")
```


```
plt.axis([0.5, 4.5, 1e-2, 1e1])
for k in range(1,5):
    #print error_SOR[k][0]
    plt.semilogy(k, error_SOR[k][0], 'g*')    
plt.title("SOR error at point 1")
```


```

```

d720e96a36cc56429eeed2d553137e54631f0b50
ipynb
Jupyter Notebook
12-iterative-solvers/IterationExample.ipynb
rachelslaybaugh/NE155
5a08229eb11eebdd60e5ec1b4c0d41a541e7d82f
CC-BY-3.0
2015-08-22T05:28:25.000Z
2021-01-15T02:00:39.000Z
12-iterative-solvers/IterationExample.ipynb
rachelslaybaugh/NE155
5a08229eb11eebdd60e5ec1b4c0d41a541e7d82f
CC-BY-3.0
2015-04-01T00:18:04.000Z
2016-10-31T20:14:58.000Z
12-iterative-solvers/IterationExample.ipynb
rachelslaybaugh/NE155
5a08229eb11eebdd60e5ec1b4c0d41a541e7d82f
CC-BY-3.0
2015-01-21T20:12:08.000Z
2020-09-20T08:01:10.000Z
Qwen/Qwen-72B
1. YES
2. YES