Как корректно дифференцировать интерполятор solve_ivp в SciPy

Почему derivative в SciPy падает на OdeSolution? Даем 2 решения: np.gradient по t_eval и векторизацию интерполятора solve_ivp для стабильной производной.

Когда вы решаете ОДУ в SciPy, а затем пытаетесь продифференцировать интерполированное решение с помощью функции derivative из SciPy, легко столкнуться с проблемами форм массивов и бродкастинга. Типичный симптом — ошибка несоответствия размеров, если передать интерполятор OdeSolution напрямую в дифференциатор или обернуть его в простую функцию. Исправление вполне прямолинейно, стоит лишь понять, как derivative из SciPy ожидает вызывать вашу функцию.

Постановка задачи

Рассмотрим модель экспоненциального роста и её численное решение. После интегрирования нужно вычислить производную решения в наборе моментов времени.

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)

Наивная попытка дифференцировать интерполятор напрямую приводит к сбою:

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

Обёртка интерполятора в функцию, возвращающую скаляр, но без векторизации, тоже не помогает:

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

Что происходит

Объект, возвращаемый solve_ivp с dense_output=True, предоставляет интерполятор через soln.sol. Этот интерполятор умеет вычислять решение в произвольные моменты времени, возвращая вектор для систем или, в скалярном случае, массив из одного элемента. Функция численного дифференцирования в SciPy ожидает опрашивать вашу функцию на массиве входных значений и получать ответы в форме, совместимой с векторизованной обработкой. Если вызываемый объект не векторизован или возвращает формы, которые не бродкастятся как ожидается, возникнут ошибки несоответствия размеров.

Есть два надёжных варианта. Можно посчитать численный градиент по дискретным значениям решения средствами NumPy. Либо, если хотите остаться с scipy’s derivative, явно векторизовать интерполятор, чтобы он вёл себя как скалярная функция, принимающая массивы.

Решение 1: использовать градиент из NumPy на дискретных значениях

Решить на заранее заданной сетке и продифференцировать эти выборки — самый прямой путь. В этом подходе производная получается численно.

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)

Здесь ans.y[0] — выборочные значения решения, а ans.t совпадает с grid. Результат num_grad — численная производная в каждой точке времени.

Решение 2: векторизовать интерполятор и использовать scipy’s derivative

Если вам удобнее scipy’s differentiate.derivative, нужно заставить интерполятор вести себя как векторизованная скалярная функция. Векторизация обеспечивает корректную работу на массивах входов и избегает конфликтов форм. Полученная таким образом производная тоже является численной.

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

Вызов возвращает объект _RichResult. Значения производной находятся в поле df, а также доступны дополнительные метаданные — например, оценки ошибок и число итераций. Так вы получаете дифференцированный сигнал, непосредственно связанный с внутренним интерполянтом интегратора.

Зачем это нужно

Послеработка с решениями ОДУ часто включает вычисление невязок, чувствительностей или проверок согласованности — всё это опирается на производные. Важно понимать, что оба показанных пути дают численные производные. NumPy.gradient дифференцирует дискретные значения решения. SciPy’s derivative дифференцирует интерполированное решение, которое формирует solve_ivp. На практических задачах оба подхода хорошо согласуются для гладких проблем, как в показанной модели экспоненциального роста.

Выводы

Если нужны производные в узлах сетки решения, решайте с t_eval и применяйте np.gradient — это кратко и надёжно. Если хотите дифференцировать непрерывный интерполянт, который предоставляет solve_ivp, убедитесь, что вызываемая функция векторизована, прежде чем передавать её в scipy’s differentiate.derivative. В обоих случаях вы избегаете ошибок несоответствия размеров и получаете согласованную, полностью численную оценку временной производной, подходящую для расчёта невязок и дальнейшего анализа.