ńņš. 11 |

are all contained within the starting interval [x(ā’1) , x(0) ].

In Figure 6.3 (right), the ļ¬rst two iterations of both the secant and Regula

Falsi methods are shown, starting from the same initial data x(ā’1) and x(0) .

Notice that the iterate x(1) computed by the secant method coincides with

that computed by the Regula Falsi method, while the value x(2) computed

(2)

by the former method (and denoted in the ļ¬gure by xSec ) falls outside the

searching interval [x(ā’1) , x(0) ].

In this respect, the Regula Falsi method, as well as the bisection method,

can be regarded as a globally convergent method.

y y f (x)

f (x)

(2)

xSec

x(2) x(1) x(0) x(ā’1) x(1) x(0)

(ā’1)

x

x x

x(2)

FIGURE 6.3. The ļ¬rst two steps of the Regula Falsi method for two diļ¬erent

functions

Newtonā™s method.

Assuming that f ā C 1 (I) and that f (Ī±) = 0 (i.e., Ī± is a simple root of f ),

if we let

āk ā„ 0

qk = f (x(k) ),

and assign the initial value x(0) , we obtain the so called Newtonā™s method

f (x(k) )

ā’ k ā„ 0.

(k+1) (k)

x =x , (6.16)

f (x(k) )

254 6. Rootļ¬nding for Nonlinear Equations

0

10

y

(1)

ā’5

10

f (x)

(2)

ā’10

(3)

10

(2)

x

x(0) x (4)

ā’15

10

a b

x(1) 0 5 10 15 20 25 30 35

FIGURE 6.4. The ļ¬rst two steps of Newtonā™s method (left); convergence histories

in Example 6.4 for the chord method (1), bisection method (2), secant method

(3) and Newtonā™s method (4) (right). The number of iterations and the absolute

error as a function of k are shown on the x-axis and y-axis, respectively

At the k-th iteration, Newtonā™s method requires the two functional evalu-

ations f (x(k) ) and f (x(k) ). The increasing computational cost with respect

to the methods previously considered is more than compensated for by a

higher order of convergence, Newtonā™s method being of order 2 (see Section

6.3.1).

Example 6.4 Let us compare the methods introduced so far for the approxima-

tion of the root Ī± 0.5149 of the function f (x) = cos2 (2x) ā’ x2 in the interval

(0, 1.5). The tolerance Īµ on the absolute error has been taken equal to 10ā’10 and

the convergence histories are drawn in Figure 6.4 (right). For all methods, the

initial guess x(0) has been set equal to 0.75. For the secant method we chose

x(ā’1) = 0.

The analysis of the results singles out the slow convergence of the chord

method. The error curve for the Regula Falsi method is similar to that of se-

cant method, thus it was not reported in Figure 6.4.

It is interesting to compare the performances of Newtonā™s and secant methods

(both having order p > 1), in terms of their computational eļ¬ort. It can indeed

be proven that it is more convenient to employ the secant method whenever the

number of ļ¬‚oating point operations to evaluate f are about twice those needed

for evaluating f (see [Atk89], pp. 71-73). In the example at hand, Newtonā™s

method converges to Ī± in 6 iterations, instead of 7, but the secant method takes

ā¢

94 ļ¬‚ops instead of 177 ļ¬‚ops required by Newtonā™s method.

The chord, secant, Regula Falsi and Newtonā™s methods are implemented

in Programs 47, 48, 49 and 50, respectively. Here and in the rest of the

chapter, x0 and xm1 denote the initial data x(0) and x(ā’1) . In the case of

the Regula Falsi method the stopping test checks is |f (x(k) )| < toll, while

for the other methods the test is |x(k+1) ā’ x(k) | < toll. The string dfun

contains the expression of f to be used in the Newton method.

6.2 A Geometric Approach to Rootļ¬nding 255

Program 47 - chord : The chord method

function [xvect,xdif,fx,nit]=chord(a,b,x0,nmax,toll,fun)

x=a; fa=eval(fun); x=b; fb=eval(fun); r=(fb-fa)/(b-a);

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=x-fx(nit)/r; err=abs(xn-x);

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

end;

Program 48 - secant : The secant method

function [xvect,xdif,fx,nit]=secant(xm1,x0,nmax,toll,fun)

x=xm1; fxm1=eval(fun); xvect=[x]; fx=[fxm1]; x=x0; fx0=eval(fun);

xvect=[xvect;x]; fx=[fx;fx0]; err=toll+1; nit=0; xdif=[];

while (nit < nmax & err > toll),

nit=nit+1; x=x0-fx0*(x0-xm1)/(fx0-fxm1); xvect=[xvect;x];

fnew=eval(fun); fx=[fx;fnew]; err=abs(x0-x); xdif=[xdif;err];

xm1=x0; fxm1=fx0; x0=x; fx0=fnew;

end;

Program 49 - regfalsi : The Regula Falsi method

function [xvect,xdif,fx,nit]=regfalsi(xm1,x0,toll,nmax,fun)

nit=0; x=xm1; f=eval(fun); fx=[f]; x=x0; f=eval(fun); fx=[fx, f];

xvect=[xm1,x0]; xdif=[]; f=toll+1; kprime=1;

