Methods for Solution of System of Equations: Direct and LU Decomposition Methods

Methods for Solution of the System of Equations: Ax = b Direct Methods: one obtains the exact solution (ignoring the round-off errors) in a finite number of steps. These group of methods are more efficient for dense and banded matrices. Gauss Elimination; Gauss-Jordon Elimination LU-Decomposition Thomas Algorithm (for tri-diagonal banded matrix)

2

LU-Decomposition: A general method Ax = b In most engineering problems, the matrix A remains constant while the vector b changes with time. The matrix A describes the system and the vector b describes the external forcing. e.g., all network problems (pipes, electrical, canal, road, reactors, etc.); structural frames; many financial analyses. If all b s are available together, one can solve the system by augmented matrix but in practice, they are not! Instead of performing n3 floating point operations to solve whenever a new b becomes available, it is possible to solve the system by performing n2 floating point operations if a LU Decomposition is available for matrix A LU-decomposition requires n3 floating point operations!

Consider the system: (b changes!) Ax = b Perform a decomposition of the form A = LU, where L is a lower-triangular and U is an upper-triangular matrix! LU-decomposition requires n3 floating point operations! For any given b, solve Ax = LUx = b This is equivalent to solving two triangular systems: Solve Ly = b using forward substitution to obtain y (~n2 operations) Solve Ux = y using back substitution to obtain x (~n2 operations) Most frequently used method for engineering applications! We will derive LU-decomposition from Gauss Elimination!

? ? ???= ?????? ? = min ?,? ?=1 Define the elements of matrix L as: ???= ??? Define the elements of matrix U as: ???= ??? ? Any element aij is actively modified for the first p steps ? ???= ?????? ? = min ?,? ?=1 This is equivalent to matrix multiplication: A = LU

? ???= ?????? ? = min ?,? ?=1 ?11 0 0 ?12 ?22 0 ?13 ?23 ?33 ?11 ?21 ?31 ?12 ?22 ?32 ?13 ?23 ?33 ?11 ?21 ?31 0 0 0 ?22 ?32 = ?33 ?11= ?11?11?12= ?11?12?13= ?11?13 ?21= ?21?11?22= ?21?12+ ?22?22 ?23= ?21?13+ ?22?23 ?31= ?21?11?32= ?31?12+ ?32?22 ?33= ?31?13+ ?32?23+ ?33?33 12 Unknowns and 9 equations! 3 free entries! In general, n2 equations and n2 + n unknows! n free entries!

Doolittles Algorithm: Define: ???= ???= 1 The U matrix: i j ? 1 ???= ??? ?????? ? = 1,2, ?; ? = ?,? + 1, ? ?=1 The L matrix: j < i ? ???= ?????? ?=1 ? 1?????? ??? ???=??? ?=1 ? = ? + 1, ?; ? = 1,2, ?

Crouts Algorithm: ? ???= ?????? ? = min ?,? ?=1 Define: ???= ???= 1 The L matrix: j i ? ? 1 ???= ?????? ???= ??? ?????? ?=1 ?=1 ? = 1,2, ?; ? = ?,? + 1, ? The U matrix: i < j ? ? 1?????? ??? ???=??? ?=1 ???= ?????? ?=1 ? = 1,2, ?; ? = ? + 1, ?

Doolittles Algorithm (33 example): The L matrix: j < i ?11= ?22= ?33= 1 ???=??? ?=1 ??? j = 1, i = 2, 3: ?21=?21 j = 2, i = 3: ?32=?32 ?31?12 The U matrix: i j ? 1 2 4 1 0 1 0 0 1 ?21 ?31 ?32 ? 1?????? ? = ? + 1, ?; ? = 1,2, ? 1 ?11 , ?31=?31 ?11 ?22 ?11 0 0 ?12 ?22 0 ?13 ?23 ?33 1 3 5 ???= ??? ?????? ? = 1,2, ?; ? = ?,? + 1, ? ?=1 i = 1, j = 1, 2, 3: ?11= ?11, ?12= ?12, ?13= ?13 i = 2, j = 2, 3: ?22= ?22 ?21?12, ?23= ?23 ?21?13 i = 3, j = 3: ?33= ?23 ?31?13 ?32?23

