Soft realtime simulation with scipy solve_ivp

Question

I've a system of ODEs which I so far solved via solve_ivp.

scipy.integrate.solve_ivp(fun=model, t_span=(0.0, t_end), y0=[s0])

My problem is, that I want to solve the ODEs in a running simulation, where new values are constantly injected into the simulation and results are displayed. The simulation may run for several hours. My approach was to call solve_ivp repeatedly, roughly as follows (this allows me to show intermediate results and fetch new data, which is not shown here):

t = 0.0
s = s0
while t < t_end:
    result = scipy.integrate.solve_ivp(fun=model, t_span=(t, t + t_step), y0=[s])
    s = result.y[0][-1]
    t += t_step

I wrote a few test cases, where I solved the ODEs analytically and for reasons I do not fully understand yet, repeatedly calling solve_ivp is always closer to the analytical solution (without manually tweaking the different parameters of solve_ivp). My question is more along the lines: If there is something wrong in this approach, or if there is some functionality in scipy or another package which may be better suited to my needs?

user_11423283 · Answer

A poorly documented option is to use the integrators that underly solve_ivp directly. For example, RK45. When an external event changes something, you would need to restart the integration as most (all?) of these integrators use multiple prior steps to construct the value for the next time point. If you change values or the derivative function mid-integration, you will introduce subtle hard to understand bugs.

Here is an example that changes the derivative function based on an external event to the integrator. In this case, it is simply based on the clock time. Note that this example will give different results on different machines depending on speed of execution.

import datetime
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate

def fun1(t, y):
    return -y
    
def fun2(t, y):
    return 2 * y
    
def use_fun2():
    return datetime.datetime.utcnow().time().microsecond >= 500000
    
time_bound = 20.

# max_step used to make sure this runs slow enough
# try changing this to see the difference
max_step = 0.0001
y0 = [5]
fun=fun1
rk45 = integrate.RK45(fun, 0, y0, time_bound, max_step=max_step)
t = []
y = []

while rk45.status == "running" and rk45.y[0] < 200:
    if use_fun2():
        if fun is fun1:
            fun = fun2
            rk45 = integrate.RK45(fun, rk45.t, rk45.y, time_bound, max_step=max_step)
    else:
        if fun is fun2:
            fun = fun1
            rk45 = integrate.RK45(fun, rk45.t, rk45.y, time_bound, max_step=max_step)
    
    t.append(rk45.t)
    y.append(rk45.y)
    rk45.step()
    
import matplotlib.pyplot as plt
plt.plot(t, y)
plt.show()

Soft realtime simulation with scipy solve_ivp

1 Answers