# Solving Linear Systems: Iterative Methods
<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a><br />This notebook by Xiaozhou Li is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  
All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

## General Form 
For solving the linear system
$$
    Ax = b,
$$
with the exact solution $x^{*}$.  The general form based on the fixed point interation:
\begin{equation}
     \begin{split}
          x^{(0)} & = \text{initial guess} \\
          x^{(k+1)} & = g(x^{(k)}) \quad k = 0,1,2,\ldots,
     \end{split}
\end{equation}
where
$$
    g(x) = x - C(Ax - b).
$$
Difficult: find a matrix $C$ such that 
$$
    \lim\limits_{k\rightarrow\infty}x^{(k)} = x^{*}
$$
and the algorithm needs to be converge fast and economy. 

**Example 1**
\begin{equation*}
    A = \left[\begin{array}{ccc} 9& -1 & -1 \\ -1 & 10 & -1 \\ -1 & -1& 15\end{array}\right],\quad b = \left[\begin{array}{c} 7 \\ 8 \\ 13\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[1, 1, 1]}^T$


```python
import numpy as np
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
plt.rcParams['axes.unicode_minus']=False
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from IPython.display import clear_output, display
```


```python
def IterC(A, b, C, x0, x_star, iters):
    x = np.copy(x0)
    print ('Iteration No.     Numerical Solution     Max norm error ')
    print (0, x, np.linalg.norm(x_star-x, np.inf))
    for i in range(iters):
        x = x + np.dot(C, b - np.dot(A,x))
        print (i+1, x, np.linalg.norm(x_star-x,np.inf))
```


```python
A = np.array([[9., -1., -1.],[-1.,10.,-1.],[-1.,-1.,15.]])
b = np.array([7.,8.,13.])
```

**Naive Choice** 

Choosing $C = I$, then 
$$g(x) = x - (Ax - b),$$
and the fixed-point iteration
$$x^{(k+1)}  = (I - A)x^{(k)} + b \quad k = 0,1,2,\ldots. $$
Let the intial guess $x_0 = [0, 0, 0]^T$.