Crouts Algorithm (33 example): Define: ???= ???= 1 1 5 3 ?11 ?21 ?31 0 0 0 The L matrix: j i ?22 ?32 ?33 ? 1 ???= ??? ?????? ? = 1,2, ?; ? = ?,? + 1, ? ?=1 2 4 1 0 0 ?12 1 0 ?13 ?23 1 The U matrix: i < j ?11= ?22= ?33= 1 ???=??? ?=1 ? 1?????? ??? Verify this computation sequence! ? = 1,2, ? 1; ? = ? + 1, ?

LU Theorem: Let Ak be a sequence of matrices formed by first k rows and k columns of a n n square matrix A. If det (Ak) 0 for k= 1, 2, (n-1), then there exist an upper triangular matrix U and a lower triangular matrix L such that, A = LU. Furthermore, if the diagonal elements of either L or U are unity, i.e. lii or uii = 1 for i= 1,2, n, both L and U are unique. Ak ?11 ?21 ??1 ??1 ?12 ?22 ??2 ??2 ?1? ?2? ??? ??? ?1? ?2? ??? ??? For k = 1, the theorem is trivially valid! A = b Ak-1 ?11 ?21 ?? 1,1 ??1 ?12 ?22 ?? 1,2 ??2 ?1,? 1 ?2,? 1 ?? 1,? 1 ??? ?1? ?2? ?? 1,? ??? Let s assume that, the theorem is valid for (k - 1) and prove it for k Ak = cT

A Ak k = = L Lk kU Uk k ?? ? ?? ? ?? ? ? ? ?? ? ?? ? 1 = ??? ??? = A Ak k- -1 1: exists uniquely (by assumption). Also note L Lk k- -1 1U Uk k- -1 1 = that the following is valid, det(A Ak k- -1 1) = det(L Lk k- -1 1).det(U Uk k- -1 1) 0 L Lk k- -1 1u u = = b b : Since det(L Lk k- -1 1) 0, the triangular system has a unique solution for the vector u lTUk-1= cT or Uk-1Tl = c: Since det(U Uk k- -1 1) 0, the triangular system has a unique solution for the vector l lTu + ukk = akk : Since l and u are unique, ukk is unique. Condition for existence: det(A Ak k- -1 1) = det(L Lk k- -1 1).det(U Uk k- -1 1) 0

Diagonalization (LDU theorem): Let A be a n n invertible matrix then there exists a decomposition of the form A = LDU where, L is a n n lower triangular matrix with diagonal elements as 1, U is a n n upper triangular matrix with diagonal elements as 1, and D is a n n diagonal matrix. (Can you prove it? also done in MTH 102) Example of a 3 3 matrix: ?11 0 0 ?12 ?22 0 ?13 ?23 ?33 ?11 ?21 ?31 ?12 ?22 ?32 ?13 ?23 ?33 1 0 1 0 0 1 ?21 ?31 = ?32 ?13?11 1 ?12?11 1 0 1 ?11= ?11 0 0 0 0 0 1 0 1 0 0 1 ?23?22 ?21 ?31 ?22= ?22 0 = 0 0 ?32 ?33= ?33 For symmetric matrix: U = LT and A = LDLT Note that the entries of the diagonal matrix D are the pivots!

For positive definite matrices, pivots are positive! Therefore, a diagonal matrix D containing the pivots can be factorized as: D = D1/2D1/2 Example of a 3 3 matrix ?11 0 0 0 0 0 ?22 0 ?33 ?11 0 0 0 ?11 0 0 0 = ?22 0 0 ?22 0 0 0 ?33 0 ?33 For positive definite matrices: A = LDLT = L D1/2D1/2 LT However, D1/2LT = (LD1/2)T. Denote: L D1/2 = L1 Therefore, A = L1 L1T . This is also a LU-Decomposition where one needs to evaluate only one triangular matrix L1.

Cholesky Algorithm (for symmetric positive definite matrices): ? ???= ?????? ? = min ?,? ?=1 A = LLT .where U = LT. Elements of matrix L are to be evaluated! Diagonal elements of the L matrix: j = i = p, ljk = ukj ? ? 1 2+ 2 ???= ?????? ???= ??? ??? ?=1 ?=1 12 ? 1 2 ???= ??? ??? ? = 1,2, ? ?=1 Off-diagonal elements of the L matrix: j < i, p = j, ljk = ukj ? 1 ? 1?????? ??? ???=??? ?=1 ???= ??????+ ?????? ?=1 ? = 1,2, ?; ? = ? + 1, ?

