Bisection and Newton's method

Halley's Comet is now 35 astronomical units from the Sun, farther than Neptune, and is slowly turning round: it will be back near the Sun in 2061. Its orbit is known with great precision, but to say where on that orbit the comet is on a given day you have to solve Kepler's equation: M = E − e · sin E. There is no formula for its solution, and for more than four hundred years it has been solved by trial, refining E until the equality holds. Below are five real comets and four ways of solving it. Move the time, change the comet and the method.

In the equation e is the eccentricity of the orbit: 0.017 for the Earth's almost circular orbit, 0.968 for Halley's Comet, 0.9992 for comet NEOWISE, whose ellipse is almost a parabola. M is the mean anomaly: it grows at a constant rate from perihelion to perihelion, so it is simply time in other units. E is the eccentric anomaly, the angle from which the position on the ellipse follows at once: the distance to the Sun is a(1 − e cos E). The second chart shows f(E) = E − e · sin E − M: the root is where the curve crosses the axis. Bisection always finds it in 43 steps, Newton in a handful: four for Halley's Comet today. But switch to comet Hale–Bopp or NEOWISE and look at the strip under the slider: the red marks are days on which Newton from the standard start E₀ = M never converges.

The idea

Bisection: if a continuous function has different signs at the ends of an interval, there is a root inside. Halve the interval and keep the half whose ends still have different signs. Each step halves the place where the root hides, guaranteed, whatever the function. Newton's method: near a point the graph almost coincides with its tangent. Where the tangent crosses the axis is easy to find, so step there: x ← x − f(x)/f′(x). Near the root the error is squared at every step, and the number of correct digits doubles. But only near the root.

The two methods have opposite characters, and that is the main thing to know about them. Bisection needs an interval with different signs at its ends and the ability to compute the sign of f; it needs no derivative, no smoothness, and there is no such thing as a bad start. But its speed is always the same, one bit per step, for a simple function and a complicated one alike. Newton needs the derivative and a good starting point, and in return it gets a speed bisection does not have. The school side of both, the intermediate value theorem and the proof of quadratic convergence, is in “Mathematics, the Queen of the Sciences”: the chapter on limits and the chapter on extrema.

A thousand digits

Take √2, the root of x² − 2. Bisection on the interval [1, 2] and Newton from x₀ = 1, step by step:

step   bisection           digits   Newton              digits
   1   1.250000000000000       0   1.500000000000000       1
   2   1.375000000000000       1   1.416666666666667       2
   3   1.437500000000000       1   1.414215686274510       5
   4   1.406250000000000       2   1.414213562374690      11
   5   1.421875000000000       2   1.414213562373095      16
bisection needs 50 steps for all 16 digits, Newton 5

Five Newton steps give all 16 digits that fit into a floating-point number. Bisection needs 50 steps for the same. Where does the doubling come from? Let r be the root and x = r + ε. Expand f around the root: f(x) ≈ f′(r)ε + ½f″(r)ε². Substituting into the Newton step gives the new error

For √2 the factor is 1/(2√2) ≈ 0.35, and an error of 10⁻⁵ becomes 3.5 · 10⁻¹¹: five correct digits become ten or eleven. A rigorous proof with a bound on the remainder is in “Mathematics, the Queen of the Sciences”.

In floating point the doubling has no room to grow beyond sixteen digits. In big-number arithmetic it does. Below Newton computes √2, ∛2 and 1/7 to a thousand decimal places, and next to it bisection makes as many steps.

For √2 Newton's correct digits go 1, 2, 5, 11, 24, 48, 97, 195, 391, 783, each step about twice as many, and by the eleventh step there are more than a thousand. In the same eleven steps bisection extracts three or four digits: a step gives it one binary digit, and a decimal digit costs log₂ 10 ≈ 3.32 steps. Look at 1/7: the Newton step for the equation 1/x − 7 = 0 is x ← x(2 − 7x), with no division at all. This is how numbers are divided where division is expensive: in big-number arithmetic and even in ordinary programs; the GCC compiler with the -mrecip flag replaces single-precision division by a fast approximate reciprocal plus one Newton step.

