<< . .

. 38
( : 95)



. . >>

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

<< . .

. 38
( : 95)



. . >>