while (nit < nmax & (abs(f) > toll),

nit=nit+1; dim=length(xvect);

x=xvect(dim); fxk=eval(fun); xk=x; i=dim;

while (i >= kprime), i=i-1; x=xvect(i); fxkpr=eval(fun);

if ((fxkpr*fxk) < 0), xkpr=x; kprime=i; break; end;

end;

x=xk-fxk*(xk-xkpr)/(fxk-fxkpr); xvect=[xvect, x]; f=eval(fun);

fx=[fx, f]; err=abs(x-xkpr); xdif=[xdif, err];

end;

Program 50 - newton : Newtonā™s method

function [xvect,xdif,fx,nit]=newton(x0,nmax,toll,fun,dfun)

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

while (nit < nmax & err > toll),

nit=nit+1; x=xvect(nit); dfx=eval(dfun);

if (dfx == 0), err=toll*1.e-10;

disp(ā™ Stop for vanishing dfun ā™);

else,

xn=x-fx(nit)/dfx; err=abs(xn-x); xdif=[xdif; err];

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

256 6. Rootļ¬nding for Nonlinear Equations

end;

end;

6.2.3 The Dekker-Brent Method

The Dekker-Brent method combines the bisection and secant methods, pro-

viding a synthesis of the advantages of both. This algorithm carries out an

iteration in which three abscissas a, b and c are present at each stage. Nor-

mally, b is the latest iterate and closest approximation to the zero, a is

the previous iterate and c is the previous or an older iterate so that f (b)

and f (c) have opposite signs. At all times b and c bracket the zero and

|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

and, as a consequence, can lead to inaccurate results. For the sake of stabil-

ity, it is therefore convenient to approximate ļ¬rst the root r1 of minimum

module, which is the most sensitive to ill-conditioning of the problem (see

Example 2.7, Chapter 2) and then to continue with the successive roots

r2 , . . . , until the root of maximum module is computed. To localize r1 , the

techniques described in Section 5.1 or the method of Sturm sequences can

be used (see [IK66], p. 126).

A further increase in accuracy can be obtained, once an approximation rj

of the root rj is available, by going back to the original polynomial pn and

generating through the Newton-Horner method (6.29) a new approximation

(0)

to rj , taking as initial guess rj = rj . This combination of deļ¬‚ation and

successive correction of the root is called the Newton-Horner method with

reļ¬nement.

Example 6.7 Let us examine the performance of the Newton-Horner method in

two cases: in the ļ¬rst one, the polynomial admits real roots, while in the second

one there are two pairs of complex conjugate roots. To single out the importance

of reļ¬nement, we have implemented (6.29) both switching it on and oļ¬ (methods

NwtRef and Nwt, respectively). The approximate roots obtained using method Nwt

are denoted by rj , while sj are those computed by method NwtRef. As for the

numerical experiments, the computations have been done in complex arithmetic,

with x(0) = 0 + i 0, i being the imaginary unit, nmax = 100 and toll = 10ā’5 . The

tolerance for the stopping test in the reļ¬nement cycle has been set to 10ā’3 toll.

1) p5 (x) = x5 + x4 ā’ 9x3 ā’ x2 + 20x ā’ 12 = (x ā’ 1)2 (x ā’ 2)(x + 2)(x + 3).

6.4 Zeros of Algebraic Equations 265

We report in Tables 6.3(a) and 6.3(b) the approximate roots rj (j = 1, . . . , 5)

and the number of Newton iterations (Nit) needed to get each of them; in the

case of method NwtRef we also show the number of extra Newton iterations for

the reļ¬nement (Extra).

rj sj

Nit Nit Extra

0.99999348047830 17 0.9999999899210124 17 10

1 ā’ i3.56 Ā· 10ā’25 1 ā’ i2.40 Ā· 10ā’28

6 6 10

2 ā’ i2.24 Ā· 10ā’13 2 + i1.12 Ā· 10ā’22

9 9 1

ā’2 ā’ i1.70 Ā· 10ā’10 ā’2 + i8.18 Ā· 10ā’22

7 7 1

ā’3 + i5.62 Ā· 10ā’6 ā’3 ā’ i7.06 Ā· 10ā’21

1 1 2

(a) (b)

TABLE 6.3. Roots of the polynomial p5 . Roots computed by the Newton-Horner

method without reļ¬nement (left), and with reļ¬nement (right)

Notice a neat increase in the accuracy of rootļ¬nding due to reļ¬nement, even with

few extra iterations.

2) p6 (x) = x6 ā’ 2x5 + 5x4 ā’ 6x3 + 2x2 + 8x ā’ 8.

The zeros of p6 are the complex numbers {1, ā’1, 1 Ā± i, Ā±2i}. We report below,

denoting them by rj , (j = 1, . . . , 6), the approximations to the roots of p6 ob-

tained using method Nwt, with a number of iterations equal to 2, 1, 1, 7, 7 and 1,

respectively. Beside, we also show the corresponding approximations sj computed

by method NwtRef and obtained with a maximum number of 2 extra iterations.

ā¢

rj sj

Nwt NwtRef

r1 1 s1 1

ā’0.99 ā’ i9.54 Ā· 10ā’17 ā’1 + i1.23 Ā· 10ā’32

r2 s2

r3 1+i s3 1+i

r4 1-i s4 1-i

-1.31 Ā· 10ā’8 + i2 ā’5.66 Ā· 10ā’17 + i2

r5 s5

r6 -i2 s6 -i2

TABLE 6.4. Roots of the polynomial p6 obtained using the Newton-Horner

method without (left) and with (right) reļ¬nement

A coding of the Newton-Horner algorithm is provided in Program 53. The

input parameters are A (a vector containing the polynomial coeļ¬cients), n

(the degree of the polynomial), toll (tolerance on the maximum variation

between successive iterates in Newtonā™s method), x0 (initial value, with

x(0) ā R), nmax (maximum number of admissible iterations for Newtonā™s

266 6. Rootļ¬nding for Nonlinear Equations

method) and iref (if iref = 1, then the reļ¬nement procedure is activated).

For dealing with the general case of complex roots, the initial datum is

automatically converted into the complex number z = x(0) + ix(0) , where

ā

i = ā’1.

The program returns as output the variables xn (a vector containing the

sequence of iterates for each zero of pn (x)), iter (a vector containing the

number of iterations needed to approximate each root), itrefin (a vector

containing the Newton iterations required to reļ¬ne each estimate of the

computed root) and root (vector containing the computed roots).

Program 53 - newthorn : Newton-Horner method with reļ¬nement

function [xn,iter,root,itreļ¬n]=newthorn(A,n,toll,x0,nmax,iref)

apoly=A;

for i=1:n, it=1; xn(it,i)=x0+sqrt(-1)*x0; err=toll+1; Ndeg=n-i+1;

if (Ndeg == 1), it=it+1; xn(it,i)=-A(2)/A(1);

else

while (it < nmax & err > toll),

[px,B]=horner(A,Ndeg,xn(it,i));

[pdx,C]=horner(B,Ndeg-1,xn(it,i));

it=it+1; if (pdx Ė=0), xn(it,i)=xn(it-1,i)-px/pdx;

err=max(abs(xn(it,i)-xn(it-1,i)),abs(px));

else,

disp(ā™ Stop due to a vanishing pā™ā™ ā™);

err=0; xn(it,i)=xn(it-1,i);

end

end

end

A=B;

