Finding the Roots of f (x) = 0

Gerald W. Recktenwald
Department of Mechanical Engineering These slides are a supplement to the book Numerical Methods with Matlab:
Portland State University Implementations and Applications, by Gerald W. Recktenwald, c 2000–2006, Prentice-Hall,
c 2000–2006 Gerald W. Recktenwald.
Upper Saddle River, NJ. These slides are copyright
[email protected] The PDF version of these slides may be downloaded or stored or printed for noncommercial,
educational use. The repackaging or sale of these slides in any form, without written
consent of the author, is prohibited.

The latest version of this PDF file, along with other supplemental material for the book,
can be found at www.prenhall.com/recktenwald or web.cecs.pdx.edu/~gerry/nmm/.

Version 1.12 August 11, 2006

Overview Example: Picnic Table Leg

Topics covered in this chapter Computing the dimensions of a picnic table leg involves a root-finding problem.
Leg assembly Detail of one leg
• Preliminary considerations and bracketing.
• Fixed Point Iteration d1
b c a
• Bisection
b2
• Newton’s Method
h
• The Secant Method d2 α
b θ α
θ
• Hybrid Methods: the built in fzero function 2α
d1

• Roots of Polynomials d2 d2

NMM: Finding the Roots of f (x) = 0 page 2 NMM: Finding the Roots of f (x) = 0 page 3
 Example: Picnic Table Leg Roots of f (x) = 0

Dimensions of a the picnic table leg satisfy Any function of one variable can be put in the form f (x) = 0.

w sin θ = h cos θ + b Example: 1.5

y=x
y = cos(x)

To find the x that satisfies f = cos(x) - x

Given overall dimensions w and h, and the material dimension, b, what is the value of θ ? 1

cos(x) = x,
An analytical solution for θ = f (w, h, b) exists, but is not obvious. 0.5

find the zero crossing of

Use a numerical root-finding procedure to find the value of θ that satisfies 0

f (x) = cos(x) − x = 0 Solution

f (θ) = w sin θ − h cos θ − b = 0 -0.5

-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4

NMM: Finding the Roots of f (x) = 0 page 4 NMM: Finding the Roots of f (x) = 0 page 5

General Considerations Root-Finding Procedure

The basic strategy is

• Is this a special function that will be evaluated often?
• How much precision is needed?
1. Plot the function.
• How fast and robust must the method be?
➣ The plot provides an initial guess, and
• Is the function a polynomial? an indication of potential problems.
• Does the function have singularities? 2. Select an initial guess.
3. Iteratively refine the initial guess
There is no single root-finding method that is best for all situations. with a root-finding algorithm.

NMM: Finding the Roots of f (x) = 0 page 6 NMM: Finding the Roots of f (x) = 0 page 7
 Bracketing Bracketing Algorithm (1)

Algorithm 6.1 Bracket Roots

A root is bracketed on the interval [a, b] if f (a) and f (b) have opposite sign. A sign
given: f (x), xmin, xmax , n
change occurs for singularities as well as roots
dx = (xmax − xmin)/n
xleft = xmin
f(a) f(a) i=0
while i < n
0 0
a b a b i←i+1
f(b) f(b) xright = xleft + dx
if f (x) changes sign in [xleft , xright]
save [xleft , xright ] for further root-finding
end
Bracketing is used to make initial guesses at the roots, not to accurately estimate the xleft = xright
values of the roots. end

NMM: Finding the Roots of f (x) = 0 page 8 NMM: Finding the Roots of f (x) = 0 page 9

Bracketing Algorithm (2) Bracketing Algorithm (3)

A simple test for sign change: f (a) × f (b) < 0 ? A better test uses the built-in sign function

or in Matlab fa = ...
fb = ...

if if sign(fa)~=sign(fb)
fa = ... save bracket
fb = ... end

if fa*fb < 0
save bracket
end
See implementation in the brackPlot function
but this test is susceptible to underflow.

NMM: Finding the Roots of f (x) = 0 page 10 NMM: Finding the Roots of f (x) = 0 page 11
 The brackPlot Function Apply brackPlot Function to sin(x) (1)

