[Python] Differentiating OdeSolution object

I would like to compute the residual of the numerical solution to an ordinary differential equation.

Here is the code

import numpy as np
from scipy.integrate import solve_ivp

def f(t, x):
return 0.5 * x

t_start = 0
t_end = 2

n = 50

t = np.linspace(t_start, t_end, 50)

x_init = 1

solution = solve_ivp(f, [t_start, t_end], [x_init], dense_output = True)


Now I would like to numerically differentiate the result at the points t. I tried naively

from scipy.differentiate import derivative

derivative(solution.sol, t)


And get the error that the sizes do not match.

I tried to overcome this by defining

def g(t):
return solution.sol(t)[0]

derivative(g, t)


and this also does not work.

How to fix these errors?

Continue reading...