Root-Finding Methods

Solving the fundamental equation \(f(x) = 0\): from the rock-solid bisection method to the blazingly fast Newton-Raphson iteration and its relatives.

The Root-Finding Problem

Given a continuous function \(f: \mathbb{R} \to \mathbb{R}\), find \(x^*\) such that \(f(x^*) = 0\). This seemingly simple problem arises everywhere:

- Finding eigenvalues: \(\det(A - \lambda I) = 0\)
- Implicit equations of state in thermodynamics
- Intersection points of curves and surfaces
- Inverse problems: given \(y\), find \(x\) such that \(g(x) = y\) (i.e., \(f(x) = g(x) - y = 0\))

The key question is always: how fast does the method converge?

1. The Bisection Method

The simplest and most robust root-finding method. By the Intermediate Value Theorem, if \(f(a)\) and \(f(b)\) have opposite signs and \(f\) is continuous on \([a,b]\), then there exists at least one root in \((a,b)\).

Algorithm

Start with \([a, b]\) where \(f(a) \cdot f(b) < 0\)
Compute midpoint \(c = (a + b)/2\)
If \(f(a) \cdot f(c) < 0\), set \(b = c\); otherwise set \(a = c\)
Repeat until \(|b - a| < \epsilon\)

Convergence: Linear. After \(n\) iterations, the interval has width \(|b-a|/2^n\). To get \(p\) digits of accuracy, need \(n \approx 3.32 \, p\) iterations. Slow but guaranteed.

Bisection Method Implementation

Python

script.py38 lines

import numpy as np

def bisection(f, a, b, tol=1e-12, max_iter=100):
    """Bisection method with convergence tracking."""
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")

history = []
    for i in range(max_iter):
        c = (a + b) / 2.0
        fc = f(c)
        history.append((i, c, abs(b - a), fc))

if abs(fc) < tol or abs(b - a) < tol:
            return c, history

if f(a) * fc < 0:
            b = c
        else:
            a = c

return c, history

# Example: Find root of x^3 - x - 2 = 0 (real root near 1.52)
f = lambda x: x**3 - x - 2

root, history = bisection(f, 1.0, 2.0)

print("Bisection Method: x^3 - x - 2 = 0")
print(f"{'Iter':>4}  {'c (midpoint)':>18}  {'interval width':>15}  {'f(c)':>12}")
print("-" * 55)
for i, c, width, fc in history:
    print(f"{i:4d}  {c:18.14f}  {width:15.2e}  {fc:12.4e}")

print(f"\nRoot: {root:.15f}")
print(f"Verification: f(root) = {f(root):.2e}")
print(f"Iterations needed: {len(history)}")
print(f"Linear convergence: interval halves each step")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. Newton-Raphson Method

The crown jewel of root finding. Linearize \(f\) at the current guess using the tangent line and solve for the next approximation:

\(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\)

Newton's iteration: the tangent line at \((x_n, f(x_n))\) intersects the x-axis at \(x_{n+1}\)

Convergence Analysis

Taylor expand \(f(x^*)\) around \(x_n\) where \(x^* = x_n + e_n\):

