Appearance
10.2 — Finding Roots Numerically
Solve x^5 - x - 1 = 0.
You cannot. Chapter 2.3 noted that Abel and Galois proved no formula in radicals exists for a general fifth-degree equation. And even for equations that do have formulas, most real problems are not polynomials at all — x = \cos x, or the equation that gives the internal rate of return on an investment, or the temperature at which two rates balance.
The answer is to compute the root numerically, to whatever accuracy you need. This chapter covers the three methods worth knowing, and the honest comparison between them.
1. Bisection: slow and never fails
Chapter 5.1's Intermediate Value Theorem: a continuous function that changes sign must cross zero in between.
The algorithm. Find a and b with f(a) and f(b) of opposite signs. Take the midpoint. Whichever half still has a sign change, keep it. Repeat.
Worked example. f(x) = x^3-x-2, with f(1) = -2 and f(2) = 4.
| Step | a | b | mid | f(\text{mid}) | keep |
|---|---|---|---|---|---|
| 1 | 1 | 2 | 1.5 | -0.125 | right |
| 2 | 1.5 | 2 | 1.75 | 1.61 | left |
| 3 | 1.5 | 1.75 | 1.625 | 0.666 | left |
| 4 | 1.5 | 1.625 | 1.5625 | 0.252 | left |
Converging on 1.5214. Each step halves the interval, so the error falls by a factor of 2 per step — linear convergence, gaining about one decimal digit every 3.3 steps.
Its virtue is that it cannot fail. Given a valid starting bracket and a continuous function, it converges, always, with a guaranteed error bound at every step. No other method here can say that.
Its vice is that it is slow, and it needs a bracket, which means you must already know roughly where the root is and it cannot find a root the function touches without crossing.
2. Newton's method: fast and occasionally reckless
Use the tangent line. From a guess x_n, follow the tangent to where it hits the axis, and take that as the next guess.
The tangent at x_n is y = f(x_n) + f'(x_n)(x-x_n). Setting y = 0:
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
Same example, f(x)=x^3-x-2, f'(x)=3x^2-1, starting at x_0=2:
x_1 = 2 - \frac{4}{11} = 1.6364
x_2 = 1.6364 - \frac{0.7963}{7.033} = 1.5232
x_3 = 1.5214
x_4 = 1.52138
Four steps to six digits, where bisection needed about twenty.
The convergence is quadratic: the number of correct digits roughly doubles each step. Two digits, then four, then eight, then sixteen. That is why Newton's method dominates numerical practice.
Why it doubles. The error in a tangent-line approximation is second order in the step, by Chapter 5.8's Taylor series — the linear term is exactly what the tangent captures, so what is left starts at the quadratic term.
How Newton's method fails
It is fast when it works and it does not always work.
Zero derivative. A flat spot sends the tangent off to infinity. Division by zero, or a wild jump.
Bad starting point. With f(x) = x^3-2x+2 and x_0 = 0, the iteration cycles between 0 and 1 forever.
Slow near a repeated root. If the root is a double root, the derivative vanishes there too, and convergence degrades from quadratic to linear.
Wandering off. Nothing keeps the iterates near the root, so a poor start can send it somewhere else entirely, or to a different root than the one you wanted.
Which root you land on can be fractal. Applying Newton's method to z^3 = 1 in the complex plane and colouring each starting point by which of the three roots it reaches produces a Newton fractal — the boundaries between the three regions are infinitely intricate, and arbitrarily close to any point of one colour are points of both others.
What real software does: hybrid methods. Brent's method combines bisection's guaranteed bracket with faster interpolation steps, using the fast step when it stays inside the bracket and falling back on bisection when it does not. It gets Newton-like speed with bisection's reliability, and it is what scipy.optimize.brentq and most library root-finders actually run.
3. The secant method
Newton's method needs the derivative, which you may not have — the function might come from a simulation or a measurement.
Replace the derivative with the slope through the last two points:
x_{n+1} = x_n - f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}
No derivative required. Convergence is slower than Newton but still better than linear, with an order of about 1.618 — the golden ratio again, arising from the recurrence relating successive errors, which is Fibonacci's from Chapter 8.5.
Since each step needs only one new function evaluation while Newton needs both a function and a derivative evaluation, the secant method is often faster in wall-clock time when derivatives are expensive.
4. Newton's method as a practical tool
Computing square roots. To find \sqrt a, solve x^2-a = 0:
x_{n+1} = x_n - \frac{x_n^2-a}{2x_n} = \frac12\left(x_n+\frac{a}{x_n}\right)
Average your guess with a divided by your guess. For \sqrt2 from x_0=1: 1.5, then 1.41667, then 1.414216, then 1.4142136 — seven digits in four steps.
This is the Babylonian method, known around 1700 BCE, and it is Newton's method three and a half thousand years early. It is also what many processors' square root instructions refine internally.
Computing a reciprocal without dividing. Solve \frac1x - a = 0:
x_{n+1} = x_n(2-ax_n)
Only multiplication and subtraction. Division is far more expensive than multiplication in hardware, so processors and GPUs compute 1/a this way from a table lookup plus one or two Newton steps.
The famous fast inverse square root in the Quake III source code — with its notorious constant 0x5f3759df and comment "what the fuck?" — is a bit-level trick to get an initial guess for 1/\sqrt x followed by one Newton iteration. It was several times faster than the library call on 1999 hardware and mattered enormously for lighting calculations.
Internal rate of return. The IRR of a cash flow is the discount rate making the net present value zero, and there is no closed form. Every spreadsheet's IRR function is running a root-finder, which is why it occasionally returns an error — the equation may have multiple roots or none.
5. Systems of nonlinear equations
For several equations in several unknowns, Newton's method generalises. Replace the derivative with the Jacobian matrix of all partial derivatives (Chapter 5.7), and the division with solving a linear system (Chapter 4.3):
J(\mathbf{x}_n)\,\Delta\mathbf{x} = -\mathbf{F}(\mathbf{x}_n), \qquad \mathbf{x}_{n+1} = \mathbf{x}_n+\Delta\mathbf{x}
Each step is a linear solve, which is why Chapter 4.3's efficiency matters so much: a nonlinear problem is a sequence of linear ones.
This is how every circuit simulator finds the operating point of a nonlinear circuit, how every power grid computes its load flow, and how every structural analysis with material nonlinearity converges.
6. Knowing when to stop
Three stopping criteria, and using the wrong one causes real problems.
The step is small: |x_{n+1}-x_n| \lt \varepsilon. Usually the right choice.
The function value is small: |f(x_n)| \lt \varepsilon. Misleading when the function is very flat near the root — you can be far from the root with a tiny function value.
A maximum iteration count. Always include this. Without it, a method that fails to converge loops forever, and in production that is a hung process rather than an error message.
Use a relative tolerance for large values. An absolute tolerance of 10^{-9} is meaningless when the root is near 10^{12}, because Chapter 10.1 says the spacing between representable doubles is already larger than that.
Every formula above, built from scratch
None of the results in this chapter are worth memorising, because each one can be rebuilt in under a minute from something simpler. What follows is that rebuilding, one result at a time, so the formula and the reason for it sit on the same page as the explanation that needed them.
Finding roots
Bisection
\text{If } f(a)f(b)<0, \text{ a root lies between } a \text{ and } b
\text{error after } n \text{ steps} \le \frac{b-a}{2^n}
Why it always works. If f is continuous and changes sign across an interval, it must cross zero somewhere inside — that is the intermediate value theorem. Halve the interval, keep the half where the sign still changes, repeat. The trap never opens.
How many steps for a given accuracy? Solve \frac{b-a}{2^n}<\epsilon:
n > \log_2\frac{b-a}{\epsilon}
To get from an interval of width 1 to 10^{-10} takes \log_2(10^{10}) = 33 steps. Slow but utterly reliable — it gains exactly one bit per step and can never diverge.
Newton's method
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
|e_{n+1}| \approx \frac{|f''|}{2|f'|}e_n^2 \qquad \text{(quadratic convergence)}
Where the error law comes from. Expand f about the root r using Taylor's series from 5.8 — series and taylor, with e_n = x_n - r:
0 = f(r) = f(x_n) - e_nf'(x_n)+\tfrac12e_n^2f''(x_n)+\cdots
Divide by f'(x_n) and rearrange:
e_n - \frac{f(x_n)}{f'(x_n)} = \frac{f''}{2f'}e_n^2 + \cdots
The left side is exactly x_{n+1}-r = e_{n+1}. So the new error is proportional to the square of the old one.
What quadratic convergence means in practice. The number of correct digits doubles each step: 10^{-2}\to10^{-4}\to10^{-8}\to10^{-16}. Four steps from a decent starting guess is usually enough, which is why this is inside the square-root routine of every calculator and standard library.
When it fails. If f'(x_n) is near zero the step is enormous and the iteration flies off. Near a repeated root, f' vanishes at the root itself and convergence drops from quadratic to linear. And a poor starting point can send it into a cycle. Production code combines bisection's safety with Newton's speed, taking a Newton step when it stays inside the bracket and bisecting when it does not.
The secant method
x_{n+1} = x_n - f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}
Newton's method with the derivative replaced by the slope through the last two points. Convergence order is \varphi = 1.618, the golden ratio of 8.1 — sets — slower than Newton per step, but each step costs one function evaluation instead of two, so in practice it is often faster overall.
Fixed-point iteration
x_{n+1} = g(x_n) \quad \text{converges to a fixed point if } |g'(x)|<1 \text{ near it}
Why the condition. Near the fixed point r, an error e_n becomes g(r+e_n)-g(r)\approx g'(r)e_n. So each step multiplies the error by g'(r). Less than 1 in size and the error shrinks; more and it grows. This one condition explains why x_{n+1}=\cos x_n converges (its derivative is small) while x_{n+1}=e^{x_n} does not.
7. Where this shows up in your life
Every spreadsheet's IRR, RATE, XIRR and Goal Seek.
Every square root and reciprocal your processor computes.
Every circuit and power-system simulation.
Every implicit step in a physics engine, where the next state satisfies an equation that must be solved rather than evaluated.
Every calibration curve inverted to get a reading from a sensor.
Every optimisation that sets a gradient to zero — which is the next chapter, and is root-finding on the derivative.
Root-finding solves f(x) = 0. The related and larger problem is finding where a function is smallest, which is what every machine learning system and every engineering design does. But first, the two operations of calculus need their numerical versions.