|f (b)| ¤ |f (c)|.

Once an interval [a, b] containing at least one root ± of the function

y = f (x) is found with f (a)f (b) < 0, the algorithm generates a sequence of

values a, b and c such that ± always lies between b and c and, at convergence,

the half-length |c ’ b|/2 is less than a ¬xed tolerance. If the function f is

su¬ciently smooth around the desired root, then the order of convergence

of the algorithm is more than linear (see [Dek69], [Bre73] Chapter 4 and

[Atk89], pp. 91-93).

In the following we describe the main lines of the algorithm as imple-

mented in the MATLAB function fzero. Throughout the parameter d will

be a correction to the point b since it is best to arrange formulae so that

they express the desired quantity as a small correction to a good approx-

imation. For example, if the new value of b were computed as (b + c)/2

(bisection step) a numerical cancellation might occur, while computing b

as b + (c ’ b)/2 gives a more stable formula.

Denote by µ a suitable tolerance (usually the machine precision) and let

c = b; then, the Dekker-Brent method proceeds as follows:

First, check if f (b) = 0. Should this be the case, the algorithm terminates

and returns b as the approximate zero of f . Otherwise, the following steps

are executed:

1. if f (b)f (c) > 0, set c = a, d = b ’ a and e = d.

2. If |f (c)| < |f (b)|, perform the exchanges a = b, b = c and c = a.

3. Set δ = 2µ max {|b|, 1} and m = (c ’ b)/2. If |m| ¤ δ or f (b) = 0 then

the algorithm terminates and returns b as the approximate zero of f .

4. Choose bisection or interpolation.

(a) If |e| < δ or |f (a)| ¤ |f (b)| then a bisection step is taken, i.e., set

d = m and e = m; otherwise, the interpolation step is executed.

(b) if a = c execute linear interpolation, i.e., compute the zero of the

straight line passing through the points (b, f (b)) and (c, f (c)) as

6.3 Fixed-point Iterations for Nonlinear Equations 257

a correction δb to the point b. This amounts to taking a step of

the secant method on the interval having b and c as end points.

If a = c execute inverse quadratic interpolation, i.e., construct

the second-degree polynomial with respect to y, that interpo-

lates at the points (f (a), a), (f (b), b) and (f (c), c) and its value

at y = 0 is computed as a correction δb to the point b. Notice

that at this stage the values f (a), f (b) and f (c) are di¬erent one

from the others, being |f (a)| > |f (b)|, f (b)f (c) < 0 and a = c.

Then the algorithm checks whether the point b + δb can be ac-

cepted. This is a rather technical issue but essentially it amounts

to ascertaining if the point is inside the current interval and not

too close to the end points. This guarantees that the length of

the interval decreases by a large factor when the function is well

behaved. If the point is accepted then e = d and d = δb, i.e.,

the interpolation is actually carried out, else a bisection step is

executed by setting d = m and e = m.

5. The algorithm now updates the current iterate. Set a = b and if

|d| > δ then b = b + d else b = b + δsign(m) and go back to step 1.

Example 6.5 Let us consider the ¬nding of roots of the function f considered in

Example 6.4, taking µ equal to the roundo¬ unit. The MATLAB function fzero

has been employed. It automatically determines the values a and b, starting from

a given initial guess ξ provided by the user. Starting from ξ = 1.5, the algorithm

¬nds the values a = 0.3 and b = 2.1; convergence is achieved in 5 iterations and

the sequences of the values a, b, c and f (b) are reported in Table 6.2.

Notice that the tabulated values refer to the state of the algorithm before step

•

3., and thus, in particular, after possible exchanges between a and b.

k a b c f (b)

0 2.1 0.3 2.1 0.5912

’2.39 · 10’2

1 0.3 0.5235 0.3

3.11 · 10’4

2 0.5235 0.5148 0.5235

’8.8 · 10’7

3 0.5148 0.5149 0.5148

’3.07 · 10’11

4 0.5149 0.5149 0.5148

