在 Python 中模拟 MATLAB 的 ode15s

Question

我正在致力于将模型从 MATLAB 转换为 Python。该模型的关键在于 MATLAB 的 ode15。在 MATLAB 执行中，ode15s 有标准选项：

options = odeset()
[t P] = ode15s(@MODELfun, tspan, y0, options, params)

作为参考，y0 是一个向量（大小为 98），MODELfun 也是如此。

我的 Python 等效尝试如下：

ode15s = scipy.integrate.ode(Model.fun)
ode15s.set_integrator('vode', method = 'bdf', order = 15)
ode15s.set_initial_value(y0).set_f_params(params)
dt = 1 
while ode15s.successful() and ode15s.t < duration:
     ode15s.integrate(ode15s.t+dt)

但这似乎不起作用。有什么建议或替代方案吗？

编辑： 查看输出后，我从 Python 得到的结果要么是 y0 的某些元素随着时间的推移没有变化，要么是 y0 的其余部分在每一步都不断变化。有类似的经历吗？

Answer 1

根据 SciPy wiki for Matlab Users，使用 ode15s 的正确方法是

scipy.integrate.ode(f).set_integrator('vode', method='bdf', order=15)

Answer 2

需要明确的一点是，与 Matlab 的 ode15s 不同，scipy 积分器“vode”不支持具有质量矩阵的模型。因此任何建议都应包含此警告。

Answer 3

对于任何遇到这个问题的人，integrate.solve_ivp 允许在求解器选项中包含质量矩阵。 BDF 或 LSODA 应该相当于求解器 ode15s。就我而言，我的实现如下：

 def dae_with_params(t, states_solved):
        return dae(t, states_solved, u_k)
 
solver_options = {'method': 'LSODA', 'rtol': 1e-3, 'atol': 1e-3, 'vectorized': False, 'jac': None, 'mass': M}

sol = solve_ivp(dae_with_params, [0, step_size], initial_states, **solver_options)