Linear systems: iterative methods¶
import numpy as np
And we set the printing format
np.set_printoptions(precision = 8)
np.set_printoptions(suppress = True)
In this lab we tackle the problem of solving linear systems of equations
$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+\ldots+a_{1n}x_{n} & = & b_{1},\\ a_{21}x_{1}+a_{22}x_{2}+\ldots+a_{2n}x_{n} & = & b_{2},\\ \vdots\\ a_{n1}x_{1}+a_{n2}x_{2}+\ldots+a_{nn}x_{n} & = & b_{n}. \end{array} $$We assume that the matrix $A=(a_{ij})$ is non-singular, i.e. det$(A)$ does not vanish, and therefore, the system has a unique solution.
Iterative methods¶
Iterative methods are, in general, more efficient than direct methods for solving sparse large systems of linear equations.
We rewrite the problem as
$$ A\mathbf{x}=\mathbf{b}\Longleftrightarrow\mathbf{x}=B\mathbf{x}+\mathbf{c} $$where $B$ is a matrix of order $n$ and $c$ is a column vector of dimension $n$.
Then, if $\mathbf{x}^{(0)}$ is an initial guess, for $k=1,2,\ldots$
$$ \mathbf{x}^{(k)}=B\mathbf{x}^{(k-1)}+\mathbf{c} $$we generate a sequence of approximate solutions $\left\{ \mathbf{x}^{(m)} \right\}$ that converges to the system solution $\mathbf{x}=\left(x_1,x_2,\ldots,x_n\right)^T$ (when the sequence converges).
Method of Jacobi¶
If we have a system of $4$ equations and $4$ unknowns
$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$and we solve the first equation for the first unknown, the second equation for the second unknown and so on
$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$And we get the values of the first iteration from the initial values as follows
$$ \begin{array}{cll} x_{1}^{{\color{red}{(1)}}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{{\color{red}{(1)}}} & = & \dfrac{b_{2}-a_{21}x_{1}^{(0)}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{{\color{red}{(1)}}} & = & \dfrac{b_{3}-a_{31}x_{1}^{(0)}-a_{32}x_{2}^{(0)}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{{\color{red}{(1)}}}& = & \dfrac{b_{4}-a_{41}x_{1}^{(0)}-a_{42}x_{2}^{(0)}-a_{43}x_{3}^{(0)}}{a_{44}} \end{array} $$If we generalize
Algorithm¶
- Given an initial guess $\mathbf{x}^{(0)}$
- For $k=1,2,\ldots,$MaxIter
- For $i=1,2,\ldots,n$, compute $$x_{i}^{{\color{red}{(k)}}}=\dfrac{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-1)}\right)$$
- If the stopping condition is met, take $\mathbf{x}^{(k)}$ as approximation of the solution.
Write a function jacobi(A,b,tol,maxiter=1000) to solve $Ax=b$
by Jacobi's method.
- Solve
$$
\begin{cases}
5x_{1}+x_{2}-x_{3}-x_{4} & =1\\
x_{1}+4x_{2}-x_{3}+x_{4} & =1\\
x_{1}+x_{2}-5x_{3}-x_{4} & =1\\
x_{1}+x_{2}+x_{3}-4x_{4} & =1
\end{cases}
$$
by Jacobi's method, for an absolute tolerance of
tol = 1.e-6using norm 2. That is, we will use as stopping criteriumnp.linalg.norm(x-xp) < tol, beingxthe approximate solution in the actual iteration andxpthe approximate solution in the previous iteration. Take as initial guess $\mathbf{x^{(0)}=0}.$A = np.array([[5.,1,-1,-1],[1,4,-1,1],[1,1,-5,-1],[1,1,1,-4]]) b = np.array([1.,1,1,1])
- As output, give the solution
xand the number of iterationsnum_iter. - Check that the solution is good comparing it with the solution
xs = np.linalg.solve(A,b). - Check that the function also works well with the $Ax=b$ system with matrices $A$ and $b$ generated as follows
n = 9 A1 = np.diag(np.ones(n))*2 A2 = np.diag(np.ones(n-1),1) A = A1 + A2 + A2.T b = np.concatenate((np.arange(1.,6),np.arange(4,0,-1)))
Hint
To update the xp vector
xp = np.copy(x)
%run Exercise1.py
Method of Gauss-Seidel¶
Algorithm¶
If we have a system of $4$ equations and $4$ unknowns
$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$and we solve the first equation for the first unknown, the second equation for the second unknown and so on
$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$And we get the values of the first iteration from the initial values, but as we get new results, we use them
$$ \begin{array}{cll} x_{1}^{\color{red}{(1)}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{\color{red}{(1)}} & = & \dfrac{b_{2}-a_{21}x_{1}^{\color{red}{(1)}}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{\color{red}{(1)}} & = & \dfrac{b_{3}-a_{31}x_{1}^{\color{red}{(1)}}-a_{32}x_{2}^{\color{red}{(1)}}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{\color{red}{(1)}}& = & \dfrac{b_{4}-a_{41}x_{1}^{\color{red}{(1)}}-a_{42}x_{2}^{\color{red}{(1)}}-a_{43}x_{3}^{\color{red}{(1)}}}{a_{44}} \end{array} $$If we generalize
Algorithm¶
- Given an initial guess $\mathbf{x}^{(0)}$
- For $k=1,2,\ldots,$MaxIter
- For $i=1,2,\ldots,n$, compute $$x_{i}^{\color{red}{(k)}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
- If the stopping condition is met, take $\mathbf{x}^{(k)}$ as approximation of the solution.
Write a function gauss_seidel(A,b,tol,maxiter=1000) to solve $Ax=b$
by Gauss-Seidel's method.
- Solve the same systems as in the previous exercise using Gauss-Seidel's method with the same tolerance and stopping criterium.
- Print the same outputs as in the previous exercise and compare the number of iterations.
%run Exercise2.py
Successive Over Relaxation (SOR)¶
Algorithm¶
- Choose an initial guess $\mathbf{x}^{(0)}$
- For $k=1,2,\ldots,MaxIter$
- For $i=1,2,\ldots,n$ compute $$ x_{i}^{\color{red}{(k)}}=\frac{\omega}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)+(1-\omega) \, x_i^{(k-1)} $$
- If the stopping criterium is satisfied, take $\mathbf{x}^{(k)}$ as approximated solution
If it converges, then $\omega \in [0,2].$
Write a function SOR(A,b,w,tol,maxiter=1000) to solve $Ax=b$
by SOR method.
- Solve the second system as in the previous exercise using the SOR method with with $\omega = 1.5$ and the same tolerance and stopping criterium.
- Print the same outputs as in the previous exercise and compare the number of iterations.
%run Exercise3.py
More exercises¶
SOR: study of $\omega$¶
Write a script that finds the optimum value of $\omega$ to solve the system of Exercise 3 with SOR method. Use maxiter = 1000 and tol = 1.e-12. Find and plot the number of iterations for the values of $\omega$ w = np.arange(0.1,2.1,0.01).
Hint
- Use
k = np.argmin(num_iter)to find the index of the optimum value.
%run Exercise4.py
Tridiagonal systems¶
In the second system, the coefficient matrix is tridiagonal. This kind of sparse matrices (with many zeros) are common in physics problems (heat equation, wave equation) that are solved with finite differences method.
For this reason, the algorithms seen in this lab are adapted to this particular case. For instance, instead of storing the full matrix $A$ of the second system as a square matrix, we store only the diagonals in a rectangular matrix $A_r$
n = 9
Ar = np.zeros((n,3))
Ar[:,0] = np.concatenate((np.array([0]),np.ones((n-1),)))
Ar[:,1] = np.ones((n),)*2
Ar[:,2] = np.concatenate((np.ones((n-1),),np.array([0])))
b = np.concatenate((np.arange(1,6),np.arange(4,0,-1)))*1.
Modify the scripts of the previous functions for them to have as coefficient matrix, the matrix $A$ as $A_r$.
- Solve the second systems of the previous exercise with the same tolerance and stopping criterium.
- Print the same outputs as in the previous exercises.
- Use $\omega = 1.5$ for the SOR method.
%run Exercise5.py