```python
C = np.eye(3)
x0 = np.zeros(3)
x_star = np.array([1.,1.,1.])

w = interactive(IterC, A=fixed(A), b=fixed(b), C=fixed(C), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Best Choice (theoretically)** 

Choosing $C = A^{-1}$, then 
$$g(x) = x - A^{-1}(Ax - b),$$
and the fixed-point iteration
$$x^{(k+1)}  = A^{-1}b \quad k = 0,1,2,\ldots. $$

* It equals to solve $Ax = b$ directly.
* However, it gives a hint that $C$ should be close to $A^{-1}$

## Some Classic Iterative Methods

Let $D$ denote the main diagonal of $A$, $L$ denote the lower triangle of $A$ (entries below the main diagonal), and $U$ denote the upper triangle (entries above the main diagonal). Then, we have $A = L + D + U$
\begin{equation}
A = \begin{pmatrix}{a_{11}} & {a_{12}} & {\cdots} & {a_{1 n}} \\ {a_{21}} & {a_{22}} & {\cdots} & {a_{2 n}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {a_{n 1}} & {a_{n 2}} & {\cdots} & {a_{n n}}\end{pmatrix},\,\,D = \begin{pmatrix}{a_{11}} & {} & {} & {} \\ {} & {a_{22}} & {} & {} \\ {} & {} & {\ddots} & {} \\ {} & {} & {} & {a_{n n}}\end{pmatrix},\,\, L= \begin{pmatrix}{0} & {} & {} & {} \\ {a_{21}} & {\ddots} & {} & {} \\ {\vdots} & {\ddots} & {0} & {} \\ {a_{n 1}} & {\cdots} & {a_{n, n-1}} & {0}\end{pmatrix},\,\, U= \begin{pmatrix}{0} & {a_{12}} & {\cdots} & {a_{1n}} \\ {} & {0} & {\ddots} & {\vdots} \\ {} & {} & {\ddots} & {a_{n-1,n}} \\ {} & {} & {} & {0}\end{pmatrix}
\end{equation}

**First Approach** 

Choosing $C = \text{diag}(A)^{-1} = D^{-1}$, then 
$$g(x) = x - D^{-1}(Ax - b),$$
and the fixed-point iteration
$$Dx^{(k+1)}  = -(L + U)x^{(k)} + b \quad k = 0,1,2,\ldots. $$


```python
C = np.diag(1./np.diag(A))
x0 = np.zeros(np.size(b))
x_star = np.array([1.,1.,1.])
```


```python
w = interactive(IterC, A=fixed(A), b=fixed(b), C=fixed(C), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


## Jacobi Method
### Matrix Form:
$$
    x^{(k+1)} = x^{(k)} - D^{-1}(Ax^{(k)} - b)
$$
or
$$
    Dx^{(k+1)} = b - (L+U)x^{(k)}
$$
### Algorithm
$$
    x^{(k+1)}_i = \frac{b_i - \sum\limits_{j < i}a_{ij}x^{(k)}_j - \sum\limits_{j > i}a_{ij}x^{(k)}_j}{a_{ii}}
$$


```python
def Jacobi(A, b, x0, x_star, iters):
    x_old = np.copy(x0)
    x_new = np.zeros(np.size(x0))
    print (0, x_old, np.linalg.norm(x_star-x_old,np.inf))
    for k in range(iters):
        for i in range(np.size(x0)): 
            x_new[i] = (b[i] - np.dot(A[i,:i],x_old[:i]) - np.dot(A[i,i+1:],x_old[i+1:]))/A[i,i]
        print (k+1, x_new, np.linalg.norm(x_star-x_new,np.inf))
        x_old = np.copy(x_new)
```


```python
w = interactive(Jacobi, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Second Approach** 

Let $D$ denote the main diagonal of $A$, $L$ denote the lower triangle of $A$ (entries below the main diagonal), and $U$ denote the upper triangle (entries above the main diagonal). Then $A = L + D + U$

Choosing $C = (L + D)^{-1}$, then 
$$g(x) = x - (L + D)^{-1}(Ax - b),$$
and the fixed-point iteration
$$(L + D)x^{(k+1)}  = Ux^{(k)} + b \quad k = 0,1,2,\ldots. $$


```python
def GS(A, b, x0, x_star, iters):
    x = np.copy(x0)
    print (0, x, np.linalg.norm(x_star-x,np.inf))
    for k in range(iters):
        for i in range(np.size(x0)): 
            x[i] = (b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i+1:],x[i+1:]))/A[i,i]
        print (k+1, x, np.linalg.norm(x_star-x,np.inf))
```


```python
w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


## Gauss-Seidel Method
### Algorithm
$$
    x^{(k+1)}_i = \frac{b_i - \sum\limits_{j < i}a_{ij}x^{(k+1)}_j - \sum\limits_{j > i}a_{ij}x^{(k)}_j}{a_{ii}}
$$
### Matrix Form:
$$
    x^{(k+1)} = x^{(k)} - (L+D)^{-1}(Ax^{(k)} - b)
$$
or
$$
    (L+D)x^{(k+1)} = b - Ux^{(k)}
$$

**Example 2**
\begin{equation*}
    A = \left[\begin{array}{ccc} 3& 1 & -1 \\ 2 & 4 & 1 \\ -1 & 2& 5\end{array}\right],\quad b = \left[\begin{array}{c} 4 \\ 1 \\ 1\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[2, -1, 1]}^T$


```python
A = np.array([[3, 1, -1],[2,4,1],[-1,2,5]])
b = np.array([4,1,1])

x0 = np.zeros(np.size(b))
x_star = np.array([2.,-1.,1.])
```


```python
w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=40,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=40), Output()), _dom_classes=('widget-intera…


**Example 3**
\begin{equation*}
    A = \left[\begin{array}{ccc} 1& 2 & -2 \\ 1 & 1 & 1 \\ 2 & 2& 1\end{array}\right],\quad b = \left[\begin{array}{c} 7 \\ 8 \\ 13\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[-3, 8, 3]}^T$


```python
A = np.array([[1, 2, -2],[1,1,1],[2,2,1]])
b = np.array([7,8,13])

#x0 = np.zeros(np.size(b))
x0 = np.array([-1,1,1])
x_star = np.array([-3.,8.,3.])

