python - 求解隐式 ODE(微分代数方程 DAE)

Question

我正在尝试使用 scipy 中的 odeint 求解二阶 ODE。我遇到的问题是该函数隐式耦合到二阶项，如简化片段所示(请忽略示例的假物理):

import numpy as np
from scipy.integrate import odeint

def integral(y,t,F_l,mass):
    dydt = np.zeros_like(y)
    x, v = y
    F_r =  (((1-a)/3)**2 + (2*(1+a)/3)**2) * v # 'a' implicit 
    a  = (F_l - F_r)/mass

    dydt = [v, a]
return dydt


y0 = [0,5]
time = np.linspace(0.,10.,21)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon,mass))

在这种情况下，我意识到可以代数求解隐式变量，但是在我的实际场景中，F_r 和 a 的求值之间有很多逻辑> 并且代数运算失败。

我相信可以使用 MATLAB 的 ode15i 解决 DAE功能，但我尽可能避免这种情况。

我的问题是 - 有没有办法解决 python 中的隐式 ODE 函数 (DAE)(最好是 scipy)？是否有更好的方法来解决上述问题？

作为最后的手段，从上一个时间步传递 a 可能是可以接受的。我如何在每个时间步后将 dydt[1] 传递回函数？

Answer 1

很旧，但值得更新，因此它可能对偶然发现此问题的任何人都有用。目前 python 中可以解决隐式 ODE 的包非常少。 GEKKO ( https://github.com/BYU-PRISM/GEKKO ) 是软件包之一，专门针对混合整数、非线性优化问题进行动态优化，但也可以用作通用 DAE 求解器。

上述“假装物理学”问题可以在 GEKKO 中按如下方式解决。

m= GEKKO()
m.time = np.linspace(0,100,101)
F_l = m.Param(value=1000)
mass = m.Param(value =1000)
m.options.IMODE=4
m.options.NODES=3
F_r = m.Var(value=0)
x = m.Var(value=0)
v = m.Var(value=0,lb=0)
a = m.Var(value=5,lb=0)
m.Equation(x.dt() == v)
m.Equation(v.dt() == a)
m.Equation (F_r ==  (((1-a)/3)**2 + (2*(1+a)/3)**2 * v)) 
m.Equation (a == (1000 - F_l)/mass)
m.solve(disp=False)
plt.plot(x)