if (iref==1), alfa=xn(it,i); itr=1; err=toll+1;

while ((err > toll*1e-3) & (itr < nmax))

[px,B]=horner(apoly,n,alfa);

[pdx,C]=horner(B,n-1,alfa); itr=itr+1;

if (pdxĖ=0)

alfa2=alfa-px/pdx;

err=max(abs(alfa2-alfa),abs(px)); alfa=alfa2;

else,

disp(ā™ Stop due to a vanishing pā™ā™ ā™); err=0;

end

end; itreļ¬n(i)=itr-1; xn(it,i)=alfa;

end

iter(i)=it-1; root(i)=xn(it,i); x0=root(i);

end

6.4 Zeros of Algebraic Equations 267

6.4.3 The Muller Method

A second example of deļ¬‚ation employs Mullerā™s method for ļ¬nding an

approximation to the root r at step 1. of the procedure described in Section

6.4.1 (see [Mul56]). Unlike Newtonā™s or secant methods, Mullerā™s method is

able to compute complex zeros of a given function f , even starting from a

real initial datum; moreover, its order of convergence is almost quadratic.

The action of Mullerā™s method is drawn in Figure 6.5. The scheme ex-

tends the secant method, substituting the linear polynomial introduced in

(6.13) with a second-degree polynomial as follows. Given three distinct

values x(0) , x(1) and x(2) , the new point x(3) is determined by setting

p2 (x(3) ) = 0, where p2 ā P2 is the unique polynomial that interpolates

f at the points x(i) , i = 0, 1, 2, that is, p2 (x(i) ) = f (x(i) ) for i = 0, 1, 2.

Therefore,

p2

x(0) x(1) x(2)

x(3)

f

FIGURE 6.5. The ļ¬rst step of Mullerā™s method

p2 (x) = f (x(2) ) + (x ā’ x(2) )f [x(2) , x(1) ] + (x ā’ x(2) )(x ā’ x(1) )f [x(2) , x(1) , x(0) ]

where

f (Ī·) ā’ f (Ī¾) f [Ī·, Ļ„ ] ā’ f [Ī¾, Ī·]

f [Ī¾, Ī·] = , f [Ī¾, Ī·, Ļ„ ] =

Ī·ā’Ī¾ Ļ„ ā’Ī¾

are the divided diļ¬erences of order 1 and 2 associated with the points Ī¾, Ī·

and Ļ„ (see Section 8.2.1). Noticing that x ā’ x(1) = (x ā’ x(2) ) + (x(2) ā’ x(1) ),

we get

p2 (x) = f (x(2) ) + w(x ā’ x(2) ) + f [x(2) , x(1) , x(0) ](x ā’ x(2) )2

having deļ¬ned

= f [x(2) , x(1) ] + (x(2) ā’ x(1) )f [x(2) , x(1) , x(0) ]

w

= f [x(2) , x(1) ] + f [x(2) , x(0) ] ā’ f [x(0) , x(1) ].

268 6. Rootļ¬nding for Nonlinear Equations

Requiring that p2 (x(3) ) = 0 it follows that

1/2

ā’w Ā± w2 ā’ 4f (x(2) )f [x(2) , x(1) , x(0) ]

x(3) (2)

=x + .

2f [x(2) , x(1) , x(0) ]

Similar computations must be done for getting x(4) starting from x(1) , x(2)

and x(3) and, more generally, to ļ¬nd x(k+1) starting from x(kā’2) , x(kā’1)

and x(k) , with k ā„ 2, according with the following formula (notice that the

numerator has been rationalized)

2f (x(k) )

ā’

(k+1) (k)

x =x . (6.30)

1/2

w ā“ w2 ā’ 4f (x(k) )f [x(k) , x(kā’1) , x(kā’2) ]

The sign in (6.30) is chosen in such a way that the module of the denomina-

tor is maximized. Assuming that f ā C 3 (J ) in a suitable neighborhood J

of the root Ī±, with f (Ī±) = 0, the order of convergence is almost quadratic.

Precisely, the error e(k) = Ī± ā’ x(k) obeys the following relation (see for the

proof [Hil87])

|e(k+1) | 1 f (Ī±)

lim = , p 1.84.

kā’ā |e(k) |p 6 f (Ī±)

Example 6.8 Let us employ Mullerā™s method to approximate the roots of the

polynomial p6 examined in Example 6.7. The tolerance on the stopping test

is toll = 10ā’6 , while x(0) = ā’5, x(1) = 0 and x(2) = 5 are the inputs to

(6.30). We report in Table 6.5 the approximate roots of p6 , denoted by sj and

rj (j = 1, . . . , 5), where, as in Example 6.7, sj and rj have been obtained by

switching the reļ¬nement procedure on and oļ¬, respectively. To compute the roots

rj , 12, 11, 9, 9, 2 and 1 iterations are needed, respectively, while only one extra

iteration is taken to reļ¬ne all the roots.

rj sj

ā’15

1 + i9.9 Ā· 10ā’18

1 + i2.2 Ā· 10

r1 s1

ā’1 ā’ i8.4 Ā· 10ā’16

r2 s2 -1

r3 0.99 + i s3 1+i

0.99 ā’ i 1ā’i

r4 s4

ā’1.1 Ā· 10ā’15 + i1.99

r5 s5 i2

ā’1.0 Ā· 10ā’15 ā’ i2

r6 s6 -i2

TABLE 6.5. Roots of polynomial p6 with Mullerā™s method without (rj ) and with

(sj ) reļ¬nement

Even in this example, one can notice the eļ¬ectiveness of the reļ¬nement procedure,

based on Newtonā™s method, on the accuracy of the solution yielded by (6.30). ā¢

6.5 Stopping Criteria 269

The Muller method is implemented in Program 54, in the special case

where f is a polynomial of degree n. The deļ¬‚ation process also includes a

reļ¬nement phase; the evaluation of f (x(kā’2) ), f (x(kā’1) ) and f (x(k) ), with

k ā„ 2, is carried out using Program 52. The input/output parameters are

analogous to those described in Program 53.

Program 54 - mulldeļ¬‚ : Mullerā™s method with reļ¬nement

function [xn,iter,root,itreļ¬n]=mulldeļ¬‚(A,n,toll,x0,x1,x2,nmax,iref)

apoly=A;

for i=1:n

