<< . .

. 37
( : 95)



. . >>

and/or its derivative.


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)

<< . .

. 37
( : 95)



. . >>