How can I solve ODEs for a set number of time steps using Python (SciPy)?

Question

I am trying to solve a set of ODEs using SciPy. The task I am given asks me to solve the differential equation for 500 time steps. How can I achieve this using SciPy?

So far, I have tried using scipy.integrate.solve_ivp, which gives me a correct solution, but I cannot control the number of time steps that it runs for. The t_span argument lets me configure what the initial and final values of t are, but I'm not actually interested in that -- I am only interested in how many times I integrate. (For example, when I run my equations with t_span = (0, 500) the solver integrates 907 times.)

Below is a simplified example of my code:

from scipy.integrate import solve_ivp

def simple_diff(t, z) :
    x, y = z
    return [1 - 2*x*y, 2*x - y]

t_range = (0, 500)
xy_init = [0, 0]

sol = solve_ivp(simple_diff, t_range, xy_init)

I am also fine with using something other than SciPy, but solutions with SciPy are preferable.

Answer 1

You can use the t_eval argument to solve_ivp to evaluate at particular time points:

import numpy as np

t_eval = np.arange(501)
sol = solve_ivp(simple_diff, t_range, xy_init, t_eval=t_eval)

However, note that this will not cause the solver to limit the number of integration steps - that is determined by error metrics.

If you absolutely must evaluate the function exactly 500 times to obtain 500 integration steps, you are describing Euler integration, which will be less accurate than the algorithm that solve_ivp uses.

Looking at the solutions to your equation, it feels like you probably want to integrate only up to t=5.

Here's what the result looks like when integrating with the above settings:

And here's the result for

t_eval = np.linspace(0, 5)
t_range = (0, 5)
sol = solve_ivp(simple_diff, t_range, xy_init, t_eval=t_eval)