6.2.1 The Bisection Method

The bisection method is based on the following property.

6.2 A Geometric Approach to Root¬nding 249

Property 6.1 (theorem of zeros for continuous functions) Given a

continuous function f : [a, b] ’ R, such that f (a)f (b) < 0, then ∃ ± ∈ (a, b)

such that f (±) = 0.

Starting from I0 = [a, b], the bisection method generates a sequence of

subintervals Ik = [a(k) , b(k) ], k ≥ 0, with Ik ‚ Ik’1 , k ≥ 1, and enjoys the

property that f (a(k) )f (b(k) ) < 0. Precisely, we set a(0) = a, b(0) = b and

x(0) = (a(0) + b(0) )/2; then, for k ≥ 0:

set a(k+1) = a(k) , b(k+1) = x(k) if f (x(k) )f (a(k) ) < 0;

set a(k+1) = x(k) , b(k+1) = b(k) if f (x(k) )f (b(k) ) < 0;

¬nally, set x(k+1) = (a(k+1) + b(k+1) )/2.

y 0

10

’2

10

f (x) ’4

10

± ’6

x 10

a b

x(1)

x(0)

’8

10

’10

10

I1 ’12

10

0 5 10 15 20 25 30

I0

FIGURE 6.1. The bisection method. The ¬rst two steps (left); convergence history

for the Example 6.3 (right). The number of iterations and the absolute error as

a function of k are reported on the x- and y-axis, respectively

The bisection iteration terminates at the m-th step for which |x(m) ’±| ¤

|Im | ¤ µ, where µ is a ¬xed tolerance and |Im | is the length of Im . As for

the speed of convergence of the bisection method, notice that |I0 | = b ’ a,

while

|Ik | = |I0 |/2k = (b ’ a)/2k , k ≥ 0. (6.8)

Denoting by e(k) = x(k) ’ ± the absolute error at step k, from (6.8) it

follows that |e(k) | ¤ (b ’ a)/2k , k ≥ 0, which implies limk’∞ |e(k) | = 0.

The bisection method is therefore globally convergent. Moreover, to get

|x ’ ±| ¤ µ we must take

(m)

log((b ’ a)/µ) log((b ’ a)/µ)

m ≥ log2 (b ’ a) ’ log2 (µ) = . (6.9)

log(2) 0.6931

In particular, to gain a signi¬cant ¬gure in the accuracy of the approxi-

mation of the root (that is, to have |x(k) ’ ±| = |x(j) ’ ±|/10), one needs

250 6. Root¬nding for Nonlinear Equations

k ’ j = log2 (10) 3.32 bisections. This singles out the bisection method as

an algorithm of certain, but slow, convergence. We must also point out that

the bisection method does not generally guarantee a monotone reduction

of the absolute error between two successive iterations, that is, we cannot

ensure a priori that

|e(k+1) | ¤ Mk |e(k) |, for any k ≥ 0, (6.10)

with Mk < 1. For this purpose, consider the situation depicted in Figure

6.1 (left), where clearly |e(1) | > |e(0) |. Failure to satisfy (6.10) does not

allow for qualifying the bisection method as a method of order 1, in the

sense of De¬nition 6.1.

Example 6.3 Let us check the convergence properties of the bisection method

in the approximation of the root ± 0.9062 of the Legendre polynomial of degree

5

x

L5 (x) = (63x4 ’ 70x2 + 15),

8

whose roots lie within the interval (’1, 1) (see Section 10.1.2). Program 46 has

been run taking a = 0.6, b = 1 (whence, L5 (a) · L5 (b) < 0), nmax = 100,

toll = 10’10 and has reached convergence in 32 iterations, this agrees with

the theoretical estimate (6.9) (indeed, m ≥ 31.8974). The convergence history is

reported in Figure 6.1 (right) and shows an (average) reduction of the error by a

factor of two, with an oscillating behavior of the sequence {x(k) }. •

The slow reduction of the error suggests employing the bisection method

as an “approaching” technique to the root. Indeed, taking few bisection

steps, a reasonable approximation to ± is obtained, starting from which a

higher order method can be successfully used for a rapid convergence to

the solution within the ¬xed tolerance. An example of such a procedure

will be addressed in Section 6.7.1.

The bisection algorithm is implemented in Program 46. The input pa-

rameters, here and in the remainder of this chapter, have the following

meaning: a and b denote the end points of the search interval, fun is the

variable containing the expression of the function f , toll is a ¬xed toler-

ance and nmax is the maximum admissible number of steps for the iterative

process.

In the output vectors xvect, xdif and fx the sequences {x(k) }, {|x(k+1) ’

x(k) |} and {f (x(k) )}, for k ≥ 0, are respectively stored, while nit denotes

the number of iterations needed to satisfy the stopping criteria. In the case

