Transcript LU Decomposition

LU Decomposition
Civil Engineering Majors
Authors: Autar Kaw
http://numericalmethods.eng.usf.edu
Transforming Numerical Methods Education for STEM
Undergraduates
7/17/2015
http://numericalmethods.eng.usf.edu
1
LU Decomposition
http://numericalmethods.eng.usf.edu
LU Decomposition
LU Decomposition is another method to solve a set of
simultaneous linear equations
Which is better, Gauss Elimination or LU Decomposition?
To answer this, a closer look at LU decomposition is
needed.
http://numericalmethods.eng.usf.edu
LU Decomposition
Method
For most non-singular matrix [A] that one could conduct Naïve Gauss
Elimination forward elimination steps, one can always write it as
[A] = [L][U]
where
[L] = lower triangular matrix
[U] = upper triangular matrix
http://numericalmethods.eng.usf.edu
How does LU Decomposition work?
If solving a set of linear equations
If [A] = [L][U] then
Multiply by
Which gives
Remember [L]-1[L] = [I] which leads to
Now, if [I][U] = [U] then
Now, let
Which ends with
and
[A][X] = [C]
[L][U][X] = [C]
[L]-1
[L]-1[L][U][X] = [L]-1[C]
[I][U][X] = [L]-1[C]
[U][X] = [L]-1[C]
[L]-1[C]=[Z]
[L][Z] = [C] (1)
[U][X] = [Z] (2)
http://numericalmethods.eng.usf.edu
LU Decomposition
How can this be used?
Given [A][X] = [C]
1. Decompose [A] into [L] and [U]
2. Solve [L][Z] = [C] for [Z]
3. Solve [U][X] = [Z] for [X]
http://numericalmethods.eng.usf.edu
When is LU Decomposition better
than Gaussian Elimination?
To solve [A][X] = [B]
Table. Time taken by methods
Gaussian Elimination
 8n 3
4n 
T 
 12 n 2  
3 
 3
LU Decomposition
 8n 3
4n 
T 
 12 n 2  
3 
 3
where T = clock cycle time and n = size of the matrix
So both methods are equally efficient.
http://numericalmethods.eng.usf.edu
To find inverse of [A]
Time taken by Gaussian Elimination
Time taken by LU Decomposition
 nCT |FE CT |BS 
 CT |LU  n  CT |FS  n  CT |BS
 8n 4
4n 2 
3

 T 
 12n 
3
3


 32n 3
20n 

 T 
 12n 2 
3 
 3
Table 1 Comparing computational times of finding inverse of a matrix using
LU decomposition and Gaussian elimination.
n
10
100
1000
10000
CT|inverse GE / CT|inverse LU
3.28
25.83
250.8
2501
http://numericalmethods.eng.usf.edu
Method: [A] Decompose to [L] and [U]
1
A  LU    21
 31
0
1
 32
0 u11
0  0
1  0
u12
u 22
0
u13 
u 23 
u 33 
[U] is the same as the coefficient matrix at the end of the forward
elimination step.
[L] is obtained using the multipliers that were used in the forward
elimination process
http://numericalmethods.eng.usf.edu
Finding the [U] matrix
Using the Forward Elimination Procedure of Gauss Elimination
 25 5 1
 64 8 1


144 12 1
Step 1:
5
1 
 25
64
 2.56; Row2  Row12.56   0  4.8  1.56
25
144 12
1 
5
1 
25
144
 5.76; Row3  Row15.76   0  4.8  1.56
25
 0  16.8  4.76
http://numericalmethods.eng.usf.edu
Finding the [U] Matrix
5
1 
25
Matrix after Step 1:  0  4.8  1.56


 0  16.8  4.76
1 
25 5
 16.8
 3.5; Row3  Row23.5   0  4.8  1.56
Step 2:
 4.8
 0
0
0.7 
1 
25 5
U    0  4.8  1.56
 0
0
0.7 
http://numericalmethods.eng.usf.edu
Finding the [L] matrix
1

 21
 31
0
1
 32
0
0
1
Using the multipliers used during the Forward Elimination Procedure
a
64
 21  21 
 2.56