\(0 = f(x^*) = f(x_n) + f'(x_n) e_n + \frac{1}{2} f''(x_n) e_n^2 + \cdots\)

The Newton step gives \(e_{n+1} = e_n + f(x_n)/f'(x_n)\), which yields:

\(e_{n+1} \approx -\frac{f''(x^*)}{2 f'(x^*)} e_n^2\) — quadratic convergence!

The number of correct digits approximately doubles each iteration. But: requires \(f'(x_n) \neq 0\) and a good initial guess.

Newton-Raphson with Convergence Comparison

Python

script.py45 lines

import numpy as np

def newton(f, df, x0, tol=1e-14, max_iter=50):
    """Newton-Raphson with error tracking."""
    x = x0
    history = [(0, x, abs(f(x)))]
    for i in range(1, max_iter + 1):
        fx, dfx = f(x), df(x)
        if abs(dfx) < 1e-15:
            print("Derivative too small!")
            break
        x = x - fx / dfx
        history.append((i, x, abs(f(x))))
        if abs(f(x)) < tol:
            break
    return x, history

# Solve x^3 - x - 2 = 0
f = lambda x: x**3 - x - 2
df = lambda x: 3*x**2 - 1

root_newton, hist_newton = newton(f, df, x0=1.5)

print("Newton-Raphson: x^3 - x - 2 = 0")
print(f"{'Iter':>4}  {'x_n':>20}  {'|f(x_n)|':>15}")
print("-" * 45)
for i, x, fx in hist_newton:
    print(f"{i:4d}  {x:20.15f}  {fx:15.2e}")

# Show quadratic convergence
print("\n--- Quadratic Convergence Verification ---")
exact = root_newton
print(f"{'Iter':>4}  {'|e_n|':>15}  {'|e_{n+1}|/|e_n|^2':>20}")
print("-" * 45)
for i in range(len(hist_newton) - 1):
    en = abs(hist_newton[i][1] - exact)
    en1 = abs(hist_newton[i+1][1] - exact)
    if en > 1e-15:
        ratio = en1 / en**2
        print(f"{i:4d}  {en:15.2e}  {ratio:20.6f}")
    else:
        print(f"{i:4d}  {en:15.2e}  {'converged':>20}")

print(f"\nRoot: {root_newton:.16f}")
print(f"|e_(n+1)| / |e_n|^2 converges to |f''/(2f')| at root")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. The Secant Method

What if we don't know \(f'(x)\)? Replace the derivative with a finite difference approximation using the two most recent iterates:

\(x_{n+1} = x_n - f(x_n) \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}\)

Uses a secant line through \((x_{n-1}, f(x_{n-1}))\) and \((x_n, f(x_n))\)

Convergence Order

The secant method has superlinear convergence of order \(p = \phi = (1 + \sqrt{5})/2 \approx 1.618\) (the golden ratio!). Not as fast as Newton (\(p=2\)), but each step requires only one function evaluation instead of two (no derivative needed). Per function evaluation, secant can be more efficient than Newton.

Secant Method Implementation

Python

script.py37 lines

import numpy as np

def secant(f, x0, x1, tol=1e-14, max_iter=50):
    """Secant method: derivative-free Newton variant."""
    history = [(0, x0, abs(f(x0))), (1, x1, abs(f(x1)))]
    for i in range(2, max_iter + 2):
        f0, f1 = f(x0), f(x1)
        if abs(f1 - f0) < 1e-15:
            break
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        history.append((i, x2, abs(f(x2))))
        if abs(f(x2)) < tol:
            break
        x0, x1 = x1, x2
    return x2, history

f = lambda x: x**3 - x - 2

root, hist = secant(f, 1.0, 2.0)

print("Secant Method: x^3 - x - 2 = 0")
print(f"{'Iter':>4}  {'x_n':>20}  {'|f(x_n)|':>15}")
print("-" * 45)
for i, x, fx in hist:
    print(f"{i:4d}  {x:20.15f}  {fx:15.2e}")

# Verify superlinear convergence with order phi = 1.618
exact = root
print("\n--- Convergence Order Estimation ---")
errors = [abs(h[1] - exact) for h in hist if abs(h[1] - exact) > 1e-15]
for i in range(2, len(errors)):
    if errors[i-1] > 1e-15 and errors[i-2] > 1e-15:
        # Estimate order: log(e_n) / log(e_{n-1})
        p = np.log(errors[i] / errors[i-1]) / np.log(errors[i-1] / errors[i-2])
        print(f"Step {i}: estimated order p = {p:.3f} (golden ratio = {(1+np.sqrt(5))/2:.3f})")

print(f"\nRoot: {root:.16f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Fixed-Point Iteration

Rewrite \(f(x) = 0\) as \(x = g(x)\) and iterate \(x_{n+1} = g(x_n)\). The Banach Fixed-Point Theorem guarantees convergence if \(|g'(x)| < 1\) in a neighborhood of the fixed point.

Convergence Criterion

If \(|g'(x^*)| = L < 1\), then the error satisfies:

\(|e_{n+1}| \leq L |e_n|\) — linear convergence with rate \(L\)

The smaller \(L\), the faster convergence. If \(g'(x^*) = 0\), we get at least quadratic convergence (Newton's method is a special case!).

Fixed-Point Iteration: Convergent vs Divergent

Python

script.py39 lines

import numpy as np

def fixed_point(g, x0, tol=1e-12, max_iter=100):
    """Fixed-point iteration x_{n+1} = g(x_n)."""
    x = x0
    history = [(0, x)]
    for i in range(1, max_iter + 1):
        x_new = g(x)
        history.append((i, x_new))
        if abs(x_new - x) < tol:
            break
        if abs(x_new) > 1e10:
            print(f"  Diverging at iteration {i}!")
            break
        x = x_new
    return x, history

# Solve x^3 - x - 2 = 0, rewritten as x = (x + 2)^(1/3)
# g(x) = (x+2)^(1/3), g'(x) = 1/(3(x+2)^(2/3))
# At root ~1.52: g'(1.52) ~ 0.13 < 1 => converges!
print("=== Convergent: x = (x + 2)^(1/3) ===")
g1 = lambda x: (x + 2)**(1.0/3.0)
root1, hist1 = fixed_point(g1, 1.0)
for i, x in hist1[:12]:
    print(f"  Iter {i:3d}: x = {x:.14f}")
print(f"  Root: {root1:.15f}\n")

# Same equation, BAD rearrangement: x = x^3 - 2
# g(x) = x^3 - 2, g'(x) = 3x^2
# At root ~1.52: g'(1.52) ~ 6.93 >> 1 => diverges!
print("=== Divergent: x = x^3 - 2 ===")
g2 = lambda x: x**3 - 2
root2, hist2 = fixed_point(g2, 1.5, max_iter=10)
for i, x in hist2[:8]:
    print(f"  Iter {i:3d}: x = {x:.6f}")

print("\n=== Lesson ===")
print("The SAME equation can converge or diverge depending on rearrangement!")
print("|g'(x*)| < 1 is essential for convergence.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Method Comparison

Method	Order	f evals/step	Requires	Guaranteed?
Bisection	\(O(1/2^n)\)	1	Bracket	Yes
Newton-Raphson	\(p = 2\)	2 (f, f')	Derivative, good \(x_0\)	Locally
Secant	\(p \approx 1.618\)	1	Two initial points	Locally
Fixed-Point	\(p = 1\) (usually)	1	\(\|g'(x^*)\| < 1\)	If contraction

Head-to-Head: Convergence Race

Python

script.py48 lines

import numpy as np

f = lambda x: x**3 - x - 2
df = lambda x: 3*x**2 - 1
exact_root = 1.5213797068045676  # pre-computed

# Bisection
a, b = 1.0, 2.0
bisect_errors = []
for _ in range(50):
    c = (a + b) / 2
    bisect_errors.append(abs(c - exact_root))
    if f(a) * f(c) < 0: b = c
    else: a = c
    if abs(b - a) < 1e-15: break

# Newton
x = 1.5
newton_errors = [abs(x - exact_root)]
for _ in range(15):
    x = x - f(x) / df(x)
    newton_errors.append(abs(x - exact_root))
    if abs(f(x)) < 1e-15: break

# Secant
x0, x1 = 1.0, 2.0
secant_errors = [abs(x0 - exact_root), abs(x1 - exact_root)]
for _ in range(20):
    f0, f1 = f(x0), f(x1)
    if abs(f1 - f0) < 1e-15: break
    x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
    secant_errors.append(abs(x2 - exact_root))
    x0, x1 = x1, x2
    if abs(f(x2)) < 1e-15: break

print("Convergence Race: |error| vs iteration")
print(f"{'Iter':>4}  {'Bisection':>14}  {'Newton':>14}  {'Secant':>14}")
print("-" * 52)
max_show = max(len(bisect_errors), len(newton_errors), len(secant_errors))
for i in range(min(max_show, 15)):
    b_err = f"{bisect_errors[i]:.6e}" if i < len(bisect_errors) else "converged"
    n_err = f"{newton_errors[i]:.6e}" if i < len(newton_errors) else "converged"
    s_err = f"{secant_errors[i]:.6e}" if i < len(secant_errors) else "converged"
    print(f"{i:4d}  {b_err:>14}  {n_err:>14}  {s_err:>14}")

print("\nNewton reaches machine precision in ~6 iterations")
print("Secant takes ~8 iterations (but no derivative needed)")
print("Bisection takes ~50 iterations (but always converges)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Practical Pitfalls

Multiple Roots

If \(f(x^*) = f'(x^*) = 0\) (double root), Newton degrades to linear convergence. Fix: use modified Newton \(x_{n+1} = x_n - m \frac{f(x_n)}{f'(x_n)}\) where\(m\) is the multiplicity.

Cycling & Divergence

Newton can cycle between points or diverge to infinity with a bad initial guess. Practical codes use hybrid methods: start with bisection to get close, then switch to Newton for fast convergence.

Near-Zero Derivatives

If \(f'(x_n) \approx 0\), the Newton step becomes huge. Always check for small derivatives and fall back to a safer method.

Catastrophic Cancellation

In the secant method, \(f(x_n) - f(x_{n-1})\) can suffer from cancellation when iterates are close. Use extended precision or reformulate.

7. Brent's Method: The Gold Standard

In practice, most scientific libraries use Brent's method (1973), which combines bisection, secant, and inverse quadratic interpolation. It has the reliability of bisection with the speed of higher-order methods. This is what scipy.optimize.brentq uses.

Using SciPy's Brent Method

Python

script.py32 lines

from scipy import optimize
import numpy as np

f = lambda x: x**3 - x - 2

# Brent's method (the standard in scientific computing)
result = optimize.brentq(f, 1.0, 2.0, full_output=True)
root = result[0]
info = result[1]

print("Brent's Method (scipy.optimize.brentq)")
print(f"Root: {root:.16f}")
print(f"Function evaluations: {info.function_calls}")
print(f"Iterations: {info.iterations}")
print(f"f(root) = {f(root):.2e}")

# Compare with Newton using scipy
print("\n--- scipy.optimize.newton ---")
root_newton = optimize.newton(f, x0=1.5, fprime=lambda x: 3*x**2 - 1)
print(f"Root: {root_newton:.16f}")

# Find ALL roots of a more complex function
print("\n--- Finding multiple roots ---")
g = lambda x: np.sin(x) - x/5
print(f"Roots of sin(x) - x/5 = 0:")
# We know roots are near 0, +/- 2.5
for x0 in [-3, 0, 3]:
    try:
        r = optimize.newton(g, x0)
        print(f"  x0 = {x0}: root = {r:.10f}, f(root) = {g(r):.2e}")
    except RuntimeError:
        print(f"  x0 = {x0}: failed to converge")

Click Run to execute the Python code

Code will be executed with Python 3 on the server