1 / 36

Scientific Computing

Scientific Computing. Linear Systems – Gaussian Elimination. Linear Systems. Linear Systems. Solve Ax=b, where A is an n  n matrix and b is an n  1 column vector We can also talk about non-square systems where A is m  n , b is m  1, and x is n  1

alaura
Download Presentation

Scientific Computing

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. Scientific Computing Linear Systems – Gaussian Elimination

  2. Linear Systems

  3. Linear Systems • Solve Ax=b, where A is an nnmatrix andb is an n1 column vector • We can also talk about non-square systems whereA is mn, b is m1, and x is n1 • Overdetermined if m>n:“more equations than unknowns” • Underdetermined if n>m:“more unknowns than equations”

  4. Singular Systems • A is singular if some row islinear combination of other rows • Singular systems can be underdetermined:or inconsistent:

  5. Using Matrix Inverses • To solve Ax = b it would be nice to use the inverse to A, that is, A-1 • However, it is usually not a good idea to compute x=A-1b • Inefficient • Prone to roundoff error

  6. Gaussian Elimination • Fundamental operations: • Replace one equation with linear combinationof other equations • Interchange two equations • Multiply one equation by a scalar • These are called elementary row operations. • Do these operations again and again to reduce the system to a “trivial” system

  7. Triangular Form • Two special forms of matrices are especially nice for solving Ax=b: • In both cases, successive substitution leads to a solution

  8. Triangular Form • A is lower triangular

  9. Triangular Form • Solve by forward substitution:

  10. Triangular Form • Solve by forward substitution:

  11. Triangular Form • Solve by forward substitution: Etc

  12. Triangular Form • If A is upper triangular, solve by back-substitution:

  13. Triangular Form • Solve by back-substitution: Etc

  14. Gaussian Elimination Algorithm • Do elementary row operations on the augmented system [A|b] to reduce the system to upper triangular form. • Then, use back-substitution to find the answer.

  15. Gaussian Elimination • Example: • Augmented Matrix form:

  16. Gaussian Elimination • Row Ops: What do we do to zero out first column under first pivot? • Zero out below second pivot:

  17. Gaussian Elimination • Back-substitute

  18. Gaussian Elimination

  19. Matlab Implementation • Task: Implement Gaussian Elimination (without pivoting) in a Matlab M-file. • Notes • Input = Coefficient matrix A, rhs b • Output = solution vector x

  20. Matlab Implementation • Class Exercise: We will go through this code line by line, using the example in Pav section 3.3.2 to see how the code works.

  21. Matlab Implementation • Matlab: • >> A=[2 1 1 3; 4 4 0 7; 6 5 4 17; 2 -1 0 7] • A = • 2 1 1 3 • 4 4 0 7 • 6 5 4 17 • 2 -1 0 7 • >> b = [7 11 31 15]' • b = • 7 • 11 • 31 • 15 >> gauss_no_pivot(A,b) ans = 1.5417 -1.4167 0.8333 1.5000

  22. Matlab Implementation • Class Exercise: How would we change the Matlab function gauss_no_pivot so we could see the result of each step of the row reduction?

  23. Gaussian Elimination - Pivoting • Consider this system: • We immediately run into a problem:we cannot zero out below pivot, or back-substitute! • More subtle problem:

  24. Gaussian Elimination - Pivoting • Conclusion: small diagonal elements are bad! • Remedy: swap the row with the small diagonal element with a row below, this is called pivoting

  25. Gaussian Elimination - Pivoting • Our Example: • Swap rows 1 and 2: • Now continue:

  26. Gaussian Elimination - Pivoting • Two strategies for pivoting: • Partial Pivoting • Scaled Partial Pivoting

  27. Partial Pivoting

  28. Matlab – Partial Pivoting • Partial Pivoting: • At step k, we are working on kth row, pivot = Akk . Search for largest Aik in kth column below (and including) Akk . Let p = index of row containing largest entry. • If p ≠ k, swap rows p and k. • Continue with Gaussian Elimination.

  29. Matlab – Partial Pivoting • Finding largest entry in a column: • >> A • A = • 2 1 1 3 • 4 4 0 7 • 6 5 4 17 • 2 -1 0 7 • >> [r,m] = max(A(2:4,2)) • r = • 5 • m = • 2 • Why isn’t m = 3?

  30. Matlab – Partial Pivoting • Swapping rows m and k: • BB = • 2 1 1 3 • 4 4 0 7 • 6 5 4 17 • 2 -1 0 7 • >> BB([1 3],:) = BB([3 1], :) • BB = • 6 5 4 17 • 4 4 0 7 • 2 1 1 3 • 2 -1 0 7

  31. Matlab – Partial Pivoting • Code change? • All that is needed is in main loop (going from rows k ->n-1) we add • % Find max value (M) and index (m) of max entry below AAkk • [M,m] = max(abs(AA(k:n,k))); • m = m + (k-1); % row offset • % swap rows, if needed • if m ~= k, AA([k m],:) = AA([m k],:); • end

  32. Matlab – Partial Pivoting • Class Exercise • Review example in Moler Section 2.6

  33. Scaled Partial Pivoting Pav, Section 3.2

  34. Practice • Class Exercise: We will work through the example in Pav, Section 3.2.2

  35. Practice • Class Exercise: Do Pav Ex 3.7 for partial pivoting (“naïve Gaussian Elimination”). • Do Pav Ex 3.7 for scaled partial pivoting.

More Related