5 1
From the first step  25
a11 25
of forward
elimination
 64 8 1


144 12 1
 31 
a31 144

 5.76
a11
25
http://numericalmethods.eng.usf.edu
Finding the [L] Matrix
From the second
step of forward
elimination
5
1 
25
 0  4.8  1.56


 0  16.8  4.76
 32 
a32  16.8

 3.5
a22
 4.8
0 0
 1
L  2.56 1 0
5.76 3.5 1
http://numericalmethods.eng.usf.edu
Does [L][U] = [A]?
0 0 25 5
1 
 1
LU   2.56 1 0  0  4.8  1.56 
5.76 3.5 1  0
0
0.7 
?
http://numericalmethods.eng.usf.edu
Example: Cylinder Stresses
To find the maximum stresses in a compound cylinder, the following
four simultaneous linear equations need to solved.
4.2857107

7
4
.
2857

10


 6.5

0

 9.2307105
0
 5.4619105
 0.15384
 4.2857107
6.5
0
4.2857107
  c1   7.887103 
  

5.4619105  c2  
0


0.15384  c3   0.007 

 
5 
 3.605710  c4  
0

0
Example: Cylinder Stresses
In the compound cylinder, the inner cylinder has an internal radius of
a = 5”, and outer radius c = 6.5”, while the outer cylinder has an
internal radius of c = 6.5” and outer radius, b=8”. Given E = 30×106
psi, ν = 0.3, and that the hoop stress in outer cylinder is given by
E
 
1  2

 1 
c3 1    c4  r 2





find the stress on the inside radius of the outer cylinder.
Find the values of c1, c2, c3 and c4 using LU decomposition.
Example: Cylinder Stresses
Use Forward Elimination Procedure of Gauss Elimination to find [U]
4.2857 107

7
4
.
2857

10


 6.5

0

Step 1
 9.2307 105
0
 5.4619 105
 4.2857 107
 0.15384
6.5
0
4.2857 107


5.4619 105 
0.15384 
5
 3.6057 10 
0
4.2857107
 1; Row2  Row11 
7
4.285710
4.2857 107

0


 6.5

0

 9.2307 105
0
3.7688 105
 0.15384
 4.2857 107
6.5
0
4.2857 107


5.4619 105 
0.15384 

 3.6057 105 
0
Example: Cylinder Stresses
Step 1 cont.
 6.5
7
7


1
.
5167

10
;
Row
3

Row
1

1
.
5167

10

7
4.2857 10

4.2857107

0


0

0

 9.2307105
0
3.7688105
 0.29384
 4.2857107
6.5
0
4.2857107



5.4619105 
0.15384 
5
 3.605710 
0
Example: Cylinder Stresses
Step 1 cont.
0
 0; Row 4  Row10  
7
4.2857 10
4.2857107

0


0

0

 9.2307105
0
3.7688105
 0.29384
 4.2857107
6.5
0
4.2857107
This is the matrix after Step 1.


5.4619105 
0.15384 
5
 3.605710 
0
Example: Cylinder Stresses
Step 2
 0.29384
7
7


7
.
7966

10
;
Row
3

Row
2

7
.
7966

10

5
3.7688 10

4.2857107

0


0

0

 9.2307105
0
3.7688105
 4.2857107
0
 26.914
0
4.2857107



5.4619 105 
0.57968 
5
 3.605710 
0
Example: Cylinder Stresses
Step 2 cont.
0
 0; Row 4  Row 20 
5
3.7688 10
4.2857107

0


0

0

 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
4.2857107
This is the matrix after Step 2.


5.4619105 
0.57968 
5
 3.605710 
0
Example: Cylinder Stresses
Step 3
4.2857107
 1.5294106 ; Row4  Row3  1.5923106 
 26.914

4.2857107

0


0

0

 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
0
This is the matrix after Step 3.



5.4619105 
0.57968 
5
5.625010 
0
Example: Cylinder Stresses
4.2857107

0
[U ]  

0

0

 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
0


5.4619105 
0.57968 
5
5.625010 
0
Example: Cylinder Stresses
Use the multipliers from Forward Elimination to find [L]
4.2857 107

