Numerical Methods in Engineering (ENGR20005)
Lecture 03 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Lecture Recap
Root Finding Overview
Numerical algorithms for finding roots of equations
Bracketing methods
2
One point iteration
Open methods
Method of false position
Bisection method
Newton-Raphson
Secant method
To find the root of
f(x)=0,
(2.6)
Lecture Recap
Fixed point iteration
rearrange the function so that x is on the left hand side (LHS) of the equation. i.e.
x=g(x).
We will rewrite the above equation as
xi+1=g(xi). (2.7)
Eq. (2.7) above provides us with a iterative formula to find the roots of Eq. (2.6). Note that the functional form of
Eq. (2.6) is not changed. Eq. (2.7) is just a rearrangement of Eq. (2.6).
New guess of x
3
Current guess of x
4
5
L3.1: Sometimes Fixed Point iteration does not work!
Example L03.1: You have seen this function before in Example L02.3 and L02.5.
f(x) = ex − 10x = 0
In those examples, you found that the root to that equation which is at approximately xr=0.1119. There are actually more than 1 root for this equation. Plot f(x) for -1
solution process will blow up. Use this condition only as a good guide. Sometimes you can still get a converged solution if you start have |dg (ξ)/dx| > 1. It really depends on the problem you
are trying to solve.
15
Difference
between old
guess, xi , and actual
root, xr
Example L03.2: Go back to Example L03.1 and see if you can explain why you could only get the root for xr1 but not xr2 using the fixed point iteration scheme.
f(x) = ex − 10x = 0 x=ex =g(x)
f (x)
Convergence criteria
10
dg(ξ) <1
dx
x
g(x)
At x=xr2
dg(ξ) >1
So you cannot get convergence using fixed point iteration scheme
dx
xr1
xr2
16
At x=xr1
dg(ξ) <1
dx
So you can get convergence using fixed point iteration scheme
End of Example L03.2
Fixed point iteration
Exercise: Let
a) Find the two roots for f(x) analytically.
f(x) = x2 − 2x − 3 = 0
b) Show that f(x) can be rearranged into the following forms
x=g1(x)= (2x+3) x=g2(x)= 3
x = g3(x) = 18
x−2 x2 − 3
2
Fixed point iteration
c) With these three expressions for gi(x), use
the fixed point iteration method to find the approximate root for f(x). Start with the initial guess of xo=4. Compare your solutions
to the answer you get in part a).
d) Sketch all the gi(x) and show what is happening graphically.
19
L3.3 Newton-Raphson Method
f (x) f(xold)
xnew
df (xold) = dx
xold
x
Slope of purple line
f(xold) xold − xnew
The value
of x=xr that you want to find
xnew = xold −
dx
f(xold)
df (xold)
tangent to the point xold
f (x)
f(xold)
The value
of x=xr that you want to find
xnew = xold − f(xold)
x new
x
x n
x
ol
ew d
df (x ) dx old
f (x)
f(xold)
The value
of x=xr that you want to find
df (x ) dx old
x
x n
x
ol
ew d
xnew = xold − f(xold)
Newton-Raphson Method
• Most widely used root finding method.
• Converges quickly.
• Only need one initial guess.
f(x)
f’(xi) 1
x
From the graph on the left, it is easy to express xi in terms of
f(xi), f’(xi) and xi+1.
f(xi) xi − xi+1
f(xi) f′(xi)
f′(xi) = xi − xi+1 =
xi+1 = xi −
f(xi) f′(xi)
xi+1
xi
(3.7)
24
Newton-Raphson Method
So to find a root using Newton-Raphson method, do the following:
1) Let the initial guess be xi
2) Find xi+1 by using Eq. (3.7)
3) Let xi= xi+1 repeat steps 2 and 3 until you feel your answer is accurate enough.
(3.7)
xi+1 = xi −
f(xi) f′(xi)
25
Example L03.3: Find the roots of the following equation
f(x) = ex − 10x = 0
Using the Newton Raphson method (Eq. (3.7)). There are actually more than 1 root for this equation. Plot f(x) for -1
f=(50./c)*(1-exp(-9*c/5))-10; df=(90./c)*exp(-9*c/5)-(50./c^2)*(1-exp(-9*c/5)); c=c-f/df
end
43
Calculate the value of f with guess value of c. Remember that you want (to find c such that) f to be zero
Check to see if absolute value of f is close to zero. I have arbitrary chosen threshold value of 1.0e-6. You can try smaller value.
Exactly the same
Bit of code as the for loop
End of Example L03.6
L3.6 Secant Method
Secant Method
In the secant method, the derivative of a function is approximated by a backward finite divided difference
f′(xi) ≈ Substitute Eq. (3.11) into Eq. (3.7)
f(xi) − f(xi−1) xi − xi−1
(3.11)
(3.7)
to obtain the following iterative formula for the Secant Method
(3.12)
Notice that the above formula requires two initial estimates of x, i.e.
xi and xi-1.
xi+1 = xi −
f(xi) f′(xi)
xi+1 = xi −
(xi − xi−1)f(xi) f(xi) − f(xi−1)
46
f (x) f (xi−1)
f(xi)
xi+1,SEC
xi+1,NR
The value
of x=xr that you want to find
xi xi−1 xi+1,NR = xi − f(xi)
xi+1,SEC = xi − (xi − xi−1)f(xi) f(xi) − f(xi−1)
x
f′(xi)
Newton-Raphson method
Secant method
Difference between Secant method and Newton-Raphson method
48
Example L03.7: Write a MATLAB program that uses the Secant method to find the root of
f(c) = 50 (1 − e− 9c ) − 10 = 0 c5
You might like to start with Lecture03M06.m.
Example L03.6 (Newton-Raphson method)
50 ( 9c ) f(c)= c 1−e−5 −10
df(c)=90e−9c −50(1−e−9c) dcc5c2 5
ci+1 = ci − f(ci) f′(ci)
Example L03.7 (Secant method)
f(c) = 50 (1 − e− 9c ) − 10 c 5
Lecture03M06.m
Lecture03M07.m
ci ci−1
f(ci) f(ci−1)
49
ci+1 = ci − (ci − ci−1)f(ci) f(ci) − f(ci−1)
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;
c=4 f=(50./c)*(1-exp(-9*c/5))-10; while abs(f)>1.0e-6
f=(50./c)*(1-exp(-9*c/5))-10; df=(90./c)*exp(-9*c/5)-(50./c^2)*(1-exp(-9*c/5)); c=c-f/df
end
end
fim1=(50./cim1)*(1-exp(-9*cim1/5))-10;
cip1=ci-(ci-cim1)*fi/(fi-fim1)
cim1=ci; ci−1 = ci
Need to explain ci=cip1; ci = ci+1 a little bit
Lecture03M07.m
ci=4
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
cim1=3
ci−1 = ci
Need to explain ci = ci+1 a little bit
i
ci-1
ci
ci+1
c0=3
1
c1=4
c2
2
c1
c2
c3
3
c2
c3
c4
4
c3
c4
c5
ci+1 = ci − (ci − ci−1)f(ci) f(ci) − f(ci−1)
Secant method
50
ci−1 = ci
ci = ci+1
End of Example L03.7