Numerical Methods in Engineering (ENGR20005)
Lecture 04 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L4.1: MATLAB: functions
MATLAB functions are elements of MATLAB program
that accepts input, xi, perform some calculations and produce outputs, yi.
x1 y1 x2 y2 x3 . .y3
MATLAB function
xN
yN
.. ..
The structure of MATLAB functions looks something like below []=NameOfFunction(x1,x2,x3,….,xN)
functiony1,y2,y3,…,yN Command A
Command B
Command C
Command D
Command E
end
Example L04.1: Write a MATLAB function that takes in two values, a and b, multiply them together and return their product
Lecture04M01.m
one output product
two inputs a and b
c=mymult(1,2) d=mymult(4,5)
function
product
product = a*b ;
=mymult(a,b)
MATLAB function mymult()
end
a=1 b=2
product=1*2=2
c=product=2
a=4 b=5
product=4*5=20
d=product=20
End of Example L04.1
Example L04.2: Convert your Lecture03M07.m into a MATLAB function that takes as an input different guesses and
outputs the root of the function
f(c) = 50 (1 − e− 9c ) − 10 = 0 c5
Lecture04M02.m
root1=SecantMethod(4) root2=SecantMethod(2)
function cr=SecantMethod(ci) cim1=ci-1 fi=(50./ci)*(1-exp(-9*ci/5))-10; while abs(fi)>1.0e-6
fi=(50./ci)*(1-exp(-9*ci/5))-10; fim1=(50./cim1)*(1-exp(-9*cim1/5))-10; cip1=ci-(ci-cim1)*fi/(fi-fim1)
cim1=ci;
ci=cip1; end
cr=cip1; end
I need a 2nd guess for ci. So I have arbitrarily said that that 2nd guess, ci-1=ci-1.
Lecture03M07.m
ci=4 cim1=3
fi=(50./ci)*(1-exp(-9*ci/5))-10; while abs(fi)>1.0e-6
fi=(50./ci)*(1-exp(-9*ci/5))-10; fim1=(50./cim1)*(1-exp(-9*cim1/5))-10; cip1=ci-(ci-cim1)*fi/(fi-fim1)
cim1=ci;
ci=cip1; end
End of Example L04.2
Example L04.3: Convert your Lecture04M02.m into a MATLAB function that takes as an input different guesses and different functions and outputs the root of any function, f(c). Test that your program works by by find the roots of the functions 50 ( 9c ) x
f(c)= c 1−e−5 −10=0 g(x)=e −10x=0
Lecture04M03.m
f=@(c) (50/c)*(1-exp(-9*c/5))-10;
g=@(x) exp(x)-10*x;
root1=SecantMethod(f,4) root2=SecantMethod(g,5)
function cr=SecantMethod(func,ci) cim1=ci-1
fi=func(ci);
while abs(fi)>1.0e-6
fi=func(ci); fim1=func(cim1); cip1=ci-(ci-cim1)*fi/(fi-fim1); cim1=ci;
ci=cip1;
end
cr=cip1;
end
Lecture04M02.m
root1=SecantMethod(4) root2=SecantMethod(2)
function cr=SecantMethod(ci) cim1=ci-1 fi=(50./ci)*(1-exp(-9*ci/5))-10; while abs(fi)>1.0e-6
fi=(50./ci)*(1-exp(-9*ci/5))-10; fim1=(50./cim1)*(1-exp(-9*cim1/5))-10; cip1=ci-(ci-cim1)*fi/(fi-fim1)
cim1=ci;
ci=cip1; end
cr=cip1; end
In class example:
Find the value of x in the domain -1
x=x*5
else
x=20
end
x=1
if x<4
x=x+1 end
x=5
if x<4
x=x+1 else
x=10 end
x=5
Final value of x
x=10
x=7
x=x*5 x=7*5 x=35
Final value of x
End of Example L04.4
L4.3
Root finding: Bracketing Method
Root Finding Overview
Numerical algorithms for finding roots of equations
Bracketing methods
Method of false position
17
One point iteration
Open methods
Bisection method
Newton-Raphson
Secant method
Bracketing Method
There are two different Bracketing Methods. They are:-
• •
These methods are known as bracketing methods because they rely on having two initial guesses.
- xl - lower bound and
- xu - upper bound.
*The guesses xl and xu must bracket (be either side of) the root.
This requirement in many cases may be hard to meet.
Bisection method Method of False position
f(x)
f(xl)
xl f(xu)
xR
x
xu
f(x)
Bracketing Method
• If f(xu) and f(xl) are of different sign, then there must be at least
xu xl xR
x
i.e.xl
fxl=exp(xl)-10*xl; fxu=exp(xu)-10*xu; fxr=exp(xr)-10*xr;
if fxr*fxl < 0 xu= xr;
elseif fxr*fxu < 0 xl= xr;
end
xr=(xl+xu)/2.0 end
Bisection method
32
Bisection method
xr =
xl + xu 2
iter
xl
xu
xr = xl + xu
x2 r
1
3.0
4.0
3.5
2
3.5
4.0
3.75
3
3.5
3.75
3.625
4
3.5
3.625
3.5625
33
20
15
10
5
0
-5
-10
-15
f(x) = ex
− 10x
v
v
v
xu
v
xl
xr
3 3.2 3.4 3.6 3.8 4 x
f(x)
Bisection method
xr =
xl + xu 2
iter
xl
xu
xr = xl + xu
x2 r
1
3.0
4.0
3.5
2
3.5
4.0
3.75
3
3.5
3.75
3.625
4
3.5
3.625
3.5625
34
20
15
10
5
0
-5
-10
-15
f(x) = ex
− 10x
xl
v
v
xr
x
u
3 3.2 3.4 3.6 3.8 4 x
f(x)
Bisection method
xr =
xl + xu 2
iter
xl
xu
xr = xl + xu
x2 r
1
3.0
4.0
3.5
2
3.5
4.0
3.75
3
3.5
3.75
3.625
4
3.5
3.625
3.5625
35
20
15
10
5
0
-5
-10
-15
f(x) = ex
− 10x
xl
v
xr xu
3 3.2 3.4 3.6 3.8 4 x
f(x)
v
Bisection method
xr =
xl + xu 2
iter
xl
xu
xr = xl + xu
x2 r
1
3.0
4.0
3.5
2
3.5
4.0
3.75
3
3.5
3.75
3.625
4
3.5
3.625
3.5625
36
20
15
10
5
0
-5
-10
-15
f(x) = ex
− 10x
xl
xv r
xu
3 3.2 3.4 3.6 3.8 4 x
f(x)
v v
End of Example L04.5
Example L04.6: Write a MATLAB program using the Bisection method to find the roots for
f(c)= 50(1−e−9c)−10 c5
Lecture04M05.m Lecture04M06.m
f (c)
xl=3.0;
xu=4.0; xr=(xl+xu)/2.0; fxr=exp(xr)-10*xr;
while abs(fxr)>1.0e-6
fxl=exp(xl)-10*xl; fxu=exp(xu)-10*xu; fxr=exp(xr)-10*xr;
if fxr*fxl < 0 xu= xr;
elseif fxr*fxu < 0 xl= xr;
end
xr=(xl+xu)/2.0 end
cl=4.0;
cu=6.0;
cr=(cl+cu)/2.0; fcr=(50/cr)*(1-exp(-9*cr/5))-10;
while abs(fcr)>1.0e-6
fcl=(50/cl)*(1-exp(-9*cl/5))-10; fcu=(50/cu)*(1-exp(-9*cu/5))-10; fcr=(50/cr)*(1-exp(-9*cr/5))-10;
if fcr*fcl < 0 cu= cr;
elseif fcr*fcu < 0 cl= cr;
end
cr=(cl+cu)/2.0 end
38
End of Example L04.6
L4.5:
Bracketing Method: Method of False Position (see Section 2.2.2 of “Book”)
Method of False Position
The method of false position is based on the following observations.
If the value of f(xl) is closer to zero than f(xu), then the root is likely to be closer to xl than xu.
See the figure below.
f(x)
f(xu)
xl xr
f(xl)
xR
xu
Using similar triangles, we can obtain the following expression
=
− f(xl) xr − xl
f(xu) xu − xr
If the value of f(xu) is closer to zero than f(xl),
then the root is likely to be closer to xu than xl.
See the figure below.
f(x) f(xl)
xr xu
f(xu)
xl xR
Using similar triangles, we can obtain the following expression
=
f(xl)
xr − xl
− f(xu) xu − xr
These two expressions
are the same!
Method of False Position
− f(xl) = f(xu) xr − xl xu − xr
You can show that the above expression reduces to
or
xr =
xr = xu − f(xu)(xl − xu)
xu f(xl) − xl f(xu) f(xl) − f(xu)
f(xl) − f(xu) 42
(4.6)
Method of False Position
So, to use the method of False Position to find roots, use Steps 1-3 of the Bisection method but replace Step 2 with the formula given in Eq. (4.6) and reproduced below.
xr = xu −
f(xu)(xl − xu) f(xl) − f(xu)
43
Method of False Position
Example L04.7: Write a MATLAB program using the method of False Position to find the roots for
f(x) = ex − 10x = 0 You might like to start with MATLAB program Lecture04M05.m
Lecture04M05.m
xl=3.0 xu=4.0
fxr=exp(xr)-10*xr
while abs(fxr)>1.0e-6
fxl=exp(xl)-10*xl; fxu=exp(xu)-10*xu; fxr=exp(xr)-10*xr;
Lecture04M07.m
xl=3.0;
xu=6.0; fxl=exp(xl)-10*xl; fxu=exp(xu)-10*xu;
fxr=exp(xr)-10*xr; while abs(fxr)>1.0e-6
fxl=exp(xl)-10*xl; fxu=exp(xu)-10*xu; fxr=exp(xr)-10*xr;
xr=(xl+xu)/2.0
if fxr*fxl < 0 xu= xr;
elseif fxr*fxu < 0 xl= xr;
end
xr=(xl+xu)/2.0
xr = xu − f(xu)(xl − xu) f(xl) − f(xu)
xr=xu-fxu*(xl-xu)/(fxl-fxu);
if fxr*fxl < 0 xu= xr;
elseif fxr*fxu < 0 xl= xr;
end
xr=xu-fxu*(xl-xu)/(fxl-fxu)
end
end
xr = xl + xu 2
xr = xu − f(xu)(xl − xu) f(xl) − f(xu)
Ensures that
xl and xu always brackets the root
xr = xl + xu 2
44
Ensures that
xl and xu always brackets the root
End of Example L04.7
Method of False Position
Comparison between Bisection method and method of False Position
Note that:
• there is an monotonic
convergence to the exact solution using the method of False position.
• The Bisection method oscillates towards the exact answer.
• Method of False position converges much faster to the exact answer.
46
Root Finding Overview
Numerical algorithms for finding roots of equations
Bracketing methods
Open methods
Method of false position
One point iteration
Bisection method
Newton-Raphson
Secant method
47