Description
1 Introduction
The bisection method is very reliable, but it can be relatively slow, and it does not generalize easily to more
than one dimension. The secant and Muller’s methods are faster, but still do not generalize easily to multiple
dimensions. In this lab we will look at Newton’s method for finding roots of functions. Newton’s method
naturally generalizes to multiple dimensions and can be much faster than bisection. On the negative side,
it requires a formula for the derivative as well as the function, and it can easily fail. Nonetheless, it is a
workhorse method in numerical analysis. We will be taking three sessions to complete this lab.
2 Stopping Tests
Root finding routines check after each step to see whether the current result is good enough. The tests that
are made are called “termination conditions” or “stopping tests”. Common tests include:
Residual size |f(x)| <
Increment size |xnew − xold| <
Number of iterations: itCount > ITMAX
The size of the residual looks like a good choice because the residual is zero at the solution; however, it turns
out to be a poor choice because the residual can be small even if the iterate is far from the true solution.
You saw this when you found the root of the function f5(x) = (x − 1)5
in the previous lab.
The size of the increment is a reasonable choice because Newton’s method (usually) converges quadratically, and when it does the increment is an excellent approximation of the true error. The third stopping
criterion, when the number of iterations exceeds a maximum, is a safety measure to assure that the iteration
will always terminate in finite time.
It should be emphasized that the stopping criteria are estimated errors. In this lab, in class, and in
later labs, you will see many other expressions for estimating errors. You should not confuse the estimated
error with the true error, |xapprox − xexact|. The true error is not included among the stopping tests
because you would need to know the exact solution to use it.
3 Failure
You know that the bisection method is very reliable and rarely fails but always takes a (sometimes large)
fixed number of steps. Newton’s method works more rapidly a good deal of the time, but does fail. And
Newton’s method works in more than one dimension. One objective of this lab will be to see how Newton
can fail and get a flavor of what might be done about it.
Although we will seem to spend a lot of time looking at failures, you should still expect that Newton’s
method will converge most of the time. It is just that when it converges there is nothing special to say, but
when it fails to converge there is always the interesting questions of, “Why?” and “How can I remedy it?”
1
4 Introduction to Newton’s Method
The idea of Newton’s method is that, starting from a guessed point x0, we find the equation of the straight
line that passes through the point (x0, f(x0)) and has slope f
0
(x0). The next iterate, x1, is simply the root
of this linear equation, i.e., the location that the straight line intersects the x-axis.
A “quick and dirty” derivation of the formula can be taken from the definition of the derivative of a
function. Assuming you are at the point x
(k)
, (I am using superscripts in parentheses to denote iteration to
reduce confusion between iteration number and vector components.)
f(x
(k) + ∆x) − f(x
(k)
)
∆x
≈ f
0
(x
(k)
).
If f were linear, this approximate equality would be true equality and the next iterate, x
(k+1) would be the
exact solution, satisfying f(x
(k+1)) = 0. In other words
f(x
(k+1)) = f(x
(k) + ∆x) = 0.
This yields
−f(x
(k)
)
∆x
≈ f
0
(x
(k)
),
or
∆x = x
(k+1) − x
(k) = −
f(x
(k)
)
f
0(x
(k))
. (1)
As you know, Newton’s method also will work for vectors, so long as the derivative is properly handled.
Assume that x and f are n-vectors. Then the Jacobian matrix is the matrix J whose elements are
Jij =
∂fi
∂xj
.
(Here, i = 1 for the first row of J, i = 2 for the second row of J, etc., so that the matrix subscripts correspond
to the usual ones.) Thus, Newton’s method for vector functions can be written
x
(k+1) − x
(k) = −J(x
(k)
)
−1
f(x
(k)
). (2)
5 Writing Python code for functions
Newton’s method requires both the function value and its derivative, unlike the bisection method that requires
only the function value. To our existing files that defined a function, we can add extra code that returns the
derivative. For instance, consider the file f1.py from the previous lab. Here’s how we can modify it so it can
be used for either bisection or Newton’s method:
def f 1 ( x ) :
”””Computes y=x ˆ2−9. ”””
y=x∗∗2
return y
def fp 1 ( x ) :
”””Computes t h e d e r i v a t i v e o f y=x ˆ2−9. ”””
yprime=2∗x
return yprime
Now, when we set up a problem, if we just need the function value, we will use a statement like
from f 1 import f 1
2
but for Newton’s method, where we need both function and derivative, we will ask for both definitions.
from f 1 import f1 , f 1p
Note that from f1 is telling us to look in the file f1.py, but import f1, f1p is telling us to look inside that
file for definitions of those names.
6 Exercise 1
In this exercise, you will modify function files from the previous lab to include the definition of the derivative.
1. Copy the files f0.py, f1.py, f2.py, f3.py and f4.py from the previous lab. Modify each file to include
derivative definitions f0p, f1p and so on. Recall that the functions are:
f0: y=x-1 yprime=???
f1: y=x^2-9 yprime=???
f2: y=x^5-x-1 yprime=???
f3: y=x*exp(-x) yprime=???
f4: y=2*cos(3*x)-exp(x) yprime=???
2. Use your new codes to compute the values of the function and derivative at x=-1 for each of these
functions.
There are several other programming strategies for computing derivatives.
• You could compute the derivative value in a new file with a new name such as df0, df1, etc.
• You could compute the derivative using some sort of divided difference formula. This method is useful
for very complicated functions.
• You could use the symbolic capabilities of the Python package sympy to symbolically differentiate the
functions. This method seems attractive, but it would greatly increase the time spent in function
evaluation.
• There are automated ways to discover the formula for a derivative from the file or symbolic formula
defining the function, and these can be used to generate the file for the derivative automatically.
7 Exercise 2
In this exercise, you will get down to the task of writing Newton’s method as a function file. In this file, you
will see how to use a variable number of arguments in a function to simplify later calls. The size of the error
and the maximum number of iterations will be optional arguments. When these arguments are omitted,
they take on their default values.
The function you are about to write is shorter than bisect.py and the instructions take you through it
step by step.
1. Create a new, empty file named newton1.py.
2. Start off the function with its signature line and follow that with comments. The comments should
repeat the signature as a comment, contain a line or two explaining what the function does, explain
what each of the parameters mean, and include your name and the date.
3
def newton1 ( func , d func , x0 , maxIts =100) :
””” Uses Newtons method t o f i n d x so t h a t func ( x ) = 0 .
Args :
func : Func t ion t o f i n d t h e r o o t o f .
d func : d e f i n e s t h e d e r i v a t i v e o f t h e f u n c t i o n .
x0 : The i n i t i a l g u e s s
maxI ts : The maximum number o f i t e r a t i o n s . D e f a ul t i s 1 0 0 .
Re turns :
x : ???
numIts : ???
”””
import numpy a s np
(You may have to complete this exercise before coming back to explain all of the parameters.)
3. From a programming standpoint, the iteration should be limited to a fixed (large) number of steps.
The reason for this is that if it does not meet some termination criterion it will run forever.
4. Define the convergence criterion
EPSILON = 5. 0 e−5
5. Initialize the solution estimate:
x = x0
6. Some people like to use while loops so that the stopping criterion appears at the beginning of the
loop. My opinion is that if you are going to count the iterations you may as well use a for statement.
Start off the loop with
fo r numIts in range ( 1 : maxIts+1) :
7. Evaluate the function and its derivative at the current point x, much as was done in bisect.py. You
may choose the variable names as you please.
8. Define a Python variable increment as the negative of the function value divided by the value of its
derivative. This is the right side of Equation (1).
9. Complete equation (1) with the statement
x += inc remen t
10. Finish the loop and the file with the following statements
e r r o r E s tim a t e = np . abs ( inc remen t ) ;
print ( f ”{numIts } , x= {x } , e r r o r e s tim a t e = { e r r o r E s tim a t e }” )
i f e r r o r E s tim a t e <EPSILON:
return x , numIts
# i f g e t here , t h e i t e r a t i o n has f a i l e d !
ra is e Excep ti on ( ” newton1 : maximum number o f i t e r a t i o n s exceeded . ” )
4
The return statement causes the function to return to its caller before all maxIts iterations are
complete. If the error estimate is not satisfied within maxIts iterations, than an exception will cause
the calculation to terminate with an error message. It is a good idea to include the name of the function
as part of the error message so you can find where the error occurred.
The print statement prints out some useful numbers during the iteration. For the first few exercises in
this lab, leave this printing in the file, but when the function is being used as part of a calculation, it
should not print extraneous information.
11. Complete the comments at the beginning of the function by replacing the “???” symbols.
12. Recall the Newton iteration should take only a single iteration when the function is linear. Test your
newton1() function on the linear function f0 that you wrote in the previous exercise. Start from an
initial guess of x=10. The default value of maxIts (100) is satisfactory, so you can omit it and use the
command
from f 0 import f0 , f 0p
from newton1 import newton1
x , numIts = newton1 ( f0 , f0p , 1 0 )
The correct solution, of course, is x=1 and it should take a single iteration. It actually takes a second
iteration in order to recognize that it has converged, so numIts should be 2. If either the solution or
the number of iterations is wrong, you have a bug: fix your code. Hint: The mistake might be in the
derivative in f0.py or in the file newton1.py.
13. Test your newton1() function on the quadratic function f1 that you wrote in the previous exercise.
Start from an initial guess of x=0.1. The default value of maxIts (100) is satsifactory, so you can omit
it. The correct solution, of course, is x=3. How many iterations are needed?
14. As you know, the theoretical convergence rate for Newton is quadratic. This means that if you compute
the quantity
r
(k)
2 =
|∆x
(k)
|
|∆x
(k−1)|
2
(the superscripts indicate iteration number) then r
(k)
2
should approach a nonzero constant as k → ∞.
You have only taken a few iterations, but if you look at the sequence of values of errorEstimate values
does it appear the the ratio r
(k)
2
is approaching a nonzero limit for this case of Newton applied to f1?
(If not, you have a bug: fix your code. The mistake could be in f1.py or in newton1.py.) Based on the
iterations you have, what do you judge is the value of the constant?
15. What is the true error in your approximate solution, |x − 3|? Is it roughly the same size as EPSILON?
16. Try again, starting at x=-0.1, you should get x=-3.
17. Fill in the following table.
Name Formula guess approxRoot No. Steps
f0 x-1 0.1 _________ _________
f1 x^2-9 0.1 _________ _________
f2 x^5-x-1 10 _________ _________
f3 x*exp(-x) 0.1 _________ _________
f4 2*cos(3*x)-exp(x) 0.1 _________ _________
f4 2*cos(3*x)-exp(x) 1.5 _________ _________
5
The theoretical convergence rate of Newton’s method is quadratic. In your work above, you used this fact
to help ensure that your code is correct. This rate depends on an accurate computation of the derivative.
Even a seemingly minor error in the derivative can yield a method that does not converge very quickly or
that diverges. In the above table, none of the cases should have taken as many as 25 iterations. If they
did, check your derivative formula. (As a general rule, if you apply Newton’s method and it takes hundreds
or thousands of iterations but it does converge, check your derivative calculation very carefully. In the
overwhelming majority of cases you will find a bug there.)
8 Exercise 3
The proofs of convergence of Newton’s method show that quadratic convergence depends on the ratio f
00
2f
0
being finite at the solution. In most cases, this means that f
0 6= 0, but it can also mean that both f
00 and f
0
are zero with their ratio remaining finite. When f
00/f0
is not finite at the desired solution, the convergence
rate deteriorates to linear.
1. Write a file f6.py defining the function and derivative for f6=(x-4)^2.
2. You recall using newton1() to find a root of f1=x^2-9 starting from x=0.1. How many iterations did
it take?
3. Now try finding the root (x=4) of f6 using Newton, again starting at x=0.1. How many iterations
does it take? Since the exact answer is x=4, the true error is abs(x-4). Is the true error larger than
EPSILON or smaller?
4. Write a file f7.py defining the function and derivative of f7=(x-4)^20. Now try finding the root (x=4)
using Newton, again starting at x=0.1. For this case, increase maxIts to some large number such
as 1000 (recall maxIts is the optional final argument to newton1()). How many iterations does it
take? Is the true error larger than EPSILON or smaller? In this case, you see that convergence rate can
deteriorate substantially.
5. Look at the final two iterations and compute the values of the ratios
r1 =
|∆x
(k+1)|
|∆x
(k)
|
, and (3)
r2 =
|∆x
(k+1)|
|∆x
(k)
|
2
. (4)
In this case, you should notice that r1 is not nearly zero, but instead is a number not far from 1.0, and
r2 has become large. This is characteristic of a linear convergence rate, and not a quadratic convergence
rate.
In practice, if you don’t know the root, how can you know whether or not the derivative is zero at
the root? You can’t, really, but you can tell if the convergence rate is deteriorating. If the convergence is
quadratic, r2 will remain bounded and r1 will approach zero. If the convergence deteriorates to linear, r2
will become unbounded and r1 will remain larger than zero.
9 Exercise 4
1. We are now going to make a new version of our Newton code. Because newton1.py is working code,
it is a good idea to keep it around so if we make errors in our changes then we can start over. Make
a copy of newton1.py called newton2.py. In the new file, change the name of the function inside to
newton2().
6
The following changes should be made to the newton2.py file
2. Just before the beginning of the loop in newton2.py, insert the following statement
inc remen t=1 # t h i s i s an a r b i t r a r y v al u e
Just before setting the value of increment inside the loop, save increment in oldIncrement:
ol d I n c r em e n t=inc remen t # t h i s i s an a r b i t r a r y v al u e
3. After the statement x+=increment;, compute r1 and r2 according to (3) and (4) respectively.
4. Comment out the existing print statement and add another print statement to display the values of
r1 and r2.
5. Use newton2() to fill in the following table, using the final values at convergence, starting from x=0.1.
Enter your estimate of the limiting values of r1 and r2 and use the term “unbounded” when you think
that the ratio has no bound. As in the previous exercise, use maxIts=1000 for these tests.
Function numIts r1 r2
f1=x^2-9 ______ ______ ______
f6=(x-4)^2 ______ ______ ______
f7=(x-4)^20 ______ ______ ______
This exercise shows that quadratic convergence sometimes fails, usually resulting in a linear convergence
rate, and you can estimate the rate. This is not always a bad situation–it is converging after all–but now
the stopping criterion is not so good. In practice, it is usually as important to know that you are reliably
close to the answer as it is to get the answer in the first place.
In the cases of (x − 4)2 and (x − 4)20, you saw that r1 turned out to be approaching a constant, not
merely bounded above and below. If r1 were exactly a constant, then
x
(∞) = x
(n) +
X∞
k=n
∆x
(k)
because it converges and the sum collapses. I have denoted the exact root by x
(∞)
. Because r1 is constant,
for k > n, ∆x
(k) = r
k−n
1 ∆x
(n)
. Hence,
x
(∞) = x
(n) + ∆x
(n) X∞
k=n
r
k−n
1
= x
(n) +
∆x
(n)
1 − r1
This equation indicates that a good error estimate would be
|x
(∞) − x
(n)
| ≈ |∆x
(n)
|
1 − r1
10 Exercise 5
We are now going to make yet another version of the code, so that we can use the new error estimate idea.
Make a copy of newton2.py called newton3.py.
7
1. Inside the file, change the name of the function from newton2() to newton3().
2. Comment out the print statement displaying r1 and r2, since it will become a distraction when large
numbers of iterations are needed.
3. Conveniently, r1 either goes to zero or remains bounded. If the sequence converges, r1 should remain
below 1, or at least its average should remain below 1. Replace the if-test for stopping to
i f e r r o r E s tim a t e < EPSILON∗(1− r 1 ) :
return x , numIts
Note: This code is mathematically equivalent to errorEstimate=abs(increment)/(1-r1), but I
have multiplied through by (1 − r1) to avoid the problems that occur when r1 ≥ 1. Convergence will
never be indicated when r1 ≥ 1 because errorEstimate is non-negative.
4. Use newton3() to fill in the following table, where you can compute the absolute value of the true
error because you can easily guess the exact solution. Use maxIts=1000 for this exercise. You have
already done these cases using the original convergence criterion in Exercise 3, so you can get those
values from that exercise.
newton1 newton1 newton3 newton3
Function start numIts true err numIts true error
f1=x^2-9 0.1 ______ ______ _______ _______
f6=(x-4)^2 0.1 ______ ______ _______ _______
f7=(x-4)^20 0.1 ______ ______ _______ _______
You should see that the modified convergence does not harm quadratic convergence (about the same
number of iterations required) and greatly improves the estimated error and stopping criterion in the
linear case. A reliable error estimate is almost as important as the correct solution–if you don’t know
how accurate a solution is, what can you do with it?
This modification of the stopping criterion is very nice when r1 settles down to a constant value quickly.
In real problems, a great deal of care must be taken because r1 can cycle among several values, some larger
than 1, or it can take a long time to settle down.
11 Choice of initial guess
The theorems about Newton’s method generally start off with the assumption that the initial guess is “close
enough” to the solution. Since you don’t know the solution when you start, how do you know when it is
“close enough?” In one-dimensional problems, the answer is basically that if you stay away from places
where the derivative is zero, then any initial guess is OK. More to the point, if you know that the solution
lies in some interval and f
0
(x) 6= 0 on that interval, then the Newton iteration will converge to the solution,
starting from any point in the interval. When there are zeros of the derivative nearby, Newton’s method can
display highly erratic behavior and may or may not converge. In the last part of the previous exercise, you
saw a case where there are several roots, with zeros of the derivative between them, and moving the initial
guess to the right moved the chosen root to the left.
12 Exercise 6
In this and the following exercise, you will be interested in the sequence of iterates, not just the final result.
Re-enable the print statement displaying the values of the iterates in newton3().
8
1. Copy the file cosmx.py from the previous lab, which defines cosmx as the function (f(x) = cos x − x).
Modify the file to add a definition cosmxp which is the derivative of the function.
2. Use newton3() to find the root of cosmx starting from the initial value x=0.5. What is the solution
and how many iterations did it take? (If it took more than ten iterations, go back and be sure your
formula for the derivative is correct.)
3. Again, use newton3() to find the root of cosmx, but start from the initial value x=12. Note that
3π < 12 < 4π, so there are several zeros of the derivative between the initial guess and the root. You
should observe that it takes the maximum number of iterations and seems not to converge.
4. Try the same initial value x=12, but also use maxIts=5000. (This is going to cause a large number of
lines of printed information, but you are going to look at some of those lines.) Does it locate a solution
in fewer than 5000 iterations? How many iterations does it take? Does it get the same root as before?
5. Look at the sequence of values of x that Newton’s method chose when starting from x=12. There is no
real pattern to the numbers, and it is pure chance that finally put the iteration near the root. Once it
is “near enough,” of course, it finds the root quickly as you can see from the estimated errors. Is the
final estimated error smaller than the square of the immediately preceeding estimated error?
You have just observed a common behavior: that the iterations seem to jump about without any real
pattern until, seemingly by chance, an iterate lands inside the circle of convergence and they converge rapidly.
This has been described as “wandering around looking for a good initial guess.” It is even more striking in
multidimensional problems where the chance of eventually landing inside the ball of convergence can be very
small.
13 Exercise 7
Another possible behavior is simple divergence to infinity. The following exercise presents a case of divergence
to infinity.
Attempt to find the root of f3 = x*exp(-x), starting from x=2. You should find it diverges in a monotone
manner, so it is clear that the iterates are unbounded. What are values of iterates 95 through 100? This
behavior can be proved, but the proof is not required for this lab.
14 Exercise 8
Sometimes people become so confident of their computational tools that they attempt the impossible. What
would happen if you attempted to find a root of a function that had no roots? Basically, the same kind of
behavior occurs as when there are zeros of the derivative between the initial guess and the root.
1. Write a file for f8=x^2+9, which defines f8 as the function value, and f8p as its derivative.
2. Apply newton3() to f8, starting from x=0.1. Describe what happens.
15 Exercise 9
But the function x
2 + 9 does have roots! The roots are complex, but Python knows how to do complex
arithmetic. (Actually the roots are imaginary, but it is all the same to Python.) All Python needs is to be
reminded to use complex arithmetic.
If you have used the letter j as a variable in your Python code, its special value as the square root of -1
has been obscured.
9
For complex constants, Python will accept the syntax 2+3j to mean the complex number whose real part
is 2 and whose imaginary part is 3.
1. Apply newton3() to f8=x^2+9, starting from the initial guess x=1+1j. (Python interprets “1j” and
“3j” as √
−1 and 3√
−1 respectively.) You should get a result close to 0+3j and it should take fewer
than 10 iterations.
2. Fill in the following table by using newton3() to find a root of f8 from various initial guesses. The
exact root, as you can easily see, is ±3j, so you can compute the true error (typically a small number),
as you did in Exercise 1.
Initial guess numIts true error
1+1j _______ _________
1-1j _______ _________
10+5j _______ _________
10+eps*1j _______ _________
Recall that eps is the smallest number you can add to 1.0 and still change it. Define eps by the
command
import numpy a s np
ep s = np . f i n f o ( np . f l o a t ) . ep s
What happened with that last case? The derivative is 2x and is not zero near x = (10 + (eps)j) ≈ 10. In
fact, the derivative is zero only at the origin. The origin is not near the initial guess nor either root. It turns
out that the complex plane is just that: a plane. While Newton’s method works in two or more dimensions,
it is harder to see when it is going to have problems and when not. We will elaborate a bit in a later section
and also return to this issue in the next lab.
16 Unpredictable convergence
The earlier, one-dimensional cases presented in this lab might lead you to think that there is some theory
relating the initial guess and the final root found using Newton’s method. For example, it is natural to
expect that Newton’s method will converge to the root nearest the initial guess. This is not true in general!
In the exercise below, you will see that it is not possible, in general, to predict which of several roots will
arise starting from a particular initial guess.
Consider the function f9(z) = z
3 + 1, where z = x + yj is a complex number with real part x and
imaginary part y. This function has the following three roots.
ω1 = −1
ω2 = (1 + √
3)/2j
ω3 = (1 −
√
3)/2j
17 Exercise 10
In the following exercise, you will choose a number of regularly-spaced points in the square given by |x| ≤ 2
and |y| ≤ 2 and then start your newton3() iteration at each of these points. You will then construct an
image by coloring each starting point in the square according to which of the four roots was found when
starting from the given point. The surprise comes in seeing that nearby initial guesses do not necessarily
result in the same root at convergence.
10
1. Be sure your newton3() does not have any active print statements in it.
2. Create a file f9.py that defines the function f9 and derivative f9p.
3. Download the file exercise10.py, which uses Newton’s method to analyze the function, and to plot the
results. Describe the information that the plot is exhibiting.
import numpy a s np
import m a t pl o tli b . p y pl o t a s p l t
from f 9 import f9 , f 9p
from newton3 import newton3
NPTS = 100
x = np . l i n s p a c e ( −2, 2 , NPTS)
y = np . l i n s p a c e ( −2, 2 , NPTS)
omega = np . z e r o s ( 3 , dtype=np . complex )
omega [ 0 ] = −1 # red
omega [ 1 ] = ( 1 + np . s q r t ( 3 ) ∗ 1 j ) / 2 # orange
omega [ 2 ] = ( 1 − np . s q r t ( 3 ) ∗ 1 j ) / 2 # p u r pl e
z = np . z e r o s ( (NPTS, NPTS) , dtype=np . complex )
fo r k in range ( 0 , NPTS) :
fo r m in range ( 0 , NPTS) :
z [ k , m] = x [ k ] + y [m] ∗ 1 j
z = z . f l a t t e n ( )
which = [ ]
fo r k in range ( 0 , len ( z ) ) :
x , n um i t e r s = newton3 ( f9 , f9p , z [ k ] , maxIts = 500 )
whichRoot = np . argmin ( np . abs ( x − omega ) )
i f whichRoot == 0 :
which += [ ’ red ’ ]
e l i f whichRoot == 1 :
which += [ ’ o r an ge ’ ]
e l i f whichRoot == 2 :
which += [ ’ p u r pl e ’ ]
e l s e :
print ( ’You have a bug ’ )
p l t . s c a t t e r ( z . r e al , z . imag , c o l o r=which )
p l t . s c a t t e r ( omega . r e al , omega . imag )
p l t . y l a b e l ( ’ Imaginary ’ )
p l t . x l a b e l ( ’ Real ’ )
p l t . show ( )
4. Execute your modified file. The position (z=x+yj) on this image corresponds to the initial guess and
the color values red, orange and purple correspond to roots 1, 2 and 3 to which the Newton iteration
converged starting from that initial guess. This image illustrates that, for example, some initial guesses
with large positive real parts converge to the root ω1 = −1 despite being closer to both of the other
roots. Please include this plot with your summary file.
It turns out that the set you have plotted is a Julia set (fractal). You can increase the detail of the plot
by increasing NPTS and waiting longer. If you wish to learn more about “Newton Fractals,” I suggest a
Wikipedia article https://en.wikipedia.org/wiki/Newton fractal or a very comprehensive paper
by J. H. Hubbard, D. Schleicher, S. Sutherland. “How to Find All Roots of Complex Polynomials by
Newton’s Method,” Inventiones Mathematicæ 146 (2001)1-33. In addition, the chatty but informative
article https://www.chiark.greenend.org.uk/∼esgtatham/newton contains interesting information
and plots.
11
18 Newton’s method revisited
One disadvantage of Newton’s method is that we have to supply not only the function, but also a derivative.
In the following exercise, we will try to make life a little easier by numerically approximating the derivative
of the function instead of finding its formula. One way would be to consider a nearby point, evaluate the
function there, and then approximate the derivative by
y
0 ≈
f(x + ∆x) − f(x)
∆x
(5)
There are many possible ways to pick step sizes ∆x that change as the iteration proceeds. One way is
considered in Exercise 13 below, but we will look at simply choosing a fixed value in the following two
exercises.
19 Exercise 11
Make a copy of newton3.py and call it newtonfd.py. The convergence criterion should involve the factor
(1-r1) in newtonfd.py because using an approximate derivative usually leads to slower linear convergence.
1. Change the function definition to be
1 def newton fd ( func , x0 , dx , maxIts ) :
and change the comment lines to reflect the new function.
2. Previously, your code probably had a line like
1 d e r i v a t i v e = d func ( x )
Now, we don’t have the dfunc() function, and instead, are going to approximate the derivative using
Equation (5) using the constant stepsize dx. Write this new approximation:
1 d e r i v a t i v e = ? ? ?
3. Check out your code using the function f1(x)=x^2-9, starting from the point x = 0.1, using an
increment dx=0.00001. Since newtonfd() doesn’t need the a user-defined derivative function, your
import statement can be simplified to only requesting the function definition:
1 from f 1 import f 1
You should observe essentially the same converged value in the same number (±1) of iterations as using
newton3().
20 Exercise 12
It turns out that the function f1(x)=x^2-9 is an “easy” function to solve using Newton iterations. More
“difficult” functions require a more delicate choice for dx.
Now we will compare the performance of the Newton method with user-supplied or approximated derivative values. We have f2(x)=x^5-x-1 and f6(x)=(x-4)^2. Use the starting point x = 5.0, compute the
number of steps taken to converge using both newton3() once, and newtonfd() using a sequence of different
stepsizes dx. You may have to increase the value of maxIts by a very great deal to get f6 to converge for
the cases of large dx.
12
using newton3 Steps for f2 + f2p steps for f6 + f6p
__________________ __________________
using newtonfd Steps for f2 steps for f6
dx = 0.00001 __________________ __________________
dx = 0.0001 __________________ __________________
dx = 0.001 __________________ __________________
dx = 0.01 __________________ __________________
dx = 0.1 __________________ __________________
dx = 0.5 __________________ __________________
As you can see, the choice of dx is critical. It is also quite difficult in more complicated problems.
As mentioned earlier, Newton’s method generally converges quite quickly. If you write a Newton solver
and observe poor or no convergence, the first thing you should look for is an error in the derivative. One way
of finding such an error is to apply a method like newtonfd() and print both the estimated derivative and
the analytical derivative computed inside your function. They should be close. This trick is especially useful
when you are using Newton’s method in many dimensions, where the Jacobian can have some correctlycomputed components and some components with errors in them.
You saw the secant method in Lab 3, and it had the update formula
x = b −
b − a
f(b) − f(a)
f(b). (6)
21 Exercise 13
1. Compare (6) with (1) and show that by properly identifying the variables x, a, and b then (6) can be
interpreted as a Newton update using an approximation for f
0
(x
(k)
). Explain your reasoning.
2. Starting from your newtonfd.py, copy it to a new file named secantfd.py and modify it to carry out
the secant method. The secant method requires an extra value for the iterations to start, but there is
only the value x0 in the signature line. Take b=x and a=x-dx only for the first iteration. Fill in the
following table.
Newtonfd Newtonfd secantfd secantfd
Function start numIts solution numIts solution
f1=x^2-9 0.1 ______ ______ ______ ______
f3=x*exp(-x) 0.1 ______ ______ ______ ______
f6=(x-4)^2 0.1 ______ ______ ______ ______
f7=(x-4)^20 0.1 ______ ______ ______ ______
You should observe that both methods converge to the same solution but that the secant method takes
more iterations, but not nearly so many as newtonfd() with a larger fixed dx. This is because the
theoretical convergence rate is quadratic for Newton and (1 + √
5)/2 ≈ 1.62 for secant because of the
approximation to the derivative.
13