How to Differentiate solve_ivp Solutions in SciPy Without Shape Mismatch: vectorize or use np.gradient

Learn to avoid SciPy solve_ivp derivative shape-mismatch errors: vectorize the OdeSolution interpolator or differentiate samples with NumPy gradient. Guide.

When you solve an ODE with SciPy and then try to differentiate the interpolated solution with scipy’s derivative, it’s easy to run into shape and broadcasting issues. The typical symptom is a size-mismatch error when you pass the OdeSolution interpolator directly to the differentiator or wrap it in a simple function. The fix is straightforward once you know how scipy’s derivative expects to call your function.

Problem setup

Consider the exponential growth model and its numerical solution. After integrating, the goal is to compute the derivative of the solution at a set of time points.

import numpy as np
from scipy.integrate import solve_ivp
def drift(t, u):
    return 0.5 * u
t0 = 0.0
t1 = 2.0
m = 50
tau = np.linspace(t0, t1, m)
u0 = 1.0
soln = solve_ivp(drift, [t0, t1], [u0], dense_output=True)

A naive attempt to differentiate the interpolator directly fails:

from scipy import differentiate
differentiate.derivative(soln.sol, tau)

Wrapping the interpolator into a scalar-returning function, but without vectorization, doesn’t solve it either:

def wrapped(s):
    return soln.sol(s)[0]
differentiate.derivative(wrapped, tau)

What’s going on

The object returned by solve_ivp with dense_output=True exposes an interpolator via soln.sol. This interpolator can evaluate the solution at arbitrary times, returning a vector for systems or, in this scalar case, a one-element array. The derivative routine from scipy expects to probe your function over an array of inputs and receive outputs in a way that is compatible with vectorized evaluation. If your callable isn’t vectorized, or if it returns shapes that don’t broadcast as expected, you’ll get size-mismatch errors.

There are two reliable ways forward. You can compute the numerical gradient on the sampled solution values with NumPy. Or, if you want to stick with scipy’s derivative, you can explicitly vectorize the interpolator so it behaves like a scalar-valued function that accepts array inputs.

Solution 1: use NumPy’s gradient on the sampled solution

Solving on a predefined grid and differentiating those samples is the most direct route. The derivative is numerical in this approach.

import numpy as np
from scipy import integrate
def rhs(t, y):
    return 0.5 * y
grid = np.linspace(0.0, 2.0, 500)
ans = integrate.solve_ivp(rhs, [grid.min(), grid.max()], y0=[1.0], t_eval=grid)
num_grad = np.gradient(ans.y[0], ans.t)

Here ans.y[0] are the sampled solution values and ans.t equals grid. The result num_grad is the numerical derivative at each time point.

Solution 2: vectorize the interpolator and use scipy’s derivative

If you prefer scipy’s differentiate.derivative, you need to make the interpolator act like a vectorized scalar function. Vectorizing ensures it can be evaluated consistently over arrays of inputs, avoiding shape mismatches. The derivative computed this way is also numerical.

import numpy as np
from scipy import integrate, differentiate
def rhs(t, y):
    return 0.5 * y
grid = np.linspace(0.0, 2.0, 500)
ans = integrate.solve_ivp(rhs, [grid.min(), grid.max()], y0=[1.0], t_eval=grid)
@np.vectorize
def interp_eval(s):
    return ans.sol(s)[0]
rich_out = differentiate.derivative(interp_eval, grid)
deriv_vals = rich_out.df

The call returns a _RichResult object. The derivative values are in the df field, along with additional metadata such as error estimates and iteration counts. This gives you a differentiated signal tied directly to the integrator’s internal interpolant.

Why this matters

Post-processing ODE solutions often includes computing residuals, sensitivities, or consistency checks, all of which rely on derivatives. It’s important to recognize that both paths shown here produce numerical derivatives. Using NumPy’s gradient differentiates the sampled solution values. Using scipy’s derivative differentiates the interpolated solution produced by solve_ivp. In practical workflows, both approaches agree well for smooth problems like the exponential model shown here.

Takeaways

If you need derivatives on the solution grid, solving with t_eval and applying np.gradient is concise and robust. If you want to differentiate the continuous interpolant that solve_ivp provides, ensure the callable is vectorized before passing it to scipy’s differentiate.derivative. In both cases you avoid shape-mismatch errors and get a consistent, fully numerical estimate of the time derivative suitable for residual evaluation and further analysis.