brackPlot is a NMM toolbox function that

>> Xb = brackPlot(’sin’,-4*pi,4*pi)
• Looks for brackets of a user-defined f (x) Xb =
1

• Plots the brackets and f (x) -12.5664 -11.2436

f(x) defined in sin.m

-9.9208 -8.5980 0.5
• Returns brackets in a two-column matrix -7.2753 -5.9525
-3.3069 -1.9842
-0.6614 0.6614 0
Syntax: 1.9842 3.3069
5.9525 7.2753
brackPlot(’myFun’,xmin,xmax) 8.5980 9.9208 -0.5
brackPlot(’myFun’,xmin,xmax,nx) 11.2436 12.5664

-1
where
-10 -5 0 5 10
x
myFun is the name of an m-file that evaluates f (x)
xmin, xmax define range of x axis to search
nx is the number of subintervals on [xmin,xmax] used to
check for sign changes of f (x). Default: nx= 20

NMM: Finding the Roots of f (x) = 0 page 12 NMM: Finding the Roots of f (x) = 0 page 13

Apply brackPlot to a user-defined Function (1) Apply brackPlot to a user-defined Function (2)

To solve 1.5
1/3
f (x) = x − x − 2 = 0 >> Xb = brackPlot(’fx3’,0,5)
Xb = 1
we need an m-file function to evaluate f (x) for any scalar or vector of x values. 3.4211 3.6842
0.5

f(x) defined in fx3.m

0

File fx3.m: Note the use of the array operator. -0.5

-1
function f = fx3(x) -1.5
% fx3 Evaluates f(x) = x - x^(1/3) - 2
f = x - x.^(1/3) - 2; -2

-2.5
0 1 2 3 4 5
x
Run brackPlot with fx3 as the input function
>> brackPlot(’fx3’,0,5)
ans =
3.4000 3.6000

NMM: Finding the Roots of f (x) = 0 page 14 NMM: Finding the Roots of f (x) = 0 page 15
 Apply brackPlot to a user-defined Function (3) Root-Finding Algorithms

Instead of creating a separate m-file, we can use an in-line function object. We now proceed to develop the following root-finding algorithms:

>> f = inline(’x - x.^(1/3) - 2’)

f = • Fixed point iteration
Inline function:
f(x) = x - x.^(1/3) - 2 • Bisection
• Newton’s method
>> brackPlot(f,0,5)
ans = • Secant method
3.4000 3.6000

These algorithms are applied after initial guesses at the root(s) are identified with
Note: When an inline function object is supplied to brackPlot, the name of the bracketing (or guesswork).
object is not surrounded in quotes:

brackPlot(f,0,5) instead of brackPlot(’fun’,0,5)

NMM: Finding the Roots of f (x) = 0 page 16 NMM: Finding the Roots of f (x) = 0 page 17

Fixed Point Iteration Convergence Criteria

Fixed point iteration is a simple method. It only works when the iteration function is An automatic root-finding procedure needs to monitor progress toward the root and stop
convergent. when current guess is close enough to the desired root.

Given f (x) = 0, rewrite as xnew = g(xold)

• Convergence checking will avoid searching to unnecessary accuracy.
• Convergence checking can consider whether two successive approximations to the root
Algorithm 6.2 Fixed Point Iteration are close enough to be considered equal.
initialize: x0 = . . . • Convergence checking can examine whether f (x) is sufficiently close to zero at the
for k = 1, 2, . . . current guess.
xk = g(xk−1)
if converged, stop
end More on this later . . .

NMM: Finding the Roots of f (x) = 0 page 18 NMM: Finding the Roots of f (x) = 0 page 19
 Fixed Point Iteration Example (1) Fixed Point Iteration Example (2)

k g1 (xk−1 ) g2 (xk−1 ) g3 (xk−1 )