TABLE 6.2. Solution of the equation cos2 (2x) ’ x2 = 0 using the Dekker-Brent

algorithm. The integer k denotes the current iteration

6.3 Fixed-point Iterations for Nonlinear Equations

In this section a completely general framework for ¬nding the roots of a

nonlinear function is provided. The method is based on the fact that, for a

given f : [a, b] ’ R, it is always possible to transform the problem f (x) = 0

258 6. Root¬nding for Nonlinear Equations

into an equivalent problem x ’ φ(x) = 0, where the auxiliary function

φ : [a, b] ’ R has to be chosen in such a way that φ(±) = ± whenever

f (±) = 0. Approximating the zeros of a function has thus become the

problem of ¬nding the ¬xed points of the mapping φ, which is done by the

following iterative algorithm:

given x(0) , let

k ≥ 0.

x(k+1) = φ(x(k) ), (6.17)

We say that (6.17) is a ¬xed-point iteration and φ is its associated iteration

function. Sometimes, (6.17) is also referred to as Picard iteration or func-

tional iteration for the solution of f (x) = 0. Notice that by construction

the methods of the form (6.17) are strongly consistent in the sense of the

de¬nition given in Section 2.2.

The choice of φ is not unique. For instance, any function of the form

φ(x) = x + F (f (x)), where F is a continuous function such that F (0) = 0,

is an admissible iteration function.

The next two results provide su¬cient conditions in order for the ¬xed-

point method (6.17) to converge to the root ± of problem (6.1). These

conditions are stated precisely in the following theorem.

Theorem 6.1 (convergence of ¬xed-point iterations) Consider the se-

quence x(k+1) = φ(x(k) ), for k ≥ 0, being x(0) given. Assume that:

1. φ : [a, b] ’ [a, b];

2. φ ∈ C 1 ([a, b]);

3. ∃K < 1 : |φ (x)| ¤ K ∀x ∈ [a, b].

Then, φ has a unique ¬xed point ± in [a, b] and the sequence {x(k) } con-

verges to ± for any choice of x(0) ∈ [a, b]. Moreover, we have

x(k+1) ’ ±

lim = φ (±). (6.18)

k’∞ x(k) ’ ±

Proof. The assumption 1. and the continuity of φ ensure that the iteration

function φ has at least one ¬xed point in [a, b]. Assumption 3. states that φ is

a contraction mapping and ensures the uniqueness of the ¬xed point. Indeed,

suppose that there exist two distinct values ±1 , ±2 ∈ [a, b] such that φ(±1 ) = ±1

and φ(±2 ) = ±2 . Expanding φ in a Taylor series around ±1 and truncating it at

¬rst order, it follows that

|±2 ’ ±1 | = |φ(±2 ) ’ φ(±1 )| = |φ (·)(±2 ’ ±1 )| ¤ K|±2 ’ ±1 | < |±2 ’ ±1 |,

for · ∈ (±1 , ±2 ), from which it must necessarily be that ±2 = ±1 .

The convergence analysis for the sequence {x(k) } is again based on a Taylor

series expansion. Indeed, for any k ≥ 0 there exists a value · (k) between ± and

x(k) such that

x(k+1) ’ ± = φ(x(k) ) ’ φ(±) = φ (· (k) )(x(k) ’ ±) (6.19)

6.3 Fixed-point Iterations for Nonlinear Equations 259

from which |x(k+1) ’ ±| ¤ K|x(k) ’ ±| ¤ K k+1 |x(0) ’ ±| ’ 0 for k ’ ∞. Thus,

x(k) converges to ± and (6.19) implies that

x(k+1) ’ ±

= lim φ (· (k) ) = φ (±),

lim (k) ’ ±

k’∞ x k’∞

3

that is (6.18).

The quantity |φ (±)| is called the asymptotic convergence factor and, in

analogy with the case of iterative methods for linear systems, the asymp-

totic convergence rate can be de¬ned as

1

R = ’ log . (6.20)

|φ (±)|