xn(1,i)=x0; xn(2,i)=x1; xn(3,i)=x2; it=0; err=toll+1; k=2; Ndeg=n-i+1;

if (Ndeg == 1), it=it+1; k=0; xn(it,i)=-A(2)/A(1);

else

while ((err > toll) & (it < nmax)),

k=k+1; it=it+1; [f0,B]=horner(A,Ndeg,xn(k-2,i));

[f1,B]=horner(A,Ndeg,xn(k-1,i)); [f2,B]=horner(A,Ndeg,xn(k,i));

f01=(f1-f0)/(xn(k-1,i)-xn(k-2,i)); f12=(f2-f1)/(xn(k,i)-xn(k-1,i));

f012=(f12-f01)/(xn(k,i)-xn(k-2,i)); w=f12+(xn(k,i)-xn(k-1,i))*f012;

arg=wĖ2-4*f2*f012; d1=w-sqrt(arg); d2=w+sqrt(arg); den=max(d1,d2);

if (denĖ=0); xn(k+1,i)=xn(k,i)-(2*f2)/den;

err=abs(xn(k+1,i)-xn(k,i));

else

disp(ā™ Vanishing denominator ā™); return; end;

end; end; radix=xn(k+1,i);

if (iref==1),

alfa=radix; itr=1; err=toll+1;

while ((err > toll*1e-3) & (itr < nmax)),

[px,B]=horner(apoly,n,alfa); [pdx,C]=horner(B,n-1,alfa);

if (pdx == 0), disp(ā™ Vanishing derivative ā™); err=0; end;

itr=itr+1; if (pdxĖ=0), alfa2=alfa-px/pdx;

err=abs(alfa2-alfa); alfa=alfa2; end;

end; itreļ¬n(i)=itr-1; xn(k+1,i)=alfa; radix=alfa;

end

iter(i)=it; root(i)=radix; [px,B]=horner(A,Ndeg-1,xn(k+1,i)); A=B;

end

6.5 Stopping Criteria

Suppose that {x(k) } is a sequence converging to a zero Ī± of the function

f . In this section we provide some stopping criteria for terminating the

iterative process that approximates Ī±. Analogous to Section 4.6, where

the case of iterative methods for linear systems has been examined, there

are two possible criteria: a stopping test based on the residual and on the

increment. Below, Īµ is a ļ¬xed tolerance on the approximate calculation of

270 6. Rootļ¬nding for Nonlinear Equations

Ī± and e(k) = Ī± ā’ x(k) denotes the absolute error. We shall moreover assume

that f is continuously diļ¬erentiable in a suitable neighborhood of the root.

1. Control of the residual: the iterative process terminates at the ļ¬rst

step k such that |f (x(k) )| < Īµ.

Situations can arise where the test turns out to be either too restrictive or

excessively optimistic (see Figure 6.6). Applying the estimate (6.6) to the

case at hand yields

1/m

|e(k) | m!

|f (x(k) )|1/m .

|Ī±| |f (m) (Ī±)||Ī±|m

In particular, in the case of simple roots, the error is bound to the residual

by the factor 1/|f (Ī±)| so that the following conclusions can be drawn:

1. if |f (Ī±)| 1, then |e(k) | Īµ; therefore, the test provides a satisfac-

tory indication of the error;

2. if |f (Ī±)| 1, the test is not reliable since |e(k) | could be quite large

with respect to Īµ;

3. if, ļ¬nally, |f (Ī±)| 1, we get |e(k) | Īµ and the test is too restrictive.

We refer to Figure 6.6 for an illustration of the last two cases.

f (x)

f (x)

x(k)

Ī±

x(k)

Ī±

FIGURE 6.6. Two situations where the stopping test based on the residual

is either too restrictive (when |e(k) | |f (x(k) )|, left) or too optimistic (when

|e(k) | |f (x(k) )|, right)

The conclusions that we have drawn agree with those in Example 2.4.

Indeed, when f (Ī±) 0, the condition number of the problem f (x) = 0 is

very high and, as a consequence, the residual does not provide a signiļ¬cant

indication of the error.

2. Control of the increment: the iterative process terminates as soon as

|x(k+1) ā’ x(k) | < Īµ.

6.5 Stopping Criteria 271

Let x(k) be generated by the ļ¬xed-point method x(k+1) = Ļ(x(k) ). Using

the mean value theorem, we get

e(k+1) = Ļ(Ī±) ā’ Ļ(x(k) ) = Ļ (Ī¾ (k) )e(k) ,

where Ī¾ (k) lies between x(k) and Ī±. Then,

x(k+1) ā’ x(k) = e(k) ā’ e(k+1) = 1 ā’ Ļ (Ī¾ (k) ) e(k)

so that, assuming that we can replace Ļ (Ī¾ (k) ) with Ļ (Ī±), it follows that

1

(x(k+1) ā’ x(k) ).

e(k) (6.31)

1 ā’ Ļ (Ī±)

Ī³

1

1

2

-1 0 1 Ļ (Ī±)

FIGURE 6.7. Behavior of Ī³ = 1/(1 ā’ Ļ (Ī±)) as a function of Ļ (Ī±)

As shown in Figure 6.7, we can conclude that the test:

- is unsatisfactory if Ļ (Ī±) is close to 1;

- provides an optimal balancing between increment and error in the case

of methods of order 2 for which Ļ (Ī±) = 0 as is the case for Newtonā™s

method;

- is still satisfactory if ā’1 < Ļ (Ī±) < 0.

Example 6.9 The zero of the function f (x) = eā’x ā’ Ī· is given by Ī± = ā’ log(Ī·).

For Ī· = 10ā’9 , Ī± 20.723 and f (Ī±) = ā’eā’Ī± ā’10ā’9 . We are thus in the case

where |f (Ī±)| 1 and we wish to examine the behaviour of Newtonā™s method in

the approximation of Ī± when the two stopping criteria above are adopted in the

computations.

We show in Tables 6.6 and 6.7 the results obtained using the test based on the

control of the residual (1) and of the increment (2), respectively. We have taken

x(0) = 0 and used two diļ¬erent values of the tolerance. The number of iterations

required by the method is denoted by nit.

According to (6.31), since Ļ (Ī±) = 0, the stopping test based on the increment

reveals to be reliable for both the values (which are quite diļ¬ering) of the stop

