ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida 24Sep2018
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Cmtrx}{\boldsymbol{\mathsf{C}}} \newcommand{\Dmtrx}{\boldsymbol{\mathsf{D}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Qmtrx}{\boldsymbol{\mathsf{Q}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \DeclareMathOperator{\rank}{rank} $
The course notes (OneNote ChEn-3170-linalg) cover basic elements of lineary system of algebraic equations. Their representation in row and column formats and the insight gained from them when searching for a solution of the system.
Some basic theoretical aspects of solving for $\xvec$ in the matrix equation $\Amtrx\,\xvec = \bvec$ are covered. $\Amtrx$ is a matrix, $\bvec$ and $\xvec$ are vectors.
The following operations between vectors and matrices are obtained directly from the buil-in functions in the numpy
package.
'''Import the NumPy package as usual'''
import numpy as np
Scalar product of two vectors: $\avec \cdot \bvec$.
'''Vector scalar product or dot product of vectors'''
a_vec = np.array( np.random.random(3) )
b_vec = np.array( np.random.random(3) )
a_vec_dot_b_vec = np.dot( a_vec, b_vec ) # clear linear algebra operation
print('a.b =', a_vec_dot_b_vec)
a_vec_x_b_vec = a_vec @ b_vec # consistent linear algebra multiplication
print('a@b =', a_vec_x_b_vec )
Matrix vector product: $\Amtrx\,\bvec$.
'''Matrix-vector product'''
a_mtrx = np.array( [ [ 2., 1., 1.], # per course notes (p 04)
[ 4., -6., 0.],
[-2., 7., 2.]])
b_vec = np.array( [5., -2., 9.]) # per course notes
a_mtrx_x_b_vec = a_mtrx @ b_vec # linear algebra matrix-vector product
print('A x b =', a_mtrx_x_b_vec)
Matrix-vector product: $\Imtrx\,\bvec = \bvec$. Note: $\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \, \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} = \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} $.
'''Identity-matrix vector product'''
i_mtrx = np.eye(3)
i_mtrx_x_b_vec = i_mtrx @ b_vec # linear algebra matrix-vector product
print('I x b =', i_mtrx_x_b_vec)
print('b =', b_vec)
Matrix-matrix product: $\Imtrx\,\Amtrx = \Amtrx$. Note: $\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \, \begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} = \begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} $.
'''Matrix-matrix product IA = A'''
i_mtrx_x_a_mtrx = i_mtrx @ a_mtrx # linear algebra matrix-matrix product
print('I x A =\n', i_mtrx_x_a_mtrx)
print('A =\n', a_mtrx)
Matrix-matrix product: $\Amtrx\,\Bmtrx = \Cmtrx$. Note: $\begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} \, \begin{pmatrix} B_{1,1} & B_{1,2} & B_{1,3} \\ B_{2,1} & B_{2,2} & B_{2,3} \\ B_{3,1} & B_{3,2} & B_{3,3} \end{pmatrix} = \begin{pmatrix} C_{1,1} & C_{1,2} & C_{1,3} \\ C_{2,1} & C_{2,2} & C_{2,3} \\ C_{3,1} & C_{3,2} & C_{3,3} \end{pmatrix} $ where each $C_{i,j}$ is a vector product of the $i$th row of $\Amtrx$ and the $j$th column of $\Bmtrx$, i.e. $C_{i,j} = \sum\limits_{k=1}^3 A_{i,k}\, B_{k,j}$.
'''Matrix-matrix product AB = C'''
b_mtrx = np.array( [[5. , 5. , 5.],
[-2., -2., -2.],
[9. , 9. , 9.]]
)
c_mtrx = a_mtrx @ b_mtrx # linear algebra matrix-matrix product
print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('C =\n', c_mtrx)
'''Matrix-matrix product BA = D'''
b_mtrx = np.array( [[5. , 5. , 5.],
[-2., -2., -2.],
[9. , 9. , 9.]]
)
d_mtrx = b_mtrx @ a_mtrx # linear algebra matrix-matrix product
print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('D =\n', d_mtrx)
NumPy has extensive support for linear algebra arrays. We collect here the relevant operations for this course. However additional resources are instead added to SciPy for general scientific computing including linear algebra.
Linear algebra operations are obtained from the linalg
sub-package of the numpy
package, and the linalg
sub-package of scipy
.
'''Import the NumPy linear algebra sub-package as usual'''
#import numpy.linalg as linalg
#from numpy import linalg # often used alternative
'''or leave it commented since the usage of np.linalg is self-documenting'''
The 2-norm or norm (magnitude) of a vector $\bvec$ is indicated as $\norm{\bvec}$ and computed as follows:
'''Vector norm (or magnitude)'''
norm_b_vec = np.linalg.norm( b_vec ) # default norm is the 2-norm
print('||b|| =', norm_b_vec) # same as magnitude
Solve for $\xvec$ in the matrix equation $\Amtrx\,\xvec = \bvec$, where $\Amtrx = \begin{pmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{pmatrix} $ and $\bvec = \begin{pmatrix} 5\\ -2\\ 9 \end{pmatrix}$.
'''Matrix solver (this is short for solution of a linear algebraic system of equations)'''
x_vec = np.linalg.solve( a_mtrx, b_vec ) # solve linear system for A, b; per course notes 04
print('solution x =', x_vec)
The residual vector defined as $\rvec = \bvec - \Amtrx\,\xvec$ is of importance. So is its norm $\norm{\rvec}$.
'''Verify the accuracy of the solution'''
res_vec = b_vec - a_mtrx @ x_vec
print('b - A x =',res_vec)
print('||b - A x|| =',np.linalg.norm( res_vec ))
The rank of a matrix of coefficients, $\rank(\Amtrx)$, of a linear algebraic system of equations determines weather the solution is unique or singular.
'''Matrix rank'''
k = np.linalg.matrix_rank( a_mtrx ) # rank; per course notes 14
print('rank(A) =',k)
print('shape(A) =',a_mtrx.shape)
if k == a_mtrx.shape[0] and k == a_mtrx.shape[1]: # flow control
print('A is non-singular; solution is unique ')
else:
print('A is singular')
Why is this matrix $\Bmtrx$ singular?
b_mtrx = np.array( [ [ 2., 1., 3.], # singular
[ 4., -6., -2.],
[-2., 7., 5.]])
k = np.linalg.matrix_rank( b_mtrx ) # rank
print('rank(B) =',k)
print('shape(B) =',b_mtrx.shape)
if k == b_mtrx.shape[0] and k == b_mtrx.shape[1]: # flow control
print('B is non-singular; solution is unique ')
else:
print('B is singular')
'''Matrix determinant'''
det_a_mtrx = np.linalg.det( a_mtrx ) # determinant; course notes 16
print('det(A) =', det_a_mtrx)
The inverse matrix is denoted as $\Amtrx^{-1}$ and is computed as the matrix that multiplies $\bvec$ and produces the solution $\xvec$, that is, $\xvec = \Amtrx^{-1}\,\bvec$.
'''Matrix inverse'''
a_mtrx_inv = np.linalg.inv( a_mtrx ) # matrix inverse; per course notes 17
print('A^-1 =\n', a_mtrx_inv)
Recall $\Amtrx^{-1}\,\Amtrx = \Imtrx$ where $\Imtrx$ is the identity matrix.
'''Identity matrix'''
i_mtrx = a_mtrx_inv @ a_mtrx # identity matrix; per course notes 17
print('A^-1 A =\n',i_mtrx)
Using the inverse, the same solution will be found: $\xvec = \Amtrx^{-1}\,\bvec$.
'''Solution using the inverse'''
x_vec_again = a_mtrx_inv @ b_vec # matrix-vector multiply; per course notes 17
print('solution x =', x_vec_again)
This is the element-by-element reciprocal of the matrix $(\Amtrx)^{-1}$, which is very different than the inverse.
'''Inverse power of a matrix'''
#a_mtrx_to_negative_1 = a_mtrx**(-1) # this will cause an error (division by zero)
Let's look at the determinant of a larger matrix, say $\Mmtrx$.
'''Generate a larger matrix from an image'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output
mtrx = plt.imread('https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cermet.png',format='png')
m_mtrx = mtrx[:,:]
plt.figure(1) # create a figure placeholder
plt.imshow( m_mtrx,cmap='gray')
plt.title('Cermet',fontsize=14)
plt.xlabel('x pixels',fontsize=12)
plt.ylabel('y pixels',fontsize=12)
print('M shape =', m_mtrx.shape)
'''Larger matrix determinant'''
det_m_mtrx = np.linalg.det( m_mtrx ) # determinant
print('max(M) =',m_mtrx.max())
print('min(M) =',m_mtrx.min())
print('det(M) = %10.3e (not an insightful number)'%det_m_mtrx) # formatting numeric output
print('rank(M) = ',np.linalg.matrix_rank( m_mtrx, tol=1e-5 ) )
Let's solve for this matrix with $\cvec$ as the right side vector, that is, $\Mmtrx\,\xvec = \cvec$.
'''Solve M x = c and plot x'''
c_vec = np.random.random(mtrx.shape[0]) # any c will do it
sol = np.linalg.solve( m_mtrx, c_vec ) # solve linear system for A, b
plt.figure(2)
plt.plot(range(c_vec.size),sol,'k')
plt.title('M x = c',fontsize=20)
plt.xlabel('n',fontsize=18)
plt.ylabel('$c_j$',fontsize=18)
print('')
res_vec = c_vec - m_mtrx @ sol
#print('c - M x =',res_vec)
print('||c - M x|| =%12.4e'%np.linalg.norm( res_vec ))
A lower triangular matrix allows for a forward solve.
'''L forward solve'''
l_mtrx = np.array( [[1., 0., 0.], # per course notes
[2., 3., 0.],
[4., 5., 6.]] )
b_vec = np.array( [1.,2.,3.] )
x_vec = np.linalg.solve( l_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
An upper triangular matrix allows for a backward solve.
'''U backward solve'''
u_mtrx = np.array( [[1., 2., 3.], # per course notes
[0, 4., 5.],
[0., 0., 6.]] )
b_vec = np.array( [1.,2.,3.] )
x_vec = np.linalg.solve( u_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Pmtrx\,\Lmtrx\,\Umtrx = \Amtrx$ can be obtained from the SciPy linear algebra package. $\Pmtrx$ is a permutation matrix if the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
'''Import only the linear algebra package'''
import scipy.linalg
import numpy as np
'''P L U factors of A'''
a_mtrx = np.array( [[1, 2, 3],
[4, 5, 6],
[7, 8, 10]] )
(p_mtrx, l_mtrx, u_mtrx) = scipy.linalg.lu( a_mtrx )
print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('Checking...')
print('PLU - A =\n', p_mtrx @ l_mtrx @ u_mtrx - a_mtrx)
'''P^-1 = P^T (i.e. the transpose of a permutation matrix is its inverse)'''
pinv_mtrx = np.linalg.inv(p_mtrx)
print('P^-1 =\n', pinv_mtrx)
print('Checking...')
print('P^-1 - P^T =\n', pinv_mtrx - p_mtrx.transpose())
'''PLU x = b; that is: Forward: L y = P^-1 b, Backward: U x = y '''
b_vec = np.array([1.,2.,3.])
y_vec = scipy.linalg.solve(l_mtrx, p_mtrx.transpose() @ b_vec) # L y = P^T b
x_vec = scipy.linalg.solve(u_mtrx, y_vec) # U x = y
print('x =', x_vec)
x_vec_gold = scipy.linalg.solve( a_mtrx, b_vec ) # solution using A x = b
print('||x - x_gold|| =',scipy.linalg.norm(x_vec-x_vec_gold))
'''Deterninant of U or L: product of the diagonal'''
det_u = np.linalg.det(u_mtrx)
print('det(U) = %8.3e'%det_u)
diag_vec = np.diagonal(u_mtrx)
prod = np.prod(diag_vec)
print('diag(U) product = %8.3e'%prod )
'''Determinant of P (always +1 or -1)'''
det_p = np.linalg.det(p_mtrx)
print('det(P) = %8.3e'%det_p)
'''Determinant of A = det(PLU)'''
det_l = np.prod( np.diagonal(l_mtrx) )
det_plu = det_p * det_l * det_u # last term is det of L
print('det(PLU) = %8.3e'%det_plu)
print('det(A) = %8.3e'%np.linalg.det(a_mtrx))
In this course various algorithms need to be programmed. These should be compared to SciPy
and/or NumPy
.
A lower triangular matrix allows for a forward solve. The algorithm for $\Lmtrx\,\xvec=\bvec$ is as follows:
\begin{equation*} x_i = \Bigl(b_i - \sum\limits_{j=1}^{i-1} L_{i,j}\,x_j \Bigr)\,L^{-1}_{i,i} \quad\ \forall \quad\ i=1,\ldots,m \end{equation*}for $i$ and $j$ with offset 1. Recall that NumPy
and Python
have offset 0 for their sequence data types.
'''L forward solve'''
l_mtrx = np.array( [[1., 0., 0.], # per course notes
[2., 3., 0.],
[4., 5., 6.]] )
b_vec = np.array( [1.,2.,3.] )
from chen_3170.help import forward_solve
x_vec = forward_solve( l_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
A upper triangular matrix allows for a backward solve. The algorithm for $\Umtrx\,\xvec=\bvec$ is as follows:
\begin{equation*} x_i = \Bigl(b_i - \sum\limits_{j=i+1}^{m} U_{i,j}\,x_j \Bigr)\,U^{-1}_{i,i} \quad\ \forall \quad\ i=m,\ldots,1 \end{equation*}for $i$ and $j$ with offset 1. Recall that NumPy
and Python
have offset 0 for their sequence data types.
'''U backward solve'''
u_mtrx = np.array( [[1., 2., 3.], # per course notes
[0, 4., 5.],
[0., 0., 6.]] )
b_vec = np.array( [1.,2.,3.] )
try:
from chen_3170.toolkit import backward_solve
except ModuleNotFoundError:
assert False, 'You need to provide your own backward_solve function here. Bailing out.'
x_vec = backward_solve( u_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
$\Lmtrx\,\Umtrx$ factorization algorithm (without using pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ computes the $\Lmtrx\,\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that
\begin{equation*} A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \quad\ \forall\ i=k+1,\ldots,m \quad\ \text{and} \quad\ j=k,\ldots,m \end{equation*}where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$.
'''L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx )
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - A =\n', l_mtrx @ u_mtrx - a_mtrx)
The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx$ can be obtained from partial pivoting strategy to the $\Lmtrx\,\Umtrx$ factorization algorithm shown above. $\Pmtrx$ is a row permutation matrix obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using partial pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,\Lmtrx\,\Umtrx$ factors.
The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row exchanges recorded in the permutation matrix $\Pmtrx$.
'''P L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(p_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='partial' )
print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PA =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx)
The factors: $\Pmtrx$, $\Qmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx\,\Qmtrx$ can be obtained from a user-developed algorithm with complete pivoting. $\Pmtrx$ is a row permutation matrix and $\Qmtrx$ is a column permutation matrix, obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using complete pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,\Qmtrx\,\Lmtrx\,\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row exchanges recorded in the permutation matrix $\Pmtrx$ and column exchanges recorded in the permutation matrix $\Qmtrx$.
'''P Q L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(p_mtrx, q_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='complete' )
print('P =\n',p_mtrx)
print('Q =\n',q_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PAQ =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx @ q_mtrx)