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.
With from __future__ import division we force the result of division with / to be always float. If we do not include this, the division of two integers with / would be an integer.
We also import the modules numpy and pprint. The pprint module provides a capability to “pretty-print” Python data.
from __future__ import division
import numpy as np
import pprint
For the system $$ L\mathbf{x}=\mathbf{b}, $$ with $L$ of the form $$ L=\left(\begin{array}{ccccc} l_{11} & 0 & 0 & \ldots & 0\\ l_{21} & l_{22} & 0 & \ldots & 0\\ l_{31} & l_{32} & l_{33} & \ldots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ l_{n1} & l_{n2} & l_{n3} & \ldots & l_{nn} \end{array}\right) $$ the $i$-th equation is, for $i=1,2,\ldots,n$, $$ l_{i1}x_{1}+l_{i2}x_{2}+\cdots+l_{ii}x_{i}=b_{i}. $$
Therefore, for $i=1,2,\ldots,n$, if $l_{ii}$ does not vanish, then $$ x_{i}=\frac{b_{i}-l_{i1}x_{1}-\cdots-l_{ii-1}x_{i-1}}{l_{ii}}= \frac{1}{l_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}l_{ij}x_{j}\right). $$
If $l_{ii}= 0$, exit.
Else, set $$ x_{i}=\frac{1}{l_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}l_{ij}x_{j}\right). $$
Write a function subs_forward(L,b) whose output $x$ is the solution of the system $Lx=b$, being L a lower triangular matriz and b a column matrix.
Use it to solve the system
$$\left(\begin{array}{rrr} 2 & 0 & 0\\ 1 & 2 & 0\\ 1 & 1 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 2\\ -1\\ 0 \end{array}\right) $$Check the solution computing $Lx$ and comparing it with $b$.
# define matrices A and b using numpy arrays
L = np.array([[2, 0, 0], [1, 2, 0], [1, 1, 2]])
b = np.array([2, -1, 0])
Check that the function also works well with the $Lx=b$ system with matrices $L$ and $b$ generated as follows
n = 5
np.random.seed(1) # for the random numbers to be repeatable
L = np.random.random((n,n)) # generate a random nxn matrix
L = np.tril(L) # remove elements over diagonal
b = np.random.random((n,)) # generate a random vector of dimension n
print 'L'
pprint.pprint(L)
print 'b'
pprint.pprint(b)
%run Exercise1.py
Similarly, for the system
$$ U\mathbf{x}=\mathbf{b}' $$with $U$ of the form
$$ U=\left(\begin{array}{ccccc} u_{11} & u_{12} & u_{13} & \ldots & u_{1n}\\ 0 & u_{22} & u_{23} & \ldots & u_{2n}\\ 0 & 0 & u_{33} & \ldots & u_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & 0 & u_{nn} \end{array}\right) $$the $i$-th equation is, for $i=1,2,\ldots,n$,
$$ u_{ii}x_{i}+u_{ii+1}x_{i+1}+\cdots+u_{in}x_{n}=b_{i} $$Therefore, for $i=1,2,\ldots,n$, if $l_{ii}$ does not vanish then
$$ x_{i}=\frac{b_{i}-u_{ii+1}x_{i+1}-\cdots-u_{in}x_{n}}{u_{ii}}= \frac{1}{u_{ii}}\left(b_{i}-\sum_{j=i+1}^{n}u_{ij}x_{j}\right), $$
- If $u_{ii}=0$, exit
- Else, set
Write a function subs_backward(U,b) whose output $x$ is the solution of the system $Ux=b$, being U a upper triangular matriz and b a column matrix.
Use it to solve the system
$$\left(\begin{array}{rrr} 2 & 1 & 1\\ 0 & 2 & 1\\ 0 & 0 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 9\\ 4\\ 4 \end{array}\right) $$n = 5
np.random.seed(2) # for the random numbers to be repeatable
U = np.random.random((n,n)) # generate a random nxn matrix
U = np.triu(U) # remove elements under diagonal
b = np.random.random((n,)) # generate a random vector of dimension n
print 'U'
pprint.pprint(U)
print 'b'
pprint.pprint(b)
%run Exercise2.py
Idea of Gauss' method is to transform the original system into an upper triangular system with the same solution. To do this, we may use the following operations:
Once the matrix is triangular we can solve the system with backward elimination.
In general, it is necessary to allow the swapping of rows, but we will code first the algorithm without it.
- If $a_{kk}=0$, exit
- Else, for $i=k+1,\ldots,n,\quad \quad b_{i}=b_{i}-\dfrac{a_{ik}}{a_{kk}}b_{k}. $
- For $i=k+1,\ldots,n, \quad a_{ij}=a_{ij}-\dfrac{a_{ik}}{a_{kk}}a_{ik}. $
Write a function triang(A,b) whose output are the transformed upper triangular matrix U and vector c of a given matrix A and right hand side of the system b.
Use it for the matrices of the system $Ax=b$
$$\left(\begin{array}{rrr} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 2\\ 4\\ 6 \end{array}\right) $$U and c.n = 5
np.random.seed(3) # for the random numbers to be repeatable
A = np.random.random((n,n)) # generate a random nxn matrix
b = np.random.random((n,)) # generate a random vector of dimension n
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Exercise3.py
We express $A$ as a matrix multiplication of two matrices: a Lower triangular matrix, $L$, and an Upper triangular matrix, $U$,
$$ A=LU. $$Consider the system
$$ A\mathbf{x}=\mathbf{b}, $$where $A$ admits the $LU$ decomposition. Then
$$ A\mathbf{x}=\mathbf{b}\Longleftrightarrow LU\mathbf{x}=\mathbf{b} $$Solving $Ax=b$ is done in two steps
For the time being, we are not considering swapping rows.
Use the $LU$ decomposition of $A$ to solve the system $Ax=b$, with
$$ A=\left(\begin{array}{rrr} 2 & 1 & 1\\ 1 & 2 & 2\\ 1 & 1 & 2 \end{array}\right)\quad\quad b=\left(\begin{array}{c} 2\\ 4\\ 6 \end{array}\right) $$LU(A) that gives the $LU$ decomposition of $A$.subs_backward and subs_forward functions to solve the resulting systems.n = 5
np.random.seed(4) # for the random numbers to be repeatable
A = np.random.random((n,n)) # generate a random nxn matrix
b = np.random.random((n,)) # generate a random vector of dimension n
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Exercise4.py
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} $$
- For $i=1,2,\ldots,n$ compute $$ x_{i}^{(k)}=\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-1)}) $$
- If the stopping criterium is satisfied, take $\mathbf{x}^{(k)}$ as approximated solution
Write a function jacobi(A,b,tol) to solve $Ax=b$
by Jacobi's method.
tol = 1.e-6 using norm 2. That is, we will use as stopping criterium np.linalg.norm(x-xp) < tol, being x the approximate solution in the actual iteration and xp the approximate solution in the previous iteration. Take as initial guess $\mathbf{x^{(0)}=0}.$x and the number of iterations num_iter.xs = np.linalg.solve(A,b).n = 5
np.random.seed(5)
A = np.random.random((n,n)) + 2*np.eye(n)
b = np.random.random((n,))
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Exercise5.py
- For $i=1,2,\ldots,n$ compute $$ x_{i}^{(k)}=\frac{1}{a_{ii}}(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{(k)}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}) $$
- If the stopping criterium is satisfied, take $\mathbf{x}^{(k)}$ as approximated solution
Write a function gauss_seidel(A,b,tol) to solve $Ax=b$
by Gauss-Seidel's method.
%run Exercise6.py