of the bisection method, the code returns as soon as the half-length of the

search interval is less than toll.

Program 46 - bisect : Bisection method

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

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

while (nit < nmax & err > toll)

6.2 A Geometric Approach to Root¬nding 251

nit=nit+1; c=(a+b)/2; x=c; fc=eval(fun); xvect=[xvect;x];

fx=[fx;fc]; x=a; if (fc*eval(fun) > 0), a=c; else, b=c; end;

err=abs(b-a); xdif=[xdif;err];

end;

6.2.2 The Methods of Chord, Secant and Regula Falsi and

Newton™s Method

In order to devise algorithms with better convergence properties than the

bisection method, it is necessary to include information from the values

attained by f and, possibly, also by its derivative f (if f is di¬erentiable)

or by a suitable approximation.

For this purpose, let us expand f in a Taylor series around ± and truncate

the expansion at the ¬rst order. The following linearized version of problem

(6.1) is obtained

f (±) = 0 = f (x) + (± ’ x)f (ξ), (6.11)

for a suitable ξ between ± and x. Equation (6.11) prompts the following

iterative method: for any k ≥ 0, given x(k) , determine x(k+1) by solving

equation f (x(k) ) + (x(k+1) ’ x(k) )qk = 0, where qk is a suitable approxima-

tion of f (x(k) ).

The method described here amounts to ¬nding the intersection between

the x-axis and the straight line of slope qk passing through the point

(x(k) , f (x(k) )), and thus can be more conveniently set up in the form

’1

x(k+1) = x(k) ’ qk f (x(k) ), ∀k ≥ 0.

We consider below four particular choices of qk .

y

y

f (x) f (x)

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

x(0) x

x ±

x(3)

a ± b a b

FIGURE 6.2. The ¬rst step of the chord method (left) and the ¬rst three steps

of the secant method (right). For this method we set x(’1) = b and x(0) = a

252 6. Root¬nding for Nonlinear Equations

The chord method. We let

f (b) ’ f (a)

∀k ≥ 0

qk = q = ,

b’a

from which, given an initial value x(0) , the following recursive relation is

obtained

b’a

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

f (x(k) ), (6.12)

f (b) ’ f (a)

In Section 6.3.1, we shall see that the sequence {x(k) } generated by (6.12)

converges to the root ± with order of convergence p = 1.

The secant method. We let

f (x(k) ) ’ f (x(k’1) )

∀k ≥ 0

qk = , (6.13)

x(k) ’ x(k’1)

from which, giving two initial values x(’1) and x(0) , we obtain the following

relation

x(k) ’ x(k’1)

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

f (x(k) ), (6.14)

(k) ) ’ f (x(k’1) )

f (x

If compared with the chord method, the iterative process (6.14) requires

an extra initial point x(’1) and the corresponding function value f (x(’1) ),

as well as, for any k, computing the incremental ratio (6.13). The bene¬t

due to the increase in the computational cost is the higher speed of con-

vergence of the secant method, as stated in the following property which

can be regarded as a ¬rst example of the local convergence theorem (for the

proof see [IK66], pp. 99-101).

Property 6.2 Let f ∈ C 2 (J ), J being a suitable neighborhood of the root

± and assume that f (±) = 0. Then, if the initial data x(’1) and x(0) are

chosen in J su¬ciently close to ±, the sequence (6.14) converges to ± with

√

order p = (1 + 5)/2 1.63.

The Regula Falsi (or false position) method. This is a variant of

the secant method in which, instead of selecting the secant line through

the values (x(k) , f (x(k) ) and (x(k’1) , f (x(k’1) ), we take the one through

(x(k) , f (x(k) ) and (x(k ) , f (x(k ) ), k being the maximum index less than k

such that f (x(k ) ) · f (x(k) ) < 0. Precisely, once two values x(’1) and x(0)

have been found such that f (x(’1) ) · f (x(0) ) < 0, we let

x(k) ’ x(k )

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

f (x(k) ), (6.15)

(k) ) ’ f (x(k ) )

f (x

6.2 A Geometric Approach to Root¬nding 253

Having ¬xed an absolute tolerance µ, the iteration (6.15) terminates at the

m-th step such that |f (x(m) )| < µ. Notice that the sequence of indices k

is nondecreasing; therefore, in order to ¬nd at step k the new value of k ,

it is not necessary to sweep all the sequence back, but it su¬ces to stop at

the value of k that has been determined at the previous step. We show in

Figure 6.3 (left) the ¬rst two steps of (6.15) in the special case in which

x(k ) coincides with x(’1) for any k ≥ 0.

The Regula Falsi method, though of the same complexity as the secant

method, has linear convergence order (see, for example, [RR78], pp. 339-

340). However, unlike the secant method, the iterates generated by (6.15)

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)