154 Views

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. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**MA/CS 375**Fall 2002 Lecture 21 MA/CS 375 Fall 2002**This Lecture**• Introduce the idea of pivoting • This allows us to deal with the case where the algorithm we discussed broke down • i.e. what to do to avoid the diagonal term we need to invert became very small or even zero. MA/CS 375 Fall 2002**Matrices and Systems**MA/CS 375 Fall 2002**Review**Example: Upper Triangular Systems MA/CS 375 Fall 2002**Review**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 2002**Review**Backwards Substitution MA/CS 375 Fall 2002**Review**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 2002**Review**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 2002**Review**LU Factorization Sequence cont. MA/CS 375 Fall 2002**Review**LU Factorization cont. What does M look like: MA/CS 375 Fall 2002**Review**LU Factorization cont. What does M-1look like: Notice: M and its inverse are very closely related. MA/CS 375 Fall 2002**Review**LU Factorization cont. • In summary: • set L=M-1 • set U=MA • the L is upper triangular • the U is lower triangular • and A = LU, we are done. MA/CS 375 Fall 2002**Review**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 2002**Review**Matlab Routine ForLU FactorizationAnd Solution of Ax=b • Build rhs • Solve Ly=b • Solve Ux=y MA/CS 375 Fall 2002**Class Example**• Try using your LU solver on the following matrix: • The inverse is: • What L and U do you get? MA/CS 375 Fall 2002**New stuff:**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 2002**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 2002**Example of a Permutation Matrix P**What does P do to a matrix A under multiplication ? MA/CS 375 Fall 2002**Example of a Permutation Matrix P**What does P do to a matrix A under multiplication ? MA/CS 375 Fall 2002**Example of a Permutation Matrix P**In this case P has swapped thetwo bottom rows of A. MA/CS 375 Fall 2002**LU with Partial Pivoting**MA/CS 375 Fall 2002**LU with Partial Pivoting**• find pivot • swap rows • Do elimination MA/CS 375 Fall 2002**LU with Partial Pivoting**• build matrix A • call our LU routine • check that: • A(piv,:) = L*U • Looks good MA/CS 375 Fall 2002**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 2002**Notes on LU Factorization**• Big fuss over nothing, right ? • What happens if I want to solve 40000 problems with Ax=b for 40000 different b vectors?. • Cost is O(40000*N2) for the solves and O(N3) for the LU factorization. • i.e. the cost of the factorization is negligable compared to the cost of the solves !!. MA/CS 375 Fall 2002**From Now On**• When you need to solve a square non-symmetric system of equations (Ax=b) with lots of different b vectors use the Matlab built in routine lu: • [L,U,P] = lu(A); • y = L\(Pb); • x = U\y; • For a one time solve use: x=A\b MA/CS 375 Fall 2002**Interlude on Norms**• We all know that the “size” of 2 is the same as the size of -2, namely: • ||2|| = sqrt(2*2) = 2 • ||-2|| = sqrt((-2)*(-2)) = 2 • The ||..|| notation is known as the norm function. MA/CS 375 Fall 2002**Norm of A Vector**• We can generalize the scalar norm to a vector norm, however now there are more choices for a vector : MA/CS 375 Fall 2002**Norm of A Vector**• Matlab command: • norm(x,1) • norm(x,2) • norm(x,inf) MA/CS 375 Fall 2002**Norm of A Matrix**• We can generalize the vector norm to a matrix norm, however now there are more choices for a matrix : MA/CS 375 Fall 2002**Norm of A Matrix**• Matlab command: • norm(A,1) • norm(A,2) • norm(A,inf) • norm(A,fro) MA/CS 375 Fall 2002**Matrix Norm in Action**• Theorem 5: (page 185 van Loan) • If is the stored version of then where and i.e. an error of order ||A||1eps occurs when a real matrix is stored in floating point. MA/CS 375 Fall 2002**Next Lecture**• We will see matrix norms in action… • Those of you who own a digital camera, please bring it in. • Otherwise bring a cd with some picture files on (jpeg, gif, tiff,…) – I have the previously taken pictures too.. • We will start talking about interpolation MA/CS 375 Fall 2002