tolerance Īµ. The test based on the residual, instead, yields an acceptable estimate

of the root only for very small tolerances, while it is completely wrong for large

values of Īµ.

ā¢

272 6. Rootļ¬nding for Nonlinear Equations

|f (x(nit) )| |Ī± ā’ x(nit) | |Ī± ā’ x(nit) |/Ī±

Īµ nit

10ā’10 5.9 Ā· 10ā’11 5.7 Ā· 10ā’2

22 0.27

10ā’3 9.1 Ā· 10ā’4

7 13.7 66.2

TABLE 6.6. Newtonā™s method for the approximation of the root of

f (x) = eā’x ā’ Ī· = 0. The stopping test is based on the control of the residual

|x(nit) ā’ x(nitā’1) | |Ī± ā’ x(nit) | |Ī± ā’ x(nit) |/Ī±

Īµ nit

10ā’10 8.4 Ā· 10ā’13

26 0 0

10ā’3 1.3 Ā· 10ā’6 8.4 Ā· 10ā’13 4 Ā· 10ā’12

25

TABLE 6.7. Newtonā™s method for the approximation of the root of

f (x) = eā’x ā’ Ī· = 0. The stopping test is based on the control of the incre-

ment

6.6 Post-processing Techniques for Iterative

Methods

We conclude this chapter by introducing two algorithms that aim at ac-

celerating the convergence of iterative methods for ļ¬nding the roots of a

function.

6.6.1 Aitkenā™s Acceleration

We describe this technique in the case of linearly convergent ļ¬xed-point

methods, referring to [IK66], pp. 104ā“108, for the case of methods of higher

order.

Consider a ļ¬xed-point iteration that is linearly converging to a zero Ī± of

a given function f . Denoting by Ī» an approximation of Ļ (Ī±) to be suitably

determined and recalling (6.18) we have, for k ā„ 1

x(k) ā’ Ī»x(kā’1) x(k) ā’ Ī»x(k) + Ī»x(k) ā’ Ī»x(kā’1)

Ī± =

1ā’Ī» 1ā’Ī» (6.32)

Ī»

(x(k) ā’ x(kā’1) ).

= x(k) +

1ā’Ī»

Aitkenā™s method provides a simple way of computing Ī» that is able to

accelerate the convergence of the sequence {x(k) } to the root Ī±. With this

aim, let us consider for k ā„ 2 the following ratio

x(k) ā’ x(kā’1)

(k)

Ī» = (kā’1) , (6.33)

ā’ x(kā’2)

x

and check that

lim Ī»(k) = Ļ (Ī±). (6.34)

kā’ā

6.6 Post-processing Techniques for Iterative Methods 273

Indeed, for k suļ¬ciently large

x(k+2) ā’ Ī± Ļ (Ī±)(x(k+1) ā’ Ī±)

and thus, elaborating (6.33), we get

x(k) ā’ x(kā’1) (x(k) ā’ Ī±) ā’ (x(kā’1) ā’ Ī±)

lim Ī»(k) = lim = lim (kā’1)

kā’ā x(kā’1) ā’ x(kā’2) ā’ Ī±) ā’ (x(kā’2) ā’ Ī±)

