# LU Factorization with MATLAB | MATLAB Tutorial

**by** Gunjan Gupta

**LU Factorization with MATLAB**

**LU Factorization with MATLAB**

**LU Factorization**

[A]{x}={b} 1
**LU Factorization**

Equation 1 can be rearranged to give

[A]{x}-{b}=0 2Suppose that above eq. 2 could be expressed as an upper triangular system. For example, a 3×3 system:

Eq. 3 can also be expressed in matrix notation and rearranged to give

[U]{x}-{d}=0 4Now assume that there is a lower diagonal matrix with 1’s on the diagonal.

that has the property that when Eq. 4 is premultiplied by it, Eq. 2 is the result. That is,

[L]{[U]{x}-{d}}=[A]{x}-{b} 6If this equation holds, it follows from the rules for matrix multiplication that

[L][U]=[A] 7and

[L]{d}={b} 8A two step strategy(Fig.1) for obtaining solutions can be based on Eqs. 3, 7, and 8:

factorization step. [*LU*] is factored or “decomposed” into lower [*A*] and upper [*L*] triangular matrices.*U*] and [*L*] are used to determine a solution {*U*} for a right-hand side {*x*}. This step itself consists of two steps. First, Eq. 8 is used to generate an intermediate vector {*b*} by forward substitution. Then, the result is substituted into Eq. 3 which can be solved by back substitution for {*d*}.*x*

**Figure 1: The steps in LU factorization**

**Gauss Elimination as LU Factorization**

**Gauss Elimination as LU Factorization**

Gauss elimination can be used to decompose [* A*] into [

*] and [*

*L**]. This can be easily seen for [*

*U**], which is a direct product of the forward elimination. Recall that the forward-elimination step is intended to reduce the original coefficient matrix [*

*U**] to the form*

*A*which is in the desired upper triangular format, the matrix [L] is also produced during the step. This can be readily illustrated for a three-equation system,

The first step in Gauss elimination is to multiply row 1 by the factor

and subtract the result from the second row to eliminate a21 . Similarly, row 1 is multiplied by

and the result subtracted from the third row to eliminate a31 . The final step is to multiply the modified second row by

and subtract the result from the third row to eliminate .

After elimination, the [A] matrix can therefore be written as

This matrix, in fact, represents an efficient storage of the LU factorization of [A] [A] -> [L][U]

where

and

**Problem 1: ** Derive LU factorization based on the Gauss elimination,

Coefficient matrix:

After forward elimination, the following upper triangular matrix was obtained:

The elements *a*_{21} and *a*_{31} were eliminated by using the factors

and the element *a*_{32} was eliminated by using the factor

Thus, the lower triangular matrix is

Consequently, the LU factorization is

This result can be verified by performing the multiplication of [L][U] to give

After the matrix is decomposed, a solution can be generated for a particular right-hand side vector {* b*}. This is done in two steps. First, a forward-substitution step is executed by solving Eq. 8 for {

*}. At the end of this step, the right-hand side will be in the same state that it would have been had we performed forward manipulation on [*

*d**] and {*

*A**} simultaneously.*

*b*The forward-substitution step can be represented concisely as

The second step then merely amounts to implementing back substitution to solve Eq. 3.

**Problem 2: **Complete the problem 1 by generating the final solution with forward and back substitution.

the forward – elimination phase of conventional Gauss elimination resulted in

The forward-substitution phase is implemented by applying Eq. 8:

We can solve the first equation for *d*_{1} = 7.85, which can be substituted into the second equation to solve for

d2=-19.3-0.03333(7.85)=-19.5617

Both *d*_{1} and *d*_{2} can be substituted into the third equation to give

d3=71.4-0.1(7.85)+0.02713(-19.5617)=70.0843

Thus,

This can be substituted into Eq. 3, [U]{x} = {d}:

**Matlab m-file:**

function x = LUGauss(A,b) % input: % A = coefficient matrix % b = right hand side vector % output: % x = solution vector [m,n]=size(A); if m~=n error('Matrix A must be square'); end [L,U] = lu(A); % LU factorization disp('L = '); disp(L) disp('U = '); disp(U) d = L\b; x = U\d; % solution disp('x = '); disp(x)

The M-file can be used to solve the problem 1 and 2.

>> A = [3 -.1 -.2;.1 7 -.3;.3 -.2 10];

b = [7.85; -19.3; 71.4];

>> LUGauss(A,b);

L =

1.0000 0 0

0.0333 1.0000 0

0.1000 -0.0271 1.0000

U =

3.0000 -0.1000 -0.2000

0 7.0033 -0.2933

0 0 10.0120

x =

3.0000

-2.5000

7.0000

**LU Factorization with Pivoting**

Just as for standard Gauss elimination, partial pivoting is necessary to obtain reliable solutions with * LU *factorization. One way to do this involves using a permutation matrix. The approach consists of the following steps:

factorization with pivoting of a matrix [*LU*] can be represented in matrix form as*A*

*][*

*P**] = [*

*A**][*

*L**]*

*U*The upper triangular matrix, [* U*], is generated by elimination with partial pivoting, while storing the multiplier factors in [

*] and employing the permutation matrix, [*

*L**], to keep track of the row switches.*

*P*] and [*L*] are used to perform the elimination step with pivoting on {*P*} in order to generate the intermediate right-hand-side vector, {*b*}. This step can be represented concisely as the solution of the following matrix formulation:*d*

*]{*

*L**} = [*

*d**]{*

*P**}*

*b**]{*

*U**} = {*

*x**}*

*d***Problem 3: **Compute the LU factorization and find the solution for the system,

Before elimination, we set up the initial permutation matrix:

We immediately see that pivoting is necessary, so prior to elimination we switch the rows:

At the same time, we keep track of the pivot by switching the rows of the permutation matrix:

We then eliminate *a*_{21} by subtracting the factor *l*_{21} = *a*_{21}/*a*_{11 }= 0.0003/1 = 0.0003 from the second row of * A*. In so doing, we compute that the new value of

*′*

*a*_{22}= 3 – 0.0003(1) = 2.9997. Thus, the elimination step is complete with the result:

Before implementing forward substitution, the permutation matrix is used to reorder the right-hand-side vector to reflect the pivots as in

Then, forward substitution is applied as in

which can be solved for *d*_{1} = 1 and *d*_{2} = 2.0001 − 0.0003(1) = 1.9998. At this point, the system is

Applying back substitution gives the final result:

The * LU *factorization algorithm requires the same total flops as for Gauss elimination. The only difference is that a little less effort is expended in the factorization phase since the operations are not applied to the right-hand side. Conversely, the substitution phase takes a little more effort.

**Matlab m-file:**

function x = LUPivot(A,b) % input: % A = coefficient matrix % b = right hand side vector % output: % x = solution vector [m,n]=size(A); if m~=n error('Matrix A must be square'); end L = eye(n); P = L; U = A; for k=1:n [pivot m]=max(abs(U(k:n,k))); m = m+k-1; end end

The M-file can be used to solve the problem 3.

>> A = [0.0003 3.000;1.0000 1.0000];

>> b = [2.0001 ; 1.0000];

>> LUPivot(A,b);

L =

1.0000 0

0.0003 1.0000

U =

1.0000 1.0000

0 2.9997

x =

0.3333

0.6667