SciPy Root Finding: root, root_scalar and brentq

SciPy gives you three doors into root finding, and picking the wrong one is why scipy.optimize.root feels awkward for simple problems.

• brentq — one equation, and you know two points where the sign flips. It’s fast, and it always converges.
• root_scalar — one equation, with a choice of method behind one interface.
• root — a system of equations, taking and returning a vector.

Most people land here with a single equation, so it’s one of the first two they want. I reach for brentq unless there’s a reason not to.

Everything that follows was run on SciPy 1.18.1 on Python 3.12.5; the screenshots are that terminal.

The same job, three ways

Here are all three solving something you can check by eye, the square root of 2:

from scipy.optimize import brentq, root, root_scalar

# one equation, one unknown, and you know a bracket where the sign flips
print("brentq      :", brentq(lambda x: x**2 - 2, 0, 2))

# the same job through the general scalar interface
print("root_scalar :", root_scalar(lambda x: x**2 - 2, bracket=[0, 2]).root)

# a system of equations: pass a vector, get a vector back
def system(v):
    x, y = v
    return [x + y - 3, x - y - 1]

print("root        :", root(system, [0, 0]).x)

Output:

brentq      : 1.4142135623731364
root_scalar : 1.4142135623731364
root        : [2. 1.]

Notice the shapes. The scalar functions hand back a number, while root returns an object whose .x is an array, because a system has one answer per unknown.

root_scalar and the methods behind it

root_scalar is a front door to several algorithms. Which one actually runs depends on what you hand it:

from scipy.optimize import root_scalar

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

bracketing = root_scalar(f, bracket=[1, 2], method="brentq")
print("brentq  ->", f"root {bracketing.root:.10f}", "| iterations", bracketing.iterations,
      "| converged", bracketing.converged)

bisection = root_scalar(f, bracket=[1, 2], method="bisect")
print("bisect  ->", f"root {bisection.root:.10f}", "| iterations", bisection.iterations)

# newton needs a starting guess instead of a bracket, and a derivative helps
newton = root_scalar(f, x0=1.5, fprime=lambda x: 3 * x**2 - 1, method="newton")
print("newton  ->", f"root {newton.root:.10f}", "| iterations", newton.iterations)

print("f(root) =", f"{f(bracketing.root):.2e}")

Output:

brentq  -> root 1.5213797068 | iterations 8 | converged True
bisect  -> root 1.5213797068 | iterations 39
newton  -> root 1.5213797068 | iterations 4
f(root) = 0.00e+00

Pass bracket and you get a bracketing method. Pass x0 and you get an open one such as newton, which is faster when it works.

Bisection is the tortoise here, and that’s the point. It can’t fail once it has a bracket, it just takes its time.

Method	What it needs	Guaranteed?	Use when
brentq	A bracket	Yes	The default choice for one equation
bisect	A bracket	Yes	You want the simplest possible behaviour
newton	x0, ideally fprime	No	You have a derivative and a good guess
secant	x0 and x1	No	No derivative available

Why brentq asks for a sign change

A bracketing method works by shrinking an interval it knows contains a root. It knows that because the function’s positive at one end and negative at the other.

Give it two points on the same side and it refuses before doing any work:

from scipy.optimize import brentq

# brentq needs the function to change sign between the two ends
f = lambda x: x**2 - 2

print("f(0) =", f(0), " f(2) =", f(2), "-> sign changes, fine")
print("root:", brentq(f, 0, 2))

print("f(2) =", f(2), " f(3) =", f(3), "-> both positive")
brentq(f, 2, 3)

What SciPy raises:

f(0) = -2  f(2) = 2 -> sign changes, fine
root: 1.4142135623731364
f(2) = 2  f(3) = 7 -> both positive
Traceback (most recent call last):
  File "C:\pyguides\scipy_brentq_sign_error.py", line 10, in <module>
    brentq(f, 2, 3)
  File "C:\pyguides\venv\Lib\site-packages\scipy\optimize\_zeros_py.py", line 858, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: f(a) and f(b) must have different signs

This is the SciPy root-finding error I get asked about most, and it’s usually honest. Your interval really doesn’t contain a root, or it holds two and they cancel out.

Plot the function first if you’re unsure. A quick sign check across a few points is enough to find a bracket.

When the solver does not converge

Bracketing methods can’t fail. Open methods can, and they tell you rather than throwing:

from scipy.optimize import root_scalar

# newton has no bracket to fall back on, so a bad guess can wander off
flat = lambda x: x**3 - 2 * x + 2

result = root_scalar(flat, x0=0, fprime=lambda x: 3 * x**2 - 2, method="newton")
print("converged:", result.converged)
print("flag     :", result.flag)
print("iterations:", result.iterations)

# the bracketing method has no such problem
safe = root_scalar(flat, bracket=[-3, 0], method="brentq")
print("brentq   :", f"{safe.root:.10f}", "| converged", safe.converged)

Output:

converged: False
flag     : convergence error
iterations: 50
brentq   : -1.7692923542 | converged True

Always read .converged on a scalar result, or .success on a system result, before you trust the number next to it. A non-converged result still carries a .root, and it’s meaningless.

Systems of equations with root

This is what scipy.optimize.root is for. Your function takes a vector and returns one residual per equation, and you want them all at zero:

import numpy as np
from scipy.optimize import root

# two curves crossing: x^2 + y^2 = 25 and y = x - 1
def equations(v):
    x, y = v
    return [x**2 + y**2 - 25, y - x + 1]

solution = root(equations, [1, 1])

print("success :", solution.success)
print("message :", solution.message)
print("solution:", np.round(solution.x, 10))
print("residual:", np.round(equations(solution.x), 12))

other = root(equations, [-5, -5])       # a different guess finds the other crossing
print("second root:", np.round(other.x, 10))

Output:

success : True
message : The solution converged.
solution: [4. 3.]
residual: [0. 0.]
second root: [-3. -4.]

The residuals matter more than the message. Printing them is the cheapest way to confirm the answer really is a root and not just where the solver gave up.

Two crossings, two answers. Nonlinear systems don’t have one solution, so your guess is part of the question, much as it is for constrained optimisation.

Giving it a Jacobian

Supply the derivatives and the solver stops estimating them by finite differences:

import numpy as np
from scipy.optimize import root

def equations(v):
    x, y = v
    return [x**2 + y**2 - 25, y - x + 1]

def jacobian(v):
    x, y = v
    return [[2 * x, 2 * y],
            [-1, 1]]

without = root(equations, [1, 1])
with_jac = root(equations, [1, 1], jac=jacobian)

print(f"no jacobian  : {without.nfev:>3} function calls")
print(f"with jacobian: {with_jac.nfev:>3} function calls, {with_jac.njev} jacobian calls")
print("same answer  :", np.allclose(without.x, with_jac.x))

Output:

no jacobian  :  14 function calls
with jacobian:  12 function calls, 1 jacobian calls
same answer  : True

On two equations the saving is small. On twenty it’s the difference between a solve you wait for and one you don’t.

More SciPy and NumPy on this site:

• Linear programming with SciPy linprog
• NumPy data types
• Create a 2D array with NumPy
• Find the maximum value in an array
• Round numbers to 2 decimal places
• Getting started with SciPy

Frequently asked questions