Theorem 6.1 ensures convergence of the sequence {x(k) } to the root ± for

any choice of the initial value x(0) ∈ [a, b]. As such, it represents an example

of a global convergence result.

In practice, however, it is often quite di¬cult to determine a priori the

width of the interval [a, b]; in such a case the following convergence result

can be useful (see for the proof, [OR70]).

Property 6.3 (Ostrowski theorem) Let ± be a ¬xed point of a func-

tion φ, which is continuous and di¬erentiable in a neighborhood J of ±. If

|φ (±)| < 1 then there exists δ > 0 such that the sequence {x(k) } converges

to ±, for any x(0) such that |x(0) ’ ±| < δ.

Remark 6.2 If |φ (±)| > 1 it follows from (6.19) that if x(n) is su¬ciently

close to ±, so that |φ (x(n) )| > 1, then |± ’ x(n+1) | > |± ’ x(n) |, thus

no convergence is possible. In the case |φ (±)| = 1 no general conclusion

can be stated since both convergence and nonconvergence may be possible,

depending on the problem at hand.

Example 6.6 Let φ(x) = x ’ x3 , which admits ± = 0 as ¬xed point. Although

φ (±) = 1, if x(0) ∈ [’1, 1] then x(k) ∈ (’1, 1) for k ≥ 1 and it converges (very

slowly) to ± (if x(0) = ±1, we even have x(k) = ± for any k ≥ 1). Starting from

x(0) = 1/2 the absolute error after 2000 iterations is 0.0158. Let now φ(x) = x+x3

having also ± = 0 as ¬xed point. Again, φ (±) = 1 but in this case the sequence

•

x(k) diverges for any choice x(0) = 0.

We say that a ¬xed-point method has order p (p non necessarily being an

integer) if the sequence that is generated by the method converges to the

¬xed point ± with order p according to De¬nition 6.1.

260 6. Root¬nding for Nonlinear Equations

Property 6.4 If φ ∈ C p+1 (J ) for a suitable neighborhood J of ± and an

integer p ≥ 0, and if φ(i) (±) = 0 for 0 ¤ i ¤ p and φ(p+1) (±) = 0, then the

¬xed-point method with iteration function φ has order p + 1 and

x(k+1) ’ ± φ(p+1) (±)

p ≥ 0.

lim = , (6.21)

k’∞ (x(k) ’ ±)p+1 (p + 1)!

Proof. Let us expand φ in a Taylor series around x = ± obtaining

p

φ(i) (±) (k) φ(p+1) (·) (k)

’±= (x ’ ±)i + (x ’ ±)p+1 ,

(k+1)

x

i! (p + 1)!

i=0

for a certain · between x(k) and ±. Thus, we have

x(k+1) ’ ± φ(p+1) (·) φ(p+1) (±)

lim = lim = .

k’∞ (x(k) ’ ±)p+1 k’∞ (p + 1)! (p + 1)!

3

The convergence of the sequence to the root ± will be faster, for a ¬xed

order p, when the quantity at right-side in (6.21) is smaller.

The ¬xed-point method (6.17) is implemented in Program 51. The variable

phi contains the expression of the iteration function φ.

Program 51 - ¬xpoint : Fixed-point method

function [xvect,xdif,fx,nit]=¬xpoint(x0,nmax,toll,fun,phi)

err=toll+1; nit=0; xvect=x0; x=x0; fx=eval(fun); xdif=[];

while (nit < nmax & err > toll),

nit=nit+1; x=xvect(nit); xn=eval(phi); err=abs(xn-x);

xdif=[xdif; err]; x=xn; xvect=[xvect;x]; fx=[fx;eval(fun)];

end;

6.3.1 Convergence Results for Some Fixed-point Methods

Theorem 6.1 provides a theoretical tool for analyzing some of the iterative

methods introduced in Section 6.2.2.

The chord method. Equation (6.12) is a special instance of (6.17), in

which we let φ(x) = φchord (x) = x’q ’1 f (x) = x’(b’a)/(f (b)’f (a))f (x).

