LU Decomposition Method 

? 2006 Kevin Martin, Autar Kaw, Jamie Trahan 

University of South Florida  

United States of America 

kaw@eng.usf.edu 

NOTE: This worksheet demonstrates the use of Maple to illustrate LU Decomposition method, a technique used in solving a system of simultaneous linear equations. 

Introduction 

When solving multiple sets of simultaneous linear equations with the same coefficient matrix but different right hand sides, LU Decomposition is advantageous over other numerical methods in that it proves to be numerically more efficient in computational time than other techniques. In this worksheet, the reader can choose a system of equations and see how each step of LU decomposition method is conducted. To learn more about LU Decomposition method as well as the efficiency of its computational time click here. 

 

LU Decomposition method is used to solve a set of simultaneous linear equations, [A] [X] = [C], where [A]_nxn is a non-singular square coefficient matrix, [X]_nx1 is the solution vector, and [C]_nx1 is the right hand side array. When conducting LU decomposition method, one must first decompose the coefficient matrix [A]_nxn into a lower triangular matrix [L]_nxn, and upper triangular matrix [U]_nxn. These two matrices can then be used to solve for the solution vector [X]_nx1 in the following sequence: 

 

Recall that  

[A] [X] = [C]. 

Knowing that  

[A] = [L] [U] 

then first solving with forward substitution  

[L] [Z] = [C] 

and then solving with back substitution 

[U] [X] = [Z] 

gives the solution vector [X]. 

 

 

A simulation of  LU Decomposition method follows. 

Section 1: Input Data 

Below are the input parameters to begin the simulation. This is the only section that requires user input. Once the values are entered, Maple will calculate the solution vector [X]. 

Input Parameters: 

n = number of equations 

[A] = nxn coefficient matrix 

[RHS] = nx1 right hand side array 

 

 

>	restart; n:=4; A:=Matrix([[12.,7.,3.,6.7],[1.,5.,1.,9.],[13.,12.,4.001,8],[5.6,3.,7.,1.003]]); RHS:=[22,7.,29.001,5.301];

 

 

 

(2.1)

 

Section 2: LU Decomposition Method 

 

The following sections divide LU Decomposition method into 3 steps:  

1.) Finding the LU decomposition of the coefficient matrix [A]_nxn  

2.) Forward substitution  

3.) Back substitution 

 

 

 

2.1 LU Decomposition 

How does one decompose a non-singular matrix [A], that is how do you find [L] and [U]? 

The following procedure decomposes the coefficient matrix [A] into a lower triangular matrix [L] and upper triangular matrix [U], given [A] = [L][U]. 

• For [U], the elements of the matrix are exactly the same as the coefficient matrix one obtains at the end of forward elimination steps in Na?ve Gauss Elimination.

 

• For [L], the matrix has 1 in its diagonal entries. The non-zero elements on the non-diagonal elements are multipliers that made the corresponding entries zero in the upper triangular matrix during forward elimination.

 

Parameter names: 

n = number of equations 

A = nxn coefficient matrix 

 

>

LUdecompose:=proc(n,A) local k,i,multiplier,j,sum,AA,L,U: #Defining the [L] and [U] matrices. L:=Matrix(1..n,1..n); U:=Matrix(1..n,1..n); #Initializing diagonal of [L] to be unity. for i from 1 by 1 to n do   L[i,i]:=1.0; end do: #Dumping [A] matrix into a local matrix [AA] because input variables should not be placed on LHS of an assignment in a procedure. for i from 1 by 1 to n do   for j from 1 by 1 to n do      AA[i,j]:=A[i,j];   end do: end do: #Conducting forward elimination steps of Na?ve Gaussian Elimination to obtain [L] and [U] matrices. for k from 1 by 1 to n-1 do   for i from (k+1) by 1 to n do   #Computing multiplier values.   multiplier:=AA[i,k]/AA[k,k]:   #Putting multiplier in proper row and column of [L] matrix.      L[i,k]:=multiplier:      for j from (k+1) by 1 to n do         #Eliminating (i-1) unknowns from the ith row to generate an upper triangular matrix.         AA[i,j]:=AA[i,j]-multiplier*AA[k,j]:      end do:   end do: end do: #Dumping the end of forward elimination coefficient matrix into the [U] matrix. for i from 1 by 1 to n do   for j from i by 1 to n do      U[i,j]:=AA[i,j]:   end do: end do: return(L,U): end proc:

 