Cholesky Algorithm: Example Problem A is a positive definite matrix. This implies:- A is a symmetric matrix, A = AT For any non-zero column vector x of n real numbers, xTAx > 0 Hence, it can be decomposed into LLT by Cholesky algorithm

Cholesky Algorithm: Example Problem Solve for coefficients of the L matrix, column by column (increasing j) Diagonal elements of the L matrix: ?? ? ? ? ???= ??? ??? ? = ?,?, ? ?=? Off-diagonal elements of the L matrix: j < i, ???=??? ?=? ? ??????? ??? ? = ? + ?, ? ? = ?,?, ?;

Cholesky Algorithm: Example Problem

Cholesky Algorithm: Example Problem

Cholesky Algorithm: Example Problem

If you already have a LU-decomposition available for matrix A: Can you use it to compute the inverse? Can you use it compute the determinant?

Banded Matrix Band Width = a + b - 1 A system of equation with Tri- diagonal coefficient matrix. Total number of elements = n2. Non-zero elements = 3n-2

Tri-Diagonal Matrix: Origin Solution of differential equations of the form: ( ) x ( ) x ( ) x ( ) x + + = p y q y r y s ( ) 0 ( ) l = = and y a y b Interpolation, such as by cubic splines: 0 i i+1 n 1 2 i-1 n-1 Node Numbers:- x0 xi xi+1 xn x1 x2 xi-1 x-values:- xn-1 y-values:- y0 yi yi+1 yn y1 y2 yi-1 yn-1

Thomas Algorithm (for Tridiagonal) No need to store n2 + n elements! Store only 4n elements in the form of four vectors l, d, u and b ith equation is: Notice: l1 = un = 0 + + = l x d x u x b + 1 1 i i i i i i i

Thomas Algorithm Initialize two new vectors and as 1 = d1and 1 = b1 Take the first two equations and eliminate x1 : u x d x l + + + = x u x 1 1 1 2 1 = x b 2 1 2 2 2 3 2 + = Resulting equation is: where, x u x 2 2 2 3 2 l l 2 2 = = d u b 2 2 1 2 2 1 1 1 Similarly, we can eliminate x2, x3

Thomas Algorithm At the ith step: + = + x u + x x d 1 1 1 1 i i i i i x = l x u b + 1 1 i i i i i i i + = Eliminate xi-1 to obtain: where, x u x +1 i i i i i l l i i = = d u b 1 1 i i i i i i 1 1 i i

Thomas Algorithm Last two equations are: + = x u x 1 1 1 1 n x n n n n + = l d x b 1 n n n n n = nx Eliminate xn-1 to obtain: n n l l n n = = b d u 1 1 n n n n n n 1 1 n n

Thomas Algorithm Given: four vectors l, d, u and b Generate: two vectors and as 1 = d1and 1 = b1 l l i i = = b d u 1 i i i 1 i i i 1 i 1 i i= 2, 3, .. n u x Solution: + 1 n i i i = = n x x i n i i = n-1 3, 2, 1 FP operations: 8(n-1) + 3(n-1) + 1 = 11n - 10

Thomas Algorithm: Example Thomas Algorithm: Example 2 2 1 0 0 0 x 1 0 2 1 1 0 0 x 2 = 0 1 1 1 1 x 3 l 0 1 0 x 4 l u x = = = = ; ; ; i i d d u b b = = + 1 ; n i i i x x 1 1 1 1 1 1 i i i i i i n i 1 1 i i n i 1 3 1 4 1 1 l = = = ) 1 = = ) 1 = = ) 1 = ; 2 2 ( ; 2 ( ; 1 ( 2 d u 1 2 2 1 3 4 l 2 2 / 3 2 3 4 / 3 4 1 1 1 1 3 = = = = = = = = 2 ; 0 0 0 ; 0 1 0 ; 1 0 1 b 1 2 2 1 3 4 2 / 3 2 4 / 3 4 1 3 ) 1 3 ) 1 2 ) 1 1 ( 0 ( 0 ( u x = = = = = = = = = 3 3 4 4 ; 3 ; 3 ; 2 1 x x x x 4 3 2 1 4 / 3 / 3 2 2 4 3