lecture8

# lecture8 - Numerical Methods in Engineering ENGR 391 System...

This preview shows page 1. Sign up to view the full content.

This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Numerical Methods in Engineering ENGR 391 System of linear equations Faculty of Engineering and Computer Sciences Concordia University LU decomposition Recall from last lecture • We have seen how we can decompose [A]=[L][U] • Useful to solve systems of equations Ax=b: Ax=LUx=b Write z=Ux => Lz=b Step 1: Decompose the matrix A by LU decomposition: A=LU Step 2: Compute z by forward substitution using Lz=b Step 3: Compute x by backward substitution using Ux=z. Problem • Since LU decomposition is related to Gauss elimination same pitfalls Division by zero 0 1 1 1 1 1 1 1 1 How can we obtain the LU decomposition of the above matrix? Partial pivoting • Remember that we introduced partial pivoting strategy in the Gauss algorithm to avoid numerical problems • How can we implement partial pivoting in the A=LU decomposition ? • This can be done with the help of permutation matrices. Permutation matrix Definition: A permutation matrix P is a matrix consisting of all zeros, except for a single 1 in every row and column. Equivalently, a permutation matrix is created by applying arbitrary row exchanges for the identity matrix. 10 01 01 10 001 001 010 100 010 100 010 100 001 001 100 010 100 010 100 010 001 001 Theorem: Let P be a permutation matrix formed by a particular set of row exchanges applied to the identity matrix. Then, for any matrix A, PA is the matrix obtained applying exactly the same set of row exchange to A. Example: 100ab c ab c 001d e f g h i 010g h i d e f PA=LU decomposition Example: find the PA=LU decomposition of following matrix: 21 A 44 5 4 13 1 5 1 1 4 1 2 Solution: 01021 00144 10013 4 1 0 044 1 002 2 100 8 1 2 4 Solving equations with the PA=LU decomposition • Step 1: Decompose the matrix A by PA=LU decomposition: PA=LU • Step 2: Compute z by forward substitution using Lz=Pb • Step 3: Compute x by backward substitution using Ux=z. An example 2 x 1y 5 z 5 4x 4 y 4z 0 x 3y z 21 44 13 [A] 5 6 x 5 4y 0 1 z {x} 6 = {b} Step 1: PA=LU 01021 00144 10013 5 4 1 1 1 4 1 2 0 044 1 002 2 100 8 1 2 4 Step 2: Compute z with Lz=Pb 1 1 4 1 2 0 1 1 2 0 z1 0 z2 z3 1 0105 0010 1006 z1 z2 0 6 z3 8 0 6 5 Step 3: Compute x with Ux=z 44 4x 0 02 2 y 6 00 8 z 8 x 1 y 2 z 1 Matrix inverse If a matrix [A] is square, there is another matrix [A]-1, called the inverse of [A]; such as: [A][A]-1=[A]-1[A]=[I] ; Identity matrix 10 AA 1 A 1 A 01 0 0 0 0 1 0 00 Matrix inverse To compute the inverse matrix, the first column of [A]-1 is obtained by solving the problem (for 3 3 matrix): Ax b 1 0 0 ; second column: Ax ; and third column: A x b b 0 1 0 0 0 1 The best way to implement such a calculation is to use LU decomposition. An other way to read the definition • Each column j of A-1 is the solution of : 0 0 Ax bj 1 0 0 j Computing A-1 Without Partial Pivoting • Step 1: Decompose A=LU • Step 2: Solve all equations Ax=bj - Compute z by forward substitution using Lz=bj - Compute x by backward substitution using Ux=z • Step 3: Write A-1 With Partial Pivoting • Step 1: Decompose PA=LU • Step 2: Solve all equations PAx=Pbj - Compute z by forward substitution using Lz=Pbj - Compute x by backward substitution using Ux=z • Step 3: Write A-1 An example 1 A 2 1 2 1 2 31 1 Step 1: A=LU Without Partial Pivoting 1 2 1 1 0 2 1 2 2 1 00 7 10 3 31 1 3 01 2 3 0 1 0 2 Step 2: Solve all Ax=bj 1 2 1x 1 1 2 1x 0 2 1 2y 0 2 1 2y 1 z 0 z 0 31 1 31 1 2 1x 0 2 1 2y 0 z 1 31 1 1 Step 2: Solve all Ax=bj 1 2 2 1 31 1x 2y 1 1 0 z 1 0 z1 1 z1 2 0 0 1 0 z2 7 1 z3 3 0 0 z2 z3 3 1 2 1x 0 0 3 0 0y 2z 1 2 5 3 x y z 1 2 2 3 5 6 1 2 5 3 Step 2: Solve all Ax=bj 1 2 2 1 31 1x 2y 1 0 1 z 1 0 0 1 2 3 0 0 z1 0 z1 0 2 0 0 1 0 z2 7 1 z3 3 1 0 z2 z3 1 7 3 3 1x 0 y 2z 0 x 1 7 3 y z 1 2 1 3 7 6 Step 2: Solve all Ax=bj 1 2 2 1 31 1x 2y 1 0 0 z 1 0 0 1 2 3 0 0 z1 0 2 1 0 1 0 z2 7 1 z3 3 0 1 3 1x 0 x y 0 y 2z 1 z 0 z1 z2 z3 1 2 0 1 2 0 0 1 Step 3: Write A-1 x y z 1 2 2 3 5 6 1 2 1 3 7 6 x y z 1 2 1 2 1 2 31 1 1 1 2 2 3 5 6 x y z 1 2 1 3 7 6 1 2 0 1 2 1 2 0 1 2 Gauss-Seidel Method • An iterative method Approach • The system Ax=b is reshaped by solving the first equation for x1, the second equation for x2, and the third for x3, …and nth equation for xn. Gauss-Seidel Method For a 3x3 system x1 b1 a12 x2 a13 x3 a11 x2 b2 a21 x1 a23 x3 a22 x3 b3 a31 x1 a32 x2 a33 A simple way to obtain initial guesses for all x’s to start the solution process is to assume that they are zero. These zeros can be substituted into x1 equation to calculate a new x1=b1/a11. Then we substitute this new value of x1 along with the previous guess of zero for x3 into the second equation to find x2. The process is repeated for the third equation to calculate a new estimate for x3 Gauss-Seidel Method • New x1 is obtained by using the new estimates of x2 and x3. The procedure is repeated until the convergence criterion is satisfied: a ,i xij xij xi j 1 100% s For all i, where j and j-1 are the present and previous iterations. Fig. 11.4 Gauss Seidel Jacobi method p.291 k In Gauss-Seidel method, the values x1k 1 , x 2 1 ,..., xik 1 k k k computed in the current iteration as well as xi 2 , xi 3 ,..., x n k1 are used in finding the value xi 1 This implies that always the most recent approximations are used during the computation. The general expression is: k1 i x 1 bi aii i1 j1 aij x k1 j NEW n ji1 aij xik ; i=1,2,...,n; k=1,2,3, ….. OLD In Jacobi iteration method, all the new values are computed using the values at the previous iteration. This implies that both the present and the previous set of values have to be stored. Gauss-Seidel method will improve the storage requirement as well as the convergence. Convergence Criterion for GaussSeidel Method • • The Gauss-Seidel method has two fundamental problems as many iterative method: – It is sometimes nonconvergent, and – If it converges, converges very slowly. Recalling that sufficient conditions for convergence of two linear equations, u(x,y) and v(x,y) are u x u y 1 v x v y 1 • Similarly, in case of two simultaneous equations, the Gauss-Seidel algorithm can be expressed as u ( x1 , x2 ) b1 a11 a12 x2 a11 v( x1 , x2 ) b2 a22 a21 x1 a22 u x1 v x1 u x2 0 a21 a22 v x2 a12 a11 0 • Substitution into convergence criterion of two linear equations yield: a12 a11 a21 a22 1 1 • In other words, the absolute values of the slopes must be less than unity for convergence: a11 a12 a22 a21 For n equations : n aii ai , j j1 ji Figure 11.5 ...
View Full Document

## This note was uploaded on 07/11/2011 for the course ENGR 391 taught by Professor Hoidick during the Winter '09 term at Concordia Canada.

Ask a homework question - tutors are online