Report

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] 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 (n2) M (n2) M (1) M (1) M (0) M (0) A b (0) (0) A ( n 1) xb 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 j2 j2 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 disadvantages; namely it 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 (n2) 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 (n2) 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