Why bisection never misses

The number of steps does not depend on the function, only on the width of the interval and the precision needed. The interval [0, 2π] for Kepler's equation shrinks to 10⁻¹² in ⌈log₂(2π · 10¹²)⌉ = 43 steps, the same for comet Churyumov–Gerasimenko and for NEOWISE, as you saw in the scene. It is the same binary search, only over real numbers instead of array indices, and the same technique as binary search on the answer.

Code

Both methods in general form: bisection needs a function and an interval, Newton a function, its derivative and a start. The program prints the table for √2 you saw above and two examples where Newton loses its speed or runs away; more on them below.

Stopping is a separate question. Bisection is stopped by the width of the interval, which is known exactly. Newton by the size of the step: if a correction is smaller than the tolerance, the following ones will be smaller still. Comparing x with the previous value for equality will not do: in floating point Newton can rock between two neighbouring numbers. That happens here: on the sixth step √2 differs from the correct double by one unit in the last place.

And this is Kepler's equation on real data: the orbital elements of five comets from the JPL database and four ways of solving it, for 16 October 2026.

16 Oct 2026            e        M°       E°   r, au   bisect  N(E0=M)  N(E0=pi)  Kepler
1P              0.967936  192.925  186.575   35.17      43        4         4    667
C/2020 F3       0.999178    0.333   18.473   18.75      43       14        10    466

The last column is Kepler's own method from the “Epitome of Copernican Astronomy” of 1621: substitute E into the right-hand side of E = M + e · sin E again and again. It always converges, but slowly: the error is multiplied each time by about e · |cos E|. Halley's Comet, now near aphelion where cos E ≈ −1, needs 667 steps: a factor of 0.96 eats the error very unhurriedly. Newton needs four.

Where Newton breaks

The promise “the digits double” comes with the proviso “near a simple root”. Every word in it matters.

Far from the root. In Kepler's equation the derivative 1 − e · cos E is close to zero near perihelion if e is close to one: the tangent is almost horizontal and throws the iterate very far. The classic example from the 1997 paper by Charles and Tatum is e = 0.992, M = 0.13π:

e = 0.992, M = 0.13pi, Newton from M: 0.4084 4.806 -1.13 -0.0178 50.07 -1761 202.2
from pi: 7 steps; bisection: 43 steps; E = 1.382957944863

The root is 1.383, while the iterates jump around in hundreds and thousands of radians. The program does not print the path beyond the sixth step: it depends on the last bit of the sine, and Python and Go, computing the same thing, diverge. The start E₀ = π solves the problem in seven steps. The program below checks this on all five comets: it takes 4000 evenly spaced values of M and runs Newton from both starts, for at most 50 steps.

For the first three comets both starts always converge. For Hale–Bopp (e = 0.995) the start E₀ = M failed to converge in our Python run in 58 cases out of 4000, for NEOWISE in 91; the start π never failed and never took more than 12 steps. The Go version counted 70 and 94. The code is the same, but sine and cosine in Go and in the C system library that Python uses are different programs, and they occasionally disagree in the last bit, while in this region the trajectory is so sensitive that the last bit decides whether a start converges. So this program has no reference output, and on another machine the numbers may differ slightly. Hence the practical rule of astronomical software: for large e start from π, and more reliably, back Newton up with bisection.

A multiple root. If the graph touches the axis instead of crossing it, f′(r) = 0, and the quadratic shrinking of the error is gone. For (x − 1)² from x₀ = 2 the error only halves: 0.5, 0.25, 0.125, …, the same speed as bisection, only without its guarantees.

No smoothness. For ∛x the root is simple, but the derivative there is infinite. The tangent from a point x crosses the axis at −2x, and the method runs away: −2, 4, −8, 16, …

Cycles. For x³ − 2x + 2 Newton from zero lands on one, from one back on zero, and so on forever. How many such traps there are is best seen on the complex plane.

Newton fractals

Newton's method works for complex numbers too: the same formula z ← z − p(z)/p′(z). The polynomial z³ − 1 has three roots: 1 and (−1 ± i√3)/2. Colour each point of the plane by the root the method reaches from it, brighter where it gets there fast.