To solve
1/3 1/3 0 3 3 3
x−x −2=0 g1(x) = x + 2
1 3.4422495703 1 3.5266442931
rewrite as ` ´3
2 3.5098974493 −1 3.5213801474
1/3
xnew = g1(xold) = xold + 2 g2(x) = x − 2
3 3.5197243050 −27 3.5213797068
or 6 + 2x1/3 4 3.5211412691 −24389 3.5213797068
` ´3 g3(x) =
xnew = g2(xold) = xold − 2 3 − x2/3 5 3.5213453678 −1.451 × 1013 3.5213797068
or 6 3.5213747615 −3.055 × 1039 3.5213797068
1/3
6 + 2xold 7 3.5213789946 −2.852 × 10118 3.5213797068
xnew = g3(xold) = 2/3
3− xold 8 3.5213796042 ∞ 3.5213797068
9 3.5213796920 ∞ 3.5213797068
Are these g(x) functions equally effective?

Summary: g1(x) converges, g2(x) diverges, g3(x) converges very quickly

NMM: Finding the Roots of f (x) = 0 page 20 NMM: Finding the Roots of f (x) = 0 page 21

Bisection Bisection (2)

Given a bracketed root, halve the interval while continuing to bracket the root For the bracket interval [a, b] the midpoint is

1
xm = (a + b)
f (x1) 2
f (b1) A better formula, one that is less susceptible to round-off is

b−a
xm = a +
a x2 x1 b 2
f (a1)

NMM: Finding the Roots of f (x) = 0 page 22 NMM: Finding the Roots of f (x) = 0 page 23
 Bisection Algorithm Bisection Example

Algorithm 6.3 Bisection Solve with bisection: k a b xmid f (xmid )

0 3 4
initialize: a = . . ., b = . . . x−x
1/3
−2=0 1 3 4 3.5 -0.01829449
for k = 1, 2, . . .
2 3.5 4 3.75 0.19638375
xm = a + (b − a)/2
3 3.5 3.75 3.625 0.08884159
if sign (f (xm)) = sign (f (xa))
a = xm 4 3.5 3.625 3.5625 0.03522131
else 5 3.5 3.5625 3.53125 0.00845016
b = xm 6 3.5 3.53125 3.515625 -0.00492550
end 7 3.51625 3.53125 3.5234375 0.00176150
if converged, stop 8 3.51625 3.5234375 3.51953125 -0.00158221
end 9 3.51953125 3.5234375 3.52148438 0.00008959
10 3.51953125 3.52148438 3.52050781 -0.00074632

NMM: Finding the Roots of f (x) = 0 page 24 NMM: Finding the Roots of f (x) = 0 page 25

Analysis of Bisection (1) Analysis of Bisection (2)

„ «n
Let δn be the size of the bracketing interval at the nth stage of bisection. Then δn 1 δn
„ «
−n
= =2 or n = log2
δ0 2 δ0

δ0 = b − a = initial bracketing interval δn function

n
1 δ0 evaluations
δ1 = δ0
2 5 3.1 × 10−2 7
1 1 10 9.8 × 10−4 12
δ2 = δ1 = δ0
2 4 20 9.5 × 10−7 22
...
„ «n
δn 1 −n −10
=⇒ = =2 30 9.3 × 10 32
„ «n δ0 2
1 40 9.1 × 10−13 42
δn = δ0
2 δn
„ «
or n = log2 50 8.9 × 10−16 52
δ0

NMM: Finding the Roots of f (x) = 0 page 26 NMM: Finding the Roots of f (x) = 0 page 27
 Convergence Criteria Convergence Criteria on x

An automatic root-finding procedure needs to monitor progress toward the root and stop f (x)
when current guess is close enough to the desired root. true root
˛ ˛
tolerance Absolute tolerance: ˛xk − xk−1˛ < δx
• Convergence checking will avoid searching to unnecessary accuracy. on f (x) x
˛ ˛
• Check whether successive approximations are close enough to be considered the same: ˛x − x
˛ k
˛
k−1 ˛
tolerance Relative tolerance: ˛ ˛ < δ̂x
on x ˛ b−a ˛
|xk − xk−1| < δx