w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
B = np.eye(3) - np.dot(np.diag(1./np.diag(A)),A)
print(B)
print (np.linalg.eig(B))
```

    [[ 0. -2.  2.]
     [-1.  0. -1.]
     [-2. -2.  0.]]
    (array([-6.16610918e-06+1.06799587e-05j, -6.16610918e-06-1.06799587e-05j,
            1.23322184e-05+0.00000000e+00j]), array([[-0.57734849-6.16605800e-06j, -0.57734849+6.16605800e-06j,
             0.57735383+0.00000000e+00j],
           [ 0.57735027+3.08300999e-06j,  0.57735027-3.08300999e-06j,
            -0.57735027+0.00000000e+00j],
           [ 0.57735205+0.00000000e+00j,  0.57735205-0.00000000e+00j,
            -0.57734671+0.00000000e+00j]]))


**Example 4**
\begin{equation*}
    A = \left[\begin{array}{cc} 1& 2 \\ 3 & 1 \end{array}\right],\quad b = \left[\begin{array}{c} 5 \\ 5\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[1, 2]}^T$


or
\begin{equation*}
    A = \left[\begin{array}{cc} 3& 1 \\ 1 & 2 \end{array}\right],\quad b = \left[\begin{array}{c} 5 \\ 5\end{array}\right],
\end{equation*}


```python
#A = np.array([[1, 2],[3,1]])
A = np.array([[3, 1],[1,2]])
b = np.array([5,5])

#x0 = np.zeros(np.size(b))
x0 = np.array([0,0])
x_star = np.array([1.,2.,])

w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Example 5**
Are Jacobi iteration and Gauss-Seidel iteration convergent for the following equations?
\begin{equation*}
    A_1 = \left[\begin{array}{ccc} 3& 0 & 4 \\ 7 & 4 & 2 \\ -1 & 1 & 2\end{array}\right],\quad A_2 = \left[\begin{array}{ccc} -3& 3 & -6 \\ -4 & 7 & -8 \\ 5 & 7 & -9\end{array}\right],
\end{equation*}

* Consider the **spectral radius** of the iterative matrix
* $B_J = -D^{-1}(L+U)$ and $B_{GS} = -(L+D)^{-1}U$


```python
def Is_Jacobi_Gauss(A):
    L = np.tril(A,-1)
    U = np.triu(A,1)
    D = np.diag(np.diag(A))
    
    B_J = np.dot(np.diag(1./np.diag(A)), L+U)
    B_GS = np.dot(np.linalg.inv(L+D),U)
    
    rho_J = np.linalg.norm(np.linalg.eigvals(B_J), np.inf)
    rho_GS = np.linalg.norm(np.linalg.eigvals(B_GS), np.inf)
    
    print ("Spectral Radius")
    print ("Jacobi: ", rho_J)
    print ("Gauss Sediel: ", rho_GS)
    
A1 = np.array([[3, 0, 4],[7, 4, 2], [-1,1,2]])
A2 = np.array([[-3, 3, -6], [-4, 7, -8], [5, 7, -9]])
Is_Jacobi_Gauss(A2)
```

    Spectral Radius
    Jacobi:  0.8133091054692768
    Gauss Sediel:  1.1111111111111105


## Successive Over-Relaxation (SOR)
### Algorithm
$$
    x^{(k+1)}_i = x^{(k)} + \omega \frac{b_i - \sum\limits_{j < i}a_{ij}x^{(k+1)}_j - \sum\limits_{j \geq i}a_{ij}x^{(k)}_j}{a_{ii}}
$$
### Matrix Form:
$$
    x^{(k+1)} = x^{(k)} - \omega(\omega L+D)^{-1}(Ax^{(k)} - b)
$$
or
$$
    (\omega L+D)x^{(k+1)} = ((1-\omega)D - \omega U)x^{(k)} + \omega b
$$


```python
def SOR(A, b, x0, x_star, omega, iters):
    x = np.copy(x0)
    print (0, x, np.linalg.norm(x_star-x,np.inf))
    for k in range(iters):
        for i in range(np.size(x0)): 
            x[i] = x[i] + omega*(b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i:],x[i:]))/A[i,i]
        print (k+1, x, np.linalg.norm(x_star-x,np.inf))
```


```python
def SOR2(A, b, x0, x_star, omega, iters):
    x = np.copy(x0)
    for k in range(iters):
        for i in range(np.size(x0)): 
            x[i] = x[i] + omega*(b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i:],x[i:]))/A[i,i]
    return (np.linalg.norm(x_star-x,np.inf))
```