2.2 Forward Substitution 

Now that the [L] and [U] matrices have been formed, the forward substitution step, [L] [Z] = [C], can be conducted, beginning with the first equation as it has only one unknown, 

 

 

(3.3.1)

 

Subsequent steps of forward substitution can be represented by the following formula: 

 

 

 

(3.3.2)

 

The following procedure conducts forward substitution steps to solve for [Z]. 

 

Parameter names: 

n = number of equations 

[L]= nxn lower triangular matrix 

[C] = nx1 right hand side array 

 

>

forward_substitution:=proc(n,L,C) local Z,i,sum,j: #Defining the [Z] vector. Z:=Array(1..n): #Solving for the first equation as it has only one unknown. Z[1]:=C[1]/L[1,1]; #Solving for the remaining (n-1) unknowns using Equation (3.3.2). for i from 2 by 1 to n do sum:=0;   #Generating the summation term in Equation (3.3.2).   for j from 1 by 1 to i-1 do      sum:=sum+L[i,j]*Z[j]:   end do: Z[i]:=(C[i]-sum)/L[i,i]; end do: return(Z): end proc:

 

 

 

 

2.3 Back Substitution 

Now that [Z] has been calculated, it can be used in the back substitution step, [U] [X] = [Z], to solve for solution vector [X]_nx1, where [U]_nxn is the upper triangular matrix calculated in Step 2.1, and [Z]_nx1 is the right hand side array. 

Back substitution begins with the n^th equation as it has only one unknown. 

 

 

(3.4.1)

 

 

The remaining unknowns are solved for using the following formula: 

 

 

(3.4.2)

 

 

 

The following procedure solves for [X]. 

Parameter names: 

n = number of equations 

[U] = nxn upper triangular matrix 

[Z] = nx1 solution vector from forward substitution step 

 

 

>

back_substitution:=proc(n,U,Z) local i,sum,j,X: #Defining the [X] vector. X:=Array(1..n); #Solving for the nth equation as it has only one unknown. X[n]:=Z[n]/U[n,n]; #Solving for the remaining (n-1) unknowns working backwards from the (n-1)th equation to the first equation. for i from n-1 by -1 to 1 do   #Initializing series sum to zero.   sum:=0;      #Calculating summation term in Equation (3.4.2).      for j from i+1 by 1 to n do         sum:=sum+U[i,j]*X[j];      end do:   #Using Equation (3.4.2) to calculate solution vector [X].   X[i]:=(Z[i]-sum)/U[i,i]; end do: return(X); end proc:

 

Section 3: Results 

Below, Maple returns the lower triangular matrix, [L]nxn , the upper triangular matrix, [U]nxn , the intermediate solution vector [Z]nx1, and the final solution vector [X]nx1. 

 

>	LU:=LUdecompose(n,A): L:=LU[1]; U:=LU[2]; Z:=forward_substitution(n,L,RHS); X:=back_substitution(n,U,Z);

 

 

 

 

(3.4.1.1)

 

 

 

 

References 

[1] Autar Kaw, Holistic Numerical Methods Institute, http://numericalmethods.eng.usf.edu/mws, See 

Introduction to Systems of Equations 

How does LU Decomposition method work? 

Saving of computational time for finding inverse of a matrix using LU decomposition 

 

 

Conclusion 

Maple helped us apply our knowledge of LU Decomposition method to solve a system of n simultaneous linear equations. 

 

Question 1: Solve the following set of simultaneous linear equations using LU decomposition method. 

 

(3.4.3.1)

 

 

Question 2: Use LU decomposition repeatedly to find the inverse of   

 

(3.4.3.2)

 

 

Question 3: Look at the [Z] matrix in LU decomposition step [L] [Z] = [C] of Question #1. Is it the same as the [RHS] matrix at the end of forward elimination steps in Naive Gauss Elimination method? If no, is this a coincidence?  

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.