LU Decomposition Method

Jamie Trahan, Autar Kaw, Kevin Martin
University of South Florida
United States of America
kaw@eng.usf.edu

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 coefficent 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. The user can change the values that are highlighted and Mathematica will calculate the solution vector [X].

Number of equations

n = 4

[A] nxn coefficient matrix

A = 
      Table[{{175, 25, 5, 1}, {512, 64, 8, 1}, {1728, 144, 12, 1}, {8, 4, 2, 
            1}}]; A // MatrixForm

[RHS] nx1 right hand side array

RHS = Table[{106.8, 177.2, 279.2, 79.3}]; RHS // MatrixForm

Section 2: LU Decomposition method

This section divides LU Decomposition into 3 steps:
1) Decomposition of coefficient matrix (Subscript[[A], nxn])
2) Forward Substitution
3) Back Substitution

Step 1: Finding [L] and [U]

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 Naive 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.

LU Decomposition procedure:

Variable names:
n = number of equations
U = nxn upper triangular matrix
L = nxn lower triangular matrix

\(LUdecompose[n_, A_] := 
      Module[{AA, multiplier}, 
L = 
          Array[0, {n, n}]; 
U = 
          Array[0, {n, n}]; 
AA = 
          Array[0, {n, n}]; 
For[i = 1, 
          i ≤ n, \(i++\), 
          L[\([i, i]\)] = 1; 
For[j = 1, j <= n, \(j++\), 
            AA[\([i, j]\)] = A[\([i, j]\)]; 
            If[j != i, L[\([i, j]\)] = 0]]]; 
For[k = 1, 
          k ≤ 
            n - 1, \(k++\), 
\(For[i = k + 1, 
                i ≤ n, \(i++\), 
multiplier = 
                  N[AA[\([i, k]\)]\/AA[\([k, k]\)]]; 
L[\([i, k]\)] = multiplier; 
For[j = 1, 
                  j ≤ n, \(j++\), 
                  AA[\([i, j]\)] = 
                    AA[\([i, j]\)] - 
                      multiplier*AA[\([k, j]\)]]]
For[
                i = 1, i ≤ n, \(i++\), 
For[
                  j = 1, j ≤ n, \(j++\), 
                  U[\([i, j]\)] = AA[\([i, j]\)]; 
                  If[j < i, 
                    U[\([i, j]\)] = 
                      0]]];\)]; 
Print["L", "=", 
          L // MatrixForm]; 
Print["U", "=", 
          U // MatrixForm]]\)

LUdecompose[n, A]

Notice that matrix [L] has only one unknown to be solved for in its first row, and matrix [U] has only one unknown to be solved for in its last row. This will prove useful in solving for the solution vector [X] in the following steps of LU decomposition method.

Step 2: Forward Substitution

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

(Subscript[z, 1]) = ((Subscript[c, 1])/(Subscript[l, 1, 1]))                                                           Equation (2.1)

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

\(TraditionalForm\`z\_i = \((c\_i - \([∑\_\(j = 1\)\%\(i - 1\)\((l\_\(i, j\)*z\_j)\)]\)\_\(i = 2.  . n\))\)\/l\_\(i, i\)\)                               Equation (2.2)

The procedure below conducts forward substitution steps to solve for (Subscript[[Z], nx1]).

Variable names:
n = number of equations
RHS = nx1 right hand side array
L = nxn lower triangular matrix

\(forwardsubstitution[n_, L_, C_] := 
      Module[{}, 
Z = 
          Array[z, n]; 
Z[\([1]\)] = 
          C[\([1]\)]\/L[\([1, 1]\)]; 
For[i = 2, 
          i ≤ n, \(i++\), 
          summ = 0; 
For[j = 1, 
            j ≤ i - 1, \(j++\), 
            summ = summ + 
                L[\([i, j]\)]*
                  Z[\([j]\)]]; 
Z[\([i]\)] = \((C[\([i]\)] - summ)\)\/L[\([i, i]\)]]; 
Print["Z", "=", 
          Z // MatrixForm]]\)

forwardsubstitution[n, L, RHS]

Step 3: Back Substitution

Now that [Z] has been found, it can be used in the back substitution step, [U][X] = [Z], to solve for solution vector (Subscript[[X], nx1]), where (Subscript[[U], nxn])is the upper triangular matrix calculated in Step 1, and (Subscript[[Z], nx1])is the right hand side array.

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

(Subscript[x, n]) = ((Subscript[z, n])/(Subscript[u, n, n]))                                                          Equation (2.3)

The remaining unknowns are solved for working backwards from the (((n - 1))^th)equation to the first equation using the formula below

\(x\_i = \((z\_i - \([∑\_\(j = i + 1\)\%n\((u\_\(i, j\)*x\_j)\)]\)\_\(i = n - 1.   .1\))\)\/u\_\(i, i\)\)                        Equation (2.4)

The following procedure solves for (Subscript[[X], nx1]).

Variable names:
n = number of equations
Z = nx1 right hand side array
U = nxn upper triangular matrix

\(backsubstitution[n_, U_, Z_] := 
      Module[{}, 
X = 
          Array[x, n]; 
X[\([n]\)] = 
          Z[\([n]\)]\/U[\([n, n]\)]; 
For[i = n - 1, 
          i >= 1, \(i--\), 
summ = 
            0; 
For[j = i + 1, j ≤ n, \(j++\), 
            summ = summ + 
                U[\([i, j]\)]*
                  X[\([j]\)]]; 
X[\([i]\)] = \((Z[\([i]\)] - summ)\)\/U[\([i, i]\)]]; 
Print["X", "=", 
          X // MatrixForm]]\)

The solution vector (Subscript[[X], nx1]) is:

backsubstitution[n, U, Z]

Conclusion

Mathematica helped us apply our knowledge of LU Decomposition method to solve a system of n simultaneous linear equations. The results are as follows:

Print["L", "=", L // MatrixForm]
Print["U", "=", U // MatrixForm]
forwardsubstitution[n, L, RHS]
backsubstitution[n, U, Z]

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

(({{5, 6, 2.3, 6}, {9, 2, 3.5, 7}, {3.5, 6, 2, 3}, {1.5, 2, 1.5, 6}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}}))= (({{4}, {5}, {6.7}, {7.8}}))

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

(({{5, 6, 2.3, 6}, {9, 2, 3.5, 7}, {3.5, 6, 2, 3}, {1.5, 2, 1.5, 6}}))

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

References

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