```python
def SOR3(A, b, x0, x_star, omega, iters):
    x = np.copy(x0)
    print (0, np.linalg.norm(x_star-x,np.inf))
    for k in range(iters):
        for i in range(np.size(x0)): 
            x[i] = x[i] + omega*(b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i:],x[i:]))/A[i,i]
        print (k+1, np.linalg.norm(x_star-x,np.inf))
```


```python
A = np.array([[9., -1., -1.],[-1.,10.,-1.],[-1.,-1.,15.]])
b = np.array([7.,8.,13.])
```


```python
x0 = np.array([0.,0.,0.])
x_star = np.array([1.,1.,1.])
omega = 1.01

w = interactive(SOR, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), omega=fixed(omega), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Example 6**
\begin{equation*}
    A = \left[\begin{array}{ccc} 2& -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1& 2\end{array}\right],\quad b = \left[\begin{array}{c} 1 \\ 0 \\ 1.8\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[1.2, 1.4, 1.6]}^T$


```python
A = np.array([[2, -1, 0],[-1, 2, -1], [0, -1, 2]])
b = np.array([1., 0, 1.8])

x0 = np.array([1.,1.,1.])
x_star = np.array([1.2,1.4,1.6])
omega = 1.2

w = interactive(SOR, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), omega=fixed(omega), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
num = 21
omega = np.linspace(0.8, 1.8, num)
err1 = np.zeros(num)
for i in range(num):
    err1[i] = SOR2(A, b, x0, x_star, omega[i], 10)
print (err1)
plt.plot(omega, np.log10(err1), 'o')
```

**Example 7**
\begin{equation*}
    A = \left[\begin{array}{cccc} -4& 1 & 1 & 1 \\ 1 & -4 & 1 & 1 \\ 1 & 1& -4 &1 \\ 1 & 1 &1 & -4\end{array}\right],\quad b = \left[\begin{array}{c} 1 \\ 1 \\ 1 \\ 1\end{array}\right],
\end{equation*}
has the exact solution $x^{*} = {[-1, -1, -1, -1]}^T$


```python
A = np.array([[-4, 1, 1, 1],[1, -4, 1, 1], [1, 1, -4, 1], [1, 1, 1, -4]])
b = np.array([1, 1, 1, 1])

x0 = np.zeros(np.size(b))
x_star = np.array([-1,-1,-1,-1])
omega = 1.25

