Root Finding with SciPy

A root is a value of x where a function equals zero — equivalently, it's the solution to the equation f(x) = 0. Root finding comes up constantly: solving equilibrium conditions in physics, finding where two curves cross, computing break-even points in finance, or solving polynomial equations analytically intractable. SciPy provides `brentq` for reliable 1D root finding and `root` for systems of nonlinear equations.

### Finding a Single Root with brentq

`brentq` finds the root of a function within a bracket `[a, b]` where the function changes sign — meaning f(a) and f(b) have opposite signs, so a root must lie between them. It combines bisection, secant, and inverse quadratic interpolation to converge reliably and quickly.

from scipy.optimize import brentq
import numpy as np
import matplotlib.pyplot as plt

f = lambda x: np.cos(x) - x  # root of cos(x) = x (fixed point)

root = brentq(f, 0, 2)
print(f"Root at x = {root:.6f}")
print(f"Verify: cos({root:.4f}) = {np.cos(root):.6f}")

x = np.linspace(0, 2, 300)
plt.figure(figsize=(7, 4))
plt.plot(x, f(x), color="steelblue", linewidth=2, label="cos(x) − x")
plt.axhline(0, color="gray", linewidth=0.8, linestyle="--")
plt.plot(root, 0, "o", color="crimson", markersize=10, zorder=5,
         label=f"Root at x = {root:.4f}")
plt.title("Root of cos(x) = x")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.tight_layout()
plt.show()

Root at x = 0.739085
Verify: cos(0.7391) = 0.739085

- `brentq(f, 0, 2)` requires f(0) and f(2) to have opposite signs — `f(0) = cos(0)−0 = 1 > 0` and `f(2) = cos(2)−2 ≈ −2.4 < 0`, so the bracket is valid.
- The root of `cos(x) = x` is known as the Dottie number (≈ 0.7391) — it's the unique value where repeated application of cosine converges.
- If you pass a bracket where f doesn't change sign, `brentq` raises a `ValueError`. Always plot first when you're unsure where the root lies.

### Finding Multiple Roots

When a function has several roots, you find them by identifying multiple sign-change brackets — often by plotting first — and calling `brentq` once per bracket. This is more robust than a single call to a generic solver.

from scipy.optimize import brentq
import numpy as np
import matplotlib.pyplot as plt

f = lambda x: np.sin(x) * np.exp(-0.1 * x)

# Identify sign-change brackets by inspecting the plot
intervals = [(2, 4), (5, 7), (8, 10), (11, 13)]
roots = [brentq(f, a, b) for a, b in intervals]
print("Roots:", [f"{r:.4f}" for r in roots])

x = np.linspace(0, 14, 500)
plt.figure(figsize=(9, 3.5))
plt.plot(x, f(x), color="steelblue", linewidth=2)
plt.axhline(0, color="gray", linewidth=0.8, linestyle="--")
plt.plot(roots, [0] * len(roots), "o", color="crimson", markersize=9,
         zorder=5, label=f"{len(roots)} roots found")
plt.title("Multiple Roots of sin(x)·exp(−0.1x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.tight_layout()
plt.show()

Roots: ['3.1416', '6.2832', '9.4248', '12.5664']

- Each call to `brentq` finds exactly one root within its bracket, so you get as many roots as intervals you specify.
- The exponential envelope `exp(-0.1x)` damps the amplitude over time but doesn't shift the zero-crossings of `sin(x)`, which remain at multiples of π.
- In practice, scan for sign changes automatically: iterate over consecutive sample pairs and collect any `(x[i], x[i+1])` pair where `f(x[i]) * f(x[i+1]) < 0`.

### Solving a System of Nonlinear Equations

`root` generalizes root finding to multiple variables: it solves f(x) = 0 where both x and f are vectors. You provide the system as a function that takes a vector and returns a vector of residuals.

from scipy.optimize import root
import numpy as np
import matplotlib.pyplot as plt

# Find intersection of a circle (radius 2) and a line (x - y = 1)
def system(xy):
    x, y = xy
    return [
        x**2 + y**2 - 4,   # circle: x² + y² = 4
        x - y - 1           # line: x - y = 1
    ]

sol = root(system, x0=[1.5, 0.5])
print(f"Solution: x={sol.x[0]:.4f}, y={sol.x[1]:.4f}")
print(f"Residual: {[f'{r:.2e}' for r in system(sol.x)]}")

theta = np.linspace(0, 2 * np.pi, 300)
x_line = np.linspace(-2, 2, 200)
plt.figure(figsize=(5, 5))
plt.plot(2 * np.cos(theta), 2 * np.sin(theta), color="steelblue", linewidth=2, label="Circle x²+y²=4")
plt.plot(x_line, x_line - 1, color="orange", linewidth=2, label="Line x−y=1")
plt.plot(sol.x[0], sol.x[1], "o", color="crimson", markersize=10, zorder=5, label="Root")
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.axhline(0, color="gray", linewidth=0.5)
plt.axvline(0, color="gray", linewidth=0.5)
plt.gca().set_aspect("equal")
plt.title("Intersection of Circle and Line")
plt.legend()
plt.tight_layout()
plt.show()

Solution: x=1.8229, y=0.8229
Residual: ['0.00e+00', '0.00e+00']

- `root(system, x0=[1.5, 0.5])` takes the function and an initial guess vector; it solves the full system simultaneously.
- The residual `system(sol.x)` should be close to `[0, 0]` — checking it confirms the solution is accurate, not just converged to a local failure.
- The system has two intersections; `x0` steers the solver toward the one nearby. To find the other, try `x0=[-0.8, -1.8]`.

`brentq` is the workhorse for 1D root finding — reliable whenever you can bracket the root. For optimization (finding minima rather than zeros), see [SciPy optimize basics](/tutorials/scipy-optimize-basics).