<<

. 11
( 26)



>>

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)
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
( 26)



>>