w = interactive(SOR, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), omega=fixed(omega), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
w = interactive(GS, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=100,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters'), Output()), _dom_classes=('widget-interact',))



```python
num = 21
omega = np.linspace(0.8, 1.8, num)
err1 = np.zeros(num)
for i in range(num):
    err1[i] = SOR2(A, b, x0, x_star, omega[i], 10)
print (err1)
plt.plot(omega, np.log10(err1), 'o')
```

**Example 8**
\begin{equation*}
    A=\begin{pmatrix}{3} & {-1} & {0} & 0 & 0 & \frac{1}{2}  \\ {-1} & {3} & {-1} & {0} & \frac{1}{2} & 0\\ {0} & {-1} & {3} & {-1} & {0} & 0 \\ 0& {0} & {-1} & {3} & {-1} & {0} \\ {0} & \frac{1}{2} & {0} & {-1} & {3} & {-1} \\ \frac{1}{2} & {0} & 0 & 0 & {-1} & {3}\end{pmatrix},\,\,b=\begin{pmatrix}\frac{5}{2} \\ \frac{3}{2} \\ 1 \\ 1 \\ \frac{3}{2}    \\ \frac{5}{2}  \end{pmatrix}
\end{equation*}
has the exact solution $x^{*} = {[1, 1, 1, 1, 1, 1]}^T$


```python
n0 = 6
A = 3*np.eye(n0) - np.diag(np.ones(n0-1),-1) - np.diag(np.ones(n0-1),+1)
for i in range(n0):
    if (abs(n0-1 - 2*i) > 1):
        A[i, n0-1-i] = - 1/2
print (A)
x_star = np.ones(n0)
b = np.dot(A, x_star)

x0 = np.zeros(np.size(b))
omega = 1.25

w = interactive(SOR, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), omega=fixed(omega), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```

    [[ 3.  -1.   0.   0.   0.  -0.5]
     [-1.   3.  -1.   0.  -0.5  0. ]
     [ 0.  -1.   3.  -1.   0.   0. ]
     [ 0.   0.  -1.   3.  -1.   0. ]
     [ 0.  -0.5  0.  -1.   3.  -1. ]
     [-0.5  0.   0.   0.  -1.   3. ]]



    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
num = 21
omega = np.linspace(0.8, 1.8, num)
err1 = np.zeros(num)
for i in range(num):
    err1[i] = SOR2(A, b, x0, x_star, omega[i], 10)
print (err1)
plt.plot(omega, np.log10(err1), 'o')
```


```python
w = interactive(Jacobi, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), iters=widgets.IntSlider(min=0,max=100,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters'), Output()), _dom_classes=('widget-interact',))


## Sparse Matrix Computations
A coefficient matrix is called sparse if many of the matrix entries are known to be zero. Often, of the $n^2$ eligible entries in a sparse matrix, only $\mathcal{O}(n)$ of them are nonzero. A full matrix is the opposite, where few entries may be assumed to be zero. 

**Example 9**
Consider the $n$-equation version of
\begin{equation*}
    A=\begin{pmatrix}{3} & {-1} & {0} & 0 & 0 & \frac{1}{2}  \\ {-1} & {3} & {-1} & {0} & \frac{1}{2} & 0\\ {0} & {-1} & {3} & {-1} & {0} & 0 \\ 0& {0} & {-1} & {3} & {-1} & {0} \\ {0} & \frac{1}{2} & {0} & {-1} & {3} & {-1} \\ \frac{1}{2} & {0} & 0 & 0 & {-1} & {3}\end{pmatrix}, 
\end{equation*}
has the exact solution $x^{*} = {[1, 1,\ldots, 1]}^T$ and $b = A x^{*}$

* First, let us have a look about the matrix $A$


```python
n0 = 10
A = 3*np.eye(n0) - np.diag(np.ones(n0-1),-1) - np.diag(np.ones(n0-1),+1)
for i in range(n0):
    if (abs(n0-1 - 2*i) > 1):
        A[i, n0-1-i] = - 1/2
plt.spy(A)
#plt.show()
```

* How about the $PA = LU$ for the above matrix $A$?
* Are the $L$ and $U$ matrices still sparse? 


```python
import scipy.linalg
#P, L, U = scipy.linalg.lu(A)

#plt.spy(L)
#plt.show()
```

Gaussian elimination applied to a sparse matrix usually causes **fill-in**, where the coefficient matrix changes from sparse to full due to the necessary row operations. For this reason, the efficiency of Gaussian elimination and its $PA = LU$ implementation become questionable for sparse matrices, leaving iterative methods as a feasible alternative.

* Let us solve it with SOR method


```python
x_star = np.ones(n0)
b = np.dot(A, x_star)

x0 = np.zeros(np.size(b))
omega = 1.25
w = interactive(SOR3, A=fixed(A), b=fixed(b), x0=fixed(x0), x_star=fixed(x_star), omega=fixed(omega), iters=widgets.IntSlider(min=0,max=200,value=0, step=10))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=200, step=10), Output()), _dom_classes=('wid…


## Application for Solving Laplace's Equation

### Laplace's equation
Consider the Laplace's equation given as 
$$
    \nabla^2 u = 0,\quad\quad (x,y) \in D,
$$
where $\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$, and the boundary conditions are given as


### Finite Difference Approximation
Here, we use a rectangular grid $(x_i,y_j)$, where
$$
    x_i = i\Delta x, \,\,\text{for }\, i = 0,1,\ldots,N+1;\quad y_j = j\Delta y,\,\,\text{for }\, j = 0,1,\ldots,M+1.
$$
Five-points scheme:
$$
    -\lambda^2 u_{i+1,j} + 2(1+\lambda^2)u_{i,j} - \lambda^2u_{i-1,j} - u_{i,j+1} - u_{i,j-1} = 0,\quad\text{for}\,\, i = 1,\ldots,N,\,\, j = 1,\ldots,M，
$$
where $\lambda = \frac{\Delta y}{\Delta x}$. The boundary conditions are 
- $x = 0: u_{0,j} = g_L(y_j), \quad\text{for }\, j = 1,\ldots,M$,
- $x = a: u_{N+1,j} = g_R(y_j), \quad\text{for }\, j = 1,\ldots,M$,
- $y = 0: u_{i,0} = g_B(x_i), \quad\text{for }\, i = 1,\ldots,N$,
- $y = b: u_{i,M+1} = g_T(x_i), \quad\text{for }\, i = 1,\ldots,N$.


```python
def generate_TD(N, dx, dy):
    T = np.zeros([N,N])
    a = - (dy/dx)**2
    b = 2*(1 - a)
    for i in range(N):
        T[i,i] += b
        if (i < N-1):
            T[i,i+1] += a
        if (i > 0):
            T[i,i-1] += a
    D = -np.identity(N)
    return T, D

def assemble_matrix_A(dx, dy, N, M):
    T, D = generate_TD(N, dx, dy)
    A = np.zeros([N*M, N*M])
    for j in range(M):
        A[j*N:(j+1)*N,j*N:(j+1)*N] += T
        if (j < M-1):
            A[j*N:(j+1)*N,(j+1)*N:(j+2)*N] += D
        if (j > 0):
            A[j*N:(j+1)*N,(j-1)*N:j*N] += D
    return A
```


```python
N = 4
M = 4
dx = 1./(N+1)
dy = 1./(M+1)
T, D = generate_TD(N, dx, dy)
#print (T)
A = assemble_matrix_A(dx, dy, N, M)
#print (A)
plt.spy(A)
plt.show()
```


```python
# Set boundary conditions
def gL(y):
    return 0.

def gR(y):
    return 0.

def gB(x):
    return 0.

def gT(x):
    return 1.
    #return x*(1-x)*(4./5-x)*np.exp(6*x)
    
def assemble_vector_b(x, y, dx, dy, N, M, gL, gR, gB, gT):
    b = np.zeros(N*M)
    # Left BCs
    for j in range(M):
        b[(j-1)*N] += (dy/dx)**2*gL(y[j+1]) 
    
    # Right BCs
    # b +=
    
    # Bottom BCs
    # b +=
    
    # Top BCs:
    for i in range(N):
        b[(M-1)*N+i] += gT(x[i+1])
    return b
```


```python
from mpl_toolkits import mplot3d
from mpl_toolkits.mplot3d import axes3d
```


```python
def Laplace_solver(a, b, N, M, gL, gR, gB, gT):
    dx = b/(M+1)
    dy = a/(N+1)
    x = np.linspace(0, a, N+2)
    y = np.linspace(0, b, M+2)
    
    A = assemble_matrix_A(dx, dy, N, M)
    b = assemble_vector_b(x, y, dx, dy, N, M, gL, gR, gB, gT)
    
    v = np.linalg.solve(A,b)
    
    # add boundary points + plotting
    u = np.zeros([(N+2),(M+2)])
    #u[1:(N+1),1:(M+1)] = np.reshape(v, (N, M))
    # Top BCs
    for i in range(N+2):
        u[i,M+1] = gT(x[i])
    u = np.transpose(u)
    u[1:(M+1),1:(N+1)] = np.reshape(v, (M, N))

    
    X, Y = np.meshgrid(x, y)
    #Z = np.sin(2*np.pi*X)*np.sin(2*np.pi*Y)

    fig = plt.figure()
    #ax = plt.axes(projection='3d')
    ax = fig.add_subplot(1, 1, 1, projection='3d')

    ax.plot_surface(X, Y, u, rstride=1, cstride=1,
                    cmap='viridis', edgecolor='none')
    ax.set_title('surface')
    plt.show()
    
Laplace_solver(1, 1, 40, 40, gL, gR, gB, gT)
```


```python
def Jacobi_tol(A, b, x0, tol):
    x_old = np.copy(x0)
    x_new = np.zeros(np.size(x0))
    for i in range(np.size(x0)): 
        x_new[i] = (b[i] - np.dot(A[i,:i],x_old[:i]) - np.dot(A[i,i+1:],x_old[i+1:]))/A[i,i]
    iters = 1
    while ((np.linalg.norm(x_new-x_old,np.inf)) > tol):
        x_old = np.copy(x_new)
        for i in range(np.size(x0)): 
            x_new[i] = (b[i] - np.dot(A[i,:i],x_old[:i]) - np.dot(A[i,i+1:],x_old[i+1:]))/A[i,i]
        iters += 1
    return x_new, iters
```


```python
def GS_tol(A, b, x0, tol):
    x_old = np.copy(x0)
    x = np.copy(x0)
    for i in range(np.size(x0)): 
        x[i] = (b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i+1:],x[i+1:]))/A[i,i]
    iters = 1
    while ((np.linalg.norm(x-x_old,np.inf)) > tol):
        x_old = np.copy(x)
        for i in range(np.size(x0)): 
            x[i] = (b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i+1:],x[i+1:]))/A[i,i]
        iters += 1
    return x, iters