If f (±) = 0, φchord (±) = 1 and the method is not guaranteed to converge.

Otherwise, the condition |φchord (±)| < 1 is equivalent to requiring that

0 < q ’1 f (±) < 2.

Therefore, the slope q of the chord must have the same sign as f (±),

and the search interval [a, b] has to satisfy the constraint

f (b) ’ f (a)

(b ’ a) < 2 .

f (±)

6.4 Zeros of Algebraic Equations 261

The chord method converges in one iteration if f is a straight line, otherwise

it converges linearly, apart the (lucky) case when f (±) = (f (b)’f (a))/(b’

a), for which φchord (±) = 0.

Newton™s method. Equation (6.16) can be cast in the general framework

(6.17) letting

f (x)

φN ewt (x) = x ’ .

f (x)

Assuming f (±) = 0 (that is, ± is a simple root)

f (±)

φN ewt (±) = 0, φN ewt (±) = .

f (±)

If the root ± has multiplicity m > 1, then the method (6.16) is no longer

second-order convergent. Indeed we have (see Exercise 2)

1

φN ewt (±) = 1 ’ . (6.22)

m

If the value of m is known a priori, then the quadratic convergence of

Newton™s method can be recovered by resorting to the so-called modi¬ed

Newton™s method

f (x(k) )

x(k+1) = x(k) ’ m k ≥ 0.

, (6.23)

f (x(k) )

To check the convergence order of the iteration (6.23), see Exercise 2.

6.4 Zeros of Algebraic Equations

In this section we address the special case in which f is a polynomial of

degree n ≥ 0, i.e., a function of the form

n

ak xk ,

pn (x) = (6.24)

k=0

where ak ∈ R are given coe¬cients.

The above representation of pn is not the only one possible. Actually,

one can also write

k

pn (x) = an (x ’ ±1 ) ...(x ’ ±k )

m1 mk

, ml = n

l=1

where ±i and mi denote the i-th root of pn and its multiplicity, respectively.

Other representations are available as well, see Section 6.4.1.

262 6. Root¬nding for Nonlinear Equations

Notice that, since the coe¬cients ak are real, if ± is a zero of pn , then

its complex conjugate ± is a zero of pn too.

¯

Abel™s theorem states that for n ≥ 5 there does not exist an explicit

formula for the zeros of pn (see, for instance, [MM71], Theorem 10.1). This,

in turn, motivates numerical solutions of the nonlinear equation pn (x) = 0.

Since the methods introduced so far must be provided by a suitable search

interval [a, b] or an initial guess x(0) , we recall two results that can be useful

to localize the zeros of a polynomial.

Property 6.5 (Descartes™ rule of signs) Let pn ∈ Pn . Denote by ν the

number of sign changes in the set of coe¬cients {aj } and by k the number

of real positive roots of pn (each counted with its multiplicity). Then, k ¤ ν

and ν ’ k is an even number.

Property 6.6 (Cauchy™s Theorem) All zeros of pn are contained in the

circle “ in the complex plane

“ = {z ∈ C : |z| ¤ 1 + ·k } , max |ak /an |.

where ·k =

0¤k¤n’1

This second property is of little use if ·k 1. In such an event, it is con-

venient to perform a translation through a suitable change of coordinates.

6.4.1 The Horner Method and De¬‚ation

In this section we describe the Horner method for e¬ciently evaluating a

polynomial (and its derivative) at a given point z. The algorithm allows for

generating automatically a procedure, called de¬‚ation, for the sequential

approximation of all the roots of a polynomial.

Horner™s method is based on the observation that any polynomial pn ∈

Pn can be written as

pn (x) = a0 + x(a1 + x(a2 + . . . + x(an’1 + an x) . . . )). (6.25)

Formulae (6.24) and (6.25) are completely equivalent from an algebraic

standpoint; nevertheless, (6.24) requires n sums and 2n ’ 1 multiplications

to evaluate pn (x), while (6.25) requires n sums and n multiplications. The

