1 / 86

Engineering Computation

Engineering Computation. Part 3. Learning Objectives for Lecture. 1. Motivate Study of Systems of Equations and particularly Systems of Linear Equations 2. Review steps of Gaussian Elimination 3. Examine how roundoff error can enter and be magnified in Gaussian Elimination

gur
Download Presentation

Engineering Computation

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. EngineeringComputation Part 3 E. T. S. I. Caminos, Canales y Puertos

  2. Learning Objectives for Lecture • 1. Motivate Study of Systems of Equations and particularly Systems of Linear Equations • 2. Review steps of Gaussian Elimination • 3. Examine how roundoff error can enter and • be magnified in Gaussian Elimination • 4. Introduce Pivoting and Scaling as defenses against roundoff. • 5. Consider what an engineer can do to generate well formulated problems. E. T. S. I. Caminos, Canales y Puertos

  3. Systems of Equations • In Part 2 we have tried to determine the value x, satisfying f(x)=0. In this part we try to obtain the values x1,x2, xn, satisfying the system of equations: • These systems can be linear or nonlinear, but in this part we deal with linear systems: E. T. S. I. Caminos, Canales y Puertos

  4. Systems of Equations • where a and b are constant coefficients, and n is the number of equations. • Many of the engineering fundamental equations are based on conservation laws. In mathematical terms, these principles lead to balance or continuity equations relating the system behavior with respect to the amount of the magnitude being modelled and the extrenal stimuli acting on the system. E. T. S. I. Caminos, Canales y Puertos

  5. Systems of Equations • Matrices are rectangular sets of elements represented by a single symbol. If the set if horizontal it is called row, and if it is vertical, it is called column. Column 3 Row 2 Column vector Row vector E. T. S. I. Caminos, Canales y Puertos

  6. Systems of Equations • There are some special types of matrices: Identity matrix Symmetric matrix Diagonal matrix Upper triangular matrix E. T. S. I. Caminos, Canales y Puertos

  7. Systems of Equations Lower triangular matrix Banded matrix Half band width All elements are null with the exception of thoise in a band centered around the main diagonal. This matrix has a band width of 3 and has the name of tridiagonal. E. T. S. I. Caminos, Canales y Puertos

  8. Systems of Equations • Linear Algebraic Equations • a11x1 + a12x2 + a13x3+ … + a1nxn = b1 • a21x1 + a22x2 + a23x3+ … + a2nxn = b2 • ….. • an1x1 + an2x2 + an3x3+ … + anxn = bn • where all aij's and bi'sare constants. • In matrix form: n x n n x 1 n x 1 or simply [A]{x} = {b} E. T. S. I. Caminos, Canales y Puertos

  9. Systems of Equations • Matrix representation of a system Matrix product: Resulting dimensions E. T. S. I. Caminos, Canales y Puertos

  10. Systems of Equations • Graphic Solution: Systems of equations are hyperplanes (straight lines, planes, etc.). The solution of a system is the intersection of these hyperplanes. Compatible and determined system. Vectors are linearly independent. Unique solution. Determinant of A is non-null. E. T. S. I. Caminos, Canales y Puertos

  11. Systems of Equations Incompatible system, Linearly dependent vectors. Null determinant of A. There is no solution. Compatible but undetermined system. Linearly dependent vectors. Null determinant of A. There exists an infinite number of solutions. E. T. S. I. Caminos, Canales y Puertos

  12. Systems of Equations Compatible and determined system. Linearly independent vectors. Nonnull determinant of A, but close to zero. There exists a solution but it is difficult to find precisely. It is an ill conditioned system leading to numerical errors. E. T. S. I. Caminos, Canales y Puertos

  13. Gauss elimination • Naive Gauss elimination method: The Gauss’ method has two phases: Forward elimination and backsustitution. In the first, the system is reduced to an upper triangular system: • First, the unknown x1 is eliminated. To this end, the first row is multiplied by -a21/a11 and added to the second row. The same is done with all other succesive rows (n-1 times) until only the first equation contains the first unknown x1. Pivot equation substract pivot E. T. S. I. Caminos, Canales y Puertos

  14. Gauss elimination • This operation is repeated with all variables xi, until an upper triangular matrix is obtained. • Next, the system is solved by backsustitution. • The number of operations (FLOPS) used in the Gauss method is: Pass 1 Pass2 E. T. S. I. Caminos, Canales y Puertos

  15. Gauss elimination • 1. Forward Elimination (Row Manipulation): • a. Form augmented matrix [A|b]: b. By elementary row manipulations, reduce [A|b] to [U|b'] where U is an upper triangular matrix: DO i = 1 to n-1 DO k = i+1 to n Row(k) = Row(k) - (aki/aii)*Row(i) ENDDO ENDDO E. T. S. I. Caminos, Canales y Puertos

  16. Gauss elimination • 2. Back Substitution • Solve the upper triangular system [U]{x} = {b´} xn = b'n / unn DO i = n-1 to 1 by (-1) END E. T. S. I. Caminos, Canales y Puertos

  17. Gauss elimination (example) Consider the system of equations To 2 significant figures, the exact solution is: We will use 2 decimal digit arithmetic with rounding. E. T. S. I. Caminos, Canales y Puertos

  18. Gauss elimination (example) Start with the augmented matrix: Multiply the first row by –1/50 and add to second row. Multiply the first row by –2/50 and add to third row: Multiply the second row by –6/40 and add to third row: E. T. S. I. Caminos, Canales y Puertos

  19. Gauss elimination (example) Now backsolve: (vs. 0.091, et = 2.2%) (vs. 0.041, et = 2.5%) (vs. 0.016, et = 0%) E. T. S. I. Caminos, Canales y Puertos

  20. Gauss elimination (example) • Consider an alternative solution • interchanging rows: • After forward elimination, we obtain: • Now backsolve: • x3 = 0.095 (vs. 0.091, et = 4.4%) • x2 = 0.020 (vs. 0.041, et = 50%) • x1 = 0.000 (vs. 0.016, et = 100%) • Apparently, the order of the equations matters! E. T. S. I. Caminos, Canales y Puertos

  21. Gauss elimination (example) • WHAT HAPPENED? • When we used 50 x1 + 1 x2 + 2 x3 = 1 to solve for x1, there was little change in other equations. • When we used 2 x1 + 6 x2 + 30 x3 = 3 to solve for x1 it made BIG changes in the other equations. Some coefficients for other equations were lost! • The second equation has little to do with x1. • It has mainly to do with x3. • As a result we obtained LARGE numbers in the table, significant roundoff error occurred and information was lost. • Things didn't go well! • If scaling factors | aji / aii | are  1 then the effect of roundoff errors is diminished. E. T. S. I. Caminos, Canales y Puertos

  22. Gauss elimination (example) Effect of diagonal dominance: As a first approximation roots are: xi bi / aii Consider the previous examples: E. T. S. I. Caminos, Canales y Puertos

  23. Gauss elimination (example) • Goals:1. Best accuracy (i.e. minimize error) • 2. Parsimony (i.e. minimize effort) • Possible Problems: • A. Zero on diagonal term  ÷ by zero. • B. Many floating point operations (flops) cause numerical precision problems and propagation of errors. • C. System may be ill-conditioned: det[A]  0. • D. No solution or an infinite # of solutions: det[A] = 0. • Possible Remedies: • A. Carry more significant figures (double precision). • B. Pivot when the diagonal is close to zero. • C. Scale to reduce round-off error. E. T. S. I. Caminos, Canales y Puertos

  24. Gauss elimination (pivoting) • PIVOTING • A. Row pivoting(Partial Pivoting) - • In any good routine, at each step i, find • maxk | aki | for k = i, i+1, i+2, ..., n • Move corresponding row to pivot position. • (i) Avoids zero aii • (ii) Keeps numbers small & minimizes round-off, • (iii) Uses an equation with large | aki | to find xi • Maintains diagonal dominance. • Row pivoting does not affect the order of the variables. • Included in any good Gaussian Elimination routine. E. T. S. I. Caminos, Canales y Puertos

  25. Gauss elimination (pivoting) • B. Column pivoting - • Reorder remaining variables xj • for j = i, . . . ,n so get largest | aji | • Column pivoting changes the order of the unknowns, xi, and thus leads to complexity in the algorithm. Not usually done. • C. Complete or Full pivoting • Performing both row pivoting and column pivoting. • (If [A] is symmetric, needed to preserve symmetry.) E. T. S. I. Caminos, Canales y Puertos

  26. Gauss elimination (pivoting) • How to fool pivoting: • Multiply the third equation by 100 and then performing pivoting will yield: • Forward elimination then yields (2-digit arithmetic): Backsolution yields: x3 = 0.095 (vs. 0.091, et = 4.4%) x2 = 0.020 (vs. 0.041, et = 50.0%) x1 = 0.000 (vs. 0.016, et = 100%) The order of the rows is still poor!! E. T. S. I. Caminos, Canales y Puertos

  27. Gauss elimination (scaling) • SCALING • A. Express all equations (and variables) in comparable units so all elements of [A] are about the same size. • B. If that fails, and maxj |aij| varies widely across the rows, replace each row i by: • aij • This makes the largest coefficient |aij| of each equation equal to 1 and the largest element of [A] equal to 1 or -1 • NOTE: Routines generally do not scale automatically; scaling can cause round-off error too! • SOLUTIONS • • Don't actually scale, but use hypothetical scaling factors to determine what pivoting is necessary. • • Scale only by powers of 2: no roundoff or division required. E. T. S. I. Caminos, Canales y Puertos

  28. Gauss elimination (scaling) • How to fool scaling: • A poor choice of units can undermine the value of scaling. • Begin with our original example: • If the units of x1 were expressed in µg instead of mg the matrix might read: • Scaling then yields: • Which equation is used to determine x1? Why bother to scale ? E. T. S. I. Caminos, Canales y Puertos

  29. Gauss elimination (operation counting) • OPERATION COUNTING • In numerical scientific calculations, the number of multiplies & divides often determines CPU time. (This represents the numerical effort!) • One floating point multiply or divide (plus any associated adds or subtracts) is called a FLOP. (The adds/subtracts use little time compared to the multiplies/divides.) • FLOP = FLoating point OPeration. • Examples: a * x + b • a / x – b E. T. S. I. Caminos, Canales y Puertos

  30. Gauss elimination (operation counting) Useful identities in counting FLOPS: O(mn) means that there are terms of order mn and lower. E. T. S. I. Caminos, Canales y Puertos

  31. Gauss elimination (operation counting) • Simple Example of Operation Counting: • DO i = 1 to n • Y(i) = X(i)/i – 1 • ENDDO • X(i) and Y(i) are arrays whose values change when i changes. In each iteration • X(i)/i – 1 • represents one FLOP because it requires one division (& one subtraction). • The DO loop extends over i from 1 to n iterations: E. T. S. I. Caminos, Canales y Puertos

  32. Gauss elimination (operation counting) Another Example of Operation Counting: DO i = 1 to n Y(i) = X(i)X(i) + 1 DO j = i to n Z(j) = [ Y(j) / X(i) ]Y(j) + X(i) ENDDO ENDDO With nested loops, always start from the innermost loop. [Y(j)/X(i)] * Y(j) + X(i) represents 2 FLOPS E. T. S. I. Caminos, Canales y Puertos

  33. Gauss elimination (operation counting) For the outer i-loop:X(i)• X(i) + 1 represents 1 FLOP = 3n +2n2 - n2 - n = n2 + 2n = n2 + O(n) E. T. S. I. Caminos, Canales y Puertos

  34. Gauss elimination (operation counting) Forward Elimination: DO k = 1 to n–1 DO i = k+1 to n r = A(i,k)/A(k,k) DO j = k+1 to n A(i,j)=A(i,j) – r*A(k,j) ENDDO B(i) = B(i) – r*B(k) ENDDO ENDDO E. T. S. I. Caminos, Canales y Puertos

  35. Gauss elimination (operation counting) Operation Counting for Gaussian Elimination Back Substitution: X(n) = B(n)/A(n,n) DO i = n–1 to 1 by –1 SUM = 0 DO j = i+1 to n SUM = SUM + A(i,j)*X(j) ENDDO X(i) = [B(i) – SUM]/A(i,i) ENDDO E. T. S. I. Caminos, Canales y Puertos

  36. Gauss elimination (operation counting) Operation Counting for Gaussian Elimination Forward Elimination Inner loop: Second loop: = (n2 + 2n) – 2(n + 1)k + k2 E. T. S. I. Caminos, Canales y Puertos

  37. Gauss elimination (operation counting) Operation Counting for Gaussian Elimination Forward Elimination (cont'd) Outer loop = E. T. S. I. Caminos, Canales y Puertos

  38. Gauss elimination (operation counting) Operation Counting for Gaussian Elimination Back Substitution Inner Loop: Outer Loop: E. T. S. I. Caminos, Canales y Puertos

  39. Gauss elimination (operation counting) Total flops = Forward Elimination + Back Substitution = n3/3 + O (n2) + n2/2 + O (n)  n3/3 + O (n2) To convert (A,b) to (U,b') requires n3/3, plus terms of order n2 and smaller, flops. To back solve requires: 1 + 2 + 3 + 4 + . . . + n = n (n+1) / 2 flops; Grand Total: the entire effort requires n3/3 + O(n2) flops altogether. E. T. S. I. Caminos, Canales y Puertos

  40. Gauss-Jordan Elimination • Diagonalization by both forward and backward elimination in each column. • Perform elimination both backwards and forwards until: • Operation count for Gauss-Jordan is: • (slower than Gauss elimination) E. T. S. I. Caminos, Canales y Puertos

  41. Gauss-Jordan Elimination Example (two-digit arithmetic): x1 = 0.015 (vs. 0.016, et = 6.3%) x2 = 0.041 (vs. 0.041, et = 0%) x3 = 0.093 (vs. 0.091, et = 2.2%) E. T. S. I. Caminos, Canales y Puertos

  42. Gauss-Jordan Matrix Inversion The solution of: [A]{x} = {b} is: {x} = [A]-1{b} where [A]-1 is the inverse matrix of [A] Consider: [A] [A]-1 = [ I ] 1) Create the augmented matrix: [ A | I ] 2) Apply Gauss-Jordan elimination: ==> [ I | A-1 ] E. T. S. I. Caminos, Canales y Puertos

  43. Gauss-Jordan Matrix Inversion Gauss-Jordan Matrix Inversion (with 2 digit arithmetic): MATRIX INVERSE [A-1] E. T. S. I. Caminos, Canales y Puertos

  44. Gauss-Jordan Matrix Inversion CHECK: [ A ] [ A]-1 = [ I ] [ A]-1 { b } = { x } Gaussian Elimination E. T. S. I. Caminos, Canales y Puertos

  45. LU decomposition • LU decomposition - The LU decomposition is a method that uses the elimination techniques to transform the matrix A in a product of triangular matrices. This is specially useful to solve systems with different vectors b, because the same decomposition of matrix A can be used to evaluate in an efficient form, by forward and backward sustitution, all cases. E. T. S. I. Caminos, Canales y Puertos

  46. Initial system Decomposition Transformed system 1 Substitution Transformed system 2 Backward sustitution Forward sustitution LU decomposition E. T. S. I. Caminos, Canales y Puertos

  47. LU decomposition • LU decomposition is very much related to Gauss method, because the upper triangular matrix is also looked for in the LU decomposition. Thus, only the lower triangular matrix is needed. • Surprisingly, during the Gauss elimination procedure, this matrix L is obtained, but one is not aware of this fact. The factors we use to get zeroes below the main diagonal are the elements of this matrix L. Substract E. T. S. I. Caminos, Canales y Puertos

  48. LU decomposition resto resto E. T. S. I. Caminos, Canales y Puertos

  49. LU decomposition (Complexity) Basic Approach Consider [A]{x} = {b} a) Gauss-type "decomposition" of [A] into [L][U] n3/3 flops [A]{x} = {b} becomes [L] [U]{x} = {b}; let [U]{x}{d} b) First solve [L] {d} = {b} for {d} by forward subst. n2/2 flops c) Then solve [U]{x} = {d} for {x} by back substitution n2/2 flops E. T. S. I. Caminos, Canales y Puertos

  50. LU Decompostion: notation [A] = [L] + [U0] [A] = [L0] + [U] [A] = [L0] + [U0] + [D] [A] = [L1] [U] [A] = [L] [U1] E. T. S. I. Caminos, Canales y Puertos

More Related