29. Newton's method
§1.8 gave you bisection: guaranteed to work, one bit per step. Newton's method trades the guarantee for speed, and the speed is dramatic — it roughly doubles the number of correct digits every iteration.
It is also the algorithm your computer actually runs when you ask for a square root, and the direct ancestor of gradient descent (§14.0).
The idea
You want to solve f(x) = 0 and you have a guess x_n.
Replace f by its tangent line at x_n — the linearization from the last lesson — and solve that instead. Lines are easy to zero.
The tangent is y = f(x_n) + f'(x_n)(x - x_n). Setting y=0:
x = x_n - \frac{f(x_n)}{f'(x_n)}
Call that x_{n+1} and repeat:
\boxed{x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}}
Slide down the tangent line to the axis, and start again from there. That's the entire method.
Why it's so fast
Bisection halves the error each step: e_{n+1} \approx \frac{e_n}{2} — linear convergence, one bit per iteration.
Newton squares it:
e_{n+1} \approx \frac{|f''(r)|}{2|f'(r)|}\,e_n^2
quadratic convergence. Error 10^{-3} becomes 10^{-6}, then 10^{-12}, then below machine precision. In practice: correct digits double per step, and you're done in four or five iterations.
The reason is exactly the error law from §3.1. The tangent line approximates f with error O(e_n^2), so the root of the tangent misses the true root by O(e_n^2). Newton's method is fast because linear approximation is second-order accurate.
Computing square roots
To find \sqrt a, solve f(x) = x^2 - a = 0. With f'(x) = 2x:
x_{n+1} = x_n - \frac{x_n^2 - a}{2x_n} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)
Average your guess with a divided by your guess. If x_n is too small then a/x_n is too big, and the average lands between — closer than either.
This is the Babylonian method, in use around 1700 BC, roughly 3,400 years before Newton. For a = 2 starting from x_0 = 1:
| n | x_n | correct digits |
|---|---|---|
| 0 | 1 | 0 |
| 1 | 1.5 | 1 |
| 2 | 1.41666... | 3 |
| 3 | 1.4142156... | 6 |
| 4 | 1.41421356237469... | 12 |
| 5 | 1.41421356237309... | 16 (machine limit) |
Digits: 0, 1, 3, 6, 12, 16. Doubling until floating point runs out.
When it fails
The guarantee is gone, and these are the ways you lose it.
f'(x_n) = 0. Division by zero — the tangent is horizontal and never meets the axis. Even a near-zero derivative throws the next iterate far away.
Cycling. f(x) = x^3 - 2x + 2 starting at x_0 = 0 gives x_1 = 1, x_2 = 0, x_3 = 1, forever. A perfect two-cycle that never converges.
Divergence. f(x) = \sqrt[3]{x} from any nonzero start doubles its distance from the root each step and runs to infinity.
Convergence to the wrong root. With several roots, which one you land on depends on the starting point in a way that is genuinely chaotic — the boundaries between basins of attraction are fractals. Newton fractals, from applying this to z^3 = 1 in the complex plane, are exactly that picture.
Halved rate at multiple roots. If the root is a double root, f' vanishes there too, and convergence degrades from quadratic to linear.
Practical implementations therefore hedge: bracket the root first, use Newton,
and fall back to bisection whenever a step would leave the bracket. That hybrid
(Brent's method) is what scipy.optimize.brentq runs, and it has bisection's
guarantee with something close to Newton's speed.
The connection to optimization
To minimize g rather than solve g = 0, apply Newton to g' = 0:
x_{n+1} = x_n - \frac{g'(x_n)}{g''(x_n)}
In several variables that becomes
\mathbf{x}_{n+1} = \mathbf{x}_n - H^{-1}\nabla g
with H the Hessian — Newton's method for optimization, §14.2. Replacing H^{-1} with a constant \alpha gives plain gradient descent: cheaper per step, far more steps. That trade-off is the central design question in numerical optimization, and it starts here.
Doing it in Python
The Babylonian square root, digits doubling:
from math import sqrt
a, x = 2.0, 1.0
true = sqrt(a)
print(f"{'n':>3} {'x_n':>20} {'error':>12} {'correct digits':>16}")
for n in range(1, 7):
x = 0.5 * (x + a / x)
err = abs(x - true)
digits = 16 if err == 0 else -__import__('math').log10(err)
print(f"{n:>3} {x:>20.15f} {err:>12.2e} {min(digits, 16):>16.1f}")
print("\ndigits roughly double each step until floating point runs out")
Newton versus bisection on the same problem:
from math import log10
def f(x):
return x**3 - x - 2
def fp(x):
return 3*x**2 - 1
# Newton
x = 2.0
newton = []
for _ in range(6):
x = x - f(x) / fp(x)
newton.append(x)
# bisection
lo, hi = 1.0, 2.0
bisect = []
for _ in range(6):
mid = (lo + hi) / 2
lo, hi = (lo, mid) if f(lo) * f(mid) <= 0 else (mid, hi)
bisect.append((lo + hi) / 2)
root = newton[-1]
print(f"{'step':>5} {'newton':>20} {'bisection':>20}")
for i, (n, b) in enumerate(zip(newton, bisect), 1):
print(f"{i:>5} {n:>20.15f} {b:>20.15f}")
print(f"\nafter 6 steps: newton is exact to machine precision,")
print(f"bisection has the root to about {log10(2**6):.1f} digits")
The failure modes, each on purpose:
def newton(f, fp, x0, steps=8):
xs, x = [x0], x0
for _ in range(steps):
d = fp(x)
if d == 0:
xs.append(float('nan'))
break
x = x - f(x) / d
xs.append(x)
return xs
print("cycling: f(x) = x^3 - 2x + 2 from x0 = 0")
print(" ", [round(v, 6) for v in newton(lambda x: x**3 - 2*x + 2,
lambda x: 3*x**2 - 2, 0.0, 6)])
print("\ndiverging: f(x) = cbrt(x) from x0 = 0.1")
cbrt = lambda x: x**(1/3) if x >= 0 else -((-x)**(1/3))
dcbrt = lambda x: (1/3) * abs(x)**(-2/3)
print(" ", [f"{v:.3e}" for v in newton(cbrt, dcbrt, 0.1, 6)])
print("\nthe first cycles forever, the second doubles its distance every step")
Which root you land on, and how sensitive that is:
def newton(f, fp, x0, steps=60):
x = x0
for _ in range(steps):
d = fp(x)
if d == 0:
return float('nan')
x = x - f(x) / d
return x
# x^3 - x has roots at -1, 0, 1
f, fp = lambda x: x**3 - x, lambda x: 3*x**2 - 1
print(f"{'start':>10} {'converges to':>14}")
for x0 in (0.5, 0.55, 0.57, 0.5773, 0.578, 0.6, 1.5):
r = newton(f, fp, x0)
print(f"{x0:>10} {round(r, 6) if r == r else 'diverged':>14}")
print("\nnear x0 = 1/sqrt(3) = 0.5774 the derivative vanishes and the")
print("outcome flips unpredictably. those boundaries are fractal.")
Worked example
Solve \cos x = x to 10 decimal places.
§1.8 proved a solution exists in [0,1] by the IVT, but bisection would need about 34 iterations to reach 10 digits. Newton needs four.
Set f(x) = \cos x - x, so f'(x) = -\sin x - 1:
x_{n+1} = x_n - \frac{\cos x_n - x_n}{-\sin x_n - 1} = x_n + \frac{\cos x_n - x_n}{\sin x_n + 1}
Starting from x_0 = 1 (anywhere sensible works — f' \approx -1.8, comfortably nonzero):
| n | x_n |
|---|---|
| 0 | 1 |
| 1 | 0.7503638678 |
| 2 | 0.7391128909 |
| 3 | 0.7390851334 |
| 4 | 0.7390851332 |
Converged: x = 0.7390851332.
This is the Dottie number — the unique real fixed point of cosine, and what
you converge to by pressing cos on a calculator repeatedly from any starting
value. (That repeated pressing is fixed-point iteration, which converges
linearly at about 0.67 per step; Newton gets there in a fifth of the iterations.)
Why f' mattered here. f'(x) = -\sin x - 1 is at most 0 and equals zero only at x = \frac\pi2 + 2k\pi — nowhere near our root. So no division-by-zero risk, and convergence from essentially any start. Checking that the derivative stays away from zero on your region of interest is the standard due diligence before trusting Newton.
Your turn
1. Use Newton's method to find \sqrt{10} starting from x_0 = 3. Do three iterations.
2. Set up Newton's iteration for \sqrt[3]{a}.
3. Why does Newton's method fail for f(x) = x^{1/3} at any x_0 \neq 0?
4. How many bisection steps match 4 Newton steps starting from an error of 10^{-1}?
Solutions
1. f(x) = x^2 - 10, so the iteration is x_{n+1} = \frac12\left(x_n + \frac{10}{x_n}\right):
x_1 = \tfrac12\left(3 + \tfrac{10}{3}\right) = \tfrac12(6.3333) = 3.166667
x_2 = \tfrac12\left(3.166667 + 3.157895\right) = 3.162281
x_3 = \tfrac12\left(3.162281 + 3.162276\right) = 3.1622777
True: \sqrt{10} = 3.16227766\ldots — 8 correct digits in 3 steps from a start with 1.
2. f(x) = x^3 - a, f'(x) = 3x^2:
x_{n+1} = x_n - \frac{x_n^3 - a}{3x_n^2} = \frac{2x_n}{3} + \frac{a}{3x_n^2} = \frac13\left(2x_n + \frac{a}{x_n^2}\right)
Same shape as the square-root case: a weighted average of the guess and a correction. The general n-th root iteration is \frac1n\left((n-1)x_n + \frac{a}{x_n^{n-1}}\right).
3. With f(x) = x^{1/3} and f'(x) = \frac13x^{-2/3}:
x_{n+1} = x_n - \frac{x_n^{1/3}}{\frac13x_n^{-2/3}} = x_n - 3x_n^{1/3}x_n^{2/3} = x_n - 3x_n = -2x_n
Each step doubles the distance from the root and flips the sign. From $x_0 = 0.1$ you get -0.2, 0.4, -0.8, 1.6, \ldots — divergence, guaranteed, from every nonzero start.
The cause is the vertical tangent at 0 (§2.1): f' blows up near the root, so the tangent line is nearly vertical and meets the axis far past it, every time. Newton assumes the function looks linear near the root, and x^{1/3} emphatically doesn't.
4. Newton: 10^{-1} \to 10^{-2} \to 10^{-4} \to 10^{-8} \to 10^{-16}. Four steps take you from 1 digit to 16.
Bisection gains \log_{10}2 \approx 0.301 digits per step, so reaching 15 more digits needs
\frac{15}{0.301} \approx \boxed{50 \text{ steps}}
Roughly 12 bisection steps per Newton step at this accuracy — and the gap widens as you demand more digits, because Newton's advantage is exponential in the digit count.
Check yourself in code
Run Newton's method for \sqrt2 and watch the error square each step.
Using x_{n+1} = x_n - \frac{x_n^2-2}{2x_n} from x_0 = 1, print 5 iterations: x_n to 15 decimals and |x_n - \sqrt2| in scientific notation with 2 decimals.
Print exactly this:
iter 1 x=1.500000000000000 err=8.58e-02
iter 2 x=1.416666666666667 err=2.45e-03
iter 3 x=1.414215686274510 err=2.12e-06
iter 4 x=1.414213562374690 err=1.59e-12
iter 5 x=1.414213562373095 err=0.00e+00
from math import sqrt
x = 1.0
for i in range(1, 6):
# one Newton step for f(x) = x^2 - 2, then report x and the error
print(f"iter {i} ...")
from math import sqrt
x = 1.0
for i in range(1, 6):
x = x - (x*x - 2) / (2*x)
print(f"iter {i} x={x:.15f} err={abs(x - sqrt(2)):.2e}")
Newton's method replaces f by its tangent line and solves that instead: x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}. Because linear approximation is second-order accurate, the error squares each iteration and correct digits double — four or five steps against bisection's fifty. The price is the guarantee: it can cycle, diverge, or find the wrong root, and it needs f' to stay away from zero. Applied to g' = 0 instead of f = 0, it becomes an optimizer, and that's the road to §14.
Next: a rule that makes every \frac00 limit from §1 mechanical.