second expression, known as nested multiplications algorithm, is the basic

ingredient of Horner™s method. This method e¬ciently evaluates the poly-

nomial pn at a point z through the following synthetic division algorithm

bn = an ,

(6.26)

bk = ak + bk+1 z, k = n ’ 1, n ’ 2, ..., 0,

which is implemented in Program 52. The coe¬cients aj of the polynomial

are stored in vector a ordered from an back to a0 .

6.4 Zeros of Algebraic Equations 263

Program 52 - horner : Synthetic division algorithm

function [pnz,b] = horner(a,n,z)

b(1)=a(1); for j=2:n+1, b(j)=a(j)+b(j-1)*z; end; pnz=b(n+1);

All the coe¬cients bk in (6.26) depend on z and b0 = pn (z). The polynomial

n

n’1

bk xk’1

qn’1 (x; z) = b1 + b2 x + ... + bn x = (6.27)

k=1

has degree n ’ 1 in the variable x and depends on the parameter z through

the coe¬cients bk ; it is called the associated polynomial of pn .

Let us now recall the following property of polynomial division:

given two polynomials hn ∈ Pn and gm ∈ Pm with m ¤ n, there exist

an unique polynomial δ ∈ Pn’m and an unique polynomial ρ ∈ Pm’1 such

that

hn (x) = gm (x)δ(x) + ρ(x). (6.28)

Then, dividing pn by x ’ z, from (6.28) it follows that

pn (x) = b0 + (x ’ z)qn’1 (x; z),

having denoted by qn’1 the quotient and by b0 the remainder of the di-

vision. If z is a zero of pn , then b0 = pn (z) = 0 and thus pn (x) =

(x ’ z)qn’1 (x; z). In such a case, the algebraic equation qn’1 (x; z) = 0

yields the n ’ 1 remaining roots of pn (x). This observation suggests adopt-

ing the following de¬‚ation procedure for ¬nding the roots of pn . For m =

n, n ’ 1, . . . , 1:

1. ¬nd a root r of pm using a suitable approximation method;

2. evaluate qm’1 (x; r) by (6.26);

3. let pm’1 = qm’1 .

In the two forthcoming sections some de¬‚ation methods will be ad-

dressed, making a precise choice for the scheme at point 1.

6.4.2 The Newton-Horner Method

A ¬rst example of de¬‚ation employs Newton™s method for computing the

root r at step 1. of the procedure in the previous section. Implement-

ing Newton™s method fully bene¬ts from Horner™s algorithm (6.26). In-

deed, if qn’1 is the associated polynomial of pn de¬ned in (6.27), since

264 6. Root¬nding for Nonlinear Equations

pn (x) = qn’1 (x; z) + (x ’ z)qn’1 (x; z) then pn (z) = qn’1 (z; z). Thanks to

this identity, the Newton-Horner method for the approximation of a root

(real or complex) rj of pn (j = 1, . . . , n) takes the following form:

(0)

given an initial estimate rj of the root, solve for any k ≥ 0

(k) (k)

pn (rj ) pn (rj )

(k+1) (k) (k)

’ ’

rj = rj = rj . (6.29)

(k) (k) (k)

pn (rj ) qn’1 (rj ; rj )

Once convergence has been achieved for the iteration (6.29), polynomial

de¬‚ation is performed, this de¬‚ation being helped by the fact that pn (x) =

(x ’ rj )pn’1 (x). Then, the approximation of a root of pn’1 (x) is carried

out until all the roots of pn have been computed.

Denoting by nk = n ’ k the degree of the polynomial that is obtained at

each step of the de¬‚ation process, for k = 0, . . . , n ’ 1, the computational

cost of each Newton-Horner iteration (6.29) is equal to 4nk . If rj ∈ C, it

(0)

is necessary to work in complex arithmetic and take rj ∈ C; otherwise,

(k)

indeed, the Newton-Horner method (6.29) would yield a sequence {rj } of

real numbers.

The de¬‚ation procedure might be a¬ected by rounding error propagation