Numerical integration and roots

Why this matters

Two questions come up constantly in engineering and science work:

• Where is it zero? The water depth in a channel that carries a given flow, the friction factor in the Colebrook equation, the temperature at which a reaction's equilibrium constant equals 1, the operating point of a diode circuit. These are root-finding problems: solve f(x) = 0 when no formula gives x directly.
• What is the total? Energy from a power-versus-time log, work from a force-displacement curve, the area under a chromatography peak, the distance from sampled velocity. These are integrals, often of data rather than of a formula.

Both have simple, reliable numerical methods you should understand, plus SciPy functions that do them robustly.

Root finding: bisection

If f is continuous and f(a) and f(b) have opposite signs, there is a root between them. Bisection halves the bracket each step:

def bisect(f, a, b, tol=1e-10):
    fa = f(a)
    if fa * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    while b - a > tol:
        m = 0.5 * (a + b)
        fm = f(m)
        if fa * fm <= 0:
            b = m                 # root is in [a, m]
        else:
            a, fa = m, fm         # root is in [m, b]
    return 0.5 * (a + b)

It is slow (one binary digit per step, about 34 steps for 1e-10 on a unit bracket) but guaranteed to converge once you have a bracket. That reliability is why it is the fallback in every serious solver.

Root finding: Newton's method

Newton follows the tangent line to where it crosses zero:

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

def newton(f, df, x0, tol=1e-12, max_iter=50):
    x = x0
    for _ in range(max_iter):
        step = f(x) / df(x)
        x -= step
        if abs(step) < tol:
            return x
    raise RuntimeError("Newton did not converge")

Near a root it converges quadratically: the number of correct digits roughly doubles each step. For √2 (root of x² - 2) from x0 = 1: 1.5, 1.41667, 1.4142157, 1.41421356237. But it needs the derivative, and from a bad start it can jump away or cycle. A good practice is to bracket first, then polish with Newton.

Worked example: open-channel flow

A rectangular channel 3 m wide must carry Q = 6 m³/s. With Manning's equation, the flow depth y satisfies f(y) = Q_manning(y) - 6 = 0. You cannot solve that for y by algebra, but f(0.1) < 0 and f(5) > 0, so bisection on [0.1, 5] finds the depth to millimetre precision in about 13 halvings. The same pattern solves the Colebrook friction factor, bubble-point temperatures and diode operating points.

Integration: the trapezoid rule

For samples y_i at points x_i, join neighbours with straight lines and add up the trapezoids:

∫ y dx ≈ Σ (x_{i+1} - x_i) (y_i + y_{i+1}) / 2

Vectorised, with no loop:

import numpy as np

def trapezoid(y, x):
    y = np.asarray(y, float); x = np.asarray(x, float)
    return float(np.sum(np.diff(x) * (y[1:] + y[:-1]) / 2))

It works for unevenly spaced data, which is what loggers often produce. Its error falls like h². NumPy has it built in: np.trapezoid(y, x) in NumPy 2.0 and later (older versions call it np.trapz).

Worked example: energy from a power log

A motor's power is logged every 10 s for a minute: [0, 1.2, 2.5, 2.5, 2.4, 1.0, 0] kW. The energy is trapezoid(P, t) = 96 kJ (with t in seconds and P in kW, the result is in kJ). Units of an integral are the units of y times the units of x.

Integration: Simpson's rule

Simpson fits parabolas through each pair of intervals. With an even number n of equal intervals of width h:

∫ f dx ≈ h/3 [ f0 + 4 f1 + 2 f2 + 4 f3 + ... + 4 f_{n-1} + f_n ]

Its error falls like h⁴. For ∫₀^π sin x dx = 2 with 10 intervals, the trapezoid rule gives 1.9835 and Simpson gives 2.00011.

The SciPy toolbox

from scipy.optimize import brentq, root_scalar, fsolve
from scipy.integrate import quad, simpson, trapezoid

root = brentq(lambda y: f(y), 0.1, 5)          # bracketing, fast and safe
value, err = quad(np.sin, 0, np.pi)              # adaptive, for functions you can call
area = simpson(y, x=x)                           # for sampled data

• brentq combines bisection's safety with faster interpolation; it is the default choice for a bracketed scalar root.
• fsolve and root handle systems of non-linear equations.
• quad adapts its sampling to the function and returns an error estimate. Use it when you can evaluate the function; use trapezoid or simpson when you only have samples.

Common mistakes

• No bracket check: bisection on an interval without a sign change silently returns garbage. Raise an error instead.
• Newton from a bad start or with a derivative near zero: it diverges. Bracket first.
• Odd number of intervals for Simpson's rule.
• Forgetting units: the integral of kW over seconds is kJ, not kWh.
• Integrating noisy data with high-order rules: Simpson amplifies noise less than differentiation does, but smoothing or a fitted model may be better.

Quick recap

• Bisection: needs a sign change, always converges, slowly. Newton: fast near a root, needs f' and a good start.
• Trapezoid: sum(diff(x) * (y[1:] + y[:-1]) / 2), any spacing, error ~ h². Simpson: even intervals, error ~ h⁴.
• In real work: brentq for roots, quad for functions, trapezoid and simpson for data.

Try it

The Do tab asks you to write bisection, a vectorised trapezoid rule and Simpson's rule, and test them on known answers.