7
4
.
2857

10


 6.5

0

From the 1st step of
forward elimination
1

 21
 31

 41
0
1
 32
 42
0
0
1
 43
0
0
0

1
 9.2307 105
0
 5.4619 105
 0.15384
 4.2857 107
6.5
0
4.2857 107


5.4619 105 
0.15384 

 3.6057 105 
0
a21 4.2857107
 21 

1
7
a11 4.285710
a31
 6.5
7
 31 



1
.
5167

10
a11 4.2857107
a41
0
 41 

0
7
a11 4.285710
Example: Cylinder Stresses
4.2857107

0
From the 2nd step of 
forward elimination 
0

0

1

 21
 31

 41
0
1
 32
 42
0
0
1
 43
0
0
0

1
 9.2307105
0
3.7688105
 0.29384
 4.2857107
6.5
0
4.2857107


5.4619105 
0.15384 
5
 3.605710 
0
 32
a32  0.29384
7



7
.
7966

10
a22 3.7688105
 42
a42
0


0
5
a22 3.768810
Example: Cylinder Stresses
4.2857107

0
From the 3rd step of 
forward elimination 
0

0

1

 21
 31

 41
0
1
 32
 42
0
0
1
 43
0
0
0

1
 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
4.2857107


5.4619105 
0.57968 
5
 3.605710 
0
a43 4.2857107
 43 

 1.5294106
a33
 26.914
Example: Cylinder Stresses
1


1
L  
 1.5167107

0

0
1
 7.7966107
0
0
0
1
 1.5924106
0
0
0

1
Example: Cylinder Stresses
Does [L][U] = [A]?
LU  
1
0
0


1
1
0

 1.5167107  7.7966107
1

0
0
 1.5924106

0 4.2857107  9.2307105

0
0


0 
0
3.7688105  4.2857107 5.4619105 
?

0 
0
0
 26.914
0.57968


1 
0
0
0
5.6250105 
Example: Cylinder Stresses
Set [L][Z] = [C]
1
0
0


1
1
0

 1.5167107  7.7966107
1

0
0
 1.5924106

Solve for [Z]
0  z1   7.887103 


0  z 2  
0


0  z3   0.007 

  
1  z 4  
0

z1  7.887103

z1  z 2  0

 1.51667107 z1   7.79662107 z 2  z 3  0.007
 1.5924106 z 3  z 4  0
Example: Cylinder Stresses
Solve for [Z]
z1  7.887103
z2   z1

   7.887103

 7.887103

 

 0.007   1.5167 10   7.887 10    7.7966 10  7.887 10 
z3  0.007   1.5167 107 z1   7.7966107 z 2
7
 1.195310 2
3
7
3
Example: Cylinder Stresses
Solving for [Z] cont.


  1.592410  1.195310 
z4    1.5924106 z3
6
 19034
 z1    7.887103 
z  
3 
7.88710 
2


Z 

 z 3  1.195310 2 

  
z 4   19034 
2
Example: Cylinder Stresses
Set [U][C] = [Z]
4.2857107  9.2307105

0
0

5
7
5
0
3
.
7688

10

4
.
2857

10
5
.
4619

10



0
0
 26.914
0.57968 

5
0
0
0
5
.
6250

10


 c1   7.887103 
c  
3 
7
.
887

10
2

 

2
c3  1.195310 

  
c
19034

 4  
The four equations become:


4.2857107 c1   9.2307105 c2  0c3  0c4  7.887103


3.7668105 c2   4.2857107 c3  5.4619105 c4  7.887103
 26.914c3  0.57968c4  1.1953102
5.6250105 c4  19034
Example: Cylinder Stresses
Solve for [C]
19034
5.6250 105
 3.3837 10 2
c4 
1.1953102  0.57968c4
c3 
 26.914
1.195310 2  0.57968c4

 26.914
1.195310 2  0.579683.383710 2

 26.914
 2.846910 4


Example: Cylinder Stresses
Solve for [C] cont.



 
7.887103   4.2857107 c3  5.4619105 c4
c2 
3.7688105
7.887103   4.2857107  2.846910 4  5.4619105  3.383710 2