```


```python
def SOR_tol(A, b, x0, omega, tol):
    x_old = np.copy(x0)
    x = np.copy(x0)
    for i in range(np.size(x0)): 
            x[i] = x[i] + omega*(b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i:],x[i:]))/A[i,i]
    iters = 1
    while ((np.linalg.norm(x-x_old,np.inf)) > tol):
        x_old = np.copy(x)
        for i in range(np.size(x0)): 
            x[i] = x[i] + omega*(b[i] - np.dot(A[i,:i],x[:i]) - np.dot(A[i,i:],x[i:]))/A[i,i]
        iters += 1
    return x, iters
```


```python
def CG_tol(A, b, x0, x_star, tol):
    r_new = b - np.dot(A, x0) 
    r_old = np.copy(np.size(x0))
    d_old = np.zeros(np.size(x0))
    x = np.copy(x0)
    iters = 0
    while ((np.linalg.norm(x-x_star,np.inf)) > tol):
        if (iters == 0):
            d_new = np.copy(r_new)
        else:
            beta = np.dot(r_new,r_new)/np.dot(r_old,r_old)
            d_new = r_new + beta*d_old
        Ad = np.dot(A, d_new)
        alpha = np.dot(r_new,r_new)/np.dot(d_new,Ad)
        x += alpha*d_new
        d_old = d_new
        r_old = r_new
        r_new = r_old - alpha*Ad
        iters += 1
    return x, iters
