使用 scipy solve_ivp 进行软实时模拟答案

【问题标题】：Soft realtime simulation with scipy solve_ivp使用 scipy solve_ivp 进行软实时模拟
【发布时间】：2021-12-01 01:10:37
【问题描述】：

我有一个 ODE 系统，到目前为止我通过 solve_ivp 解决了它。

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

我的问题是，我想在运行模拟中求解 ODE，其中新值不断注入模拟并显示结果。模拟可能会运行几个小时。我的做法是反复调用solve_ivp，大致如下（这样可以展示中间结果，获取新数据，这里没有展示）：

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

我写了一些测试用例，在其中我解析求解了 ODE，由于我还不完全理解的原因，重复调用 solve_ivp 总是更接近解析解（无需手动调整 solve_ivp 的不同参数）。我的问题更接近于：如果这种方法有问题，或者 scipy 或其他包中的某些功能可能更适合我的需求？

【问题讨论】：

标签： python scipy

【解决方案1】：

一个记录不充分的选项是直接使用solve_ivp 下的集成器。例如，RK45。当外部事件发生变化时，您需要重新启动集成，因为这些集成器中的大多数（全部？）都使用多个先前的步骤来构造下一个时间点的值。如果您在集成过程中更改值或导数函数，您将引入难以理解的细微错误。

这里是一个基于外部事件改变微分函数到积分器的例子。在这种情况下，它只是基于时钟时间。请注意，此示例将根据执行速度在不同机器上给出不同的结果。

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()

【讨论】：