Numerical Methods in Engineering (ENGR20005)
Lecture 05 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Preallocating Memory
Preallocating Memory
close all
clear all
tic
n = 10000000;
a = zeros(n,1); % Preallocate Memory
a(1) = 0;
for i = 1:n-1
a(i+1) = a(i)+i;
end
toc
Preallocating Memory
Preallocating Memory
Credits:Wilson
L5.1:
Examples of applications of Systems of linear algebraic equations
This example is taken from “Introduction to MATLAB 6 for Engineers” by W. J. Palm III
Example L05.1: Solve the circuit below for currents i1, i2,…,i5. You are given that R1=5, R2=150, R3=100, R4=250 and R5=200 kilo ohms. V1=100 and V2=50 volts.
i1 R1 i3 R3 R5 i5
++
V1 R2 R4 V2 –
i4
–
i2
7
This example is taken from “Introduction to MATLAB 6 for Engineers” by W. J. Palm III
i1 R1 i3 R3 R5 i5
++
V1 R2 R4 V2 –
–
i4 Applying Kirchoff’s voltage law gives
We also know that
i2
R1i1 + R2i2 = V1 −R2i2 + R3i3 + R4i4 = 0
R4i4 − R5i5 = V2 i1 − i2 − i3 = 0
i3 − i4 − i5 = 0 8
Applying Kirchoff’s voltage law gives
R1i1 + R2i2 = V1 −R2i2 + R3i3 + R4i4 = 0
We also know that
i1 − i2 − i3 = 0 i3 − i4 − i5 = 0
R4i4 − R5i5 = V2
Hence, in order to find the currents, i1, i2,…,i5, you need to solve
where
R1 R2 0 0 0 i1
0 −R2 R3 R4 0 i2
V1 V0
0 0 0 R4 −R5 i3
=
0
1 −1 −1 0 0 i4 0 0 1 −1 −1 i5
R1 R2 0 0 0 0 −R2 R3 R4 0
001−1−1 9
02 [A] {X} = {C}
i1
i2 {X}= i3
V
[A]= 0 0 0 R4 −R5 , 1 −1 −1 0 0
,
01
i4 i5
{C}= V2
0 0
Lecture05M01.m
R1=5; R2=150; R3=100; R4=250; R5=200;
V1=100; V2=50;
A=[R1R20 00; 0 -R2 R3 R4 0; 0 0 0 R4 -R5;
1 -1 -1 0 0;
0 0 1 -1 -1];
C=[V1;0;V2;0;0];
x1=A\C x2=linsolve(A,C)
R1 R2 0 0 0 0 −R2 R3 R4 0
[A]= 0 0 0 R4 −R5
1 −1 −1 0 0 0 0 1 −1 −1
x1 and x2 gives the same results
V1 {C}= V0
02 0
Two different MATLAB commands to solve
[A] {X} = {C}
MATLAB output
10
End of Example L05.1
Example L05.2
Use method of joints to analyse the truss below and determine the forces in its bars and the reaction forces applied due to the fixed and roller supports.
A
7m
B
10N
24m
25m
C
0 378> 9> 8> 10 9> 0 7>< >= >< 0 >=
26 (7/25) (24/25) 0
6 (24/25) (7/25) 0
6 0 (24/25) 1 0 0 07 =0
Tension and compression forces inside beams
NAB NAC NBC
RCy RBx RBy
0 0 0 0
64 0 (7/25) 0 1 0 0 75> > > 0 > (7/25) 0 1 0 1 0 >: >; >: 0 >;
(24/25) 0 0001 0 12
Reaction forces at support
26 (7/25) (24/25) 0 0 0 0 37 6 (24/25) (7/25) 0 0 0 0 7 6 0 (24/25) 1 0 0 07 640 (7/25)010075
(7/25) 0 1010 (24/25) 0 0 001
8> NAB 9> 8> 9> >
NBC =
> RCy > > > >:RBx>;>: >;
RBy
10 0 0 0 0 0
Lecture05M02.m
clear all
A=[-7/25 24/25 0 0 0 0;
-24/25 -7/25 0 0 0 0;
0
0 7/25
-24/25 -1000;
C=[-10; 0;
0; 0; 0; 0];
24/25 0
7/25 0
0 1 0 0; 1 0 1 0; 0001];
res=A\C;
NAB=res(1)
NAC=res(2)
NBC=res(3)
RCy=res(4)
RBx=res(5)
RBy=res(6)
NAB = 2.8000 NAC = -9.6000 NBC = 9.2160 RCy = 2.6880 RBx = -10
RBy = -2.6880
MATLAB command to solve
[A] {X} = {C}
13
End of Example L05.2
How do we solve all
these
problems?
Mathematical equations that govern fluid flow
∂ρ + ∂(ρu) + ∂(ρv) + ∂(ρw) = 0 ∂t ∂x ∂y ∂z
∂(ρu)+∂(ρu2)+∂p−∂τxx +∂(ρuv)−∂τxy +∂(ρuw)−∂τxz =0 ∂t ∂x ∂x ∂x ∂y ∂y ∂z ∂z
∂(ρv)+∂(ρuv)−∂τxy +∂(ρv2)+∂p−∂τyy +∂(ρvw)−∂τyz =0 ∂t ∂x ∂x ∂y ∂y ∂y ∂z ∂z
∂(ρw) ∂(ρuw) ∂τ ∂(ρvw) ∂τ ∂(ρw2) ∂
∂t
+ − xz+ − yz+ + Navier-Stokes Equations
∂x ∂x ∂y ∂y ∂z ∂
+∂(Et +p)u−∂(uτxx)−∂(vτxy)−∂(wτxz)+∂(qx) ∂t ∂x ∂x ∂x ∂x ∂x
+∂(Et +p)v−∂(uτxy)−∂(vτyy)−∂(wτyz)+∂(qy) ∂y ∂y ∂y ∂y ∂y
+∂(Et +p)w−∂(uτxz)−∂(vτyz)−∂(wτzz)+∂(qz)=0 ∂z ∂z ∂z ∂z ∂z
p−∂τzz =0 z ∂z
∂(Et )
Conservation of mass Conservation of momentum Conservation of energy
Navier-Stokes Equations
∂ρ + ∂(ρu) + ∂(ρv) + ∂(ρw) = 0 ∂t ∂x ∂y ∂z
∂(ρu)+∂(ρu2)+∂p−∂τxx +∂(ρuv)−∂τxy +∂(ρuw)−∂τxz =0 ∂t ∂x ∂x ∂x ∂y ∂y ∂z ∂z
∂(ρv)+∂(ρuv)−∂τxy +∂(ρv2)+∂p−∂τyy +∂(ρvw)−∂τyz =0 ∂t ∂x ∂x ∂y ∂y ∂y ∂z ∂z
∂(ρw)+∂(ρuw)−∂τxz +∂(ρvw)−∂τyz +∂(ρw2)+∂p−∂τzz =0 ∂t ∂x ∂x ∂y ∂y ∂z ∂z ∂z
∂(Et)+∂(Et +p)u−∂(uτxx)−∂(vτxy)−∂(wτxz)+∂(qx) ∂t ∂x ∂x ∂x ∂x ∂x
+∂(Et +p)v−∂(uτxy)−∂(vτyy)−∂(wτyz)+∂(qy) ∂y ∂y ∂y ∂y ∂y
+∂(Et +p)w−∂(uτxz)−∂(vτyz)−∂(wτzz)+∂(qz)=0 ∂z ∂z ∂z ∂z ∂z
Navier-Stokes Equations
∂ρ + ∂(ρu) + ∂(ρv) + ∂(ρw) = 0 ∂t ∂x ∂y ∂z
∂(ρu)+∂(ρu2)+∂p−∂τxx +∂(ρuv)−∂τxy +∂(ρuw)−∂τxz =0 ∂t ∂x ∂x ∂x ∂y ∂y ∂z ∂z
∂(ρv)+∂(ρuv)−∂τxy +∂(ρv2)+∂p−∂τyy +∂(ρvw)−∂τyz =0 ∂t ∂x ∂x ∂y ∂y ∂y ∂z ∂z
∂(ρw)+∂(ρuw)−∂τxz +∂(ρvw)−∂τyz +∂(ρw2)+∂p−∂τzz =0 ∂t ∂x ∂x ∂y ∂y ∂z ∂z ∂z
∂(Et)+∂(Et +p)u−∂(uτxx)−∂(vτxy)−∂(wτxz)+∂(qx) ∂t ∂x ∂x ∂x ∂x ∂x
+∂(Et +p)v−∂(uτxy)−∂(vτyy)−∂(wτyz)+∂(qy) ∂y ∂y ∂y ∂y ∂y
+∂(Et +p)w−∂(uτxz)−∂(vτyz)−∂(wτzz)+∂(qz)=0 ∂z ∂z ∂z ∂z ∂z
http://www.news.navy.mil/view_single.asp?id=9813
http://commons.wikimedia.org/wiki/McLaren
For real life engineering applications, how
many grid points do you need to solve?
http://commons.wikimedia.org/wiki/Qantas
~1016
L5.2:
Solving
Systems of linear algebraic equations
Systems of Linear Algebraic Equations
In the two examples above, we look at some examples of systems of algebraic equations that represents real engineering problems. These problems are relatively small. 5 unknowns For the problem on the Example L05.1 and 6 unknowns for L05.2.
R1 R2 0 0 0 i1 V1 i
2 (7/25) 6 (24/25)
0 38> NAB 0 7> NAC
9> 8> 10 9> >>0>
(24/25) 0 (7/25) 0 (24/25) 1
6
6 0 (7/25) 0
a11x1+a12x2+a13x3+……+a1nxn =c1 a21x1+a22x2+.a23x3+……+a2nxn =c2
an1x1+an2x2+.an3x3+……+annxn =cn
where all the a’s and c’s are constants. We want find all the xi’s such that all the above equations are simultaneously satisfied.
0 −R2 R3 R4 0 i2 = V0 0 0 0 R −R 3 2
L05.1
6 0
0 0
0 0
0 0
1 0
0 1
0 0
07>
L05.2
45 i
1 −1 −1 0 0 4 0
4 (7/25) 0 1 (24/25) 0 0
07>RCy 05>RBx
> > 0 > >>0> >; >: 0 >;
1 >: RBy
For many engineering problems, you need to deal with systems of equations that have n unknowns. Sometimes, n
0
0 0 1 −1 −1 i5
can be very large. In general a system of n linear algebraic equations can be written as
24
where
Systems of Linear Algebraic Equations
Sometimes, the matrix notation will be used
a11 a12 . . a1n [] a21 a22 . . a2n
a11x1+a12x2+a13x3+……+a1nxn =c1 a21x1+a22x2+a23x3+……+a2nxn =c2
. an1x1+an2x2+an3x3+……+annxn =cn
[A] {X} = {C}
x1 c1 {X}= x2 , {C}= c2
(5.1)
A = a31 a32 . . a3n ,
….. x. .
x3
an1 an2 . . ann n cn
In the following slides, sometimes to save save space, and for brevity, I might augment the [A] and {C} matrices
25
c3
a11x1+a12x2+a13x3+……+a1nxn =c1 a21x1+a22x2+a23x3+……+a2nxn =c2
. an1x1+an2x2+an3x3+……+annxn =cn
a11 a12 . . a1n x1 c1
a21 a22 . . a2n x2 = c2
a31 a32 . . a3n x3 c3 …….
an1 an2 . . ann xn cn
a11 a12 a13 . . . . a1n | c1 a21 a22 a23 . . . . a2n | c2 a31 a32 a33 . . . . a3n | c3
……..|. vv
……..|. ……..|. ……..|.
an1 an2 an3 . . . . ann | cn 26
Systems of Linear Algebraic Equations
Equation (5.1) can be solved by inverting the matrix [A] and then
multiplying [C] by[A]-1.
{X} = [A]−1 {C}
[A] {X} = {C}
But this might be difficult and time consuming, sometimes
impossible to do if matrix [A] is large. Can we solve {X} numerically?
Turns out that you can….and sometimes without needing to find [A]-1. This is the topic for upcoming lectures.
27
The above set of equations is often very difficult to solve, especially when n is large. For practical engineering problems, sometimes n can be 109 (or more!)
Before we go on and solve the general problem, it is appropriate now to point out that there are some linear problems that are easy to solve…..even when n is very large.
Let’s now consider the case when [A] is diagonal, upper triangular and lower triangular.
Upper triangular
Lower triangular
a 0 0 0 . . . 0 11 0011 0
Diagonal
a 0 0 . . . . 0
a11 a12 a13 a14 . . . a1n
0 a22 0 . . . . 0
0 a22 a23 a24 . . . a2n
a21 a22 0 0 . . . 0 a a a 0 . . . 0
0 0 a33 a34 . . . a3n A=…….. A=000a…a A=
0 0 a . . . . 0
[] 033 31 32 33
a41 a42 a43 a44 . . . 0 …….. ……..
. . . . . . . . 44
[] []
4n …….. ……..
…….0 …….. …….. 0 ……..
00…0a 000….a aaa….a
nn
nn n1 n2 n3 nn
28
L5.3:
Systems of linear algebraic equations that are easy to solve
Diagonal systems
Consider first the case when [A] is diagonal. It is easy to solve this system of equations.
I will illustrate using a small 3×3 system.
100|3 020|4 003|9
a11 0 0 | c1 0 a22 0 | c2 0 0 a33 | c3
First row of matrix
1×1 = 3
x1 = 3/1 = 3
a11x1 = c1
x1 = c1/a11
2nd row of matrix
2×2 = 4
x2 = 4/2 = 2
a22x2 = c2
x2 = c2/a22
3rd row of matrix
3×3 = 9
x3 =9/3=3
a x = c 33 3 3
31 x3 =c3/a33
If you have a 3×3 system, the solution process is
a11 0 0 | c1
0 a22 0 | c2
0 0 a | c 33 3
1 1 11 a22x2 = c2
x2 = c2/a22
a33x3 = c3
x3 = c3/a33
So in general, if you have a nxn system,
a11 0 0 . . . . 0 | c1 0 a22 0 . . . . 0 | c2
0 0 a33 . . . . 0 | c3 ……..|.
. 0. . . . . . . | .
. . . . . . . . | . x=c/a …….0|.
a x = c 11 1 1
32
0
0 0 0 . . . 0 ann | cn The solution for xi is
(5.2)
aiixi = ci
xi = ci/aii
Upper triangular systems
Now let us consider the case when [A] is upper triangular.
1 2 3 | 14
0 −4 −7 | −26 0 0 −9/4 | −18/4
a11 a12 a13 | c1
0 a22 a23 | c2 0 0 a33 | c3
This system is easy to solve providing you know where to start. Start with bottom the move up. Starting with the last row gives
− 9 x3 = − 18 44
x3 = 2/1 = 2
a33x3 = c3
x3 = c3/a33
Next, look at 2nd last row
−4×2 − 7×3 = − 26
x2 = −26 + 7 ⋅ 2 = 3
−4
a22x2 + a23x3 = c2
x2 = c2 − a23x3
a22
Next, look at 1st row
x1 + 2×2 + 3×3 = 14
x1 = 14 − 2 ⋅ 3 − 3 ⋅ 2 = 2
1
a11x1 + a12x2 + a13x3 = c1
x1 = c1 − a12x2 − a13x3
a11
34
If you have a 3×3 system, the solution process is
a11 a12 a13 | c1 0 a22 a23 | c2
0 0 a33 | c3
This system is easy to solve providing you know
where to start. Starting with the last row gives
a33x3 = c3
x3 = c3/a33
Next, look at 2nd last row
a22x2 + a23x3 = c2 Next, look at 1st row
a11x1 + a12x2 + a13x3 = c1
x1 = c1 − a12x2 − a13x3
In general, when you have a nxn upper triangular system
a11 a12 a13 a14 . . . a1n | c1 0 a22 a23 a24 . . . a2n | c2
0 0 a33 a34 . . . a3n | c3 0 0 0 a44 . . . a4n | c4 . . . . . . . . | .
……..|.
x2 = c2 − a23x3 a22
. . . . 0 0 0 .
The solution for xi is
. . . .
. . . ann
| . | cn
xn =
cn ann
and
for i=n−1,n−2,n−3,..,2,1
xi =
c i − ∑n a i j x j j=i+1
aii
a11
(5.3)
35
Example L05.3:Write a MATLAB function to implement Eq.(5.3) for upper triangular systems and use it
to solve .
x1 x2 x3 x4
1235 0 8 −7 3 0 0 −9 7 0006
14 −26 12 3
=
Lecture05M03.
A=[1 2 3 5; 0 8 -7 3;
0 0 -9 7 0 0 0 6];
C=[14;-26;12;3]
x=UpperTriangularSolver(A,C)
%Proof that you have the right answer by calculating Ax-C=0, which should be %very close to zero if Ax=C.
A*x-C
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
(5.3)
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
36
xn = cn ann
ci − and xi=
n
aii
for i=n−1,n−2,n−3,..,2,1
∑ aijxj j=i+1
a11 a12 a13 a14 . . . a1n
0 a22 a23 a24 . . . a2n
0 0 a33 a34 . . . a3n
[A] =
……..
0 0 0 a44 . . . a4n
0
…….. ……..
0 0 0 . . . . ann
c1 c2 c3
x1 x2 x3
{X} =
n
{C} =
n
x.
c.
37
End of Example L05.3
Lower triangular systems
Now let us consider the case when [A] is lower triagular.
5 0 0 | 12 −2−40| 8
7 2 3 | −4
a11 0 0 | c1
a21 a22 0 | c2 a31 a32 a33 | c3
Start from the top and go down. The first row
gives us
5×1 = 12
x1 = 12/5 = 2.4
a11x1 = c1
x1 = c1 a11
Next, look at 2nd last row
−2×1 − 4×2 = 8
x2 = 8+2⋅(2.4) =3.2 −4
a21x1 + a22x2 = c2
x2 = c2 − a21x1
a22
Next, look at last row
7×1 + 2×2 + 3×3 = − 4
x3 = −4−2⋅(−3.2)−7⋅(2.4) =−4.8
3
a31x1 + a32x2 + a33x3 = c3
x3 = c3 − a32x2 − a31x1
a33
If you have a 3×3 system, the solution process is
a11 0 0 | c1 a21 a22 0 | c2
a31 a32 a33 | c3 The first row gives us
a11 Next, look at 2nd last row
a21x1 + a22x2 = c2
Next, look at last row
a31x1 + a32x2 + a33x3 = c3
a11x1 = c1 x1 = c1
x2 = c2 − a21x1 a22
x3 = c3 − a32x2 − a31x1 a33
In general, when you have a n x n lower triangular system
a11 0 0 0 . . . 0 | c1 a21 a22 0 0 . . . 0 | c2 a31 a32 a33 0 . . . 0 | c3 a41 a42 a43 a44 . . . 0 | c4
……..|. ……..|. ……..|.
an1 an2 an3 . . . . ann | cn The solution for xi is
x1 = c1/a11
( i−1 ) ci − ∑ aijxj
j=1 for aii
xi =
i = 2,3,4,…..,n (5.4)
41
Example L05.4:Write a MATLAB function to implement Eq.(5.4) for lower triangular system and use it
to solve
x1 x2 x3 x4
5000 −2 −4 0 0 7230 6897
12 8 −4 1
=
Lecture05M04.m
clear all
A=[5 0 0 0; -2 -4 0 0;
7 2 3 0; 6 8 9 7];
C=[12;8;-4;1]
x(1)=C(1)/A(1,1);
for i=2:n
sum=sum+A(i,j)*x(j);
x=LowerTriangularSolver(A,C)
%Check that you have the right answer by calculating Ax-C=0, which should be %very close to zero if Ax=C.
A*x-C
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
sum=0; for j=1:i-1
end
x(i)=(C(i)-sum)/A(i,i);
end end
xi =
(
i−1 ) ci − ∑ aij xj
j=1 for aii
42
(5.4)
x1 = c1/a11
i = 2,3,4,…..,n
[A] =
..0
a000…0 a11 a 0 0 . . . 0
a21 a22 a 0 .
a31 a32 a33 a . . . 0
0
41 42 43 44 …….. …….. ……..
.
c3 n c.
an1 an2 an3
. . . ann
x1 {X}= x2
c1 c2
x3 x.
{C} =
n
43
End of Example L05.4
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.
45
L5.4:
Solving linear systems of algebraic equations Direct Methods: Gauss Elimination Section 3.1.1 of “Book”
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
47
Systems of Linear Algebraic Equations
We want to solve where
[A] {X} = {C}
a11 a12 . . a1n [A]= a21 a22 . . a2n ,
x1 {X}= x2 ,
c1 {C}= c2
a31 a32 . . a3n In augmented form
a11 a12 a13 . . . . a1n | c1 a21 a22 a23 . . . . a2n | c2 a31 a32 a33 . . . . a3n | c3
……..|. ……..|. ……..|. ……..|.
an1 an2 an3 . . . . ann | cn 48
c3 an1 an2 . . ann n cn
x3
….. x. .
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′
Easy to solve See Example L05.3
Gauss Elimination
49
| c4′ . 0. . . . . . . | . ……..|. ……..|.
0 0 0 . . . . an′n | cn′
Illustrate Gauss Elimination with a simple 3×3 example
Suppose you want to solve:
[1 2 3] x1 {14} 3 2 2 x2 = 16 232×3 17
[ ] [1 2 3] {14} A = 3 2 2 {C}= 16 .
232 17
First augment [A] with {C}
123| 322| 232|
14 16 17
50
a11 a12 a13 x1 c1 a21 a22 a23 x2 = c2
a31 a32 a33 x3 [] a11a12a13
c3
.
A = a21 a22 a23 {C}= a31 332 a33
c1 c2 c3
a(1) a(1) a(1) | 11 12 13
a(1) a(1) a(1) | 21 22 23
a(1) a(1) a(1) | 31 32 33
a(1) 1,4
a (1) 2,4
a (1) 3,4
• Eliminate x1 from second and third row.
• Eliminate x1 from second and third row.
• To eliminate x1 from 2nd row, multiply row 1 by 3/1 and subtract from row 2.
• To eliminate x1 from 3rd row, multiply row 1 by 2 /1 and subtract from row 3.
• To eliminate x1 from 2nd row, multiply row 1 by a21(1) /a11(1) and subtract from row 2.
• To eliminate x1 from 3rd row, multiply row 1 by a31(1) /a11(1) and subtract from row 3.
0 0
1 2 3 | 14 322|16 2 3 2 | 17
a(1) a(1) a(1) | a(1) 11 12 13 14
a(1) a(1) a(1) |a(1) 21 22 23 24
a(1) a(1) a(1) | a(1) 31 32 33 34
1 2
3 | 14
a(1)
11 12 13 14
21
a(1) 11 22
11 11 11 11
a(2) a(2) a(2) | a(2) 11 12 13 14
a(2) a(2) | a(2) 22 23 24
a(2) a(2) | a(2) 32 33 34
a(1) a(1) | a(1) 21 21 21 21
3−311 2−312 2−313 | 16−3114
2−211 3−212 2−213 | 17−2114
a(1) − a(1) a(1) a(1) − a(1) a(1) a(1) − a(1) a(1) | a(1) − a(1) a(1)
a(1) 12 23
a(1) 13 24 a(1) 14
31 31 31 31
a(1) − a(1) a(1) a(1) − a(1) a(1) a(1) − a(1) a(1) | a(1) − a(1) a(1)
31 a(1) 11 32
11 11 11 11
1 2 3 | 14 −4 −7 | −26 −1 −4 | −11
0 0
a(1) 12 33
a(1) 13 34 a(1) 14
51
• Next, eliminate x2 from third row.
• To eliminate x2 from the 3rd row, multiply row 2 by -1/-4 and subtract from row 3.
1 2 3 | 14 0 −4 −7 | −26 0 −1 −4 | −11
1 2 3 | 14 0 −4 −7 | −26
0 −1−−14 −4−−1(−7) | −11−−1(−26) −4 −4 −4
1 2 3 | 14 −4 −7 | −26
0 −9/4 | −18/4
0 0
• Next, eliminate x2 from third row.
• To eliminate x2 from the 3rd row, multiply row 2 by a32(2)/ a22(2) and subtract from row 3.
a(2) a(2) a(2) | a(1) 11 12 13 14
0 a(2) a(2) | a(2) 22 23 24
0 a(2) a(2) | a(2) 32 33 34
a(2) a(2) a(2) | a(2) 11 12 13 14
0
a(2) a(2) | a(2) 22 23 24
a(3) a(3) a(3) | a(3) 11 12 13 14
a(3) a(3) | a(3) 22 23 24
0 a(3) | a(3) 33 34
32 32 32
0 a(2) − a(2) a(2) a(2) − a(2) a(2) | a(2) − a(2) a(2)
32 a(2) 22 33
22 22 22
0 0
a(2) 23 34 a(2) 24
This is now
an upper triangular system. Easy to solve.
Can use function from Example L05.3
52
3×3 system
a(1) a(1) a(1) | a(1) 11 12 13 1,4
a(1) a(1) a(1) | a(1) 21 22 23 2,4
a(1) a(1) a(1) | a(1) 31 32 33 3,4
• Eliminate x1 from second and third row.
• To eliminate x1 from 2nd row, multiply row 1
by a21(1) /a11(1) and subtract from row 2.
• To eliminate x1 from 3rd row, multiply row 1
by a31(1) /a11(1) and subtract from row 3.
a(2) a(2) a(2) | a(1) 11 12 13 14
a(2) a(2) | a(2) 22 23 24
a(2) a(2) | a(2) 32 33 34
0 0
nxn system
a(1) a(1) a(1) . . . . a(1) | a(1)
11 12 13 1n 1,n+1 a(1) a(1) a(1) . . . . a(1) | a(1)
21 22 23 2n 2,n+1 a(1) a(1) a(1) . . . . a(1) | a(1)
31 32 33 3n 3,n+1
……..|.
……..|.
……..|.
……..|. a(1) a(1) a(1) . . . . a(1) | a(1)
n1 n2 n3 nn n,n+1
Forward elimination to get matrix into solvable form. Eliminate x1 from second and subsequent equations.
To eliminate x1 from ith row, multiply row 1 by ai1(1)/ a11(1) and subtract from row i.
a(1)
a(2)=a(1)−(ai1 )a(1) for i=2,3..,n and j=1,2,3..,n+1
53
ij ij (1) 1j 11
a(2) a(2) a(2) . . . . a(2) | a(2)
11
12 13 1n 1,n+1 a(2) a(2) . . . . a(2) | a(2)
0
0
. . .
0
22 23 2n 2,n+1 a(2) a(2) . . . . a(2) | a(2)
32 33 3n 3,n+1
…….|.
…….|.
…….|.
a(2) a(2) . . . . a(2) | a(2) n2 n3 nn n,n+1
3×3 system
a(2) a(2) a(2) | a(1) 11 12 13 14
0 a(2) a(2) | a(2) 22 23 24
0 a(2) a(2) | a(2) 32 33 34
• Next, eliminate x2 from third row.
• To eliminate x2 from the 3rd row, multiply row 2
by a32(2)/ a22(2) and subtract from row 3.
a(3) a(3) a(3) | a(3) 11 12 13 14
a(3) a(3) | a(3) 22 23 24
0 a(3) | a(3) 33 34
0 0
ij ij (2) 2j 22
nxn system
a(2) a(2) a(2) . . . . a(2) | a(2) 11 12 13 1n 1,n+1
0 a(2) a(2) . . . . a(2) | a(2) 22 23 2n 2,n+1
0 a(2) a(2) . . . . a(2) | a(2) 32 33 3n 3,n+1
……..|.
……..|.
……..|. 0 a(2) a(2) . . . . a(2) | a(2)
n2 n3 nn n,n+1
Next, eliminate x2 from third and subsequent equations. Therefore, to eliminate x2 from the ith row, multiply row 2 by ai2(2)/a22(2) and subtract from row i.
a(2)
a(3)=a(2)−(ai2 )a(2) for i=3..,n and j=2,3..,n+1
a(3) a(3) a(3) . . . . a(3) | a(3)
11
12 13 1n 1,n+1 a(3) a(3) . . . . a(3) | a(3)
22 23 2n 2,n+1 a(3) . . . . a(3) | a(3)
54
0
0
. . .
0
0
. . .
0
33 3n 3,n+1
……|.
……|.
……|.
a(2) . . . . a(2) | a(2) n3 nn n,n+1
a(3) a(3) a(3) . . . . a(3) | a(3)
11
12 13 1n 1,n+1 a(3) a(3) . . . . a(3) | a(3)
0
0
. . .
0
22 23 2n 2,n+1 a(3) . . . . a(3) | a(3)
0
. . .
0
33 3n 3,n+1 ……|. ……|.
……|. a(2) . . . . a(2) | a(2)
n3
Keep
nn n,n+1 doing n-1
times
a(n) a(n) a(n) a(n) . . . 11 12 13 14
a(n) a(n) a(n) . . . 22 23 24
a(n) | a(n) 1n 1,n+1
a(n) | a(n) 2n 2,n+1
a(n) | a(n) 3n 3,n+1
a(n) | a(n) 4n 4,n+1
0 0 0 0
. .
0
0
0
0
. .
0
a(n) a(n) . . .
33
34
0
0
. .
0
a(n)… 44
0
. .
.
a(n) . . . | . 55
….|. ….|.
. . . a(n) | a(n) nn n,n+1
(5.5)
(a(k)) a(k+1) = a(k) − i,k a(k)
ij
ij a(k) k,j k,k
for k=1..,n−1, i=k+1..,n
and j=k,..,n+1
55
Example L05.5:Write a MATLAB function that uses Gauss Elimination to transform the following system to upper triangular form. Use the MATLAB backslash operator \ and the linsolve() function to convince yourself that the transformed upper triangular form and the original system gives the same answer.
x1 x2 x3 x4
1234 3489 8 20 4 3 5678
348 10
=
Lecture05M05.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);
k=1
1234|3 −2 −1 −3 | −5
0 0 0
[newA] =
−4 −8 −12 | −5
4 −20 −29 | −16
[n,n]=size(A);
Aprime=newA(1:n,1:n) Cprime=newA(:,n+1)
x1=A\C x2=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
1 2 3 4 3
[A]= 3 4 8 9 {C}= 4 8 20 4 3 8
5 6 7 8 10
[newA] =
1 2 3 4 | 3 3 4 8 9 | 4 8 20 4 3 | 8 5 6 7 8 | 10
for k=1..,n−1, i=k+1..,n
and j = k, . . , n + 1
a(k+1) = a(k) −
( a(k) )
i,k a(k)
a(k) k,j k,k
56
ij
ij
k=2
1234|3 −2 −1 −3 | −5
0 0 0
[newA] =
−6 −6 | 5
0 0
−22 −35 | −26
k=3
[newA] =
1234|3 −2 −1 −3 | −5
−22 −35 | −26 0 3.5455 | 12.1
0 0 0
0 0
Example L05.5:Write a MATLAB function that uses Gauss Elimination to transform the following system to upper triangular form. Use the MATLAB backslash operator and the linsolve() function to convince yourself that the transformed upper triangular form and the original system gives the same answer.
x1 x2 x3 x4
1234 3489 8 20 4 3 5678
348 10
=
Lecture05M05.m
A=[1 2 3 4; 3 4 8 9;
8 20 4 3; 5 6 7 8];
C=[3;4;8;10];
Aprime=newA(1:n,1:n)
k=3
[newA] =
1234|3 −2 −1 −3 | −5
−22 −35 | −26 0 3.5455 | 12.1
0 0 0
0 0
newA=GaussElimination(A,C); [n,n]=size(A);
Cprime=newA(:,n+1)
x1=A\C x2=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
1234 3
[Aprime] =
0 −2 −1 −3 0 0 −22 −35
{Cprime} = −5 −26
12.1
Should give the same answer!
0 0
0
3.5455
Using MATLAB \ operator to solve for original system A and C
Using MATLAB \ operator to solve for transformed system Aprime and Cprime
a(k+1) = a(k) −
( a(k) )
i,k a(k)
a(k) k,j k,k
57
ij
ij
End of Example L05.5
****Try solving using Gauss elimination
2×2 + 3×3 = 1 3×1 + 5×2 + 6×3 = 2 9×1 + 2×2 + 3×3 = 3
59
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
60
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).
61
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
63
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
64
End of Example L05.6