The Islamic University of Gaza Faculty of Engineering Civil Engineering Department

1. The Islamic University of Gaza Faculty of Engineering Civil Engineering Department Numerical Analysis ECIV 3306 Chapter10 LU Decomposition and Matrix Inversion

2. Introduction • Gauss elimination solves [A] {x} ={B} • It becomes insufficient when solving these equations for different values of {B} • LU decomposition works on the matrix [A] and the vector {B} separately. • LU decomposition is very useful when the vector of variables {x} is estimated for different parameter vectors {B} since the forward elimination process is not performed on {B}.

3. LU Decomposition If: L: lower triangular matrix U: upper triangular matrix Then, [A]{X}={B} can be decomposed into two matrices [L] and [U] such that: 1. [L][U] = [A] ([L][U]){X} = {B}

4. LU Decomposition Consider: [U]{X} = {D} So, [L]{D} = {B} 2. [L]{D} = {B} is used to generate an intermediate vector {D} by forward substitution. 3. Then, [U]{X}={D} is used to get {X} by back substitution.

5. Summary of LU Decomposition

6. LU Decomposition As in Gauss elimination, LU decomposition must employ pivoting to avoid division by zero and to minimize round off errors. The pivoting is done immediately after computing each column.

7. LU Decomposition System of linear equations [A]{x}={B} Step 1:Decomposition

8. LU Decomposition Step 2:Generate an intermediate vector {D} by forward substitution Step 3:Get {X} by back substitution.

9. LU Decomposition-Example

10. LU Decomposition-Example (cont’d) Use previous L and D matrices to solve the system: Step 2:Find the intermediate vector {D} by forward substitution

11. LU Decomposition-Example (cont’d) Step 3:Get {X} by back substitution.

12. Decomposition Step • % Decomposition Step for k=1:n-1 [a,o]= pivot(a,o,k,n); for i = k+1:n a(i,k) = a(i,k)/a(k,k); a(i,k+1:n)= a(i,k+1:n)-a(i,k)*a(k,k+1:n); end end

13. Partial Pivoting • %Partial Pivoting function [a,o] = pivot(a,o,k,n) [big piv]=max(abs(a(k:n,k))); piv=piv+(k-1); if piv ~= k temp = a(piv,:); a(piv,:)= a(k,:); a(k,:)=temp; temp = o(piv); o(piv)=o(k); o(k)=temp; end

14. Substitution Steps %Forward Substitution d(1)=bn(1); for i=2:n d(i)=bn(i)-a(i,1:i-1)*(d(1:i-1))'; end • % Back Substitution x(n)=d(n)/a(n,n); for i=n-1:-1:1 x(i)=(d(i)-a(i,i+1:n)*x(i+1:n)')/a(i,i); end

15. Matrix Inverse Using the LU Decomposition • LU decomposition can be used to obtain the inverse of the original coefficient matrix [A]. • Each column j of the inverse is determined by using a unit vector (with 1 in the jth raw ).

16. 2nd column of [A]-1 1st column of [A]-1 3rd column of [A]-1 Matrix Inverse: LU Decomposition [A][A]-1 = [A]-1[A] = I

17. Matrix inverse using LU decomposition Example • 1A. [L]{d}1 = {b}1 1st column of [A]-1 • 1B. Then, [U]{X}1={d}1

18. Matrix inverse using LU decomposition Example (cont’d) • 2A. [L]{d}2 = {b}2 2nd column of [A]-1 • 2B. Then, [U]{X}2={d}2

19. Matrix inverse using LU decomposition Example (cont’d) • 3A. [L]{d}3 = {b}3 • 3B. Then, [U]{X}3={d}3 3rd column of [A]-1

20. Matrix inverse using LU decomposition Example (cont’d)

21. Vector and Matrix Norms Norm is a real-valued function that provides a measure of size or “length” of vectors and matrices. Norms are useful in studying the error behavior of algorithms.

22. Vector and Matrix Norms (cont’d)

23. Vector and Matrix Norms (cont’d) • The length of this vector can be simply computed as Length or Euclidean norm of [F] • For an n dimensional vector Frobenius norm

24. Vector and Matrix Norms (cont’d)

25. Vector and Matrix Norms (cont’d) • Uniform vector norm • Uniform matrix norm (row sum Norm)

26. Vector and Matrix Norms (cont’d) • Matrix Condition umber Defined as: • For a matrix [A], this number will be greater than or equal to 1. • If the coefficients of [A] are known to t-digit precision (rounding errors~10-t) and Cond [A]=10c, the solution [X] may be valid to only t-c digits (rounding errors~10c-t).

27. Iterative Refinement • Round-off errors can be reduced by the following procedure: • Suppose an approximate solution vectors given by • Substitute the result in the original system ....(Eq.1)

28. Iterative Refinement (cont’d) • Now, assume the exact solution is: • Then: … …….(Eq.2)

29. Iterative Refinement (cont’d) Subtract Eq.2 from Eq.1, the result is a system of linear equations that can be solved to obtain the correction factors The factors then can be applied to improve the solution as specified by the equation:

30. Iterative Refinement- Example (cont’d) Solve: The exact solution is ……… 1- Solve the equations using [A]-1, such as {x}=[A]-1{c}

31. Iterative Refinement - Example (cont’d) 2- Substitute the result in the original system [A]{x}={c} 

32. Iterative Refinement - Example (cont’d) 3-Solve the system of linear equations using [A]-1 to find the correction factors to yield