Chaos on the line between two roots

Where Newton’s method goes instead ran Newton’s method for $z^3 - 1$ from every point of the plane and found three basins whose common boundary is a fractal, every point of it touched by all three. In passing it recorded that Arthur Cayley had solved the quadratic case completely in 1879 — two roots, two basins, the two half-planes on either side of a straight line — and had then written that the cubic “appears to present considerable difficulty”. The fractal was the difficulty, and nobody saw it for another century.

This essay goes back to the case Cayley solved, because the solution has more in it than a straight line. The line itself is where Newton’s method never settles, and there the method is not merely indecisive; it is chaotic, in the most exact sense the word has. The tool that shows it is a single change of variable, and the same change, applied to Newton’s method for $x^2 + 1$ on the real line — an equation with no real root to find — turns a root-finder into a generator of random numbers with a famous and strange distribution.

Two roots, two half-planes #

Newton’s method replaces a guess $z$ by $z - p(z)/p'(z)$. For $p(z) = z^2 - 1$ that is

$$N(z) = z - \frac{z^2 - 1}{2z} = \frac{z^2 + 1}{2z},$$

which is the average of $z$ and $1/z$. Cayley’s theorem is that a starting point in the right half-plane — real part positive — converges to $+1$, a point in the left half-plane to $-1$, and a point on the imaginary axis, the perpendicular bisector of the two roots, stays on that axis for ever. The same holds for any quadratic with two distinct roots: the basins are the two half-planes cut off by the perpendicular bisector, because a rotation, scaling and shift of the plane carries any two points to $\pm 1$ and Newton’s method respects such changes.

The opening figure shows the theorem working and shows something it does not say. Orbits starting well inside a half-plane arrive in two or three steps. The orbit starting at $0.04 + 0.9i$, close to the line, wanders for a while before committing to the right side. Near the line the method is slow, and on the line it never commits at all. That suggests the line is not just a boundary but a place with dynamics of its own, and a change of coordinates makes them visible.

A change of variable that makes it squaring #

Set

$$w = \frac{z - 1}{z + 1}.$$

This is a Möbius transformation of the kind the sphere that complex numbers live on studied: it sends $+1$ to $0$, $-1$ to infinity, and the imaginary axis — the points equally far from $+1$ and $-1$ — to the unit circle $|w| = 1$, since $|z - 1| = |z + 1|$ exactly there. The right half-plane, where $z$ is closer to $+1$, goes to the inside of the circle, and the left half-plane to the outside.

The remarkable fact is what Newton’s method becomes in the new variable. A line of algebra shows

$$\frac{N(z) - 1}{N(z) + 1} = \frac{z^2 - 2z + 1}{z^2 + 2z + 1} = \left(\frac{z - 1}{z + 1}\right)^2,$$

so in the $w$ coordinate, Newton’s method is $w \mapsto w^2$. The figure checks the identity at scattered points and then draws one orbit in both coordinates side by side.

Everything about the quadratic case now follows from squaring. A point inside the unit circle has $|w| < 1$, its squares fall to $0$, and $0$ is the root $+1$; the number of correct digits doubles at every step, the quadratic convergence for which Newton’s method is famous. A point outside flies to infinity, which is the root $-1$. And a point on the circle stays on the circle, its angle doubling at every step. The same map in different coordinates called this kind of relation a conjugacy: two maps that are one map seen through a change of variable, so that every dynamical fact about one transfers to the other. Newton’s method for a quadratic is squaring.

On the line, the doubling map #

The circle is where the interest is. Write a point of it as $w = e^{i\phi}$; squaring sends it to $e^{2i\phi}$, so on the circle the map is $\phi \mapsto 2\phi$, angle doubling. The orbit written as a word made doubling the standard example of chaos: write the angle as a binary fraction of a full turn, and doubling shifts the digits one place left and discards the whole part. Two angles that agree to ten binary places agree for ten steps and then are unrelated. Almost every angle’s orbit visits every part of the circle with even frequency. And the periodic points — angles that return after $n$ steps — are the fractions $k/(2^n - 1)$ of a turn, dense on the circle and every one of them repelling.

All of that transfers back to the dividing line. A starting point exactly on the imaginary axis never reaches a root; its orbit wanders the axis forever, sensitively dependent on its starting value, and returns periodically only from a dense countable set of special starts. Cayley’s straight line is a doubling map in disguise. For the cubic the boundary between basins is a fractal and the dynamics on it are complicated; for the quadratic the boundary is as simple as a set can be and the dynamics on it are as chaotic as dynamics can be.

The real line, where there is no root #

There is a second way to see the same map, and it is the one a calculator meets. Rotate the picture by a quarter turn: $z^2 - 1$ becomes $x^2 + 1$ with $x = iz$, and the imaginary axis of $z$ becomes the real line of $x$. Newton’s method for $x^2 + 1 = 0$ on the real numbers is