• Check whether f (x) is close enough zero. xk = current guess at the root
xk−1 = previous guess at the root
|f (xk )| < δf

NMM: Finding the Roots of f (x) = 0 page 28 NMM: Finding the Roots of f (x) = 0 page 29

Convergence Criteria on f (x) Convergence Criteria on f (x)

f (x) If f ′(x) is small near the root, it is easy If f ′(x) is large near the root, it is
˛ ˛
true root Absolute tolerance: ˛f (xk )˛ < δf to satisfy a tolerance on f (x) for a large possible to satisfy a tolerance on ∆x
range of ∆x. A tolerance on ∆x is more when |f (x)| is still large. A tolerance
tolerance
on f (x) Relative tolerance: conservative. on f (x) is more conservative.
x
n o
tolerance
|f (xk )| < δ̂f max |f (a0)|, |f (b0)| f (x) f (x)

on x

where a0 and b0 are the original brackets

x x

NMM: Finding the Roots of f (x) = 0 page 30 NMM: Finding the Roots of f (x) = 0 page 31
 Newton’s Method (1) Newton’s Method (2)

Expand f (x) in Taylor Series around xk

For a current guess xk , use f (xk ) and the slope f ′(xk ) to predict where f (x) crosses
the x axis.
˛
(∆x)2 d2f ˛˛
˛
df ˛˛
f (xk + ∆x) = f (xk ) + ∆x + + ...
dx2 ˛
˛
dx ˛xk 2
xk

f(x1)
Substitute ∆x = xk+1 − xk and neglect second order terms to get
′
f (xk+1) ≈ f (xk ) + (xk+1 − xk ) f (xk )
x2 x3 x1 where ˛
f(x2) ′ df ˛˛
f (xk ) =
dx ˛xk

NMM: Finding the Roots of f (x) = 0 page 32 NMM: Finding the Roots of f (x) = 0 page 33

Newton’s Method (3) Newton’s Method Algorithm

Goal is to find x such that f (x) = 0. Algorithm 6.4

Set f (xk+1) = 0 and solve for xk+1 initialize: x1 = . . .

for k = 2, 3, . . .
0 = f (xk ) + (xk+1 − xk ) f (xk )
′ xk = xk−1 − f (xk−1)/f ′(xk−1)
if converged, stop
or, solving for xk+1 end
f (xk )
xk+1 = xk −
f ′(xk )

NMM: Finding the Roots of f (x) = 0 page 34 NMM: Finding the Roots of f (x) = 0 page 35
 Newton’s Method Example (1) Newton’s Method Example (2)

Solve: Conclusion
1/3 1/3
x−x −2=0 xk − xk −2
xk+1 = xk − • Newton’s method converges
−2/3
First derivative is 1 − 13 xk much more quickly than
′ 1 −2/3
f (x) = 1 − x bisection
3 k xk f ′ (xk ) f (x)
The iteration formula is • Newton’s method requires an
1/3
xk − xk − 2 0 3 0.83975005 -0.44224957 analytical formula for f ′(x)
xk+1 = xk − −2/3 1 3.52664429 0.85612976 0.00450679
1 − 13 xk • The algorithm is simple as long
2 3.52138015 0.85598641 3.771 × 10−7 as f ′(x) is available.
3 3.52137971 0.85598640 2.664 × 10−15
• Iterations are not guaranteed to
4 3.52137971 0.85598640 0.0 stay inside an ordinal bracket.

NMM: Finding the Roots of f (x) = 0 page 36 NMM: Finding the Roots of f (x) = 0 page 37

Divergence of Newton’s Method Secant Method (1)

Since Given two guesses xk−1 and xk , the next guess at the root is where the line through
f '(x1) ≈ 0 f (xk−1) and f (xk ) crosses the x axis.
f (xk )
f(x1) xk+1 = xk −
f ′(xk )
f(b)
the new guess, xk+1, will be far from
the old guess whenever f ′(xk ) ≈ 0
x1
f(x1)

x2 a x1 b
f(a)

NMM: Finding the Roots of f (x) = 0 page 38 NMM: Finding the Roots of f (x) = 0 page 39
 Secant Method (2) Secant Method (3)

