Bisection Method — Midpoint Overflow at 1e308 Boundaries

Bisection is the simplest application of the Intermediate Value Theorem: if f is continuous and f(a) < 0 and f(b) > 0, then somewhere in [a,b] the function crosses zero. The algorithm just applies this fact repeatedly, halving the interval each time.

The guaranteed convergence of bisection makes it invaluable as a safety net. Bisection cannot diverge, cannot oscillate, cannot fail as long as the initial bracket is valid. This is why Brent's method — the gold standard for numerical root finding, used in scipy.optimize.brentq — combines bisection as a fallback with Newton-Raphson's quadratic convergence: the speed of Newton when near the root, the reliability of bisection everywhere else.

Algorithm and Implementation

Invariant: f(a) and f(b) have opposite signs → root in [a,b]. Each step halves the interval.

There's a subtle trap in the midpoint calculation. The naive (a+b)/2 can overflow when a and b are large and opposite signs. Production code uses a + (b-a)/2, which is safe.

The termination check uses abs(b-a)/2 < tol, not abs(f(mid)) < tol — that's a common source of early false convergence. The interval width is the true error bound.

def bisection(f, a: float, b: float, tol: float = 1e-10, max_iter: int = 100) -> float:
    """
    Find root of f in [a,b] where f(a)*f(b) < 0.
    Guaranteed to converge for continuous functions.
    """
    if f(a) * f(b) > 0:
        raise ValueError('f(a) and f(b) must have opposite signs')
    for i in range(max_iter):
        mid = a + (b - a) / 2          # overflow-safe midpoint
        if (b - a) / 2 < tol:
            print(f'Converged in {i+1} iterations, root ≈ {mid}')
            return mid
        if f(a) * f(mid) < 0:
            b = mid
        else:
            a = mid
    return a + (b - a) / 2

import math
root = bisection(lambda x: x**2 - 2, 1, 2)
print(f'sqrt(2) = {root}')

Error Analysis and Convergence

After n iterations, the interval width is (b-a)/2^n. The error in the root position is at most this width.

To achieve precision ε: n ≥ log2((b-a)/ε) iterations.

For [1,2] with ε=10^-10: n ≥ log2(1/10^-10) = 10×log2(10) ≈ 33.2 → 34 iterations.

This is linear convergence — one binary digit per iteration. Compare to Newton-Raphson which doubles correct digits each step. But bisection never fails to get those digits, Newton sometimes doesn't converged at all.

import math

def bisection_with_trace(f, a, b, n_steps=10):
    fa = f(a)
    print(f'Step  a              b              mid            f(mid)')
    for i in range(n_steps):
        mid = a + (b - a) / 2
        fmid = f(mid)
        print(f'{i+1:4}  {a:.10f}  {b:.10f}  {mid:.10f}  {fmid:+.2e}')
        if fa * fmid < 0:
            b, _ = mid, fmid
        else:
            a, fa = mid, fmid

bisection_with_trace(lambda x: x**2 - 2, 1.0, 2.0, 8)

Bisection vs Newton-Raphson

Best practice: Use bisection to get close (10^-3 precision), then switch to Newton-Raphson for fast final convergence. This is essentially Brent's method.

But there's a hidden cost: Newton-Raphson requires computing the derivative. For functions that are expensive to evaluate or have no analytic derivative, bisection remains the only option.

When Bisection Fails: Invalid Brackets and Discontinuities

Bisection's guarantee depends on two assumptions: 1. f(a) and f(b) have opposite signs 2. f is continuous on [a,b]

If the function has a vertical asymptote (like 1/x at 0), a sign change doesn't guarantee a root — you'll converge to a singularity. Similarly, if there are multiple roots inside the interval, bisection will only find one (the one that triggers the sign change at each step).

Another failure mode: the bracket is valid but the function is flat at the root, causing f(mid) to be exactly zero prematurely. Your abs(b-a)/2 < tol check saves you here, but if you mistakenly use abs(f(mid)) < tol you might converge to a point where f is tiny but the interval is still large.

Brent's Method: The Production-Grade Hybrid

Brent's method combines bisection's reliability with Newton-Raphson's speed. It uses inverse quadratic interpolation when possible, but falls back to bisection when interpolation would go out of bounds or converge too slowly.

The result: guaranteed convergence, and typically converges superlinearly (faster than bisection, slower than pure Newton). In practice it's the default choice for one-dimensional root-finding in scipy, MATLAB's fzero, and many other numerical libraries.

Implementation complexity is higher than plain bisection, but you don't write it yourself — you use a library. Understanding the internals helps when debug convergence issues.

from scipy.optimize import brentq
import math

# Use Brent's method (no manual bisection needed)
root = brentq(lambda x: x**2 - 2, 1, 2)
print(f'sqrt(2) = {root:.15f}')

# Even works near singularities if the bracket avoids them
root2 = brentq(lambda x: 1/x - 1, 0.5, 2.0)
print(f'1/x - 1 = 0 at x = {root2}')

Frequently Asked Questions