$$x \mapsto x - \frac{x^2 + 1}{2x} = \frac{1}{2}\left(x - \frac{1}{x}\right).$$

The equation has no real root. A real starting value can never reach $\pm i$, because the map sends reals to reals; the method looks for a root that is not on the line it is confined to.

The conjugacy is now a trigonometric identity. Put $x = \cot\theta$. Then

$$\frac{1}{2}\left(\cot\theta - \tan\theta\right) = \frac{\cos^2\theta - \sin^2\theta}{2\sin\theta\cos\theta} = \frac{\cos 2\theta}{\sin 2\theta} = \cot 2\theta,$$

so Newton’s method for $x^2 + 1$ is angle doubling, read through the cotangent. The cobweb in the figure bounces between the two branches of the map’s graph without ever approaching anything, and the circle beside it shows why: the angles $2\theta, 4\theta, 8\theta, \dots$ spread round the circle with no pattern a short orbit can reveal.

Where the orbits spend their time #

Doubling spreads almost every angle evenly round the circle. Carried back through $x = \cot\theta$, an even spread of angles gives a definite distribution of $x$: the chance that $\cot\theta$ lies in a short interval near $x$ is proportional to the angle that interval subtends, and that gives the density

$$\frac{1}{\pi(1 + x^2)}.$$

That curve is the Cauchy distribution, and an average that never settles made it famous for a single property: it has no mean. Its tails fall off like $1/x^2$, so slowly that the average of $n$ samples is distributed exactly like a single sample, however large $n$ is. So the long-run behaviour of Newton’s method for an equation with no real root is this: its iterates are, in distribution, Cauchy random variables, and the running average of the iterates never settles on anything either. A procedure designed to converge produces the most famous example of non-convergence in probability.

The figure’s caption hides a computational subtlety worth stating. Each doubling step destroys one binary digit of the starting angle, and a computer holds about fifty-three. After fifty or so steps a floating-point orbit is no longer an orbit of the true map from any nearby start; it has become an orbit of the computer’s rounding, and long such orbits drift measurably from the Cauchy histogram — the orbit a computer draws described the general phenomenon. The histogram is therefore built from many short orbits, each well within its precision, rather than from one long one.

There is a check of the density that needs no orbits at all. If $X$ has the Cauchy distribution, then so does $(X - 1/X)/2$: write $X = \cot\Theta$ with $\Theta$ spread evenly over half a turn, and the map sends it to $\cot 2\Theta$, where $2\Theta$ is again spread evenly once taken modulo half a turn. So the Cauchy distribution is carried to itself by one step of the method, exactly — the invariant density in the strictest sense. The same calculation gives the rate at which nearby orbits separate. Averaging the logarithm of the map’s slope, $\tfrac12(1 + 1/x^2)$, against the Cauchy density gives exactly $\ln 2$: one binary digit of the starting point lost per step, which is the doubling map’s rate, and the reason the floating-point orbits above last about fifty steps.

The orbits that come back #

Periodic points of the doubling map are angles that are fractions $k/(2^n - 1)$ of a half turn, here, since $\cot$ has period $\pi$. Carried through the cotangent they give every periodic point of Newton’s method for $x^2 + 1$.

For $n = 2$ the period-two points are $\cot(\pi/3) = 1/\sqrt3$ and $\cot(2\pi/3) = -1/\sqrt3$, and indeed $\tfrac12(1/\sqrt3 - \sqrt3) = -1/\sqrt3$: Newton’s method bounces between them forever. For $n = 3$ there are six more, the cotangents of $k\pi/7$, which form two orbits of length three. The counts $1, 3, 7, 15, \dots$ are $2^n - 1$, the fixed points of the $n$-fold doubling, and every one of these orbits is repelling: a start a hair away from one drifts off it at a rate that doubles the separation each step.

So the line carries infinitely many cycles, all unstable, dense in the line, with the typical orbit wandering among them. That is the definition of chaos used in dynamics — sensitive dependence, dense periodic points, a dense orbit — satisfied exactly and provably by the simplest root-finding iteration there is, on the simplest equation that has no real solution. How fast two orbits part measured the separation rate for maps where it must be estimated; here it is exactly a factor of two per step, the Lyapunov exponent $\ln 2$.

One step off the line #

The rule $x \mapsto \tfrac12(x - 1/x)$ is the Babylonian square-root method in disguise. For a positive number $c$, Heron’s rule $x \mapsto \tfrac12(x + c/x)$ is Newton’s method for $x^2 - c$, and it converges to $\sqrt c$ from every positive start, doubling its correct digits at each step; it is how square roots were computed by hand for two thousand years. Put $c = -1$ and the same rule is asked for $\sqrt{-1}$, and on the real line it has nowhere to go.

