1 / 29

MA/CS 375

MA/CS 375. Fall 2003 Lecture 19. Matrices and Systems. Upper Triangular Systems. Upper Triangular Systems. Volunteer: tell us how to solve this? ( hint: this is not difficult ). Example: Upper Triangular Systems. Upper Triangular Systems. 1) Solve the N th equation

plato-slater
Download Presentation

MA/CS 375

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. MA/CS 375 Fall 2003 Lecture 19 MA/CS 375 Fall 2003

  2. Matrices and Systems MA/CS 375 Fall 2003

  3. Upper Triangular Systems MA/CS 375 Fall 2003

  4. Upper Triangular Systems Volunteer: tell us how to solve this? (hint: this is not difficult) MA/CS 375 Fall 2003

  5. Example: Upper Triangular Systems MA/CS 375 Fall 2003

  6. Upper Triangular Systems 1) Solve the Nth equation 2) Use the value from step 1 to solve the (N-1)th equation 3) Use the values from steps 1 and 2 to solve the (N-2)th equation 4) keep going until you solve the whole system. MA/CS 375 Fall 2003

  7. Backwards Substitution MA/CS 375 Fall 2003

  8. Matlab Code for Backward Substitution MA/CS 375 Fall 2003

  9. Testing Backwards Substitution MA/CS 375 Fall 2003

  10. Finite Precision Effects • Notice that x(1) depends on the values of x(2),x(3),…,x(N) • Remembering that round off rears its ugly head when we deal with (small+large) numbers • Volunteer to design a U matrix and b vector, which do not give good answers when used in this algorithm. • (goal is to give you some intuition for cases which will break an algorithm) MA/CS 375 Fall 2003

  11. Team Exercise • Who can build a “reasonable”, non-singular, U matrix and b vector which has the worst answer for U\b • 10 Minutes only MA/CS 375 Fall 2003

  12. LU Factorization • By now you should be persuaded that it is not difficult to solve a problem of the form Lx=b or Ux=b using forward/backward elimination • Suppose we are given a matrix A and we are able to find a lower triangular matrix L and an upper triangular matrix U such that: A = LU • Then to solve Ax=b we can solve two simple problems: • Ly = b • Ux = y • (sanity check LUx = L(Ux) = Ly =b) MA/CS 375 Fall 2003

  13. Example LU Factorization • Suppose we are given: • Then we can write A = LU where: • Let’s check that:  MA/CS 375 Fall 2003

  14. LU Factorization (in words) • Construct a sequence of multiplier matrices (M1, M2, M3, M4,…., MN-1) such that: • (MN-1…. M4 M3 M2 M1)A is upper triangle. • A suitable candidate for M1 is the matrix ( 4 by 4 case): MA/CS 375 Fall 2003

  15. Check to See What M1 Does to A MA/CS 375 Fall 2003

  16. M1 AHas Zeros Below The Diagonal MA/CS 375 Fall 2003

  17. LU Factorization cont. A suitable candidate for M2 is the matrix ( 4 by 4 case): MA/CS 375 Fall 2003

  18. LU Factorization cont. A suitable candidate for M3 is the matrix ( 4 by 4 case): MA/CS 375 Fall 2003

  19. LU Factorization cont. • So in this case: U=(M3M2M1)A is upper triangle • Each of the M matrices is lower triangle • The inverse of M2 looks like: • Using the fact that the inverse of the each M matrix just requires us to negate the below diagonal terms • A = (M3M2M1)-1 U = (M1)-1 (M2)-1 (M3)-1 U • Define L = (M1)-1 (M2)-1 (M3)-1 and we are done since the product of lower matrices is a lower matrix  MA/CS 375 Fall 2003

  20. Matlab Routine ForLU Factorization MA/CS 375 Fall 2003

  21. Matlab Routine ForLU Factorization • Build L and U • A is modified in place • Total cost O(N3) ( we can get this by noticing the triply nested loop) MA/CS 375 Fall 2003

  22. Matlab Routine ForLU FactorizationAnd Solution of Ax=b • Build rhs • Solve Ly=b • Solve Ux=y MA/CS 375 Fall 2003

  23. Team Exercise • Code up the factorization • Test your code to make sure it gives thecorrect answer as Matlab • Make sure abs(LU – Aorig) is small • Make sure abs(x - Aorig\b) is small • Make sure the routine actually finished without the break. • Time the code for an N=400 case • Compare with the timings for [Lmat,Umat] = lu(Aorig); • ASK IF YOU DO NOT UNDERSTAND ANY PART OF THIS!! MA/CS 375 Fall 2003

  24. When Basic LU Does Not Work • Remember that we are required to compute 1/A(i,i) for all the successive diagonal terms • What happens when we get to the i’th diagonal term and this turns out to be 0 or very small ?? • At any point we can swap the i’th row with one lower down which does not have a zero (or small entry) in the i’th column. • This is known as partial pivoting. MA/CS 375 Fall 2003

  25. Partial Pivoting Remember before that we constructed a setof M matrices so that is upper triangle Now we construct a sequence of M matricesand P matrices so that:is upper triangle. (MN-1…. M4 M3 M2 M1)A (MN-1 PN-1 …. M4 P4 M3 P3 M2 P2 M1 P1)A MA/CS 375 Fall 2003

  26. LU with Partial Pivoting MA/CS 375 Fall 2003

  27. LU with Partial Pivoting • find pivot • swap rows • Do elimination MA/CS 375 Fall 2003

  28. LU with Partial Pivoting • build matrix A • call our LU routine • check that: • A(piv,:) = L*U • Looks good  MA/CS 375 Fall 2003

  29. Team Exercise • Make sure you all have a copy of the .m files on your laptops (courtesy of Coutsias) • Locate GEpiv • For N=100,200,400,600 • build a random NxN matrix A • time how long GEpiv takes to factorize A • end • Plot the cputime taken per test against N using loglog • On the same graph plot N^3 using loglog • Observation MA/CS 375 Fall 2003

More Related