Given Two versions of this formula are equivalent in exact math:

xk = current guess at the root xk − xk−1
» –
xk+1 = xk − f (xk ) (⋆)
xk−1 = previous guess at the root f (xk ) − f (xk−1)

Approximate the first derivative with and

f (xk )xk−1 − f (xk−1)xk
xk+1 = (⋆⋆)
′ f (xk ) − f (xk−1) f (xk ) − f (xk−1)
f (xk ) ≈
xk − xk−1
Equation (⋆) is better since it is of the form xk+1 = xk + ∆. Even if ∆ is inaccurate
Substitute approximate f ′(xk ) into formula for Newton’s method the change in the estimate of the root will be small at convergence because f (xk ) will
also be small.
f (xk )
xk+1 = xk − Equation (⋆⋆) is susceptible to catastrophic cancellation:
f ′(xk )
to get • f (xk ) → f (xk−1) as convergence approaches, so cancellation error in
»
xk − xk−1
– the denominator can be large.
xk+1 = xk − f (xk ) • |f (x)| → 0 as convergence approaches, so underflow is possible
f (xk ) − f (xk−1)

NMM: Finding the Roots of f (x) = 0 page 40 NMM: Finding the Roots of f (x) = 0 page 41

Secant Algorithm Secant Method Example

Algorithm 6.5 Solve: Conclusions

1/3
x−x −2=0 • Converges almost as quickly as
initialize: x1 = . . ., x2 = . . .
for k = 2, 3 . . . Newton’s method.
xk+1 = xk k xk−1 xk f (xk ) • No need to compute f ′(x).
−f (xk )(xk − xk−1)/(f (xk ) − f (xk−1)) 0 4 3 −0.44224957
• The algorithm is simple.
if converged, stop 1 3 3.51734262 −0.00345547
2 3.51734262 3.52141665 0.00003163
• Two initial guesses are necessary
end
3 3.52141665 3.52137970 −2.034 × 10−9 • Iterations are not guaranteed to
4 3.52137959 3.52137971 −1.332 × 10−15 stay inside an ordinal bracket.
5 3.52137971 3.52137971 0.0

NMM: Finding the Roots of f (x) = 0 page 42 NMM: Finding the Roots of f (x) = 0 page 43
 Divergence of Secant Method Summary of Basic Root-finding Methods

• Plot f (x) before searching for roots

Since • Bracketing finds coarse interval containing roots and singularities
f '(x) ≈ 0
xk − xk−1
» –
f(x2)
xk+1 = xk −f (xk ) • Bisection is robust, but converges slowly
f(x3)
f (xk ) − f (xk−1)
• Newton’s Method
the new guess, xk+1, will be far from the ⊲ Requires f (x) and f ′(x).
x1 x3 x2 old guess whenever f ′(xk ) ≈ f (xk−1) ⊲ Iterates are not confined to initial bracket.
and |f (x)| is not small. ⊲ Converges rapidly.
f (x1)
⊲ Diverges if f ′(x) ≈ 0 is encountered.
• Secant Method
⊲ Uses f (x) values to approximate f ′(x).
⊲ Iterates are not confined to initial bracket.
⊲ Converges almost as rapidly as Newton’s method.
⊲ Diverges if f ′(x) ≈ 0 is encountered.

NMM: Finding the Roots of f (x) = 0 page 44 NMM: Finding the Roots of f (x) = 0 page 45

fzero Function (1) Reverse Quadratic Interpolation

20
fzero is a hybrid method that combines bisection, secant and reverse quadratic
interpolation Find the point where the x
axis intersects the sideways 15
Syntax: parabola passing through
three pairs of (x, f (x))
r = fzero(’fun’,x0)
values. 10
r = fzero(’fun’,x0,options)
r = fzero(’fun’,x0,options,arg1,arg2,...)

5
x0 can be a scalar or a two element vector

0
• If x0 is a scalar, fzero tries to create its own bracket.
• If x0 is a two element vector, fzero uses the vector as a bracket.
−5
0 0.5 1 1.5 2