3.7688105
 4.2615103




 
 7.887103   9.2307105 c2  0c3-0c4
c1 
4.2857107
 7.887103   9.2307105  4.2615103

4.2857107
 9.2244105



Example: Cylinder Stresses
Solution:
The solution vector is
 c1   9.2244105 
c  
3 
 2    4.2615 10 
c3   2.846910 4 
  
2 
c4   3.3837 10 
The stress on the inside radius of the outer cylinder is
then given by
E 
 1   


 
c
1



c
3
4
2 
1  2 
r



4
 2  1  0.3  


2
.
8469

10
1

0
.
3

3
.
3837

10


2 
 6.5 

 30683psi
30106

1  0.32
Finding the inverse of a square matrix
The inverse [B] of a square matrix [A] is defined as
[A][B] = [I] = [B][A]
http://numericalmethods.eng.usf.edu
Finding the inverse of a square matrix
How can LU Decomposition be used to find the inverse?
Assume the first column of [B] to be [b11 b12 … bn1]T
Using this and the definition of matrix multiplication
First column of [B]
b11  1
b  0
A  21    
    
   
bn1  0
Second column of [B]
 b12  0
 b  1
A  22    
    
   
bn 2  0
The remaining columns in [B] can be found in the same manner
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Find the inverse of a square matrix [A]
 25 5 1
A   64 8 1
144 12 1
Using the decomposition procedure, the [L] and [U] matrices are found to be
1 
 1 0 0 25 5
A  LU   2.56 1 0  0  4.8  1.56
5.76 3.5 1  0
0
0.7 
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Solving for the each column of [B] requires two steps
1) Solve [L] [Z] = [C] for [Z]
2) Solve [U] [X] = [Z] for [X]
Step 1:
0 0  z1  1
 1
LZ   C   2.56 1 0  z2   0
5.76 3.5 1  z3  0
This generates the equations:
z1  1
2.56z1  z2  0
5.76z1  3.5z2  z3  0
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Solving for [Z]
z1  1
z 2  0  2.56z1
 0  2.561
 2.56
z3  0  5.76z1  3.5 z 2
 z1   1 
Z    z2    2.56
 z3   3.2 
 0  5.761  3.5 2.56
 3.2
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Solving [U][X] = [Z] for [X]
5
1  b11   1 
25
 0  4.8  1.56 b    2.56

  21  

 0
0
0.7  b31   3.2 
25b11  5b21  b31  1
 4.8b21  1.56b31  2.56
0.7b31  3.2
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Using Backward Substitution
3.2
 4.571
0.7
2.56  1.560b31
b21 
4.8
2.56  1.5604.571

 0.9524
4 . 8
1  5b21  b31
b11 
25
1  5 0.9524  4.571

 0.04762
25
b31 
So the first column of
the inverse of [A] is:
b11   0.04762
b    0.9524
 21  

b31   4.571 
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
Repeating for the second and third columns of the inverse
Second Column
 25 5 1 b12  0
 64 8 1 b   1

  22   
144 12 1 b32  0
b12   0.08333
b    1.417 
 22  

b32    5.000 
Third Column
 25 5 1  b13  0
 64 8 1 b   0

  23   
144 12 1 b33  1
b13   0.03571
b    0.4643
 23  

b33   1.429 
http://numericalmethods.eng.usf.edu
Example: Inverse of a Matrix
The inverse of [A] is
 0.04762  0.08333 0.03571
A1   0.9524 1.417  0.4643
 4.571
 5.000
1.429 
To check your work do the following operation
[A][A]-1 = [I] = [A]-1[A]
http://numericalmethods.eng.usf.edu
Additional Resources
For all resources on this topic such as digital audiovisual
lectures, primers, textbook chapters, multiple-choice
tests, worksheets in MATLAB, MATHEMATICA, MathCad
and MAPLE, blogs, related physical problems, please
visit
http://numericalmethods.eng.usf.edu/topics/lu_decomp
osition.html
THE END
http://numericalmethods.eng.usf.edu