Numerical Methods in Engineering (ENGR20005)
Lecture 06 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Lecture Recap
We have established that the following type of systems are easy to solve.
Diagonal
a11 0 0 . .0. . 0
Lower triangular
Upper triangular
aaaa…a 0
[]
00aa…a
033 34 3n aaa0…0
0 a22 0 . . . . 0 00a….0
11 12 13 14 1n a11000…0 0 a22 a23 a24 . . . a2n a21a2200…0
[] []
033 31 32 33
A=……..
A=000a…a A=
44 4n a41a42a43a44…0
…….. ……..
……..
…….. …….. …….. ……..
…….0
00 ……..
0…0a 000….a aaa….a nn nn n1 n2 n3 nn
In general, methods of solving general linear systems involves first “transforming” the general system to look like any of the 3 systems above and then solving the resulting “easier” systems.
2
Lecture Recap
Systems of Linear Algebraic Equations
The available methods to solve systems of linear algebraic equations can be broadly classified below
Solution of linear algebraic equations
Direct methods
Iterative methods
• Gauss Elimination • Gauss Jordan
• LU Decomposition
• Gauss-Seidel
• Point-Jacobi
• Successive over- relaxation
3
Difficult to solve
a11 a12 a13 . . . . a1n | c1 a21 a22 a23 . . . . a2n | c2 a31 a32 a33 . . . . a3n | c3
. . . . . . . . | . ……..|. ……..|. ……..|.
an1 an2 an3 . . . . ann | cn
a1′1 a1′2 a1′3 a1′4 . . . a1n 0 a2′2 a2′3 a2′4 . . . a2′n 0 0 a3′3 a3′4 . . . a3′n 0 0 0 a4′4 . . . a4′n
| c1′ | c2′ | c3′
Lecture Recap
Easy to solve See Example L05.3
Gauss Elimination
4
| c4′ . 0. . . . . . . | . ……..|. ……..|.
0 0 0 . . . . an′n | cn′
****Try solving using Gauss elimination
2×2 + 3×3 = 1 3×1 + 5×2 + 6×3 = 2 9×1 + 2×2 + 3×3 = 3
5
Gauss Elimination – Pitfalls
The pitfalls of the Gauss elimination method for solving a linear system of equations are:
• •
Very expensive, compared to other methods
The formula for forward elimination stage requires we divide by diagonal components of matrix[A].
•
by zero.
Need to rearrange equations to ensure that we are not dividing
For example the following system
2×2 + 3×3 = 1 3×1 + 5×2 + 6×3 = 2 9×1 + 2×2 + 3×3 = 3
6
Gauss Elimination – Pitfalls
For this case, a11=0, will run into problems in first forward elimination stage. Re-writing equation as
3×1 + 5×2 + 6×3 = 2 2×2 + 3×3 = 1 9×1 + 2×2 + 3×3 = 3
will eliminate problem. The process of interchanging rows so that you are not dividing by zero is called partial pivoting (also try to avoid dividing by small numbers).
7
L5.5:
Solving linear systems of algebraic equations Gauss Elimination & solving Upper triangular system Section 3.1.1 of “Book”
Example L05.6:Write a MATLAB function that uses Gauss Elimination to transform the following system to upper triangular form and then use the UpperTriangularSolver() function developed in Example L05.3 to solve the resulting upper triangular system.
1234×1 3 3 4 8 9 x2 = 4 82043 x3 8 5678×4 10
Difficult to solve
1234
[A]= 3 4 89 82043
5678
348 10
Easy to solve
1234
[Aprime]= 0−2 −1 −3 0 0 −22 −35
Gauss Elimination (Example L05.5)
UpperTriangularSolver() (Example L05.3)
0
0 0 0 3.5455
3 −5 −26 12.1
{C} =
{Cprime} =
x1 x2 x3 x4
9
Lecture05M06.m
A=[1 2 3 4; 3 4 8 9;
8 20 4 3;
5 6 7 8]; C=[3;4;8;10];
newA=GaussElimination(A,C);
[n,n]=size(A);
Aprime=newA(1:n,1:n) Cprime=newA(:,n+1)
x=UpperTriangularSolver(Aprime,Cprime)
function newA=GaussElimination(A,C) %Augment A and C
[n,n]=size(A);
newA=[A C];
for k=1:n-1 for i=k+1:n
factor=newA(i,k)/newA(k,k); for j=k:n+1
newA(i,j)=newA(i,j)-factor*newA(k,j); end
end end
end
function x=UpperTriangularSolver(A,C)
[n,n]=size(A) % Finding the size of the problem
x=zeros(n,1) % set up x vector with elements initially set to zero x(n)=C(n)/A(n,n);
for i=n-1:-1:1 sum=0;
for j=i+1:n
sum=sum+A(i,j)*x(j); end
x(i)=(C(i)-sum)/A(i,i); end
end
10
End of Example L05.6
Finding inverse of matrix
Suppose that you want to find an inverse of a matrix
a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44
b11 b21 b31 b41
b12 b13 b14 b22 b23 b24 b32 b33 b34 b42 b43 b44
[A]
a11 a12 a13 a14 b11 1 a21 a22 a23 a24 b21 = 0 a31 a32 a33 a34 b31 0 a41 a42 a43 a44 b41 0
a11 a12 a13 a14 b13 0 a21 a22 a23 a24 b23 = 0 a31 a32 a33 a34 b33 1 a41 a42 a43 a44 b43 0
[A]-1
=
1000 0100 0010 0001
a13 a14 b12 0
a11 a12
a21 a22 a23 a24 b22 = 1
0 0
a31 a32 a33 a34 b32 a41 a42 a43 a44 b42
b14
a21 a22 a23 a24 b24 = 0
a11 a12 a13 a14
0
12
a31 a32 a33 a34 b34 0 a41 a42 a43 a44 b44 1
L6.1
LU Decomposition: Context
Systems of Linear Algebraic Equations
The available methods to solve systems of linear algebraic equations can be broadly classified below
Solution of linear algebraic equations
Direct methods
Iterative methods
• Gauss Elimination • Gauss Jordan
• LU Decomposition
• Gauss-Seidel
• Point-Jacobi
• Successive over- relaxation
14
LU Decomposition
Aim is to solve
If [A] is a `full’ matrix, it is difficult to solve.
[A] {X} = {C}
Let’s see what happens if we re-express [A] such that
(6.1) (6.2)
[A] = [L][U] ,
where [L] and[U] are lower- and upper-triangular matrices.
The structure of [L] and [U] are shown below
l11 0 0 0 . . 0 l21 l22 0 0 . . 0
[L]= l31 l32 l33 0 . . 0 ……. ……. …….
ln1 ln2 ln3 ln4 . . lnn
1 u12 u12 u12 . . u1n
0 1 u23 u24 . . u2n
0
0 0 1 u34 . . u3n
15
[U]=
…….
0
…….
…….
0000..1
LU Decomposition
If we now substitute Eq. (6.2) into Eq. (6.1), then we will have
[A] Let
= {C} = {R}
(6.3) (6.4)
[L]
[U]
{X}
[U]{X}
and substitute Eq. (6.4) into Eq. (6.3) to give
[L] {R} = {C} 16
LU Decomposition
Thus we have replaced the original problem,
by two problems
[A] {X} = {C}
[L] {R} = {C} [U] {X} = {R}
(6.5) However, the two equations are much easier to solve because [L]
and [U] are lower- and upper-triangular matrices.
This is because, to solve the Eqs. (6.5), we only need backward and
forward substitution, respectively.
17
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
[U]{X}={R} {X}
2 forward substitution
18
3 backward substitution
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
2 forward substitution
[U]{X}={R} {X}
3 backward substitution
19
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
2 forward substitution
[U]{X}={R} {X}
3 backward substitution
20
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[L] [L]{R}={C}
{R}
[U]{X}={R} {X}
[A]{X}={C} [U]
2 forward substitution
21
3 backward substitution
}
}
LU Decomposition
Summary of the LU decomposition method in 1 slide. We would like to find {X} by solving
[A] {X} = {C} Decompose [A] into [L] and [U]
(6.6)
(6.7)
(6.8)
(6.9)
If we let
then Eq. (6.7) becomes
[L][U]{X} = {C} [U] {X} = {R}
[L] {R} = {C}
Often Eq. (6.6) is difficult to solve. Using LU Decomposition, we replace (6.6) by two problems Eqs. (6.8) and (6.9) that are much easier to solve (see Examples L05.3 and L05.4). We are given {C} and will assume we know [L] and [U].
1. Solve Eq. (6.9) to get {R}
2. Substitute {R} from step 1 into the right hand side of Eq. (6.8) and solve to get {X}.
22
LU Decomposition
Use the following procedure to solve the system
[A] {X} = {C}
1.Decompose[A] into[L] and [U]
2.Solve [L]{R} ={C} and obtain {R} by forward substitution. 3.Use {R} obtained in above step and substitute into [U]{X} = {R}
so that you can solve for {X} by backward substitution.
We know how to do steps 2 and 3. So now need to do step 1, decomposing [A] into [L] and [U]. How can we achieve that?
23
L6.2
LU Decomposition: How do we get [L] and [U]?
Sequence of computation for LU Decomposition
a11 a21 a31 a41
a12 a13 a14
a23 a24 a34 a44
a22 a32 a42
25
a33 a43
LU Decomposition
To illustrate the procedure of finding [L] and [U] for a more general matrix, first consider a smaller 4 x 4 [A] matrix.
[A] [L] [U]
l11 0 0 0
a12 a13 a14
a22 a23 a24 = l21 l22 0 0
u12 u13 u14 1 u23 u24
0 1 u34 001
l11 = a11 l21 = a21 l31 = a31 l41 = a41
a11
a21 a31 a41
a32 a33 a34 a42 a43 a44
l31 l32 l33 0 l41 l42 l43 l44
26
1 0 0 0
LU Decomposition
To illustrate the procedure of finding [L] and [U] for a more general
matrix, first consider a smaller 4 x 4 [A] matrix. [A] [L]
l11 0 0 0 a22 a23 a24 = l21 l22 0 0
[U]
u13 u14
1
0
0 1u34 001
l11 = a11 l21 = a21
l =a 31 31
l41 = a41
a11 a21 a31 a41
a12
a13 a14
a32 a33 a34 a42 a43 a44
l11u12 = a12
u12 = a12/l11
l11u13=a13 u13 = a13/l11
l11u14 = a14
u14 = a14/l11
l31 l32 l33 0 l41 l42 l43 l44
u23 u24
27
u12 1
0 0
LU Decomposition
To illustrate the procedure of finding [L] and [U] for a more general
matrix, first consider a smaller 4 x 4 [A] matrix. [A] [L]
l11 = a11 l21 = a21 l31 = a31 l41 = a41
l11u12 = a12
u12 = a12/l11
l11u13 = a13 u13 = a13/l11
l11u14 = a14
u14 = a14/l11
l21u12 + l22 = a22
l22 = a22 − l21u12
l31u12 + l32 = a32
l32 = a32 − l31u12
l41u12 + l42 = a42
l42 = a42 − l41u12
a11 a21 a31 a41
l11 0 0 0 a23 a24 = l21 l22 0 0
[U]
u13 u14
a12 a13 a14
1
0
0 1u34 001
a22 a32 a42
a33 a34 l31 l32 l33 0 a43 a44 l41 l42 l43 l44
u23 u24
28
u12 1
0 0
LU Decomposition
To illustrate the procedure of finding [L] and [U] for a more general
matrix, first consider a smaller 4 x 4 [A] matrix. [A] [L]
a11 a21 a31 a41
l11 0 0 0 a33 a34 = l31 l32 l33 0
[U]
001 l21u13 + l22u23 = a23
u23 = (a23 − l21u13)/l22 l u +l u =a
a12 a13 a14
1 u12
01 u24
a22 a32 a42
a23
l11 = a11
l21 = a21
l31 =a31 l =a
l11u12 = a12
u12 = a12 /l11
l11u13 =a13 u13 =a13/l11
l21u12 + l22 = a22
l22 = a22 − l21u12
l31u12+l32 =a32
l32 =a32−l31u12
l u +l =a
21 14
22 24 24
u24 =(a24−l21u14)/l22
l11u14 =a14
a43 a44 l41 l42 l43 l44
u14 00 u34
41 41 u=a/l 4112 42 42
14 14 11
l42 = a42 − l41u12
a24
29
l21 l22 0 0
u13 u23
1 0
LU Decomposition
If you keep going you get
l33 = a33 − l31u13 − l32u23, l43 = a43 − l41u13 − l42u23
a34 − l31u14 − l32u24 l33
u34 =
l44 = a44 − l41u14 − l42u24 − l43u34
30
LU Decomposition
In general, the formula for getting elements of [L] and [U] from elements of an n x n [A] matrix can be written as
Crout’s Algorithm
li1=ai1 for i=1,2,…,n u1j=a1j for j=2,3,…,n
Beginl11for loop j−1
lij =aij−∑likukj k=1
for k=j+1,j+2,…,n uii = 1.0 for i = 1,2,…,n
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
j = 2,3,…, n
for i=j,j+1,…,n
End for loop
This procedure is known as Crout’s LU decomposition algorithm.
31
L6.3
LU Decomposition: A numerical example
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
[U]{X}={R} {X}
2 forward substitution
33
3 backward substitution
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
2 forward substitution
[U]{X}={R} {X}
3 backward substitution
34
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
2 forward substitution
[U]{X}={R} {X}
3 backward substitution
35
}
}
LU Decomposition – Example 1
Let’s illustrate the LU decomposition with a simple 3×3 example
[1 2 3] x1 {14} 3 2 2 x2 = 16 232×3 17
[1 2 3] x1 {C}=[14] [A]= 3 2 2 {X}= x2 16
232 x3 17 We decompose [A] into [L] and [U] as follows
[1 2 3] l11 0 0 1 u12 u13 [A]= 3 2 2 = l21 l22 0 0 1 u23
2 3 2 l31 l32 l33 0 0 1
Note that to find [L] and [U] for this problem, we need to find 9
unknowns, l11, l21, l22, l31,l32, l33, u12, u13, u23. 36
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[A]=[1 2 3]= l l 0 322 2122
2 3 2 l31 l32 l33
Consider the elements in the first column of the [A] matrix:
u12 u13 1u23
Crout’s Algorithm
0 1
li1=ai1 for i=1,2,3 ⇒
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
l11 0 0
l =a =3 21 21
l31=a31=2
j =2,3,…,n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
37
1 0 0
l11 =a11 =1
l11
Begin for loop
l=a−j−1lu ij ij ∑k=1 ik kj
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[A]=[1 2 3]= l11 0 0 3 2 2
2 3 2 l31 l32 l33
Consider the elements in the first column of the [A] matrix:
u12 u13 1 u23
Crout’s Algorithm
0 1
li1 =ai1 for
i=1,2,3 ⇒
l11 =a11 =1 l31 =a31 =2
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
l21 l22 0
38
1 0 0
l21 =a21 =3
l11
Begin for loop
j−1
lij=aij−∑likukj k=1
j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[1 2 3] l11 0 0 u12 u13 [A]= 3 2 2 = l21 l22 0 1 u23
232 01
Consider the elements in the first column of the [A] matrix:
Crout’s Algorithm
l =a for i=1,2,3 ⇒ l =a =1
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
i1 i1
11 11
l11
Begin for loop
l=a−j−1lu
i j i j ∑k = 1 i k k j
l
=a =3 21
j =2,3,…,n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
21
l31 l32 l33
39
1 0 0
l31 =a31 =2
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[A]=[1 2 3]= l l 0 1 u12 u13 3 2 2 21 22 0 1 u23
2 3 2 l31 l32 l33 0 0 1 Consider the elements in the first column of the [A] matrix:
li1 =ai1 for i=1,2,3 ⇒ l11 =a11 =1 l21 =a21 =3 l31 =a31 =2
l11 0 0
Now consider the (unused) elements in the 1st row of the [A] matrix
a1j a
u= forj=2,3⇒u=12=2
Crout’s Algorithm
1j l 11
12 l 11
li1 = ai1 for i = 1,2,…,n u =a1j for j=2,3,…,n
u=a13=3
1j l11
Begin for loop
j−1
lij=aij−∑likukj k=1
j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n
13
l11
40
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
0 1 u13 3222122 0u23
l11 0 0
[A]=[1 2 3]= l l
232 l31l32l33 0 1
Consider the elements in the first column of the [A] matrix:
li1 =ai1 for i=1,2,3 ⇒ l11 =a11 =1 l21 =a21 =3 l31 =a31 =2
Now consider the (unused) elements in the 1st row of the [A] matrix
u =a1j for j=2,3 ⇒ 1j l
Crout’s Algorithm
li1 = ai1 for i = 1,2,…,n u =a1j for j=2,3,…,n
u12 = a12 = 2 l11
11
u=a13=3
1j l11
13
l11
41
Begin for loop j = 2,3,…, n j−1
lij=aij−∑likukj for i=j,j+1,…,n k=1
ujk= i=1 for k=j+1,j+2,…,n ljj
j−1
ajk − ∑ ljiuik
u12 1
0
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
0 1 u12 3222122 01
l11 0 0
[A]=[1 2 3]= l l
2 3 2 l31 l32 l33 0 0
Consider the elements in the first column of the [A] matrix:
100 li1 =ai1 for i=1,2,3 ⇒ l11 =a11 =1 [L]= 3 l22 0
l21=a21=3 2 l32 l33 l31 =a31 =2
Now consider the (unused) elements in the 1st row of the [A] matrix
u=a1jforj=2,3⇒u=a12=2
Crout’s Algorithm
li1 = ai1 for i = 1,2,…,n
123 0 1 u23
001
1j
12
u =a1j for j=2,3,…,n 1j
ll
11 11 l
Begin for loop j = 2,3,…, n
a13 11 [U] =
u==3 13 l
42
l=a−j−1lu fori=j,j+1,…,n i j i j ∑k = 1 i k k j
for k=j+1,j+2,…,n
11
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
u13
u23 1
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[1 2 3] l11 0 0 1 u12 u13 [A]= 3 2 = l21 l22 0 0 1 u23
2 2 l31 l32 l33 0 0 1 Now let’s look at the “unused” elements in 2nd column of [A]
li2 =ai2−∑2−1likuk2 k=1
for
i=2,3
⇒
l22 =a22−l21u12 =2−3⋅2=−4 l32 =a32−l31u12 =3−2⋅2=−1
2 3
43
ujk= i=1 ljj
Crout’s Algorithm
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
l11
Begin for loop
j−1
lij=aij−∑likukj k=1
j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
j−1
ajk − ∑ ljiuik
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[A]=[1 2 3]= l11 0 0 1 u12 u13 3 2 0 1 u23
2 2 l31 l32 l33 0 0 1 Now let’s look at the “unused” elements in 2nd column of [A]
l21 l22 0
2
3
li2 =ai2−∑2−1likuk2 k=1
for i=2,3 ⇒
l22 =a22−l21u12 =2−3⋅2=−4
l32 =a32−l31u12 =3−2⋅2=−1 Crout’s Algorithm
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
l11
Begin for loop
j−1
lij=aij−∑likukj k=1
j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
j−1
ajk − ∑ ljiuik
44
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[1 2 3] l11 0 0 1 u12 u13 [A]= 3 2 = l21 l22 0 0 1 u23
22 001
Now let’s look at the “unused” elements in 2nd column of [A]
2
3
li2 =ai2−∑2−1likuk2 k=1
for i=2,3 ⇒
l22 =a22−l21u12 =2−3⋅2=−4 Crout’s Algorithm
l31 l32 l33
l32 =a32−l31u12 =3−2⋅2=−1
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
l11
Begin for loop
j−1
lij=aij−∑likukj k=1
j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii=1.0 for i=1,2,…,n
j−1
ajk − ∑ ljiuik
45
ujk= i=1 ljj
End for loop
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[A]=[1 2 3]= l11 0 0 1 u12 3201
2−1
li2 =ai2−∑likuk2 for i=2,3 ⇒
2 3
l21 l22 0
22
l31 l32 l33
Now let’s look at the “unused” elements in 2nd column of [A]
k=1
Now let’s look at the “unused” elements in 2nd row of [A]
2−1
a2k − ∑ l2iuik
l22 =a22−l21u12 =2−3⋅2=−4 l32 =a32−l31u12 =3−2⋅2=−1
00
u2k= i=1 for k=3 ⇒ l22
u23 = a23 −l21u13 = 2−3⋅3 = 7 l22 −44
46
u13
u23 1
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[1 2 3] l11 0 0 1 u12 u13 [A]= 3 2 = l21 l22 0 0 1 u23
2 2 l31 l32 l33 0 0 1
2 3
Now let’s look at the “unused” elements in 2nd column of [A]1 2−1
0 0 li2 =ai2−∑likuk2 for i=2,3 ⇒ l22 =a22−l21u12 =2−3⋅2=−4 [L]= 3 −4 0
k=1 l32 =a32−l31u12 =3−2⋅2=−1 Now let’s look at the “unused” elements in 2nd row of [A]
2−1 a2k−∑l2iuik
a −l u 2−3⋅3 7
2 −1 l33 [U]=[12 3]
u2k= i=1 l22
fork=3 ⇒u23=23
2113= = l22 −44
0 1 7/4 001
47
LU Decomposition – Example 1
Consider the following 3×3 matrix A
[1 2 3] l11 0 0 1 u12 [A]= 3 2 = l21 l22 0 0 1
2200
Now let’s look at the last element.
2 3
l31 l32 l33
∑3−1
l =a − l u i3 i3 k=1ikk3
Crout’s Algorithm
li1 = ai1 for i = 1,2,…,n u1j=a1j for j=2,3,…,n
for i=3
11
Beginl for loop j=2,3,…,n
48
j−1
lij=aij−∑likukj k=1
u13
u23 1
for i=j,j+1,…,n for k=j+1,j+2,…,n
End for loop
uii=1.0 for i=1,2,…,n
l33 =a33−(l31u13+l32u23)=2−(2⋅3−1⋅74)=−94
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
[U]{X}={R} {X}
2 forward substitution
49
3 backward substitution
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
2 forward substitution
[U]{X}={R} {X}
3 backward substitution
50
}
}
LU Decomposition
The flow of information in the LU decomposition method.
1 decomposition
[A]{X}={C} [L] [U]
[L]{R}={C} {R}
[U]{X}={R} {X}
2 forward substitution
3 backward substitution
51
}
}
LU Decomposition – Example 1
Thus we get the upper and lower diagonal matrices
[A]=[L][U]=[1 2 3] 322
232
[L]=10 0,[U]=[123], 3−4 0 017/4
[L]{R}={C}, solve for {R} by forward substitution.
l11 0 0 l21 l22 0
r1 c1 r2 = c2
2 −1 −9/4
0 0 1
l31 l32 l33
r3 c3
52
LU Decomposition – Example 1
[L]{R}={C}, solve for {R} by forward substitution.
l31 l32 l33 From the first row,
From the 2nd row, c2 − l21r1 l21r1+l22r2=c2 ⇒ r2= l22
From the third row,
16 − 3 ⋅ 14 = −4
=6.5
{R}= 6.5
l11 0 0 l21 l22 0
r1 c1 r2 = c2
r3 c3 l11r1 = c1 ⇒ r1 = c1 = 14 = 14
r3=c3−l31r1−l32r2 =17−2⋅14+6.5=2
{14} 2
l33
−9/4
l11 1
53
LU Decomposition – Example 1
Lastly, we need to solve
{R} = 6.5 {14}
[U] {X} = {R} 1 u12 u13 x1
2
r1 0 1 u23 x2 = r2
001 x3 r3
We can now use back substitution. Start with last (third) row
x3 = r3 = 2 7
x2 = r2 − u23x3 = 6.5 − 4 ⋅ 2 = 3
x1 = r1 − u12x2 − u13x3 = 14 − 2 ⋅ 3 − 3 ⋅ 2 = 2 .
The solution is therefore, X = [232] 54
Example L06.1: Write a MATLAB function that finds the lower and upper triangular matrices, [L] and [U] for any general nxn matrix. Test your function on this [A] matrix.
1234 3489 8 20 4 3 5678
[A] =
Crout’s Algorithm
li1=ai1 for i=1,2,…,n u1j = a1j for j = 2,3,…,n
Beginl11for j−1
loop
j = 2,3,…,n
for i=j,j+1,…,n
for k=j+1,j+2,…,n uii = 1.0 for i = 1,2,…,n
lij =aij−∑likukj k=1
j−1
ajk − ∑ ljiuik
ujk= i=1 ljj
End for loop
55
Lecture06M01.m
A=[1 2 3 4; 3 4 8 9;
8 20 4 3;
5 6 7 8]; C=[3;4;8;10];
[L,U]=LUDecomposition(A)
%
% Check L*U-A must be close to zero %
L*U-A
function [L,U]=LUDecomposition(A)
[n,n]=size(A); L=zeros(n,n); U=zeros(n,n);
for j=2:n for i=j:n
sum=0; for k=1:j-1
sum=sum+L(i,k)*U(k,j) end
L(i,j)=A(i,j)-sum;
end
for k=j+1:n sum=0;
for i=1:j-1 sum=sum+L(j,i)*U(i,k)
end
U(j,k)=(A(j,k)-sum)/L(j,j);
end end
for i=1:n U(i,i)=1.0;
end end
for i=1:n L(i,1)=A(i,1)
end
for j=2:n U(1,j)=A(1,j)/L(1,1)
end
l11 = a11 l21 = a21 l31 = a31 l41 = a41
Crout’s Algorithm
for i = 1,2,…, n for j = 2,3,…, n
Begin for loop j = 2,3,…, n
for i=j,j+1,…,n
for k=j+1,j+2,…,n
End for loop
uii = 1.0 for i = 1,2,…,n
l33 =a33−l31u13−l32u23
li1 = ai1
l11u12 = a12
u12 = a12/l11
l11u13 = a13 u13 = a13/l11
l11u14 = a14
u14 = a14/l11
u1j = a1j l11
j−1
lij = aij − ∑ likukj k=1
j−1
ajk − ∑ ljiuik
ujk = i=1 ljj
56
u34 = a34 − l31u14 − l32u24 l33
Example L06.2: Use your LU decomposition function from the previous example to obtain the solution to the following system of equations.
1234×1 3 3 4 8 9 x2 = 4 82043 x3 8 5678×4 10
57
1234×1 3 3 4 8 9 x2 = 4 82043 x3 8 5678×4 10
[A] = [L][U]
[U] {X} = {R}
[L] {R} = {C} 58
Lecture05M03.m
A=[1 2 3 4; 3 4 8 9;
8 20 4 3;
5 6 7 8]; C=[3;4;8;10];
1234 3
[A]= 3 4 8 9 {C}= 4 8 20 4 3 8
10
5678
[L,U]=LUDecomposition(A)
R=LowerTriangularSolver(L,C) X=UpperTriangularSolver(U,R)
function [L,U]=LUDecomposition(A)
[n,n]=size(A); L=zeros(n,n); U=zeros(n,n);
for i=1:n L(i,1)=A(i,1)
end
for j=2:n
U(1,j)=A(1,j)/L(1,1) end
for j=2:n for i=j:n
sum=0; for k=1:j-1
sum=sum+L(i,k)*U(k,j) end
L(i,j)=A(i,j)-sum; end
for k=j+1:n sum=0;
for i=1:j-1 sum=sum+L(j,i)*U(i,k)
end
U(j,k)=(A(j,k)-sum)/L(j,j); end
end
for i=1:n U(i,i)=1.0;
end end
function x=UpperTriangularSolver(A,C)
[n,n]=size(A) % Finding the size of the problem
x=zeros(n,1) % set up x vector with elements initially set to zero x(n)=C(n)/A(n,n);
for i=n-1:-1:1 sum=0;
for j=i+1:n
sum=sum+A(i,j)*x(j); end
x(i)=(C(i)-sum)/A(i,i); end
end
function x=LowerTriangularSolver(A,C)
[n,n]=size(A); % Finding the size of the problem
x=zeros(n,1); % set up x vector with elements initially set to zero x(1)=C(1)/A(1,1);
for i=2:n
sum=0; for j=1:i-1
sum=sum+A(i,j)*x(j); end
x(i)=(C(i)-sum)/A(i,i); end
end
LU Decomposition- Cost
To solve the system [A] {X} = {C}, the number of operations is:
1. decompose [A] into [L] [U]:
2. solve [L] {R}= {C}:
3. solve [U] {X}= {C}:
2/3 n3 operations n2 operations
n2 operations
Compare that with directly inverting matrix [A]
1. compute [A]-1: 2n3 operations
2. multiply{X}= [A]-1{C}: 2n2 operations
LU method takes 2/3 of the time required compared to directly inverting matrix 59
L6.4
Systems of linear equations: Iterative method
Systems of Linear Algebraic Equations
The available methods to solve systems of linear algebraic equations can be broadly classified below
Solution of linear algebraic equations
Direct methods
Iterative methods
• Gauss Elimination • Gauss Jordan
• LU Decomposition
• Gauss-Seidel
• Point-Jacobi
• Successive over- relaxation
61
Iterative Methods
If the size of [A] becomes very large, the cost of direct methods can be prohibitive (scales as n3).
An alternative are iterative methods. Let’s rewrite [A] as
[A] = [M] − [N]
With this, our system of equations ([A] {X}= {C}) becomes
[M] {X} = [N] {X} + {C} 62
Iterative Methods
An iterative method of solving for the vector{X} would be [M] {X}(k+1) = [N] {X}(k) + {C} ,
where {X}(k) is current (k’th) guess of true solution, {X}.
We would like to choose the matrix [M] such that Eq. (6.13) is
easier to solve than the original expression [A]{X}={C}.
(6.13)
Here, we will cover two different methods:
1. Point Jacobi method:[M] consists of only the diagonal elements
of [A].
2. Gauss-Seidel method:[M] consists of the lower triangular
elements of [A].
63
Diagonal
m11 0 0 .
0 m22 0 .
. . . 0 . . . 0
Lower triangular
m11 0 0 0 . . . 0 m21 m22 0 0 . . . 0
[A]{X} = {C}
[A] = [M] − [N]
[M] {X}(k+1) = [N] {X}(k) + {C} ,
0 0 m33 . . . . 0
m31 m32 m33 0 . . . 0 mmmm…0
[M]= . 0. . . . . . . …….. ……..
[M] =
……..
…….0 0 0 0 . . . 0 mnn
Point Jacobi method
41 42 43 44 …….. …….0
0
0
64
mn1 mn2 mn3 . . . . mnn Gauss-Seidel method
Iterative Methods
The general procedure to solving linear systems of equations with an iterative method is
1. Define the [M] and [N] matrices
2. Obtain an initial guess,{X}(0)
3. Solve [M]{X}(k+1) = [N]{X}(k)+{C} for new guess
4. Define and obtain maximum absolute value of residual vector,
{r}= [A]{X}(k)-{C}
5. If maximum value of |{r}|>𝜺, repeat step 3. Else, end.
65
Point Jacobi
In the point Jacobi method, [M] consist of only the diagonal elements of [A]. The negative of all other elements is placed into [N]. This method will be illustrated with a 3 x 3 system
a11 a12 a13 x1 c1 a21 a22 a23 x2 = c2 a31 a32 a33 x3 c3
[M] and [N] are
[A]{X} = {C} 0 0 0 −a12 −a13
[A] = [M] − [N]
[M]= 0 a22 0 , [N]= −a21 0 −a23
a11
0 0 a33 −a31 −a32 0
[M] {X}(k+1) = [N] {X}(k) + {C} For this case, the procedure to solve the system will look like
a11 0 0 x1 (k+1) 0 −a12 −a13 x1 (k) c1
0 a22 0 x2 0 0 a33 x3
= −a21 0 −a23 x2 + c2
−a31 −a32 0 x3 c3 66
Point Jacobi
For this case, system that we need to solve will look like
a11 0 0 x1 (k+1) 0 −a12 −a13 x1 (k) c1
0 a22 0 x2 0 0 a33 x3
= −a21 0 −a23 x2 + c2 −a31 −a32 0 x3 c3
Solving this system of equation will give
x(k+1) = a11 −a12x(k) − a13x(k) + c1 11(23)
x(k+1) = a22 −a21x(k) − a23x(k) + c2 21(13)
x(k+1) = a33 −a31x(k) − a32x(k) + c3 31(12)
67
So for 3×3 system
a11 a12 a13 x1 c1 a21 a22 a23 x2 = c2 a31 a32 a33 x3 c3
x(k+1) = a11 −a12x(k) − a13x(k) + c1 11(23)
x(k+1) = a22 −a21x(k) − a23x(k) + c2 21(13)
x(k+1) = a33 −a31x(k) − a32x(k) + c3 31(12)
Point Jacobi
So for nxn system a11 a12 . . a1n
x1 c1 x2 = c2
x(k+1)= 1 ∑n −aijx(k)+ci
x1(k), x2(k), x3(k) are previous guesses of x1, x2, x3
x1(k+1), x2(k+1), x3(k+1) are latest (better) guesses of x1, x2, x3.
a21 a22 . . a2n
a31a32..a3nx3c3 iaj
……. an1 an2 . . ann xn cn
ii
j=1, j≠i
xi(k) are previous guesses of xi
xi(k+1) are latest (better) guesses of xi.
68
Example L06.3:Write a MATLAB function that uses the Point Jacobi method to iteratively solve a system of equation [A]{X}={C}. Test the function by finding the solution to
9 4 1 x1 {−17} 1 6 0 x2 = 4 1−2−6 x3 14
Lecture06M03.m
A=[9 4 1; 1 6 0;
1 -2 -6]; C=[-17; 4;
14];
x=PointJacobi(A,C)
while max(abs(residual)) > 1.0e-8 for i=1:n
sum=0.0; for j=1:n
if j~=i sum=sum-A(i,j)*xold(j);
end end
x(i)=(sum+C(i))/A(i,i); end
xold=x;
residual=A*x-C; end
end
1−2−60 14 −14
Set initial guess of all your xi to be zero
x1 (0) {0} x2 =0 x3 0
Initial values of residuals are quite large
function x=PointJacobi(A,C) xold=zeros(size(C));
x=xold;
[n,n]=size(A); residual=A*x-C;
9410−17 17 residual= 1 6 0 0 − 4 =− −4
{ } {
} { }
x(k+1)=1 ∑n −aijx(k)+ci i aii j=1,j≠i j
x1 x2 x3
(1)
−2.44 0.5 −2.5
Magnitude of residuals decreases with every iteration.
=
residual = 69
−5.5 −3.44 −2.44
1.88 residual = 0.61
−0.53
residual =
−0.497 −0.21 −0.0062
L6.5
Iterative method: Gauss-Seidel method
Systems of Linear Algebraic Equations
The available methods to solve systems of linear algebraic equations can be broadly classified below
Solution of linear algebraic equations
Direct methods
Iterative methods
• Gauss Elimination • Gauss Jordan
• LU Decomposition
• Gauss-Seidel
• Point-Jacobi
• Successive over- relaxation
71
Gauss-Seidel
The Gauss-Seidel method is the most commonly used iterative method. For conciseness, the method will be illustrated with a 3×3 system
a11x1 + a12x2 + a13x3 = c1 a21x1 + a22x2 + a23x3 = c2 a31x1 + a32x2 + a33x3 = c3
[A]{X} = {C}
For this method, the matrix [M] consists of only the lower triangular elements of [A].
72
Gauss-Seidel
[M] will contain the lower triangular elements of [A]
The negative of all other elements is placed in [N].
a11 0 0 a21 a22 0
[A]=[M]−[N]
[M] =
a31 a32 a33
Thus,
[N] =
0 −a12 −a13 0 0 −a23 000
a11 0 0 x1 (k+1) 0 −a12 −a13 x1 (k) c1
a21 a22 0 x2 = 0 0 −a23 x2 + c2
a31 a32 a33 x3 0 0 0 x3 c3
[M] {X}(k+1) = [N] {X}(k) + {C}
73
Gauss-Seidel
a11 0 0 x1 (k+1) 0 −a12 −a13 x1 (k) c1 a21 a22 0 x2 = 0 0 −a23 x2 + c2 a31 a32 a33 x3 0 0 0 x3 c3
This is lower triangular system. Solving this system in sequence from the top row to the bottom row will give you
c −a x(k)−a x(k)
x(k+1) = 1 12 2 13 3 (6.13)
1
x(k+1) = 2
c −a x(k+1)−a x(k+1) 3 31 1 a 32 2
33
a11
c −a x(k+1)−a x(k)
2 21 1a 23 3 22
(6.14) (6.15)
x(k+1) = 3
We can now start the iterative process. WARNING, you have to
be careful that none of your aii’s are zero !! 74
Gauss-Seidel
1. Assume initial guess values for x1, x2, and x3: x1(0), x2(0) and x3(0).
2. Substitute the guess values for x2(0), and x3(0) into Eq. (6.16) to
obtain a better estimate for x1. Call this value x1(1) .
3. Substitute x1(1) and x3(0) into Eq. (6.17) to obtain a new estimate
for x2 =x2(1).
4. Use x1(1) and x2(1) in Eq. (6.18) to obtain a better estimate for x3(1)
= x3(1).
1
c −a x(k)−a x(k)
x(k+1) = 1 12 2 13 3 (6.16)
x(k+1) = 2
2 21 1a 23 3
(6.17) (6.18)
a11
c −a x(k+1)−a x(k)
x(k+1) = 3
3 31 1 a 33
32 2
22
c −a x(k+1)−a x(k+1)
75
Gauss-Seidel
The flow of information in the Gauss-Seidel method
x= 1
c-ax-ax 1 12 2 13 3
a
11
x2 = x3 =
x= 1
c-ax-ax 2 21 1 23 3
First iteration
a22 c-ax-ax
3 31 1 32 2 a33
c-ax-ax 1 12 2 13 3
a
11
c-ax-ax
x2 = 2 21 1 23 3 Second iteration
a22 c-ax-ax
x3 =
3 31 1 32 2 a33
76
Gauss-Seidel
If [A] is an n x n matrix, then the generic formula for xi(k+1) is
x
=
+ci
1 i−1 (k+1) n (k) −∑aijx − ∑ aijx
(k+1)
i aii j=1 j j=i+1 j
Note that one can extend the point Jacobi and Gauss-Seidel methods to solve nonlinear systems too!
77
Example L06.4:Write a MATLAB function that uses the Gauss-Seidel method to iteratively solve a system of equation [A]{X}={C}. Test the function by finding the solution to
9 4 1 x1 {−17} 1 6 0 x2 = 4 1−2−6 x3 14
Lecture06M04.m
A=[9 4 1; 1 6 0;
1 -2 -6]; C=[-17; 4;
14]; x=GaussSeidel(A,C)
function x=GaussSeidel(A,C) x=zeros(size(C)); [n,n]=size(A);
residual=A*x-C
while max(abs(residual)) > 1.0e-8 for i=1:n
sum=0.0; for j=1:n if j~=i
sum=sum-A(i,j)*x(j); end
end
x(i)=(sum+C(i))/A(i,i);
end
residual=A*x-C; end
end
Very similar to Example L06.3
x
=
+ci
1 i−1 (k+1) N (k) −∑aijx − ∑ aijx
(k+1)
i aii j=1 j j=i+1 j
78