Linear systems: direct methods¶
import numpy as np
We set the format
np.set_printoptions(precision = 2)
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.
Backward substitution¶
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), $$Algorithm¶
- $x_{n}=\dfrac{b_{n}}{u_{nn}}$
- For $i=n-1,\ldots,1$
- 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) $$U = np.array([[2., 1, 1], [0, 2, 1], [0, 0, 2]])
b = np.array([9., 4, 4])
- Check the solution computing $Ux$ and comparing it with $b$.
- Check that the function also works well with the $Ux=b$ system with matrices $L$ and $b$ generated as follows
n = 5 np.random.seed(2) U = np.random.random((n,n)) U = np.triu(U) # make zero elements under the diagonal b = np.random.random(n)
Hints
- We can obtain the number of rows and columns of a matrix
Awithn_rows, n_columns = A.shape - The number of elements of
bcan be obtained withn = len(b). - We initialize a vector of zeros
xwith the same shape asbwithx = np.zeros(len(b))orx = np.zeros_like(b)
%run Exercise1.py
Gaussian elimination¶
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:
- We add a row multiplied by a constant to other $r_i\rightarrow r_i +\lambda r_k,\quad i\ne k$
- We can interchange rows.
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.
Algorithm¶
- For $k=1,\ldots,n-1$
- If $a_{kk}=0$, exit
- Else, for $i=k+1,\ldots,n,$
- $f = \dfrac{a_{ik}}{a_{kk}}$
- $b_{i}=b_{i}-f\,b_{k}$
- $a\,\mathrm{row}_i \leftarrow a\,\mathrm{row}_i - f\times a\,\mathrm{row}_k$
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) $$- Print the matrices
Uandc. - Solve the equivalent system $Ux=c$ and check the solution computing $Ax$ and comparing it with $b$.
- Check that the function also works well with the $Ax=b$ system with matrices $A$ and $b$ generated as follows
n = 5 np.random.seed(3) A = np.random.random((n,n)) b = np.random.random((n,))
Hint
- We can extract the row
ifrom matrixAwithA[i,:]or justA[i].
%run Exercise2.py
Inverse matrix with Gauss-Jordan method¶
To compute $A_{n\times n}$ by Gauss-Jordan method, we will follow the steps:
Write a $n\times 2n$ matrix $M$, built concatenating matrices $A$ and $I$ (identity matrix) $n\times n$. That is, $M = [A\,|\,I]$
Operating by rows, we transform matrix $A$ into matrix $I$. Then, the right hand matrix will be the inverse of $A$. That is, from $M = [A\,|\,I]$ we arrive to $M = [I\,|\,A^{-1}]$ with the following operations:
For $k = 1,\ldots,n$
- If $m_{kk}=0$, exit.
- Else,
- $m\,\mathrm{row}_k \leftarrow \dfrac{m\,\mathrm{row}_k}{m_{kk}}$
- For $i = 1,\ldots,n \quad i \ne k$
$m\,\mathrm{row}_i \leftarrow m\,\mathrm{row}_i - m_{ik}\times m\,\mathrm{row}_k$
Write a function gaussjordan(A) with input parameter the square matrix A and output the inverse of A and that computes this inverse using Gauss-Jordan method without pivoting.
Compute the inverse of the matrices from the previous exercise.
Hints
- The identity matrix $n \times n$ can be obtained with
I = np.eye(n). - We can concatenate two matrices
AandBhorizontally (if the number of rows in both matrices is the same) withM = np.concatenate((A,B),axis=1). - We can extract row
iof matrixMwithM[i,:]or justM[i]. - We can extract the last columns from column
nwithM[:,n:].
%run Exercise3.py
More exercises¶
Forward substitution¶
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). $$
Algorithm¶
- If $l_{11}=0$, exit
- Else, set $x_{1}=\dfrac{b_{1}}{l_{11}}$
- For $i=2,\ldots,n$
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$.
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)
L = np.random.random((n,n))
L = np.tril(L) # hace cero los elementos por encima de la diagonal
b = np.random.random(n)
%run Exercise4.py
LU decomposition¶
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
- Solve $L\mathbf{y}=\mathbf{b}$, by forward substitution, to get $\mathbf{y}$.
- Solve $U\mathbf{x}=\mathbf{y}$, by backward substitution, to get $\mathbf{x}$.
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 & 1\\ 1 & 1 & 2 \end{array}\right)\quad\quad b=\left(\begin{array}{c} 2\\ 4\\ 6 \end{array}\right) $$- Create a function
LU(A)that gives the $LU$ decomposition of $A$. - Obtain and print the matrices $L$ and $U$.
- Use the
subs_backwardandsubs_forwardfunctions to solve the resulting systems. - Check the solution computing $Ax$ and comparing it with $b$.
- You can initialize the matrices:
U = np.zeros((n,n))(matrix of zeros with the same shape as $A$, $n\times n$) orU = np.copy(A)(copy of $A$) Beware!U = Aworks differently.L = np.eye(n)(identity matrix $n\times n$)
Check that the function also works well with the $Ax=b$ system with matrices $A$ and $b$ generated as follows
n = 5
np.random.seed(4)
A = np.random.random((n,n))
b = np.random.random(n)
`
%run Exercise5.py