Around each root there is a calm region of its colour: from there the method converges fast and where you expect. But the boundaries between the regions are infinitely complicated. Zoom into any of them: wherever two colours meet, the third is there too, at every scale. This is a theorem, not a feature of the picture: all three regions share one and the same boundary, the Julia set of the map z ↦ z − p(z)/p′(z), and in any neighbourhood of any of its points there are starts leading to each of the three roots. Switch to z³ − 2z + 2: the black islands around 0 and 1 are the starts that get pulled into the cycle.

The same picture as text: the program puts into each cell the sign of the root the method reaches from it:

In 1879 Arthur Cayley asked how the complex plane is divided between the roots. For a quadratic equation the answer is simple and neat, the boundary is a straight line, while the cubic case, Cayley wrote, appeared to present considerable difficulty. How considerable became visible only a hundred years later, on computer screens.

What it costs

Bisection needs ⌈log₂((b − a)/ε)⌉ evaluations of the sign of f, and this number is known before it starts. Newton from a good start needs about log₂ d steps for d correct digits: 5 steps for 16 digits, 11 for a thousand. But a Newton step is dearer: it evaluates both f and f′. If the derivative is expensive, the secant method is used, a line through the last two points instead of the tangent; its digits grow not twofold but 1.618 times per step, with one evaluation.

Big-number arithmetic has one more trick. On the first steps there are few correct digits, and computing them with a thousand digits is pointless. So the precision is raised along with the steps: 2 digits, 4, 8, …, 1000. The last step costs more than all the previous ones together, and the whole root costs a few multiplications of thousand-digit numbers.

In practice the two characters are combined. Reliable solvers keep an interval with different signs, like bisection, but step with fast methods, and if a fast step jumps out of the interval or narrows it too slowly, they take a bisection step. This is how Brent's method of 1973 works; in SciPy it is the function brentq, which the documentation calls the best of its bracketing solvers.

Where it is used

• Roots and division in big-number arithmetic. In Go big.Int.Sqrt repeats the Newton step z ← ⌊(z + ⌊x/z⌋)/2⌋ from a start known to be too large, and big.Float.Sqrt solves the equation 1/t² − x = 0 to avoid division. In the pure-Python reference implementation of the decimal module (_pydecimal) the integer part of a square root is also found by Newton's method.
• The fast inverse square root. In the 1999 game Quake III Arena 1/√x was obtained from the special constant 0x5F3759DF and one Newton step; a second step was in the source but commented out, as the precision was enough.
• Celestial mechanics. The positions of planets, comets and satellites from orbital elements are Kepler's equation for each moment of time, solved by Newton or its relatives.
• Finance. The internal rate of return, the rate at which the payments pay back the investment, is found by Excel's IRR function by iteration: to an accuracy of 0.00001 % and in at most 20 attempts, otherwise the error #NUM!. The first problem below is about the same thing.
• Fitting a parameter. Any “at what speed, rate or temperature will the condition hold” with a monotone dependence is bisection; that is why solvers like brentq require an interval with different signs.

History

Johannes Kepler derived the equation itself in 1609 in “Astronomia nova”, and in 1621 proposed solving it by substitution, the very iteration in the last column of the table. Who invented “Newton's method” is still argued about. Newton, in a manuscript of 1669 printed only in 1711, refined a root of a polynomial by corrections; the technique first reached print in John Wallis's “Treatise of Algebra” of 1685. Joseph Raphson wrote it down in 1690 as a repeatable formula, hence the second name, the Newton–Raphson method. In its present form, with the derivative and for any function, it was published by Thomas Simpson in 1740, and the historian Nick Kollerstrom argued in 1992 that the method would more properly be called Simpson's. The way of extracting a square root itself, x ← (x + 2/x)/2, is much older: it is in Heron's “Metrica”, and the chapter on the square root of two of “Mathematics, the Queen of the Sciences” tells about it in detail.

Problems

In the first problem the root is sought for money: what rate a buyer on instalments really pays. In the second, Kepler's equation with tests where the start E₀ = M does not converge.