\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small Sixth Mississippi State
Conference on Differential Equations and Computational
Simulations, {\em Electronic Journal of Differential Equations},
Conference 15 (2007), pp. 163--173.\newline ISSN: 1072-6691. URL:
http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{163}
\title[\hfilneg EJDE-2006/Conf/15\hfil Approximations of continuous Newton's method]
{Approximations of continuous Newton's method:
An extension of Cayley's problem}
\author[J. Jacobsen, O. Lewis, B. Tennis\hfil EJDE/Conf/15 \hfilneg]
{Jon Jacobsen, Owen Lewis, Bradley Tennis} % in alphabetical order
\address{Jon Jacobsen \newline
Mathematics Department, Harvey Mudd College, 301 Platt Blvd,
Claremont, CA 91711, USA}
\email{jacobsen@math.hmc.edu}
\address{Owen Lewis \newline
4047 North Castle Ave., Portland, OR 97227, USA}
\email{owen.lewis@gmail.com}
\address{Bradley Tennis \newline
Computer Science Department, Stanford University, Stanford, CA
94305, USA}
\email{btennis@stanford.edu}
\dedicatory{Dedicated to Klaus Schmitt on his 65th birthday}
\thanks{Published February 28, 2007.}
\subjclass[2000]{34C35, 58C15, 28A80, 65H10}
\keywords{Newton's method; damping; fractal basins of attraction}
\begin{abstract}
Continuous Newton's Method refers to a certain dynamical system whose
associated flow generically tends to the roots of a given polynomial.
An Euler approximation of this system, with step size $h=1$,
yields the discrete Newton's method algorithm
for finding roots. In this note we contrast Euler approximations with
several different approximations of the continuous ODE
system and, using computer experiments, consider their impact
on the associated fractal basin boundaries of the roots.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\section{Introduction}
In many contexts and applications one is faced with the task of
finding a value $x$ such that $f(x)=0$ for a given function $f$.
Despite being a relatively easy problem to state, for certain
choices of $f$, this task may be formidable. In the late
17$^{th}$ century Newton considered an iterative method to find
zeros of real-valued functions $f:\mathbb{R} \to \mathbb{R}$. Starting with an
initial seed value $x_0$, Newton proposed the iterative method
\begin{equation}
x_{n+1} = x_{n} - \frac{f(x_n)}{f'(x_n)} \quad x_n \in \mathbb{R},
\label{newton1}
\end{equation}
to generate a sequence $\{x_n\}$ that converges to a root of $f$.
In 1870 Schr\"oder \cite{schroder:1870} and
in 1879, Cayley \cite{cayley:1879}
imitated the method for complex polynomials $p(z)$ and considered
the behavior of sequences
\begin{equation}
z_{n+1} = z_{n} - \frac{p(z_n)}{p'(z_n)} \quad z_n \in \mathbb{C}.
\label{newton-c}
\end{equation}
For quadratic polynomials Cayley found that the method did indeed
converge to a root. Moreover, he proved (via conjugation to the
map $z^2 + c$) that the complex plane could be divided into two
half-planes, such that seeds in either half converged to the
unique root lying in its own half. In other words, Cayley
characterized the global basins of attraction $A(\xi) =\{ z_0 \in
\mathbb{C} : z_n \to \xi\}$ for each root $\xi$.
Having succeeded in the
quadratic case, Cayley remarked that \lq\lq the next succeeding
case of the cubic equation appears to present considerable
difficulty.\rq\rq \, The task of identifying such basins for
complex polynomials is now known as \textit{Cayley's Problem}.
Nearly 40 years later, work of French mathematicians Julia and
Fatou \cite{julia:1918, fatou:1920} shed much light on Cayley's
difficulties in their groundbreaking study of iterating complex
rational mappings. In particular, for polynomials, the right-hand
side of (\ref{newton-c}) represents a rational map of the extended
complex plane to itself and their work implies that the basins of
the three roots share a common boundary, the so-called
\textit{Julia set} of $p(z)$. Without having seen an example
before, it is hard to imagine a ternary division of the complex
plane such that every point on the boundary of one set is also on
the boundary of the other two sets. In other words, imagine
partitioning the plane into three \lq\lq states\rq\rq\, such that
every point on the border of one state also simultaneously lies on
the border of the other two states. Basins of attraction with
this property are said to satisfy the Wada property, a concept
introduced by Japanese mathematician
Koniz\^o Yoneyama \cite{yoneyama:1917} in 1917.
In the early 1980's, and just over a century after Cayley's work,
computer simulations allowed one to visualize Cayley's basins and
their associated fractal structure. For the cubic $p(z)=z^3 - 1$
the basins are captured by the image in Figure 1.
Although a familiar image now, it is hard to imagine to what
degree Cayley, Julia, Fatou and other early researchers imagined
these basins to look like.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.5\textwidth]{figures/cubic}
\label{classic-newton}
\caption{Basins of Attraction for $p(z)=z^3-1$.}
\end{center}
\end{figure}
In hindsight, one can understand why Schr\"oder and Cayley's
approach worked. We now recognize that the Newton map
(\ref{newton-c}) is precisely an Euler approximation of the
differential equation
\begin{equation}
z'(t) = - \frac{p(z(t))}{p'(z(t))}, \quad t > 0,
\label{euler-newt}
\end{equation}
with step size $h=1$. Equation (\ref{euler-newt}) is known as
\textit{Continuous Newton's Method}. If $z(t)$ solves
(\ref{euler-newt}) with $z'(t) \to 0$ as $t \to \infty$, then it
should be that $z(t)$ converges to a root of $p$. Indeed, suppose
$z(t)$ solves (\ref{euler-newt}) with $z(0)=z_0$. Assuming $p(z(t)) \neq 0$ for $t > 0$,
then it follows $\frac{p'(z)}{p(z)} \dot{z}(t) = -1$, or equivalently,
$p(z(t))=p(z_0) e^{-t}$. Thus, solutions $z(t)$ flow to a zero of $p$ while
keeping the argument of $p(z(t))$ constant at $\arg(p(z_0))$.
In this way we can view
Newton's method (for real or complex polynomials) as an
approximation to this flow, which leads to roots in most cases.
For example, if $p(z)=z^3 -1,$ then we can explicitly solve for
$z(t) = \sqrt[3]{(z_0^3-1)e^{-t}+1},$ where the appropriate branch
of the cube root is chosen corresponding to the nearest root of
$p(z)$. In this case, the corresponding basins of attraction are
defined by the ternary division of the plane defined by the rays
$\theta=\frac{\pi}{3}, \theta=-\frac{\pi}{3},$ and $\theta=\pi$,
as suggested by the rightmost image in Figure 2. This is the
cubic analog of the binary division
Cayley found for quadratic polynomials.
From this point of view, the fractal basin boundaries that arise
in Newton's method (such as those in Figures 1 and 2) are due to
the discretization of the Continuous Newton flow
(\ref{euler-newt}). In one sense, the complex boundaries arise
from the numerical error inherent in the approximation of the
flow.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.32\textwidth]{figures/cnm-step1p3-28}
\includegraphics[width=0.32\textwidth]{figures/cnm-step0p94-frac19}
\includegraphics[width=0.32\textwidth]{figures/cnm-step0p4-frac6}
\caption{Basins for $p(z)=z^3-1$ with damped Newton method:
$h= 1.3, 0.9, 0.4$ respectively.}
\label{damp-newton}
\end{center}
\end{figure}
The first author to recognize the relevance of the continuous
Newton method to classical Newton's method appears to be Smale
\cite{smale:1976}, in his work on finding equilibrium points.
Several other authors, including Hirsch and Smale
\cite{hirsch:1979}, Braess \cite{braess:1977}, Saupe
\cite{saupe:1988}, Peitgen, Pr\"ufer, and Schmitt
\cite{peitgen:1988}, and Neuberger \cite{neuberger:1999}, have
also consider aspects of continuous methods. In particular, the
continuous Newton's method (\ref{euler-newt}) leads one to
consider the damped Newton algorithm defined by approximating
(\ref{euler-newt}) using Euler's method with step size $h$:
\begin{equation}
x_{n+1} = x_{n} -h \frac{f(x_n)}{f'(x_n)}, \quad x_n \in \mathbb{R},
\label{newton2}
\end{equation}
for $0 < h< 2$ (when $h=2$ the fixed points of
(\ref{newton2}) become repelling). Neuberger
\cite{neuberger:1999} observed that the associated Cayley basins
get larger and their fractal boundary structure shrinks away as
$h \to 0$, as indicated in Figure \ref{damp-newton}. This
is in agreement with the fact that smaller step sizes lead to
improved accuracy in the approximation of equation
(\ref{euler-newt}).
With this in mind, in this paper we consider how the fractal
boundaries behave with respect to improved numerical integration
techniques. For instance, what do the basins look like if one
uses a Runge-Kutta algorithm on (\ref{euler-newt}) instead of the
Euler method? How do the sizes of the fractal basin boundaries
compare for different step sizes and different methods? Methods
such as Runge-Kutta provide higher-order accuracy, but their
associated iterative map (for a given step size) is also more
complex (e.g., than (\ref{newton2})). With this in mind, we
consider several numerical schemes for approximating
(\ref{euler-newt}) and compare their associated Cayley basins.
We consider this in the context of complex polynomials as well as
for general vector fields on $\mathbb{R}^n$. For related work based on
altering explicitly the discrete iterative step (\ref{newton2})
see Varona \cite{varona:2002} and Epureanu et
al.~\cite{greenside:1998}.
\section{Refined Algorithms}
In this section we consider six different numerical algorithms for
approximating Continuous Newton's Method (\ref{euler-newt}). For
each method we consider the associated fractal basin boundaries
for the cubic $p(z)=z^3-1,$ for a range of step sizes. We first
define the methods in the context of solving $\frac{dy}{dt} =
f(t, y)$, using $w_i$ for the iterations and $h$ for the step
size. The six methods we consider are:
\begin{itemize}
\item Euler's Method. This basic numerical technique with step size $h$
leads to the iterative scheme known as Damped Newton's Method:
\[
w_{i+1} = w_i + h \, f(t_i, w_i).
\]
\item Refined Euler.
The refined Euler's method performs a size $\frac{h}{2}$ Euler step and uses this result to
perform midpoint integration:
\begin{align*}
w^* &= w_i + \frac{h}{2}f(t_i, w_i)\\
w_{i+1} &= w_i + hf\big(t_i + \frac{h}{2}, w^*\big).
\end{align*}
\item Heun's Method. Another refined Euler method, this method
performs a size $h$ Euler step and uses this result to perform trapezoidal integration:
\begin{align*}
w^* &= w_i + hf(t_i, w_i)\\
w_{i+1} &= w_i + \frac{h}{2}\left[f(t_i, w_i) + f(t_i + h,
w^*)\right].
\end{align*}
\item Runge-Kutta (Order 2).
Though the Refined Euler method and Heun's method are also second order Runge-Kutta methods, this method is an optimal one in the sense that its coefficients are chosen to minimize the error in the approximation of the Taylor series of $y$:
\begin{align*}
w^* &= w_i + \frac{2h}{3}f(t_i, w_i)\\
w_{i+1} &= w_i + \frac{h}{4}\big[f(t_i, w_i) + 3f\big(t_i +
\frac{2h}{3}, w^*\big)\big].
\end{align*}
\item Runge-Kutta (Order 4).
The classical fourth order Runge-Kutta method is defined by:
\begin{align*}
k_1 &= hf(t_i, w_i)\\
k_2 &= hf\big(t_i +\frac{h}{2}, w_i + \frac{k_1}{2}\big)\\
k_3 &= hf\big(t_i +\frac{h}{2}, w_i + \frac{k_2}{2}\big)\\
k_4 &= hf(t_i + h, w_i + k_3)\\
w_{i+1} &= w_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4).
\end{align*}
\item Adams-Bashforth (Order 2).
The two-step Adams-Bashforth method is a second order multi-step method. The initial iterate is assigned to $w_0$ and $w_1$ is calculated using a second order Runge-Kutta approximation:
\begin{align*}
w^* &= w_0 + \frac{2h}{3}f(t_i, w_0)\\
w_{1} &= w_i + \frac{h}{4}\big[f(t_i, w_0) + 3f\big(t_i + \frac{2h}{3},
w^*\big)\big]\\
w_{i+1} &= w_i + h\big[\frac{3}{2}f(t_i, w_i) -
\frac{1}{2}f(t_{i-1}, w_{i-1})\big].
\end{align*}
\end{itemize}
For each method we apply the associated iterative map for each
point in the region $[-1.5,1.5] \times [-1.5,1.5] \subset
\mathbb{R}^2$, using a mesh of $1024 \times 1024$. Each point is
iterated and colored according to which root the orbit converges
to. We use a tolerance of $\epsilon= 1 \times 10^{-10}$ for
determining which root, and set the maximum iterations at 500.
Figures 3 and 4 contain select computed basins for three fixed
step sizes: $h=1.4$, $0.9$, and $0.4$. The corresponding images
for Newton's method appear in Figure 2. In addition to considering the basins for fixed step size
$h$, one can also consider the dynamic flow of the basins as the step size is continuously varied. The accompanying website \cite{magifrac} contains several such movies.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.32\textwidth]{figures/RefinedNewton-fractal-3}
\includegraphics[width=0.32\textwidth]{figures/RefinedNewton-fractal-2}
\includegraphics[width=0.32\textwidth]{figures/RefinedNewton-fractal-1} \\
\vspace{0.025in}
\includegraphics[width=0.32\textwidth]{figures/Heun-fractal-3}
\includegraphics[width=0.32\textwidth]{figures/Heun-fractal-2}
\includegraphics[width=0.32\textwidth]{figures/Heun-fractal-1} \\
\label{r-newton}
\caption{Basins for $p(z)=z^3-1$ with (top): Refined Newton method; (bottom):
Heun's method. The images reflect step sizes (from left to right) of
$h=1.4$, $0.9$, and $0.4$.}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=0.32\textwidth]{figures/RK_4-fractal-3}
\includegraphics[width=0.32\textwidth]{figures/RK_4-fractal-2}
\includegraphics[width=0.32\textwidth]{figures/RK_4-fractal-1} \\
\vspace{0.025in}
\includegraphics[width=0.32\textwidth]{figures/RK_2-fractal-3}
\includegraphics[width=0.32\textwidth]{figures/RK_2-fractal-2}
\includegraphics[width=0.32\textwidth]{figures/RK_2-fractal-1} \\
\vspace{0.025in}
\includegraphics[width=0.32\textwidth]{figures/AB_2-fractal-3}
\includegraphics[width=0.32\textwidth]{figures/AB_2-fractal-2}
\includegraphics[width=0.32\textwidth]{figures/AB_2-fractal-1}
\caption{Basins for $p(z)=z^3-1$ with (top): Runge-Kutta, Order 4 method; (middle): Runge-Kutta, Order 2 method; (bottom): Adams-Bashforth, Order 2 method. The images reflect step sizes (from left to right) of
$h=1.4$, $0.9$, and $0.4$.}
\label{rk4-newton}
\end{center}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.3\textwidth]{figures/e-re}
\includegraphics[width=0.3\textwidth]{figures/e-rk4}
\includegraphics[width=0.3\textwidth]{figures/e-h} \\
\includegraphics[width=0.3\textwidth]{figures/e-rk2}
\includegraphics[width=0.3\textwidth]{figures/re-solid-ab-cross}
\includegraphics[width=0.3\textwidth]{figures/rk4-line-re-cross-h-diamond} \\
\caption{Fractal dimension vs. step size for (a) Euler (line) vs. Refined Euler (cross);
(b) Euler (line) vs. Runge-Kutta Order 4 (cross); (c) Euler (line) vs. Heun (cross);
(d) Euler (line) vs. Runge-Kutta Order 2 (cross); (e) Refined Euler (line) vs.
Adams-Bashforth (cross); (f) Runge-Kutta Order 4 (line), Refined Euler (cross), Heun (diamond).}
\label{dimension}
\end{center}
\end{figure}
Due to the Wada property, any
neighborhood of a point on a basin boundary will contain points whose associated iteration will converge to each and
every root of
the polynomial. In particular, trajectories corresponding to arbitrarily close initial data on the boundary will have
divergent orbits. In this sense the boundary
describes the unstable set for a given approximation. In order to compare the complexity
of these unstable sets we
estimate their fractal dimension using a box-counting algorithm.
With this approach
we find several interesting dynamics with respect to
the variation in fractal dimension. In particular, better algorithms did not necessarily
produce basin boundaries with smaller fractal dimension. Several sample dimension
comparisons are shown in Figure 5.
Although in many cases the higher-order approximations lead to basins with smaller fractal dimension,
the associated computation time can also be significantly larger, given the higher complexity of the algorithm.
\section{Newton's Method in $\mathbb{R}^n$}
There is a natural extension of continuous Newton's method to
higher dimensional systems. Namely, for $G:\mathbb{R}^n \to \mathbb{R}^n$ the associated continuous
flow for $x:\mathbb{R} \to \mathbb{R}^n$ is
\begin{equation}
\dot{x} = - [DG(x(t))]^{-1} G(x(t)),
\label{rn-cnm}
\end{equation}
where $DG$ is the Jacobian matrix of the vector field $G$.
The flow is defined for
$x \notin \mathcal{S}$, where
\begin{equation}
\mathcal{S} = \{ x \in \mathbb{R}^n : \det DG(x) = 0 \},
\end{equation}
is the singular set of $G$. The Euler approximation
of (\ref{rn-cnm}) with step size $h=1$ yields Newton's Method for $\mathbb{R}^n$:
\begin{equation}
x_{n+1} = x_n - [DG(x_n)]^{-1} G(x_n). \label{newtrn}
\end{equation}
In this case, we have left the realm of rational maps of the complex plane, and the theory
of iterated rational maps no longer applies.
For example,
the basin boundaries need not be common and there is no known general characterization as
precise as the Julia and Fatou set theory for complex maps.
Peitgen, Pr\"ufer, and Schmitt \cite{peitgen:1988} suggest the following candidates for
the Julia-like sets for discretizations of (\ref{rn-cnm}):
\begin{equation}
\mathcal{J}_{h} = {\left\{ x \in \mathbb{R}^n : N^k(z) \in \mathcal{S}
\text{ for some } k \in N \right\}}, \label{ss}
\end{equation}
where $N$ is an associated iterative map based on a discretization
of (\ref{rn-cnm}) with step size $h$, and $N^k$ denotes the
$k^{th}$ iteration of $N$. To illustrate the dynamics of the
generalized continuous Newton's Method we follow Peitgen et
al.~\cite{peitgen:1988} and consider vector fields $G$ that arise
from discretizing the nonlinear boundary value problem
\begin{equation}
\begin{gathered}
u'' + \lambda (u-u^2) = 0, \quad 0 < x < 1, \\
u(0)=u(1)=0.
\end{gathered}
\label{pss}
\end{equation}
A centered difference approximation
with a two point interior mesh for the system (\ref{pss}) yields the vector field $G:\mathbb{R}^2 \to \mathbb{R}^2$
defined by
\begin{equation}
G(x,y) = \begin{bmatrix} 2x - y - \mu (x - x^2) \\
2 y - x - \mu ( y - y^2) \end{bmatrix},
\label{vf}
\end{equation}
where $\mu = -\frac{\lambda}{9}$.
That is, roots of $G$ yield approximate values of solutions
to (\ref{pss}) at the two interior mesh points. To illustrate the
dynamics, we consider the system for $\mu =3.2$, for which
(\ref{pss}) has four solutions given by the trivial solution, the
positive solution, and two sign-changing solutions (one the
negative of the other). Figure 6 shows a comparison of the basin
boundaries of the four roots with the singular set, showing a good
match. Notice that in this setting the basin boundaries do not
satisfy the Wada property. In particular, the basin for the
positive solution does not meet the other basins, other than on
the border in the first quadrant.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.4\textwidth]{figures/mu3p2-E-Basins-step1p6b}
\includegraphics[width=0.4\textwidth]{figures/mu3p2-E-SS-step1p6}
\caption{Basins of attraction (left) and singular set (right) for (\ref{vf}) with Euler approximations
for $\mu=3.2$ and $h=1.6$. This singular set example is due to Peitgen et al.~\cite{peitgen:1988}. }
\label{alien-basins}
\end{center}
\end{figure}
In Figure 7 we show the corresponding images for three additional methods.
In each of these examples the associated singular set defined by (\ref{ss})
corresponded with the basin boundaries. Again, the boundary for
the positive root does not meet the other basins, other than in a
small neighborhood of the first quadrant.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.3\textwidth]{figures/mu3p2-pde2-step1p6-RK4-basins}
\includegraphics[width=0.3\textwidth]{figures/mu3p2-AB2-step1p6-basins}
\includegraphics[width=0.3\textwidth]{figures/mu3p2-step1p6-basins-Heun}
\caption{Basins of attraction for (\ref{vf}) with $\mu = 3.2,$ and
$h=1.6$ using (left) Runge-Kutta Order 4; (middle)
Adams-Bashforth; and (right) Heun's Method.}
\label{alien-basins-altmethods}
\end{center}
\end{figure}
Generically, the basins of attraction alone need not reveal the
singular set. For example, Figure 8 shows a comparison of the
basins with the singular set for $\mu = 2.1$. In this case
(\ref{pss}) has only the trivial solution and the positive
solution, and the basins do not reflect the singular set. This is
true generically for the case $1 < \mu < 3$, for which equation
(\ref{pss}) has only the the trivial and positive solution
(compare to Figure 6, with three nonpositive solutions). Further
properties of the orbits can reveal the singular set. For
example, Peitgen et al.~demonstrate that the singular set
structure can be revealed by further partitioning the basins into
level sets based on rates of convergence \cite{peitgen:1988}.
The images in this article were generated using C++ code. We
have developed an application that allows the interested reader to
generate their own images, or sequence of images, for any two
dimensional vector field of interest \cite{magifrac}.
This reference also contains several movies
that accompany the material in this paper, and give further
insight into the global dynamics for these different methods. In
particular, the movies show very clearly how the basins vary as
the step size $h$ varies. For each system, when $h$ is
sufficiently small, the images appear the same, up to a small
rescaling, such as demonstrated in Figure 2. However, this is not
generically true. For example, the approximation methods in
Figure 4 show notable differences in the basins for the same
three step size parameters used in Figure 2.
Finally, it is interesting to note the behavior of the basins as
$h$ approaches the critical value for which the roots become
repelling fixed points. To illustrate the complexity that can
arise we end with a few images near the edge of stability in
Figure 9.
\begin{figure}
\begin{center}
\includegraphics[width=0.4\textwidth]{figures/mu2p1-E-basins-step-h1p2-pde2}
\includegraphics[width=0.4\textwidth]{figures/mu2p1-SS-step1p2-pde2}
\caption{Basins of attraction (left) and singular set (right) for (\ref{vf}) with Euler approximations
(similar images result for the other methods) with $\mu = 2.1$, and $h=1.2$. With only two roots, the basin boundaries do not capture the singular set.}
\label{alien-basins2}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=0.4\textwidth]{figures/e}
\includegraphics[width=0.4\textwidth]{figures/rk2} \\
\vspace{0.02in}
\includegraphics[width=0.4\textwidth]{figures/AB2-step1p999-maxit80k-dim1p72}
\includegraphics[width=0.4\textwidth]{figures/h}
\caption{Basins on the edge of stability for:
(top left) Newton's Method, $h=1.995$; (top right) Runge-Kutta Order 2, $h=1.354$;
(bottom left) Adams-Bashforth, $h=1.999$; (bottom right) Heun's method, $h=1.995$. Reference \cite{magifrac} contains several animations of the basins as the step size is varied from the stable to the unstable regime.}
\label{final-basins-altmethods}
\end{center}
\end{figure}
\subsection*{Acknowledgments} We would like to thank John M.
Neuberger (Northern Arizona University) for many helpful
discussions on the ideas in this paper.
\begin{thebibliography}{10}
\bibitem{braess:1977}
D.~Braess, \emph{\"{U}ber die {E}inzugsbereiche der {N}ullstellen von
{P}olynomen beim {N}ewton-{V}erfahren}, Numer. Math. \textbf{29} (1977/78),
no.~1, 123--132.% \MR{MR0474765 (57 \#14398)}
\bibitem{cayley:1879}
A.~Cayley, \emph{Desiderata and suggestions. {N}o. 3.-the {N}ewton-{F}ourier
imaginary problem}, Amer. J. Math \textbf{2} (1879), no.~1, 97.
\bibitem{greenside:1998}
B.~I. Epureanu and H.~S. Greenside, \emph{Fractal basins of attraction
associated with a damped {N}ewton's method}, SIAM Rev. \textbf{40} (1998),
no.~1, 102--109 (electronic). %\MR{MR1612502 (99d:65153)}
\bibitem{fatou:1920}
P.~Fatou, \emph{Sure les equations fonctionnelles}, Bull. Soc. Math., France
\textbf{47--48} (1919--1920), 161--314.
\bibitem{hirsch:1979}
M.~W. Hirsch and S.~Smale, \emph{On algorithms for solving {$f(x)=0$}}, Comm.
Pure Appl. Math. \textbf{32} (1979), no.~3, 281--313.
\bibitem{julia:1918}
G.~Julia, \emph{Memoire sure l'iteration des fonction rationelles}, J. Math.
Pures et Appl. \textbf{81} (1918), 47--235.
\bibitem{magifrac}
Magifractificator, \emph{Available electronically},
{\tt{http://www.math.hmc.edu/magifrac/}}.
\bibitem{neuberger:1999}
J.~W. Neuberger, \emph{Continuous {N}ewton's method for polynomials}, Math.
Intelligencer \textbf{21} (1999), no.~3, 18--23.
\bibitem{peitgen:1988}
H.-O. Peitgen, M.~Pr{\"u}fer, and K.~Schmitt, \emph{Global aspects of the
continuous and discrete {N}ewton method: a case study}, Acta Appl. Math.
\textbf{13} (1988), no.~1-2, 123--202.
\bibitem{saupe:1988}
D.~Saupe, \emph{Discrete versus continuous {N}ewton's method: a case study},
Acta Appl. Math. \textbf{13} (1988), no.~1-2, 59--80.
\bibitem{schroder:1870}
E.~Schr{\"{o}}der, \emph{Ueber unendlich viele algorithmen zur aufl{\"{o}}sung
der gleichungen}, Math. Ann. \textbf{2} (1870), 317--365.
\bibitem{smale:1976}
S.~Smale, \emph{A convergent process of price adjustment and global {N}ewton
methods}, J. Math. Econom. \textbf{3} (1976), no.~2, 107--120.
\bibitem{varona:2002}
J.~L. Varona, \emph{Graphic and numerical comparison between iterative
methods}, Math. Intelligencer \textbf{24} (2002), no.~1, 37--46.
\bibitem{yoneyama:1917}
K.~Yoneyama, \emph{Theory of continuous sets of points}, Tohoku Math. J.
\textbf{12} (1917), 43--158.
\end{thebibliography}
\end{document}