kā’ā (x

kā’ā

x(k) ā’ Ī±

ā’1

Ļ (Ī±) ā’ 1

x(kā’1) ā’ Ī±

= lim = = Ļ (Ī±)

ā’Ī± 1

(kā’2)

x

kā’ā

1ā’

1 ā’ (kā’1) Ļ (Ī±)

ā’Ī±

x

which is (6.34). Substituting in (6.32) Ī» with its approximation Ī»(k) given

by (6.33), yields the updated estimate of Ī±

Ī»(k)

(x(k) ā’ x(kā’1) )

(k)

Ī± x + (6.35)

1ā’Ī» (k)

which, rigorously speaking, is signiļ¬cant only for a suļ¬ciently large k.

However, assuming that (6.35) holds for any k ā„ 2, we denote by x(k) the

new approximation of Ī± that is obtained by plugging (6.33) back into (6.35)

(x(k) ā’ x(kā’1) )2

ā’ (k) k ā„ 2.

(k) (k)

x =x , (6.36)

(x ā’ x(kā’1) ) ā’ (x(kā’1) ā’ x(kā’2) )

This relation is known as Aitkenā™s extrapolation formula.

Letting, for k ā„ 2,

x(k) = x(k) ā’ x(kā’1) , x(k+1) ā’

2 (k)

( x(k) ) = x(k) ,

x =

formula (6.36) can be written as

( x(k) )2

x(k) = x(k) ā’ k ā„ 2.

, (6.37)

2 x(kā’1)

Form (6.37) explains the reason why method (6.36) is more commonly

known as Aitkenā™s 2 method.

For the convergence analysis of Aitkenā™s method, it is useful to write (6.36)

as a ļ¬xed-point method in the form (6.17), by introducing the iteration

function

xĻ(Ļ(x)) ā’ Ļ2 (x)

Ļ (x) = . (6.38)

Ļ(Ļ(x)) ā’ 2Ļ(x) + x

This function is indeterminate at x = Ī± since Ļ(Ī±) = Ī±; however, by

applying Lā™Hospitalā™s rule one can easily check that limxā’Ī± Ļ (x) = Ī±

274 6. Rootļ¬nding for Nonlinear Equations

under the assumption that Ļ is diļ¬erentiable at Ī± and Ļ (Ī±) = 1. Thus, Ļ

is consistent and has a continuos extension at Ī±, the same being also true

if Ī± is a multiple root of f . Moreover, it can be shown that the ļ¬xed points

of (6.38) coincide with those of Ļ even in the case where Ī± is a multiple

root of f (see [IK66], pp. 104-106).

From (6.38) we conclude that Aitkenā™s method can be applied to a ļ¬xed-

point method x = Ļ(x) of arbitrary order. Actually, the following conver-

gence result holds.

Property 6.7 (convergence of Aitkenā™s method) Let x(k+1) = Ļ(x(k) )

be a ļ¬xed-point iteration of order p ā„ 1 for the approximation of a simple

zero Ī± of a function f . If p = 1, Aitkenā™s method converges to Ī± with order

2, while if p ā„ 2 the convergence order is 2p ā’ 1. In particular, if p = 1,

Aitkenā™s method is convergent even if the ļ¬xed-point method is not. If Ī±

has multiplicity m ā„ 2 and the method x(k+1) = Ļ(x(k) ) is ļ¬rst-order con-

vergent, then Aitkenā™s method converges linearly, with convergence factor

C = 1 ā’ 1/m.

Example 6.10 Consider the computation of the simple zero Ī± = 1 for the func-

tion f (x) = (x ā’ 1)ex . For this, we use three ļ¬xed-point methods whose iteration

functions are, respectively, Ļ0 (x) = log(xex ), Ļ1 (x) = (ex + x)/(ex + 1) and

Ļ2 (x) = (x2 ā’ x + 1)/x (for x = 0). Notice that, since |Ļ0 (1)| = 2, the corre-

sponding ļ¬xed-point method is not convergent, while in the other two cases the

methods have order 1 and 2, respectively.

Let us check the performance of Aitkenā™s method, running Program 55 with

x = 2, toll = 10ā’10 and working in complex arithmetic. Notice that in the case

(0)

of Ļ0 this produces complex numbers if x(k) happens to be negative. According

to Property 6.7, Aitkenā™s method applied to the iteration function Ļ0 converges

in 8 steps to the value x(8) = 1.000002 + i 0.000002. In the other two cases, the

method of order 1 converges to Ī± in 18 iterations, to be compared with the 4

iterations required by Aitkenā™s method, while in the case of the iteration function

Ļ2 convergence holds in 7 iterations against 5 iterations required by Aitkenā™s

ā¢

method.

Aitkenā™s method is implemented in Program 55. The input/output pa-

rameters are the same as those of previous programs in this chapter.

Program 55 - aitken : Aitkenā™s extrapolation

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

nit=0; xvect=[x0]; x=x0; fxn=eval(fun);

fx=[fxn]; xdif=[]; err=toll+1;

while err >= toll & nit <= nmax

nit=nit+1; xv=xvect(nit); x=xv; phix=eval(phi);

x=phix; phixx=eval(phi); den=phixx-2*phix+xv;

if den == 0, err=toll*1.e-01;

else, xn=(xv*phixx-phixĖ2)/den; xvect=[xvect; xn];

xdif=[xdif; abs(xn-xv)]; x=xn; fxn=abs(eval(fun));

6.6 Post-processing Techniques for Iterative Methods 275

fx=[fx; fxn]; err=fxn;

end

end

6.6.2 Techniques for Multiple Roots

As previously noticed in deriving Aitkenā™s acceleration, taking the incre-

mental ratios of successive iterates Ī»(k) in (6.33) provides a way to estimate

the asymptotic convergence factor Ļ (Ī±).

This information can be employed also to estimate the multiplicity of the

root of a nonlinear equation and, as a consequence, it provides a tool for

modifying Newtonā™s method in order to recover its quadratic convergence

(see (6.23)). Indeed, deļ¬ne the sequence m(k) through the relation Ī»(k) =

1 ā’ 1/m(k) , and recalling (6.22), it follows that m(k) tends to m as k ā’ ā.

If the multiplicity m is known a priori, it is clearly convenient to use the

modiļ¬ed Newton method (6.23). In other cases, the following adaptive New-

ton algorithm can be used

f (x(k) )

ā’m k ā„ 2,

(k+1) (k) (k)

x =x , (6.39)

f (x(k) )

where we have set

x(kā’1) ā’ x(kā’2)

1

(k)

m = = (kā’1) . (6.40)

1 ā’ Ī»(k) ā’ x(k) ā’ x(kā’2)

2x

Example 6.11 Let us check the performances of Newtonā™s method in its three

versions proposed so far (standard (6.16), modiļ¬ed (6.23) and adaptive (6.39)),

to approximate the multiple zero Ī± = 1 of the function f (x) = (x2 ā’ 1)p log x

(for p ā„ 1 and x > 0). The desired root has multiplicity m = p + 1. The values

p = 2, 4, 6 have been considered and x(0) = 0.8, toll=10ā’10 have always been

taken in numerical computations.

The obtained results are summarized in Table 6.8, where for each method

the number of iterations nit required to converge are reported. In the case of

the adaptive method, beside the value of nit we have also shown in braces the

ā¢

estimate m(nit ) of the multiplicity m that is yielded by Program 56.

m standard adaptive modiļ¬ed

3 51 13 (2.9860) 4

5 90 16 (4.9143) 5

7 127 18 (6.7792) 5

TABLE 6.8. Solution of problem (x2 ā’ 1)p log x = 0 in the interval [0.5, 1.5], with

p = 2, 4, 6

276 6. Rootļ¬nding for Nonlinear Equations

In Example 6.11, the adaptive Newton method converges more rapidly

than the standard method, but less rapidly than the modiļ¬ed Newton

method. It must be noticed, however, that the adaptive method yields as

a useful by-product a good estimate of the multiplicity of the root, which

can be proļ¬tably employed in a deļ¬‚ation procedure for the approximation

of the roots of a polynomial.

The algorithm 6.39, with the adaptive estimate (6.40) of the multiplicity

of the root, is implemented in Program 56. To avoid the onset of numerical

instabilities, the updating of m(k) is performed only when the variation be-

tween two consecutive iterates is suļ¬ciently diminished. The input/output

parameters are the same as those of previous programs in this chapter.

Program 56 - adptnewt : Adaptive Newtonā™s method

function [xvect,xdif,fx,nit,m] = adptnewt(x0,nmax,toll,fun,dfun)

xvect=x0; nit=0; r=[1]; err=toll+1; m=[1]; xdif=[];

while (nit < nmax) & (err > toll)

nit=nit+1; x=xvect(nit); fx(nit)=eval(fun); f1x=eval(dfun);

if (f1x == 0), disp(ā™ Stop due to vanishing derivative ā™); return; end;

x=x-m(nit)*fx(nit)/f1x; xvect=[xvect;x]; fx=[fx;eval(fun)];

rd=err; err=abs(xvect(nit+1)-xvect(nit)); xdif=[xdif;err];

ra=err/rd; r=[r;ra]; diļ¬=abs(r(nit+1)-r(nit));

if (diļ¬ < 1.e-3) & (r(nit+1) > 1.e-2),

m(nit+1)=max(m(nit),1/abs(1-r(nit+1)));

else, m(nit+1)=m(nit); end

end

6.7 Applications

We apply iterative methods for nonlinear equations considered so far in the

solution of two problems arising in the study of the thermal properties of

gases and electronics, respectively.

6.7.1 Analysis of the State Equation for a Real Gas

For a mole of a perfect gas, the state equation P v = RT establishes a re-

lation between the pressure P of the gas (in Pascals [P a]), the speciļ¬c vol-

ume v (in cubic meters per kilogram [m3 Kg ā’1 ]) and its temperature T (in

Kelvin [K]), R being the universal gas constant, expressed in [JKg ā’1 K ā’1 ]

(joules per kilogram per Kelvin).

For a real gas, the deviation from the state equation of perfect gases is

due to van der Waals and takes into account the intermolecular interaction

and the space occupied by molecules of ļ¬nite size (see [Sla63]).

6.7 Applications 277

Denoting by Ī± and Ī² the gas constants according to the van der Waals

model, in order to determine the speciļ¬c volume v of the gas, once P and

T are known, we must solve the nonlinear equation

f (v) = (P + Ī±/v 2 )(v ā’ Ī²) ā’ RT = 0. (6.41)

With this aim, let us consider Newtonā™s method (6.16) in the case of carbon

dioxide (CO2 ), at the pressure of P = 10[atm] (equal to 1013250[P a]) and

at the temperature of T = 300[K]. In such a case, Ī± = 188.33[P a m6 Kg ā’2 ]

and Ī² = 9.77 Ā· 10ā’4 [m3 Kg ā’1 ]; as a comparison, the solution computed by

assuming that the gas is perfect is v 0.056[m3 Kg ā’1 ].

Ė

We report in Table 6.9 the results obtained by running Program 50 for

diļ¬erent choices of the initial guess v (0) . We have denoted by Nit the num-

ber of iterations needed by Newtonā™s method to converge to the root v ā— of

f (v) = 0 using an absolute tolerance equal to the roundoļ¬ unit.

v (0) v (0) v (0) v (0)

Nit Nit Nit Nit

10ā’4 10ā’2 10ā’3 10ā’1

47 7 21 5

TABLE 6.9. Convergence of Newtonā™s method to the root of equation (6.41)

The computed approximation of v ā— is v Nit 0.0535. To analyze the causes

of the strong dependence of Nit on the value of v (0) , let us examine the

derivative f (v) = P ā’ Ī±v ā’2 + 2Ī±Ī²v ā’3 . For v > 0, f (v) = 0 at vM

1.99 Ā· 10ā’3 [m3 Kg ā’1 ] (relative maximum) and at vm 1.25 Ā· 10ā’2 [m3 Kg ā’1 ]

(relative minimum), as can be seen in the graph of Figure 6.8 (left).

A choice of v (0) in the interval (0, vm ) (with v (0) = vM ) thus necessarily

leads to a slow convergence of Newtonā™s method, as demonstrated in Figure

6.8 (right), where, in solid circled line, the sequence {|v (k+1) ā’ v (k) |} is

shown, for k ā„ 0.

A possible remedy consists of resorting to a polyalgorithmic approach,

based on the sequential use of the bisection method and Newtonā™s method

(see Section 6.2.1). Running the bisection-Newtonā™s method with the end-

points of the search interval equal to a = 10ā’4 [m3 Kg ā’1 ] and b = 0.1[m3 Kg ā’1 ]

and an absolute tolerance of 10ā’3 [m3 Kg ā’1 ], yields an overall convergence

of the algorithm to the root v ā— in 11 iterations, with an accuracy of the

order of the roundoļ¬ unit. The plot of the sequence {|v (k+1) ā’ v (k) |}, for

k ā„ 0, is shown in solid and starred lines in Figure 6.8 (right).

6.7.2 Analysis of a Nonlinear Electrical Circuit

Let us consider the electrical circuit in Figure 6.9 (left), where v and j

denote respectively the voltage drop across the device D (called a tunneling

diode) and the current ļ¬‚owing through D, while R and E are a resistor and

a voltage generator of given values.

278 6. Rootļ¬nding for Nonlinear Equations

4

x 10 0

6 10

ā’2

10

4

ā’4

10

ā’6

2 10

ā’8

10

0

ā’10

10

ā’12

ā’2 10

ā’14

10

ā’4

ā’16

10

ā’18

ā’6 10

0 0.02 0.04 0.06 0.08 0.1 0 10 20 30 40 50

FIGURE 6.8. Graph of the function f in (6.41) (left); increments |v (k+1) ā’ v (k) |

computed by the Newtonā™s method (circled curve) and bisection-Newtonā™s

method (starred curve)

The circuit is commonly employed as a biasing circuit for electronic de-

vices working at high frequency (see [Col66]). In such applications the pa-

rameters R and E are designed in such a way that v attains a value internal

to the interval for which g (v) < 0, where g is the function which describes

the bound between current and voltage for D and is drawn in Figure 6.9

(right). Explicitly, g = Ī±(ev/Ī² ā’ 1) ā’ Āµv(v ā’ Ī³), for suitable constants Ī±, Ī²,

Ī³ and Āµ.

j ā’5

x 10

12

g(v)

10

R 8

D v 6

4

2

+

E 0

_

ā’2

ā’4

0 0.1 0.2 0.3 0.4 0.5

FIGURE 6.9. Tunneling diode circuit (left) and working point computation

(right)

Our aim is to determine the working point of the circuit at hand, that

is, the values attained by v and j for given parameters R and E. For that,

we write Kirchhoļ¬ā™s law for the voltages across the loop, obtaining the

following nonlinear equation

1 E

ā’ Āµv 2 + Ī±(ev/Ī² ā’ 1) ā’

f (v) = v + ĀµĪ³ = 0. (6.42)

R R

From a graphical standpoint, ļ¬nding out the working point of the circuit

amounts to determining the intersection between the function g and the

6.8 Exercises 279

straight line of equation j = (E ā’ v)/R, as shown in Figure 6.9 (right).

Assume the following (real-life) values for the parameters of the problem:

E/R = 1.2Ā·10ā’4 [A], Ī± = 10ā’12 [A], Ī² ā’1 = 40 [V ā’1 ], Āµ = 10ā’3 [AV ā’2 ] and

Ī³ = 0.4 [V ]. The solution of (6.42), which is also unique for the considered

values of the parameters, is v ā— 0.3 [V ].

To approximate v ā— , we compare the main iterative methods introduced

in this chapter. We have taken v (0) = 0 [V ] for Newtonā™s method, Ī¾ = 0

for the Dekker-Brent algorithm (for the meaning of Ī¾, see Example 6.5),

while for all the other schemes the search interval has been taken equal to

[0, 0.5]. The stopping tolerance toll has been set to 10ā’10 . The obtained

results are reported in Table 6.10 where nit and f (nit) denote respectively

the number of iterations needed by the method to converge and the value

of f at the computed solution.

Notice the extremely slow convergence of the Regula Falsi method, due

to the fact that the value v (k ) always coincides with the right end-point

v = 0.5 and the function f around v ā— has derivative very close to zero. An

analogous interpretation holds for the chord method.

f (nit) f (nit)

Method Method

nit nit

ā’1.12 Ā· 10ā’15 1.09 Ā· 10ā’14

bisection 33 Dekker-Brent 11

ā’9.77 Ā· 10ā’11 2.7 Ā· 10ā’20

Regula Falsi 225 secant 11

ā’9.80 Ā· 10ā’14 ā’1.35 Ā· 10ā’20

chord 186 Newtonā™s 8

TABLE 6.10. Convergence of the methods for the approximation of the root of

equation (6.42)

6.8 Exercises

1. Derive geometrically the sequence of the ļ¬rst iterates computed by bisec-

tion, Regula Falsi, secant and Newtonā™s methods in the approximation of

the zero of the function f (x) = x2 ā’ 2 in the interval [1, 3].

2. Let f be a continuous function that is m-times diļ¬erentiable (m ā„ 1), such

that f (Ī±) = . . . = f (mā’1) (Ī±) = 0 and f (m) (Ī±) = 0. Prove (6.22) and check

that the modiļ¬ed Newton method (6.23) has order of convergence equal to

2.

[Hint: let f (x) = (x ā’ Ī±)m h(x), h being a function such that h(Ī±) = 0].

3. Let f (x) = cos2 (2x) ā’ x2 be the function in the interval 0 ā¤ x ā¤ 1.5

examined in Example 6.4. Having ļ¬xed a tolerance Īµ = 10ā’10 on the abso-

lute error, determine experimentally the subintervals for which Newtonā™s

method is convergent to the zero Ī± 0.5149.

[Solution: for 0 < x(0) ā¤ 0.02, 0.94 ā¤ x(0) ā¤ 1.13 and 1.476 ā¤ x(0) ā¤ 1.5,

the method converges to the solution ā’Ī±. For any other value of x(0) in

[0, 1.5], the method converges to Ī±].

280 6. Rootļ¬nding for Nonlinear Equations

4. Check the following properties:

(a) 0 < Ļ (Ī±) < 1: monotone convergence, that is, the error x(k) ā’ Ī±

maintains a constant sign as k varies;

(b) ā’1 < Ļ (Ī±) < 0: oscillatory convergence that is, x(k) ā’ Ī± changes sign

as k varies;

(c) |Ļ (Ī±)| > 1: divergence. More precisely, if Ļ (Ī±) > 1, the sequence

is monotonically diverging, while for Ļ (Ī±) < ā’1 it diverges with

oscillatory sign.

5. Consider for k ā„ 0 the ļ¬xed-point method, known as Steļ¬ensenā™s method

f (x(k) + f (x(k) )) ā’ f (x(k) )

f (x(k) )

x(k+1) = x(k) ā’ Ļ•(x(k) ) =

, ,

Ļ•(x(k) ) f (x(k) )

and prove that it is a second-order method. Implement the Steļ¬ensen

method in a MATLAB code and employ it to approximate the root of

the nonlinear equation eā’x ā’ sin(x) = 0.

6. Analyze the convergence of the ļ¬xed-point method x(k+1) = Ļj (x(k) ) for

computing the zeros Ī±1 = ā’1 and Ī±2 = 2 of the function f (x) = x2 ā’ x ā’ 2,

when the following iteration functions are used: Ļ1 (x) = x2 ā’ 2, Ļ2 (x) =

ā ā

2 + x Ļ3 (x) = ā’ 2 + x and Ļ4 (x) = 1 + 2/x, x = 0.

[Solution: the method is non convergent with Ļ1 , it converges only to Ī±2 ,

with Ļ2 and Ļ4 , while it converges only to Ī±1 with Ļ3 ].

7. For the approximation of the zeros of the function f (x) = (2x2 ā’ 3x ā’

2)/(x ā’ 1), consider the following ļ¬xed-point methods:

(1) x(k+1) = g(x(k) ), where g(x) = (3x2 ā’ 4x ā’ 2)/(x ā’ 1);

(2) x(k+1) = h(x(k) ), where h(x) = x ā’ 2 + x/(x ā’ 1).

Analyze the convergence properties of the two methods and determine

in particular their order. Check the behavior of the two schemes using

Program 51 and provide, for the second method, an experimental estimate

of the interval such that if x(0) is chosen in the interval then the method

converges to Ī± = 2.

[Solution: zeros: Ī±1 = ā’1/2 and Ī±2 = 2. Method (1) is not convergent,

while (2) can approximate only Ī±2 and is second-order. Convergence holds

for any x(0) > 1].

8. Propose at least two ļ¬xed-point methods for approximating the root Ī±

0.5885 of equation eā’x ā’ sin(x) = 0 and analyze their convergence.

9. Using Descartesā™s rule of signs, determine the number of real roots of the

polynomials p6 (x) = x6 ā’ x ā’ 1 and p4 (x) = x4 ā’ x3 ā’ x2 + x ā’ 1.

[Solution: both p6 and p4 have one negative and one positive real root].

ā

10. Let g : R ā’ R be deļ¬ned as g(x) = 1 + x2 . Show that the iterates of

Newtonā™s method for the equation g (x) = 0 satisfy the following proper-

ties:

|x(0) | < 1 ā’ g(x(k+1) ) < g(x(k) ), k ā„ 0, lim x(k) = 0,

(a)

kā’ā

|x(0) | > 1 ā’ g(x(k+1) ) > g(x(k) ), k ā„ 0, lim |x(k) | = +ā.

(b)

kā’ā

7

Nonlinear Systems and Numerical

Optimization

ńņš. 11 |