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.

How the Bisection Method Actually Finds Roots

The bisection method is a root-finding algorithm that repeatedly halves an interval known to contain a sign change. Given a continuous function f(x) and an interval [a, b] where f(a) and f(b) have opposite signs, it computes the midpoint m = (a + b) / 2, evaluates f(m), and replaces either a or b with m depending on the sign of f(m). This guarantees linear convergence: each iteration halves the error, so 10 iterations give about 3 decimal digits of accuracy.

In practice, the method is simple, robust, and guaranteed to converge — but it's slow. It requires O(log((b-a)/ε)) iterations for tolerance ε. Unlike Newton's method, it needs no derivative and works on any continuous function. The critical implementation detail: the midpoint calculation (a + b) / 2 can overflow when a and b are near Double.MAX_VALUE (~1.8e308). Use a + (b - a) / 2 instead.

Use bisection when you need guaranteed convergence and can tolerate linear speed — for example, calibrating sensor thresholds, finding breakpoints in piecewise models, or as a fallback when derivative-based methods fail. It's the safety net of numerical root-finding: never diverges, never needs tuning, but never fast.

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}')

Assumptions That Kill Your Root-Finding Run

The bisection method isn't magic. It's a direct application of the Intermediate Value Theorem, and that theorem comes with non-negotiable preconditions. Skip the sanity checks, and you'll be debugging silent garbage.

First, your function must be continuous on the closed interval [a, b]. Not "mostly continuous" or "continuous enough." Fully continuous. If there's a jump discontinuity at the midpoint, the method will happily converge on the wrong side of the crack.

Second, f(a) and f(b) must have opposite signs. f(a) * f(b) < 0. That's the root guarantee. If both are positive or both negative, you have zero confidence a root exists inside the bracket. The method will narrow the interval mindlessly, giving you a precise answer to a question you shouldn't have asked.

Third, you need exactly one root in that interval. Multiple roots with the same sign? The theorem is silent. The method will find one, but you won't know which, or if it even found the one you wanted.

Always validate these three before the first iteration. Production code is not a guessing game.

// io.thecodeforge — dsa tutorial

public class BracketValidator {
    public static boolean validBracket(double a, double b, Function<Double, Double> f) {
        double fa = f.apply(a);
        double fb = f.apply(b);
        
        // Check continuity is your problem - assume f is continuous
        // Minimum: opposite signs
        if (fa * fb >= 0) {
            return false; // No sign change, no guarantee
        }
        return true;
    }
    
    public static void main(String[] args) {
        Function<Double, Double> func = x -> x * x - 4;
        System.out.println(validBracket(0, 3, func)); // false
        System.out.println(validBracket(1, 3, func)); // true
    }
}

The Only Example You'll Ever Need: x³ - x - 2

Let's make this concrete. You're chasing the root of f(x) = x³ - x - 2. This polynomial has one real root near 1.5. We'll use bisection to find it to three decimal places.

Pick a = 1, b = 2. f(1) = -2, f(2) = 4. Signs are opposite. We're good. The midpoint c = 1.5. f(1.5) = -0.125. Still negative. Since f(a) and f(c) share a sign, the root must be between c and b. So we shift: a = 1.5, b stays 2.

Next iteration: c = 1.75. f(1.75) = 1.609. Positive. Now f(a) is negative, f(c) is positive → root between a and c. Set b = 1.75.

Keep halving. By iteration 10, your interval width is under 0.001. The midpoint is 1.521. f(1.521) ≈ 0.006. That's your root to three decimals.

Each iteration cuts the error in half. Linear convergence is slow but reliable. You trade speed for certainty. Use it when correctness matters more than cycles.

// io.thecodeforge — dsa tutorial

import java.util.function.Function;

public class BisectionExample {
    public static double findRoot(Function<Double, Double> f, double a, double b, double tol) {
        if (f.apply(a) * f.apply(b) >= 0) throw new IllegalArgumentException("Bad bracket");
        double c = a;
        while ((b - a) / 2 > tol) {
            c = (a + b) / 2;
            if (f.apply(c) == 0.0) break;
            if (f.apply(a) * f.apply(c) < 0) b = c;
            else a = c;
        }
        return c;
    }

    public static void main(String[] args) {
        Function<Double, Double> f = x -> x * x * x - x - 2;
        double root = findRoot(f, 1.0, 2.0, 1e-6);
        System.out.printf("Root: %.6f\n", root);
        System.out.printf("f(root): %.6f\n", f.apply(root));
    }
}