2.4 Newton-Raphson and Secant Methods: Slope Methods For Finding Roots
2.4 Newton-Raphson and Secant Methods: Slope Methods For Finding Roots
2.4 Newton-Raphson and Secant Methods: Slope Methods For Finding Roots
C HAP. 2
(1)
0 f ( p0 )
,
p1 p0
(2)
which is the slope at the point ( p0 , f ( p0 )). Equating the values of the slope m in
equations (1) and (2) and solving for p1 results in
p1 = p0
(3)
f ( p0 )
.
f ( p0 )
y = f (x)
p1
p0
x
p2
(p1, f (p1))
(p0 , f(p0))
S EC . 2.4
71
pk = g( pk1 ) = pk1
f ( pk1 )
f ( pk1 )
for k = 1, 2, . . .
g(x) = x
f (x)
f (x)
f (x) = f ( p0 ) + f ( p0 )(x p0 ) +
f (c)(x p0 )2
,
2!
where c lies somewhere between p0 and x. Substituting x = p into equation (6) and
using the fact that f ( p) = 0 produces
(7)
0 = f ( p0 ) + f ( p0 )( p p0 ) +
f (c)( p p0 )2
.
2!
If p0 is close enough to p, the last term on the right side of (7) will be small compared to the sum of the first two terms. Hence it can be neglected and we can use the
approximation
(8)
0 f ( p0 ) + f ( p0 )( p p0 ).
p1 = p0
f ( p0 )
.
f ( p0 )
When pk1 is used in place of p0 in equation (9), the general rule (4) is established. For
most applications this is all that needs to be understood. However, to fully comprehend
72
C HAP. 2
what is happening, we need to consider the fixed-point iteration function and apply
Theorem 2.2 in our situation. The key is in the analysis of g (x):
g (x) = 1
| f (x) f (x)|
<1
| f (x)|2
for all x ( p , p + ).
is a real number and let p0 > 0 be an initial approximation to A. Define the sequence
{ pk }
k=0 using the recursive rule
pk1 +
(11)
pk =
A
pk1
for k = 1, 2, . . . .
A; that is, limn pk =
A.
2
Outline of Proof. Start with the
the roots of
function f (x) = x A, and notice that
2
the equation x A = 0 are A. Now use f (x) and the derivative f (x) in formula
(5) and write down the Newton-Raphson iteration formula
(12)
g(x) = x
x2 A
f (x)
=x
.
f (x)
2x
(13)
g(x) =
x+
A
x .
2
When g(x) in (13) is used to define the recursive iteration in (4), the result is formula
(11). It can be proved that the sequence that is generated in (11) will converge for any
starting value p0 > 0. The details are left for the exercises.
An important point of Corollary 2.2 is the fact that the iteration function g(x)
involved only the arithmetic operations +, , , and /. If g(x) had involved the calculation of a square root, we would be caught in the circular reasoning that being able
to calculate the
square root would permit you to recursively define a sequence that will
converge to A. For this reason, f (x) = x 2 A was chosen, because it involved only
the arithmetic operations.
S EC . 2.4
73
Now let us turn to a familiar problem from elementary physics and see why determining the location of a root is an important task. Suppose that a projectile is fired
from the origin with an angle of elevation b0 and initial velocity v0 . In elementary
courses, air resistance is neglected and we learn that the height y = y(t) and the distance traveled x = x(t), measured in feet, obey the rules
(14)
y = v y t 16t 2
and
x = vx t,
where the horizontal and vertical components of the initial velocity are vx = v0 cos(b0 )
and v y = v0 sin(b0 ), respectively. The mathematical model expressed by the rules
in (14) is easy to work with, but tends to give too high an altitude and too long a range
for the projectiles path. If we make the additional assumption that the air resistance is
proportional to the velocity, the equations of motion become
(15)
y = f (t) = (Cv y + 32C 2 ) 1 et/C 32Ct
and
(16)
x = r (t) = Cvx 1 et/C ,
where C = m/k and k is the coefficient of air resistance and m is the mass of the
projectile. A larger value of C will result in a higher maximum altitude and a longer
range for the projectile. The graph of a flight path of a projectile when air resistance is
considered is shown in Figure 2.14. This improved model is more realistic, but requires
the use of a root-finding algorithm for solving f (t) = 0 to determine the elapsed time
until the projectile hits the ground. The elementary model in (14) does not require a
sophisticated procedure to find the elapsed time.
74
C HAP. 2
y
300
200
100
0
200
Table 2.4
k
0
1
2
3
4
400
600
800
1000
pk+1 pk
0.79773101
0.05530160
0.00025475
0.00000001
0.00000000
Height, f ( pk )
83.22097200
6.68369700
0.03050700
0.00000100
0.00000000
83.22097200
= 8.797731010.
104.3220972
Division-by-Zero Error
One obvious pitfall of the Newton-Raphson method is the possibility of division by
zero in formula (4), which would occur if f ( pk1 ) = 0. Program 2.5 has a procedure
S EC . 2.4
75
to check for this situation, but what use is the last calculated approximation pk1 in
this case? It is quite possible that f ( pk1 ) is sufficiently close to zero and that pk1
is an acceptable approximation to the root. We now investigate this situation and will
uncover an interesting fact, that is, how fast the iteration converges.
Definition 2.4. Assume that f (x) and its derivatives f (x), . . . , f (M) (x) are defined
and continuous on an interval about x = p. We say that f (x) = 0 has a root of order
M at x = p if and only if
(17) f ( p) = 0,
f ( p) = 0,
...,
f (M1) ( p) = 0,
and
f (M) ( p) = 0.
f (x) = (x p) M h(x),
where h( p) = 0.
Speed of Convergence
The distinguishing property we seek is the following. If p is a simple root of f (x) = 0,
Newtons method will converge rapidly, and the number of accurate decimal places
(roughly) doubles with each iteration. On the other hand, if p is a multiple root, the
error in each successive approximation is a fraction of the previous error. To make
this precise, we define the order of convergence. This is a measure of how rapidly a
sequence converges.
Definition 2.5. Assume that { pn }
n=0 converges to p and set E n = p pn for n 0.
If two positive constants A = 0 and R > 0 exist, and
(19)
| p pn+1 |
|E n+1 |
= lim
= A,
R
n | p pn |
n |E n | R
lim
then the sequence is said to converge to p with order of convergence R. The number A is called the asymptotic error constant. The cases R = 1, 2 are given special
76
C HAP. 2
Table 2.5
pk
pk+1 pk
E k = p pk
0
1
2
3
4
2.400000000
2.076190476
2.003596011
2.000008589
2.000000000
0.323809524
0.072594465
0.003587422
0.000008589
0.000000000
0.400000000
0.076190476
0.003596011
0.000008589
0.000000000
|E k+1 |
|E k |2
0.476190475
0.619469086
0.664202613
consideration.
(20)
If R = 1, the convergence of { pn }
n=0 is called linear.
If R = 2, the convergence of { pn }
n=0 is called quadratic.
If R is large, the sequence { pn } converges rapidly to p; that is, relation (19) implies
that for large values of n we have the approximation |E n+1 | A|E n | R . For example,
suppose that R = 2 and |E n | 102 ; then we would expect that |E n+1 | A 104 .
Some sequences converge at a rate that is not an integer,
and we will see that the
(21)
3
2
2 pk1
2
3 pk1
3
Using formula (19) with R = 2 to check for quadratic convergence, we get the values in
Table 2.5.
A detailed look at the rate of convergence in Example 2.14 will reveal that the error
in each successive iteration is proportional to the square of the error in the previous
iteration. That is,
| p pk+1 | A| p pk |2 ,
where A 2/3. To check this, we use
| p p3 | = 0.000008589
and
| p p2 |2 = |0.003596011|2 = 0.000012931
2
| p p2 |2 .
3
S EC . 2.4
Table 2.6
77
pk
0
1
2
3
4
5
..
.
1.200000000
1.103030303
1.052356420
1.026400811
1.013257730
1.006643419
..
.
pk+1 pk
E k = p pk
0.096969697
0.050673883
0.025955609
0.013143081
0.006614311
0.003318055
..
.
0.200000000
0.103030303
0.052356420
0.026400811
0.013257730
0.006643419
..
.
|E k+1 |
|E k |
0.515151515
0.508165253
0.496751115
0.509753688
0.501097775
0.500550093
..
.
Example 2.15 (Linear Convergence at a Double Root). Start with p0 = 1.2 and use
Newton-Raphson iteration to find the double root p = 1 of the polynomial f (x) = x 3
3x + 2. Using formula (20) to check for linear convergence, we get the values in Table 2.6.
Notice that the Newton-Raphson method is converging to the double root, but at
a slow rate. The values of f ( pk ) in Example 2.15 go to zero faster than the values
of f ( pk ), so the quotient f ( pk )/ f ( pk ) in formula (4) is defined when pk = p.
The sequence is converging linearly, and the error is decreasing by a factor of approximately 1/2 with each successive iteration. The following theorem summarizes the
performance of Newtons method on simple and double roots.
Theorem 2.6 (Convergence Rate for Newton-Raphson Iteration). Assume that
Newton-Raphson iteration produces a sequence { pn }
n=0 that converges to the root p
of the function f (x). If p is a simple root, convergence is quadratic and
(22)
|E n+1 |
| f ( p)|
|E n |2
2| f ( p)|
|E n+1 |
M 1
|E n |
M
Pitfalls
The division-by-zero error was easy to anticipate, but there are other difficulties that
are not so easy to spot. Suppose that the function is f (x) = x 2 4x + 5; then the
sequence { pk } of real numbers generated by formula (4) will wander back and forth
from left to right and not converge. A simple analysis of the situation reveals that
f (x) > 0 and has no real roots.
78
C HAP. 2
y = xe
-x
0.2
0.1
0.0
p0 = 2
p1 = 4
p2
x
p3
Sometimes the initial approximation p0 is too far away from the desired root and
the sequence { pk } converges to some other root. This usually happens when the slope
f ( p0 ) is small and the tangent line to the curve y = f (x) is nearly horizontal. For
example, if f (x) = cos(x) and we seek the root p = /2 and start with p0 = 3,
calculation reveals that p1 = 4.01525255, p2 = 4.85265757, . . . , and { pk } will
converge to a different root 3/2 4.71238898.
Suppose that f (x) is positive and monotone decreasing on the unbounded interval
[a, ) and p0 > a; then the sequence { pk } might diverge to +. For example, if
f (x) = xex and p0 = 2.0, then
p1 = 4.0,
p2 = 5.333333333,
...,
p15 = 19.723549434,
...,
and { pk } diverges slowly to + (see Figure 2.15(a)). This particular function has
another surprising problem. The value of f (x) goes to zero rapidly as x gets large, for
example, f ( p15 ) = 0.0000000536, and it is possible that p15 could be mistaken for a
root. For this reason we designed the stopping criterion in Program 2.5 to involve the
relative error 2| pk+1 pk |/(| pk |+106 ), and when k = 15, this value is 0.106817, so
the tolerance = 106 will help guard against reporting a false root.
Another phenomenon, cycling, occurs when the terms in the sequence { pk } tend to
repeat or almost repeat. For example, if f (x) = x 3 x 3 and the initial approximation
is p0 = 0, then the sequence is
p1 = 3.000000,
p5 = 3.000389,
p2 = 1.961538,
p6 = 1.961818,
p3 = 1.147176,
p7 = 1.147430,
p4 = 0.006579,
...
and we are stuck in a cycle where pk+4 pk for k = 0, 1, . . . (see Figure 2.15(b)).
But if the starting value p0 is sufficiently close to the root p 1.671699881, then { pk }
S EC . 2.4
p2
p3
79
p0
1
x
5
10
y = x3 x 3
15
20
25
y
1
3
p3
2
p1
p0
1
1
p2
2
1
y = arctan(x)
p2 = 1.845931751,
p3 = 2.889109054,
etc. (see Figure 2.15(c)). But if the starting value is sufficiently close to the root p = 0,
80
C HAP. 2
y = f (x)
p2
p1
p0
x
(p, 0)
(p1, f(p1))
(p0, f(p0))
p2 = 0.000335302,
p3 = 0.000000000.
The situations above point to the fact that we must be honest in reporting an answer.
Sometimes the sequence does not converge. It is not always the case that after N
iterations a solution is found. The user of a root-finding algorithm needs to be warned
of the situation when a root is not found. If there is other information concerning
the context of the problem, then it is less likely that an erroneous root will be found.
Sometimes f (x) has a definite interval in which a root is meaningful. If knowledge
of the behavior of the function or an accurate graph is available, then it is easier to
choose p0 .