Allow the starting value a little imaginary part, however small, and everything changes. A start $x_0 + i\varepsilon$ with $\varepsilon > 0$ lies in the upper half-plane, which in the rotated picture is one of Cayley’s two basins; the orbit converges to $i$, and from below to $-i$. But the convergence is slow at first in exact proportion to how small $\varepsilon$ is: in the $w$ coordinate the start sits just inside the unit circle, a distance of order $\varepsilon$ from it, and squaring has to bring that radius down to nought, which takes about $\log_2(1/\varepsilon)$ steps before the familiar doubling of digits takes over. Meanwhile the angle doubles, and the orbit’s real part wanders exactly as the chaotic real orbit would.

So the real line is a knife-edge in a precise sense. Every start off it reaches a root; the time it takes grows like the number of binary digits of closeness to the line; and on the line itself the time is infinite and the wandering never stops. That is the general picture of a basin boundary in miniature — a set of measure nought on which the dynamics are chaotic, surrounded by points whose convergence time grows without bound as they approach it — and the essay on the double root met a different face of the same slowness, where the trouble came from two roots merging rather than from a point equidistant between them.

Why the cubic is different, and why it is not #

The quadratic case looks like the trivial one, and in a sense it is: the basins are half-planes and the boundary is a line. But the conjugacy explains why the cubic is hard in a way the half-planes alone do not. For the quadratic, the dividing set is where the dynamics of squaring live, the unit circle, a smooth curve. For the cubic $z^3 - 1$ there is no Möbius transformation that turns Newton’s map into $w^3$, because the Newton map of a cubic is a rational map of degree three with three roots and a free critical point at $0$, while $w^3$ has two critical points, both fixed. The two maps have different numbers of points that matter, and no change of coordinates can make them the same.

What survives in every degree is the shape of the result: the boundary of the basins is the set where the method’s dynamics are chaotic, the so-called Julia set, and the basins are its complement. For the quadratic that set is a circle seen in the right coordinates and a line in the original ones; for the cubic it is the fractal of the first essay on Newton’s basins; and for the trapped cubics of an area that never finishes the complement includes, besides the root basins, an open region that never finishes at all. The line between two roots is the first member of the family, and the only one in which every feature can be written down.

What the drawings show and what they assume #

The orbits in the first two figures are computed in floating point and converge to the roots to within $10^{-9}$, which the figures check; that a point exactly on the dividing line stays on it is a statement about exact arithmetic, since a computed orbit near the line drifts off it within a few dozen steps as rounding errors are doubled. The figures draw orbits near the line and say what exact arithmetic would do on it; they cannot exhibit an orbit that stays on the line for a hundred steps.

The histogram is likewise a picture of the theorem rather than a proof of it. That doubling spreads almost every angle evenly is an ergodic theorem; the histogram shows its consequence in a sample of two hundred thousand points, built from short orbits for the reason given above. And the periodic points are checked by iterating each cotangent $n$ times and confirming it returns, to a relative error of $10^{-6}$ — which the large periods strain, since an orbit of period eight passes near zero, where the map divides by a small number.

Still open: when the conjugacy is only approximate #

For the quadratic, Newton’s map is exactly conjugate to squaring on the whole sphere. For polynomials of higher degree, near each root Newton’s method is conjugate to $w \mapsto w^2$ locally — that is what quadratic convergence means — but no global conjugacy exists, and the basins’ boundaries are fractal. Between the exact case and the general one sit families where the conjugacy almost holds: Newton maps of polynomials with two roots of high multiplicity, or with roots spread along a line, whose basin boundaries are close to straight lines over most of the plane and fractal only near particular points. How close, and where the fractal parts are, is described in detail for special families and not in general.

The real-line version raises its own question. Newton’s method for $x^2 + 1$ is conjugate to doubling, which is uniformly chaotic. For an equation with some real roots, the real Newton iteration typically has both converging orbits and chaotic regions, and for most real polynomials the structure of the chaotic set — its dimension, how much of the line it occupies, which periodic orbits it carries — has to be found case by case. For cubics with one real root and two complex ones it has been studied extensively; a general description for real polynomials of every degree does not exist.

A straight line with everything on it #

Cayley’s solution is usually quoted as the dull half of a story whose interesting half is the cubic. Seen through the right change of variable it is not dull at all. Newton’s method for a quadratic is squaring, its two basins are the inside and outside of a circle, and its boundary is a circle on which the method doubles angles — the textbook example of chaos, sitting on the straight line between two roots that anyone could have drawn in 1879.

The same map on the real line is what a calculator does when asked to solve $x^2 + 1 = 0$ by Newton’s method: it never gives up and never converges, its iterates have the distribution of the ratio of two independent bell curves, and their average never settles. An algorithm built to find answers, pointed at an equation with no answer it can reach, does the next most definite thing: it behaves exactly like chance.