### Linear Systems of Equations

```Systems of Linear Equations
Ax  b
Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften
ETH Hönggerberg / HCI F128 – Zürich
E-Mail: [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
1
Problem Definition
 We want to solve the problem
A ( n , n )  x ( n ,1)  b ( n ,1)
where x is the vector of unknowns, while A and b are given
 Assumptions
• A-1 exists
1. The number of equations is equal
to not
the singular
number of unknowns,
• A is
therefore A is a square matrix
• A’s columns are linearly independent
2. All components of A, b and x are real
• A’s rows are linearly independent
3. A solution exists and it is unique
• det(A) is non-zero
• rank(A) is equal to n
• Ax = 0 only if x is a null vector
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
2
Analytical Approach
 Cramer’s rule (1750):
The solution is given by
xi 
det( A i )
det( A )
b replaces the ith column
where Ai is defined as follows
 a1,1

a 2 ,1

Ai 


 a n ,1
a1, 2
a1, i 1
b1
a1, i  1
a 2 ,2
a 2 , i 1
b2
a 2 , i 1
a n ,2
a n , i 1
bn
a n ,i 1
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
a1, n 

a 2,n



a n ,n 
3
Calculation of the Determinant
 The Laplace formula (1772) allows the computation of the
determinant of a square matrix:
N
det( A ) 
a
1, i
C 1, i
i 1
 where Ci,j is the determinant of the sub-matrix obtained by removing
the ith row and the jth column of the matrix, multiplied by (-1)i+j:
C i , j  (  1)
 a
det( M i , j )  1,1

 a i 1,1
M i, j  
 a i  1,1
no ith row 

 a
 n ,1
i j
a1, j 1
a1, j  1
a i 1, j 1
a i 1, j  1
a i  1, j 1
a i  1, j  1
a n , j 1
a n , j 1
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
a1, n 

no jth column

a i 1, n 

a i  1, n 


a n , n 
4
A First Numerical Approach: Gauss Elimination
Method
 The following operations do not change the result:
1. Multiply a line by a constant
2. Substitute a line with a linear combination of multiple lines
3. Permute the order of lines
 This can be used to produce a triangular matrix, which
allows the solution to be found easily by substitution
 Example system
3

1

 1
2
2
1
1   x1   1 
    
2  x2  5
    
1   x 3   2 
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
5
Gauss Elimination Example
Triangular System
• Multiply by -3
Multiply
by-4-3
• Sum it •to• Multiply
1st
line by
• • Sum
Sumitittoto2nd
1st line
line
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
6
Generalization of the Gauss Elimination Method
 We want to find a general procedure to replace one entry in
the matrix with a zero; for this we define a multiplier l21
 a11

A  a 21

 a 31
a12
a 22
a 32
a13 

a 23

a 33 
l 21 
a 21
a11
( new )
 a 21
 l 21 a11  0
( new )
 a 22
 l 21 a12
( new )
 a 23
 l 21 a13
a 21

a 22
a 23
( old )
( old )
( old )
 Note that
 1

l
 21
 0
0
1
0
0   a11
 
0  a 21
 
1   a 31
a12
a 22
a 32
a13   a11
 
a 23  0
 
a 33   a 31
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
a12
( new )
a 22
a 32
a13 
( new ) 
a 23

a 33 
7
The Gauss Elimination in Matrix Notation
 If we put all multipliers for one column in one matrix, we get
M
(0)
A
(0)
 1
 (0)
  l 21

(0)
  l 31
0
1
0
(0)
0   a11
  (0)
0   a 21

(0)
1   a 31
(0)
a12
(0)
a 22
(0)
a 32
(0)
(0)
a13   a11

(0) 
a 23    0
(0)
a 33   0
(0)
a12
(1)
a 22
(1)
a 32
(0)
a13 
(1) 
(1)
a 23   A
(1)
a 33 
(k )
(k )
where lij 
a ij
(k )
a jj
 This way, a triangular matrix is easily obtained
A
b
( n  1)
( n  1)
M
(n2)
M
(n2)
M
(1)
M
(1)
M
(0)
M
(0)
A
b
(0)
(0)




A
( n  1)
xb
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
( n  1)
8
Example: Gauss Elimination in Matrix Notation
Total number of operations required
n 1
n 1
j2
j2
 j ( j  1)   j 
2
n
3
3
(n-1)(n-2)
operations(flops)
(flops)
n(n-1) operations
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
9
Gauss Elimination / Transformation Method
 The Gauss elimination method is relatively easy to
implement (even by hand), but has some distinct
 Changes the matrix A
 Requires and changes the coefficient vector b
 Must be rerun if the vector b changes
 If we consider the Gauss method in matrix form on the
other hand, we can see that we can use
M M
(n2)
M
(1)
M
(0)
to transform A and b; M is therefore called the Gauss
transformation matrix
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
10
The Gauss Transformation Method
 We have so far
 The final matrix A, which is an right triangular matrix
 The matrix M, which is a left triangular matrix
 The inverse of M is also a left triangular matrix
L M
1
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
 1
 (0)
 l 21

(0)
 l 31
0
1
(1)
l 31
0

0

1 
11
LR (LU) Factorization
 We define our (right triangular) solution matrix as follows
A
(n2)
 MA
(0)

R  MA
(0)
 If we multiply with L = M-1 from the left on both sides
1
LR  M M A
(0)
 A
(0)
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
12
LR Factorization (Continued)
 The starting matrix A is transformed (factorized) as
A  LR
 If we apply this to a linear equation system, we get
Ax  LRx  b
Ly  b
 
Rx  y
 This approach rids us of the disadvantages discussed
earlier, because
 For every vector b, two simple triangular systems must be solved
without factorizing again
 The matrices L and R can be stored using the elements of A
 If A is modified, it is often possible to modify L and R accordingly
without re-factorizing
 Note that Gauss elimination is still needed once to compute L and R
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
13
Pivoting
 It is easy to see from the definition of the multiplication
factors, that the diagonal elements at each step (the pivot
values) cannot be equal to zero
 This is circumvented by reordering the rows of the matrix A
by multiplication with a permutation matrix P
1
PA  PM R  LR
 This approach is referred to as LRP-factorization
(or LR-factorization with partial pivoting)
 Example:
 x2  1

 0.4 x1  0.3 x 2  0.1
a11 = 0  switch the lines  x1 = x2 = 1
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
14
Pivoting (Continued)
 Pivoting must also take into account scaling problems;
Let us consider a similar example
 2  10  20 x1  x 2  1

 0.4 x1  0.3 x 2  0.1
 Pivot elements should have large absolute values
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
15
How does Matlab do it?
 When using left divide \, Matlab chooses a procedure
depending on the properties of the problem, i.e. if A is
1. Diagonal, x is computed directly by division
2. Sparse, square and banded, then banded solvers are used;
Either Gauss elimination without pivoting or LR factorization
3. Left or right triangular, then backsubstitution is used
4. A permutation of a triangular matrix, it is permuted and 3. applies
5. Symmetric or Hermitian, then a Cholesky factorization is
attempted (A = RR*, where R* is the conjugate transpose of R);
If it fails, another indefinite symmetric factorization is attempted
6. Square but 1 through 5 do not apply, then LR factorization with
partial pivoting is applied
7. Not square, then Householder reflections are used to compute a
factorization which leads to a least-squares solution, i.e. a vector x
which minimizes the length of Ax – b
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
16
Assignment 1
1. Find online the template for an LR-factorization of a matrix
with partial pivoting. Make the function operational by
adding the lines that calculate the M matrix of the current
step and the new A matrix.
 Why does the transformation matrix T appear in the formula for M?
Explain by comparing to the definition of M in the non-pivoting case.
2. Use the function to factorize the following matrix
1

A  9

 5
2
8
0
3

1

1 
 Test if the factorization worked, i.e. if LR = A.
 Is L in the form you would expect it to be? What implications does
this have for its application in solving a linear system (see slide 13)
and how could you correct for it?
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
17
Hints
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
18
Assignment 2
1. Create a random matrix A with dimensions 4x4 and a
random column vector b with size 4x1. Solve the system
Ax = b.
2. Create a function that computes the determinant of square
matrices using the Laplace formula. Use a recursive
approach (see hints).
3. Use this function to compute the solution of the linear
system above using Cramer’s rule.
4. Do the same for linear equation systems with sizes
ranging from 5 to 9.
5. Read out the CPU time required to solve all these systems
with both methods (Cramer’s method and A\b) and
compare them.
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
19
Hints
 If you remember the Laplace formula as a sum
N
det( A ) 

i 1
N
a1, i C 1, i 

a1, i (  1)
1 i
det( M 1, i )
i 1
you can see that the calculation of the determinant requires
the calculation of a determinant. You could let the function
call itself to do that (recursive function). Remember that the
determinant of a 1x1 matrix is equal to its only element.
 There are several Matlab commands to read out timings,
however the most reliable one is
 t_start = tic;
statements;
t_elapsed = toc(t_start);
 Where the time elapsed (in second) is stored in t_elapsed
Daniel Baur / Numerical Methods for Chemical Engineers / Linear Equation Systems
20
```