1 / 39

Chapter 12 Gaussian Elimination (II)

Chapter 12 Gaussian Elimination (II). Speaker: Lung-Sheng Chien. Reference book: David Kincaid, Numerical Analysis Reference lecture note: Wen-wei Lin, chapter 2, matrix computation http://math.ntnu.edu.tw/~min/matrix_computation.html. OutLine.

madra
Download Presentation

Chapter 12 Gaussian Elimination (II)

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. Chapter 12 Gaussian Elimination (II) Speaker: Lung-Sheng Chien Reference book: David Kincaid, Numerical Analysis Reference lecture note: Wen-wei Lin, chapter 2, matrix computation http://math.ntnu.edu.tw/~min/matrix_computation.html

  2. OutLine • Weakness of A=LU- erroneous judgment: A is invertible but A=LU does not exist- unstable, A is far from LU • A=LU versus PA=LU • Pivoting strategy • Implementation of PA=LU • MATLAB usage

  3. Fail of LU: singular of leading principal minor Question: How about interchanging row 1 and row 2? and

  4. unstable of LU: near singular of leading principal minor [1] ? Theoretical: , why? Numerical: Problem: for double-precision, we only have 16 digit-accuracy, we cannot accept if then

  5. unstable of LU: near singular of leading principal minor [2] backward substitution step 1: step 2: Question: why not due to rounding error due to rounding error

  6. unstable of LU: near singular of leading principal minor [3] Case 1: Case 2: wrong solution

  7. unstable of LU: near singular of leading principal minor [4] Question: How about interchanging row 1 and row 2? LU-factorization backward substitution step 1: step 2: Key observation: without pivoting pivoting (rounding error does not occur for normal number) (rounding normal number)

  8. unstable of LU: near singular of leading principal minor [5] Case 1: Case 2: Why not 1?

  9. unstable of LU: cause and controllability define largest component of matrix without pivoting pivoting (stable) (NOT stable) Objective: control growth rate of control multiplier

  10. OutLine • Weakness of A=LU • A=LU versus PA=LU- controllability of lower triangle matrix L • Pivoting strategy • Implementation of PA=LU • MATLAB usage

  11. Recall LU example in chapter 11 [1] 6 -2 2 4 1 6 -2 2 4 12 -8 6 10 2 1 0 -4 2 2 3 -13 9 3 0.5 1 0 -12 8 1 -6 4 1 -18 -1 1 0 2 3 -14 is not maximum among since 6 -2 2 4 1 6 -2 2 4 0 -4 2 2 1 0 -4 2 2 0 -12 8 1 3 1 0 0 2 -5 0 2 3 -14 -0.5 1 0 0 4 -13 is not maximum among since

  12. Recall LU example in chapter 11 [2] 6 -2 2 4 1 6 -2 2 4 0 -4 2 2 1 0 -4 2 2 0 0 2 -5 1 0 0 2 -5 0 0 0 -3 2 1 0 0 4 -13 is not maximum among since 6 -2 2 4 6 -2 2 4 1 12 -8 6 10 -4 2 2 2 1 3 -13 9 3 2 -5 0.5 3 1 -6 4 1 -18 -3 -1 -0.5 2 1 Question: How can we control , is uniform bound possible?

  13. PA = LU in MATLAB

  14. PA = LU: assume P is given [1] 12 -8 6 10 1 12 -8 6 10 3 -13 9 3 0.25 1 0 -11 7.5 0.5 -6 4 1 -18 -0.5 1 0 0 4 -13 6 -2 2 4 0.5 1 0 2 -1 -1 since is maximum among 12 -8 6 10 1 12 -8 6 10 0 -11 7.5 0.5 1 0 -11 7.5 0.5 0 0 4 -13 0 1 0 0 4 -13 0 2 -1 -1 -2/11 1 0 0 4/11 -10/11 is maximum among since

  15. PA = LU: assume P is given [2] 12 -8 6 10 1 12 -8 6 10 0 -11 7.5 0.5 1 0 -11 7.5 0.5 0 0 4 -13 1 0 0 4 -13 0 0 4/11 -10/11 1/11 1 0 0 0 3/11 is maximum among since 12 -8 6 10 1 12 -8 6 10 3 -13 9 3 0.25 1 -11 7.5 0.5 -6 4 1 -18 -0.5 0 1 4 -13 6 -2 2 4 0.5 -2/11 1/11 1 3/11 Observation: though proper permutation, we can control Question: in fact, we cannot know permutation matrix in advance, how can we do?

  16. Sequence of matrices during LU-decomposition

  17. OutLine • Weakness of A=LU • A=LU versus PA=LU • Pivoting strategy- partial pivoting (we adopt this formulation)- complete pivoting • Implementation of PA=LU • MATLAB usage

  18. Partial pivoting and complete pivoting Partial pivoting such that (1) find a (2) swap row and row then with complete pivoting such that (1) find a (2) swap row and row (2) swap column and column then with

  19. Recall permutation matrix Let Define permutation matrix interchange row 2, 3 interchange column 2, 3 Symmetric permutation: 1 2 3 1 2 3 1 3 2 interchange row 2, 3 interchange column 2, 3 4 5 6 7 8 9 7 9 8 7 8 9 4 5 6 4 6 5

  20. Meaning of matrix coefficient 1 relationship between node and node 4 7 x1 x1 x2 x3 constraint on node 2 3 1 2 3 x3 x2 x1 node is named 4 5 6 x2 6 5 9 7 8 9 x3 8 change node 2 to node 3 and node 3 to node 2 1 write constraint for now node index 4 7 x1 x1 x2 x3 2 3 x2 x3 1 3 2 x1 6 7 9 8 x2 5 9 8 4 6 5 x3

  21. Identity matrix is invariant under symmetric permutation 1 x1 x2 x3 1 x1 x1 1 x2 1 x3 x2 1 x3 1 change node 2 to node 3 and node 3 to node 2 1 x1 x2 x3 x1 1 x1 1 1 x2 x2 x3 1 x3 1

  22. PA = LU: partial pivoting [1] 6 -2 2 4 12 -8 6 10 12 -8 6 10 6 -2 2 4 3 -13 9 3 3 -13 9 3 -6 4 1 -18 -6 4 1 -18 12 -8 6 10 1 12 -8 6 10 6 -2 2 4 0.5 1 0 2 -1 -1 3 -13 9 3 0.25 1 0 -11 7.5 0.5 -6 4 1 -18 -0.5 1 0 0 4 -13 12 -8 6 10 12 -8 6 10 0 -11 7.5 0.5 0 2 -1 -1 0 2 -1 -1 0 -11 7.5 0.5 0 0 4 -13 0 0 4 -13

  23. PA = LU: partial pivoting [2] where 1 1 1 0.25 1 0.5 1 0.5 1 0.5 1 0.25 1 Interchange columns 2, 3 0.25 1 Interchange rows 2, 3 -0.5 1 -0.5 1 -0.5 1 12 -8 6 10 1 12 -8 6 10 0 -11 7.5 0.5 0.25 1 3 -13 9 3 verify 0 2 -1 -1 0.5 1 6 -2 2 4 0 0 4 -13 -0.5 1 -6 4 1 -18

  24. PA = LU: partial pivoting [3] 12 -8 6 10 1 12 -8 6 10 0 -11 7.5 0.5 1 0 -11 7.5 0.5 0 2 -1 -1 -2/11 1 0 0 4/11 -10/11 0 0 4 -13 0 1 0 0 4 -13 12 -8 6 10 12 -8 6 10 0 -11 7.5 0.5 0 -11 7.5 0.5 0 0 4/11 -10/11 0 0 4 -13 0 0 4 -13 0 0 4/11 -10/11 where

  25. PA = LU: partial pivoting [4] 1 1 1 0.25 1 0.25 1 0.25 1 -0.5 0 1 0.5 -2/11 1 Interchange columns 3,4 0.5 -2/11 1 Interchange rows 3.4 0.5 -2/11 1 -0.5 1 -0.5 0 1 12 -8 6 10 1 12 -8 6 10 0 -11 7.5 0.5 0.25 1 3 -13 9 3 verify 0 0 4 -13 -0.5 0 1 -6 4 1 -18 0 0 4/11 -10/11 0.5 -2/11 1 6 -2 2 4

  26. PA = LU: partial pivoting [5] 12 -8 6 10 12 -8 6 10 1 0 -11 7.5 0.5 0 -11 7.5 0.5 1 0 0 4 -13 0 0 4 -13 1 0 0 4/11 -10/11 0 0 0 3/11 1/11 1 we have where 1 12 -8 6 10 0.25 1 -11 7.5 0.5 -0.5 0 1 4 -13 0.5 -2/11 1/11 1 3/11

  27. OutLine • Weakness of A=LU • A=LU versus PA=LU • Pivoting strategy • Implementation of PA=LU- commutability of GE matrix and permutation- recursive formula • MATLAB usage

  28. Implementation issue: commutability of GE matrix and permutation [1] 1 1 0.25 1 0.5 1 Interchange rows and columns 2, 3 0.5 1 0.25 1 -0.5 1 -0.5 1 1 1 0.25 1 0.25 1 Interchange rows and columns 3, 4 -0.5 0 1 0.5 -2/11 1 0.5 -2/11 1 -0.5 1

  29. Implementation issue: commutability of GE matrix and permutation [2] Observation: If we define , then where interchanges rows or columns and interchange

  30. Implementation issue: commutability of GE matrix and permutation [3] partial rows (k and p) interchange symmetric permutation on an Identity matrix

  31. Algorithm ( PA = LU ) [1] Given full matrix , construct initial lower triangle matrix use permutation vector to record permutation matrix let , and for we have compute update original matrix stores in lower triangle matrix

  32. Algorithm ( PA = LU ) [2] such that 1 find a 2 swap row and row define permutation matrix , we have then after swapping rows 3 compute for if then by swapping and compute 4 then endif

  33. Algorithm ( PA = LU ) [3] 5 decompose where by updating by updating matrix (recursion is done) then endfor

  34. Question: recursion implementation • Initialization- check algorithm holds for k=1 • Recursion formula- check algorithm holds for k, if k-1 is true • Termination condition- check algorithm holds for k=n • Breakdown of algorithm- check which condition PA=LU fails • Exception: algorithm works only for square matrix? normal abnormal Extension of algorithm

  35. Exception: algorithm works only for square matrix? [1] Case 1: Case 2:

  36. Exception: algorithm works only for square matrix? [2] Case 3: we can simplify it as Question: what is termination condition of case 3 ?

  37. OutLine • Weakness of A=LU • A=LU versus PA=LU • Pivoting strategy • Implementation of PA=LU • MATLAB usage- abs- max

  38. Command “abs”: take absolute value In PA=LU algorithm, we need to take absolute value of one column vector for pivoting abs also works for matrix

  39. Command “max”: take maximum value

More Related