```


```python
def Iterative_solver(a, b, N, M, gL, gR, gB, gT, tol):
    dx = b/(M+1)
    dy = a/(N+1)
    x = np.linspace(0, a, N+2)
    y = np.linspace(0, b, M+2)
    
    A = assemble_matrix_A(dx, dy, N, M)
    b = assemble_vector_b(x, y, dx, dy, N, M, gL, gR, gB, gT)
    
    v = np.linalg.solve(A,b)
    #tol = 1.e-8
    v0 = np.zeros(np.size(b))
    #v_J, iters = Jacobi_tol(A, b, v0, tol)
    #print ("Jacobi Method:  %4d    %7.2e" %(iters, np.linalg.norm(v - v_J, np.inf)))
    
    #v_GS, iters = GS_tol(A, b, v0, tol)
    #print ("Gauss Seidel :  %4d    %7.2e" %(iters, np.linalg.norm(v - v_GS, np.inf)))
    
    omega = 2./(1 + np.sin(np.pi*dx))
    print ("omega = ", omega)
    v_SOR, iters = SOR_tol(A, b, v0, omega, tol)
    print ("SOR Method   :  %4d    %7.2e" %(iters, np.linalg.norm(v - v_SOR, np.inf)))
    
    v_CG, iters = CG_tol(A, b, v0, v, tol)
    print ("CG Method    :  %4d    %7.2e" %(iters, np.linalg.norm(v - v_CG, np.inf)))
```


```python
Iterative_solver(1, 1, 20, 20, gL, gR, gB, gT, 1.e-4)
```

    omega =  1.7405800107385732
    SOR Method   :    44    1.03e-04
    CG Method    :    37    8.38e-05



```python

```

647703655f420fdec0b1c90aca503f6c7492687b
ipynb
Jupyter Notebook
LinearSys_IterMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
2021-11-06T02:11:33.000Z
2021-11-06T02:11:33.000Z
LinearSys_IterMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
LinearSys_IterMethod.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
Qwen/Qwen-72B
1. YES
2. YES