NMM: Finding the Roots of f (x) = 0 page 46 NMM: Finding the Roots of f (x) = 0 page 47
 fzero Function (2) fzero Function (3)

Optional parameters to control fzero are specified with the optimset function.
fzero chooses next root as
Examples:

• Result of reverse quadratic interpolation (RQI) if that result is inside the current Tell fzero to display the results of each step:
bracket. >> options = optimset(’Display’,’iter’);
• Result of secant step if RQI fails, and if the result of secant method is in inside the >> x = fzero(’myFun’,x0,options)
current bracket.
• Result of bisection step if both RQI and secant method fail to produce guesses inside Tell fzero to use a relative tolerance of 5 × 10−9:
the current bracket. >> options = optimset(’TolX’,5e-9);
>> x = fzero(’myFun’,x0,options)

Tell fzero to suppress all printed output, and use a relative tolerance of 5 × 10−4:
>> options = optimset(’Display’,’off’,’TolX’,5e-4);
>> x = fzero(’myFun’,x0,options)

NMM: Finding the Roots of f (x) = 0 page 48 NMM: Finding the Roots of f (x) = 0 page 49

fzero Function (4) Roots of Polynomials

Allowable options (specified via optimset):

3
Complications arise due to
Option type Value Effect
• Repeated roots
’Display’ ’iter’ Show results of each iteration 2
• Complex roots
’final’ Show root and original bracket
• Sensitivity of roots to small

y = f(x)
’off’ Suppress all print out
perturbations in the 1
’TolX’ tol Iterate until polynomial coefficients f1(x) f2(x) f3(x)
|∆x| < max [tol, tol ∗ a, tol ∗ b] (conditioning).
0
where ∆x = (b−a)/2, and [a, b] is the current bracket.
distinct repeated complex
real roots real roots roots
-1
The default values of ’Display’ and ’TolX’ are equivalent to 0 2 4 6 8 10
x (arbitrary units)
options = optimset(’Display’,’iter’,’TolX’,eps)

NMM: Finding the Roots of f (x) = 0 page 50 NMM: Finding the Roots of f (x) = 0 page 51
 Algorithms for Finding Polynomial Roots roots Function (1)

The built-in roots function uses the companion matrix method

• Bairstow’s method
• Müller’s method • No initial guess
• Laguerre’s method • Returns all roots of the polynomial
• Jenkin’s–Traub method • Solves eigenvalue problem for companion matrix
• Companion matrix method
Write polynomial in the form
n n−1
c1x + c2x + . . . + cnx + cn+1 = 0

Then, for a third order polynomial

>> c = [c1 c2 c3 c4];

>> r = roots(c)

NMM: Finding the Roots of f (x) = 0 page 52 NMM: Finding the Roots of f (x) = 0 page 53

roots Function (2) roots Function (3)

The eigenvalues of The statements

2 3
−c2/c1 −c3/c1 −c4/c1 −c5/c1 c = ... % vector of polynomial coefficients
6 1 0 0 0 7 r = roots(c);
A=4
6 7
0 1 0 0 5
0 0 1 0 are equivalent to

are the same as the roots of

c = ...
4 3 2 n = length(c);
c5λ + c4λ + c3λ + c2λ + c1 = 0.
A = diag(ones(1,n-2),-1); % ones on first subdiagonal
A(1,:) = -c(2:n) ./ c(1); % first row is -c(j)/c(1), j=2..n
r = eig(A);

NMM: Finding the Roots of f (x) = 0 page 54 NMM: Finding the Roots of f (x) = 0 page 55
 roots Examples

Roots of are found with

2
>> roots([1 -3 2])
f1(x) = x − 3x + 2 ans =
2
f2(x) = x − 10x + 25 2
1
2
f3(x) = x − 17x + 72.5
>> roots([1 -10 25])
ans =
5
5

>> roots([1 -17 72.5])

ans =
8.5000 + 0.5000i
8.5000 - 0.5000i

NMM: Finding the Roots of f (x) = 0 page 56