import numpy as np
Given the linear system $Ax=b$ where the coefficient matrix $A$ and the constant matrix $b$ are
n = 7
A1 = np.diag(np.ones(n))*3
A2 = np.diag(np.ones(n-1),1)
A = A1 + A2 + A2.T
b = np.arange(n,2*n)*1.
print('A')
print(A)
print('b')
print(b)
We can solve this system with the generic function
x = np.linalg.solve(A,b)
print('x')
print(x)
Let us see an specific method to solve these type of linear systems, the tridiagonal systems.
Gauss method for tridiagonal systems¶
The coefficient matrix of a tridiagonal systems is tridiagonal. This is a matrix where all the elements are zero, except the elements in the main diagonal and the adyacent ones that are over and under. For example, the following linear system
has a $7\times7$ tridiagonal coefficient matrix
$$ A = \left( \begin{matrix} a_{11} & a_{12} & 0 & 0 & 0 & 0& 0\\ a_{21} & a_{22} & a_{23} & 0 & 0 & 0& 0\\ 0 & a_{32} & a_{33} & a_{34} & 0 & 0& 0\\ 0 & 0 & a_{43} & a_{44} & a_{45} & 0& 0\\ 0 & 0 & 0 & a_{54} & a_{55} & a_{56}& 0\\ 0 & 0 & 0 & 0 & a_{65} & a_{66}& a_{67}\\ 0 & 0 & 0 & 0 & 0 & a_{76}& a_{77} \end{matrix} \right) \quad \quad b = \left( \begin{matrix} b_1\\ b_2\\ b_3\\ b_4\\ b_5\\ b_6\\ b_7 \end{matrix} \right) $$Gauss method has two steps:
- Triangularization of the augmented matrix.
- Back-substitution.
Triangularization¶
We must obtain zeros under the main diagonal using row operations. In the actual example we do
Given a tridiagonal linear system, write a function At, bt = triangular(A,b) that returns matrices At and bt that are the transformed matrices after the triangularization process. Use one and only one for loop.
- To read better the matrices, use
np.set_printoptions(precision = 2) # only two decimals
np.set_printoptions(suppress = True) # do not use exponential format
Obtain the number of rows and columns of
Awithshape. For instance,m, n = A.shape. You can also obtain the number of elements of a vector withlen. For instance,n = len(b).Initialize the triangular matrix with
At = np.copy(A)if you do not want to modifyA. Similarly withb.Triangularize $Ax=b$, where matrices
Aandbcan be generated with the code
n = 7
A1 = np.diag(np.ones(n))*3
A2 = np.diag(np.ones(n-1),1)
A = A1 + A2 + A2.T
b = np.arange(n,2*n)*1.
- Leave the function as it is and triangularize also
Aandbgenerated with the code
n = 8
np.random.seed(3)
A1 = np.diag(np.random.rand(n))
A2 = np.diag(np.random.rand(n-1),1)
A = A1 + A2 + A2.T
b = np.random.rand(n)
%run Exercise1.py
We solve the equations from the bottom to the top
Write a function x = back_subs(At,bt) that solves a triangular system with back-substitution. Use one and only one for loop. Use it with the function triangular to solve the linear systems of the previous exercise.
Hint:
- You can use a loop with a negative step. For instance where
for k in range(n1,n2,-1):
n1is the initial value,n2the next value to the final one and-1the step.
%run Exercise2.py
Storing the coefficients in a rectangular matrix of three columns¶
Let us store now the coefficients of the coefficient matrix in a rectangular matrix. The linear system was
If we store the coefficients in a rectangular matrix
$$ A_r = \left( \begin{matrix} 0 & a_{11} & a_{12}\\ a_{21} & a_{22} & a_{23}\\ a_{32} & a_{33} & a_{34}\\ a_{43} & a_{44} & a_{45}\\ a_{54} & a_{55} & a_{56}\\ a_{65} & a_{66} & a_{67}\\ a_{76} & a_{77} & 0 \end{matrix} \right)\quad = \quad \left( \begin{matrix} 0 & r_{12} & r_{13}\\ r_{21} & r_{22} & r_{23}\\ r_{31} & r_{32} & r_{33}\\ r_{41} & r_{42} & r_{43}\\ r_{51} & r_{52} & r_{53}\\ r_{61} & r_{62} & r_{63}\\ r_{71} & r_{72} & 0 \end{matrix} \right) $$Now, taking into account the new storing
And now, the triangular matrix is stored
And the back-substitution will be
Linear systems with tridiagonal matrices are quite common in physics problems (heat equation, wave equation). Thus, common algorithms to solve linear systems are adapted so these kind of probles can be solved more efficiently. For example, we can store the diagonals of the $A$ in a rectangular matrix $A_r$ with only three columns. We can store the square matrix of Exercise 1 in a rectangular matrix
n = 7
Ar = np.zeros((n,3))
Ar[:,0] = np.concatenate((np.array([0]),np.ones((n-1),)))
Ar[:,1] = np.ones((n),)*3
Ar[:,2] = np.concatenate((np.ones((n-1),),np.array([0])))
b = np.arange(n,2*n)*1.
Modify the code from the two previous exercises son we can use $A_r$ instead of $A.$
- Without modifying the functions, apply them also to the data
n = 8
np.random.seed(3)
Ar = np.zeros((n,3))
Ar[:,1] = np.random.rand(n)
Ar[:,0] = np.concatenate((np.array([0]),np.random.rand(n-1)))
Ar[0:n-1,2] = Ar[1:n,0]
b = np.random.rand(n)
%run Exercise3.py
Proposed exercises¶
Matrix product¶
We use three numpy arrays $A_{m\times n}$, $B_{n\times p}$ and $C_{m\times p}$.
As we are multipliying $A$ by $B$, the number of columns in $A$, $n$, must be the same as the number of rows in $B$. And the result matrix, $C$, must have the same number of rows as $A$ and the same number of columns as $B$.
An element of $C$ is given by
$$ c_{ij} = a_{i1}\,b_{1j}+a_{i2}\,b_{2j}+\dots+a_{in}\,b_{nj} = \sum_{k=1}^{n}a_{ik}\,b_{kj}\quad \mathrm{for} \; i = 1,\ldots,m \quad j = 1,\ldots,p \quad (1) $$
Taking into account the equation (1), write a function C = multiply_3loops(A,B) that multiplies matrices A and B and returns the result C. Use three for loops. One with index i for the row, another with index j for the column and a last index k for the cumulative addition.
- Use
shapeto get the dimensions ofAandB. For instance,m, n = A.shape - Initialize
C = np.zeros((m,p))beingmthe number ofArows andpthe number ofBcolumns. - Then, every
Celement has been initialized to zero and we get its value as a cumulative addition. That is, we start with $c_{ij}=0$, and then $$k= 1\qquad c_{ij}\leftarrow c_{ij}+a_{i1}b_{1j}$$ $$k= 2\qquad c_{ij}\leftarrow c_{ij}+a_{i2}b_{2j}$$ $$\ldots$$ $$k= n\qquad c_{ij}\leftarrow c_{ij}+a_{in}b_{nj}$$ - Use the function to multiply the matrices
A = np.array([[-3.,2],[-2,0],[-4,4],[4,-4]]) B = np.array([[4.,-3,1],[-2,1,1]])
- Without modifying the function, multiply the matrices
A1 = np.array([[-3.,2],[-2,0],[-4,4],[4,-4],[1,1]]) B1 = np.array([[4.,-3,1,1],[-2,1,1,1]])
%run Exercise4.py
Write a function C = multiply_2loops(A,B) that multiplies matrices A and B using two loops. Multiply the matrices of Exercise 4. We can remove a look if we consider that
A[i,:]is a vector of n components that contains the rowifromA.B[:,j]is a vector of n components that contains the columnjfromB.A[i,:]*B[:,j]is the element-wise product of these two vectors and it is also a vector with n components.np.sum(A[i,:]*B[:,j])is the addition of the elements of this vector.
%run Exercise5.py
We can improve the vectorization getting a column (or a row) of the matrix C with only a code line. That way, we can remove one more loop.
B[:,j]is a vector of n components that contains the columnjfromB.- In
D = A*B[:,j]every row ofDis the element-wise product of each row ofAwith the vectorB[:,j]. Thus,Dis the same size asA. C[:,j] = np.sum(D,axis=1)is a column matrix where each element is the addition of the elements of a row ofDand it is the matrix product of matrixAby column matrixB[:,j].
A = np.array([[-3.,2],[-2,0],[-4,4],[4,-4]])
B = np.array([[4.,-3,1],[-2,1,1]])
m = A.shape[0]
p = B.shape[1]
C = np.zeros((m,p))
for j in range(p):
C[:,j] = np.sum(A*B[:,j],axis=1)
print(C)
%run Exercise6.py
Compare the execution time for the three functions that multiply matrices and the numpy command np.dot(A,B) and A@B. Use matrices $A_{m\times n}$ and $B_{n\times p}$ with random elements. Use the command time.time() from the previous lab. np.random.rand(m,n) is a $m\times n$ matrix with random elements from $[0,1]$.
m = 300; n = 200; p = 400
A = np.random.rand(m,n)
B = np.random.rand(n,p)
%run Exercise7.py
If we multiply a column matrix $A_{m\times 1}$ by a row matrix $B_{1 \times n}$ the result is a matrix $C_{m\times n}$. For instance
$$ AB = \left( \begin{array}{c} a_1\\ a_2\\ a_3\\ a_4 \end{array} \right) \left( \begin{array}{ccc} b_1 & b_2 & b_3 \end{array} \right) = \left( \begin{array}{ccc} a_1b_1 & a_1b_2 & a_1b_3 \\ a_2b_1 & a_2b_2 & a_2b_3 \\ a_3b_1 & a_3b_2 & a_3b_3 \\ a_4b_1 & a_4b_2 & a_4b_3 \end{array} \right) = \left( \begin{array}{ccc} c_{11} & c_{12} & c_{13} \\ c_{21} & c_{22} & c_{23} \\ c_{31} & c_{32} & c_{33} \\ c_{41} & c_{42} & c_{43} \end{array} \right) = C $$Using two for loops, write a function C = multiply_col_row(a,b) that multiplies a column matrix by a row matrix a return the product matrix. The parameters a and b are given as 1D numpy arrays, while C will be a 2D numpy array.
Multipliy
a1 = np.array([1,3,4])
b1 = np.array([2,1])
And also
a2 = np.array([1,3])
b2 = np.array([2,1,0,1])
%run Exercise8
Given the matrices $A_{m\times n}$ and $B_{n\times p}$ where $A_1$, $A_2$, $\ldots$, $A_n$ represent the columns of the matrix $A$ and $B_1$, $B_2$, $\ldots$, $B_n$ the rows of matrix $B$, we can obtain the product of these two matrices
$$ AB = \left( \begin{array}{cccc} A_1 & A_2 & \ldots & A_n \end{array} \right) \left( \begin{array}{c} B_1\\ B_2\\ \dots\\ B_n \end{array} \right) = A_1B_1+A_2B_2+\dots+A_nB_n=C_1+C_2+\dots+C_n=C $$Where $C_1$, $C_2$, $\ldots$, $C_n$ are $m\times p$ matrices
Using a for loop, write a function C = multiply(A,B) that multiplies matrices A and B and returns matrix C. using this algorithm. To multiply column matrix by row matrix use the function Ci = multiply_col_row(A,B)
Multiply
A1 = np.array([[-3,2],[-2,0],[-4,4],[4,-4]])
B1 = np.array([[4,-3,1],[-2,1,1]])
And also
A2 = np.array([[-3,2],[-2,0],[-4,4],[4,-4],[1,1]])
B2 = np.array([[4,-3,1,1],[-2,1,1,1]])
%run Exercise9