HOME
PDF (letter size)
PDF (legal size)

Implementation of LU Decomposition and Linear Solver using Matlab

Nasser M. Abbasi

sometime in 2009   Compiled on May 24, 2020 at 5:01am

Contents

1 introduction
2 examples
 2.1 using nma_LU()
  2.1.1 example 1
  2.1.2 example 2
 2.2 using nma_LinearSolve()
  2.2.1 example 1
  2.2.2 example 2
 2.3 using nma_ForwardSub()
  2.3.1 example 1
 2.4 using nma_BackSub()
  2.4.1 example 1
3 code listing
 3.1 test script
 3.2 nma_LU.m
 3.3 nma_LinearSolve.m
 3.4 nma_ForwardSub.m
 3.5 nma_BackSub.m

1 introduction

This is MATLAB implementation for LU decomposition, forward substitution, backward substitution, and linear system solver.

The functions written are:

  1. nma_LU.m.txt LU decomposition with partial pivoting with threshold support.
  2. nma_ForwardSub.m.txt solves \(L y = b\) for \(y\)
  3. nma_BackSub.m.txt solves \(U x = y\) for \(x\)
  4. nma_LinearSolve.m.txt driver to solve \(A x = b\) for \(x\) using calling sequence \(1\rightarrow 2\rightarrow 3\)

Partial pivoting (P matrix) was added to the LU decomposition function. In addition, the LU function accepts an additional argument which allows the user more control on row exchange.

Matlab lu() function does row exchange once it encounters a pivot larger than the current pivot. This is a good thing to always try to do. But sometimes if the difference between the pivots is small, a user might not want this feature. Hence I added a ”threshold” second parameter to the nma_LU.m function to indicate how large a difference should exist for a row exchange to occur.

A row exchange will always occur if the current pivot is zero and a non-zero pivot exist to do the exchange.

To get the same exact behavior as Matlab lu() simply make this parameter zero.

Below are examples calling the nma_LU, nma_ForwardSub.m, nma_BackSub.m and nma_LinearSolve.m.

In each example below, the output is verified against Matlab own functions

2 examples

2.1 using nma_LU()

2.1.1 example 1

>> A=[1 1 2;2 -1 1;1 2 0]
A =
     1     1     2
     2    -1     1
     1     2     0

>> [L,U,P]=nma_LU(A,0)
L =
   1.00000000000000                  0                  0
   0.50000000000000   1.00000000000000                  0
   0.50000000000000   0.60000000000000   1.00000000000000

U =
   2.00000000000000  -1.00000000000000   1.00000000000000
                  0   2.50000000000000  -0.50000000000000
                  0                  0   1.80000000000000
P =
     0     1     0
     0     0     1
     1     0     0

>> [L,U,P]=lu(A)

L =

   1.00000000000000                  0                  0
   0.50000000000000   1.00000000000000                  0
   0.50000000000000   0.60000000000000   1.00000000000000

U =
   2.00000000000000  -1.00000000000000   1.00000000000000
                  0   2.50000000000000  -0.50000000000000
                  0                  0   1.80000000000000

P =
     0     1     0
     0     0     1
     1     0     0
>>

2.1.2 example 2

>> A=rand(4);
>> [L,U,P]=nma_LU(A,0)

L =
   1.00000000000000                  0                  0                  0
   0.01212703756687   1.00000000000000                  0                  0
   0.07119243718995   0.20742768803520   1.00000000000000                  0
   0.43394327408595   0.19377225100868   0.40879105345917   1.00000000000000

U =
   0.81316649730376   0.19872174266149   0.01527392702904   0.46599434167542
                  0   0.60138257315521   0.74660044907755   0.41299833684006
                  0                  0   0.11623493184110   0.32625386921423
                  0                  0                  0   0.51620218784594

P =
     0     0     1     0
     0     0     0     1
     1     0     0     0
     0     1     0     0

>> [L,U,P]=lu(A)
L =
   1.00000000000000                  0                  0                  0
   0.01212703756687   1.00000000000000                  0                  0
   0.07119243718995   0.20742768803520   1.00000000000000                  0
   0.43394327408595   0.19377225100868   0.40879105345917   1.00000000000000

U =
   0.81316649730376   0.19872174266149   0.01527392702904   0.46599434167542
                  0   0.60138257315521   0.74660044907755   0.41299833684006
                  0                  0   0.11623493184110   0.32625386921423
                  0                  0                  0   0.51620218784594

P =
     0     0     1     0
     0     0     0     1
     1     0     0     0
     0     1     0     0
>>

2.2 using nma_LinearSolve()

2.2.1 example 1

>> A=[1 1 2;2 -1 1;1 2 0]
A =
     1     1     2
     2    -1     1
     1     2     0

>> b=[1 2 1];

>> nma_LinearSolve(A,b)

ans =

     1
     0
     0

>> A\b(:)

ans =

     1
     0
     0
>>

2.2.2 example 2

>> A=rand(6);
>> b=rand(6,1);
>> nma_LinearSolve(A,b)

ans =

   0.59090034220622
  -0.56523444269280
   0.95687095978224
  -0.97248777153372
   1.00007995741472
   0.24035777097022

>> A\b(:)

ans =

   0.59090034220622
  -0.56523444269280
   0.95687095978223
  -0.97248777153372
   1.00007995741472
   0.24035777097022

>>

2.3 using nma_ForwardSub()

2.3.1 example 1

>> [L,U,P]=nma_LU(A,0);
>> nma_ForwardSub(L,b)

ans =

   0.83849604493808
   0.36727512318587
   0.12405626870025
  -0.14539724685973
   0.17813906538571
  -0.19809655526705

2.4 using nma_BackSub()

2.4.1 example 1

>> nma_BackSub(U,ans)

ans =

   0.29867870305809
  -0.84855613142087
   0.48347828223154
  -1.68311779577975
   1.49928530116874
   1.53825192677360

3 code listing

3.1 test script

3.2 nma_LU.m

3.3 nma_LinearSolve.m

3.4 nma_ForwardSub.m

3.5 nma_BackSub.m