Root finding and special functions

Bracket a root or Newton-iterate towards one, read the return value honestly, and use scipy.special instead of re-deriving the functions yourself.

Choosing a root finder

import numpy as np
from scipy import optimize

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

root_b, info = optimize.brentq(f, 2, 3, xtol=1e-12, full_output=True)

root_n = optimize.newton(f, x0=2.0, fprime=lambda x: 3 * x ** 2, tol=1e-12)

sys = lambda v: [v[0] ** 2 + v[1] ** 2 - 4, v[0] - v[1]]
root_f, out = optimize.fsolve(sys, [1.0, 1.0], full_output=True)

print(root_b, root_n, root_f, out["ier"])

Function	Needs	Guarantee
`brentq`	A bracket with a sign change	Converges if the bracket exists and f is continuous
`brenth`	Same as brentq	Similar; slightly different convergence order
`newton`	A start point, optionally a derivative	Local; can diverge or leave the domain
`fixed_point`	A contraction mapping	Requires the derivative magnitude below one near the solution
`fsolve`	A starting vector	Least-squares style; not guaranteed to find a root
`root`	A method name and start	A common interface over several algorithms

Using them honestly

# always check convergence, do not just take the number
res = optimize.root_scalar(f, bracket=[2, 3], method="brentq")
print(res.converged, res.iterations, res.root, f(res.root))

# scan for brackets when you have no idea where the roots are
xs = np.linspace(-5, 5, 200)
fs = f(xs)
candidates = [(a, b) for a, b, fa, fb in zip(xs[:-1], xs[1:], fs[:-1], fs[1:]) if fa * fb < 0]
print(len(candidates))

# polynomials: get every root at once
print(np.roots([1, 0, -2, -5]))

Pass full_output=True or use root_scalar; a returned number is not a promise that it is a root.
Sanity-check by evaluating f(root). If it is far from zero, the tolerance was not met.
For a polynomial, use numpy.roots — it finds all roots at once and is more robust than repeated bracketing.

Special functions you will meet

from scipy import special
import numpy as np

print(special.gammaln(1000))              # log(999!) without overflow
print(special.comb(50, 25, exact=True))   # exact integer binomial

print(special.erfc(3.0))                  # 1 - erf, stable in the tail
print(special.logsumexp([1000.0, 999.0, 998.0]))

j = special.jv(0, np.linspace(0, 10, 5))  # Bessel functions for wave problems
print(j)

def softmax(z):
    return np.exp(z - special.logsumexp(z))

⚠️

Computing 1 - erf(x) or exp(-x) for large x destroys precision through cancellation and underflow. erfc, expm1, log1p and logsumexp exist precisely to keep accuracy in those regimes.

FAQ

brentq raises that f(a) and f(b) must have different signs. What does that mean?

There is no root, an even number of roots, or a discontinuity inside the bracket. Sample the interval, plot it, or widen the bracket. A sign change is a necessary condition, not a guarantee of a unique root.

Why is fsolve giving an answer that is not a solution?

fsolve minimises the residual from a starting guess and can settle somewhere that merely looks flat given rounding error. Supply an analytic Jacobian, check ier in the output, and re-run from several start points to see whether the answer is stable.

Optimisation and curve fitting Interpolation, smoothing and integration

Last refreshed 2026-09-18.

Root finding and special functions

Choosing a root finder

Using them honestly

Special functions you will meet

FAQ

Related