解决一个隐式常微分方程组 (代数微分方程 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（首选scipy）中解决隐式ODE函数（DAE）？以及是否有更好的方式来表达上述问题以解决它？

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

- Shaun_M

---

你也可以看一下DAE工具。 - winkmal

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)

- Siva-Sg

---

你好，我想问一下你是否知道如何解决具有非常数系数的二阶隐式问题？例如，考虑方程$f(x)y''(x) + g(x) y'(x)*y(x)$。快速查看文档并没有给出明确的答案，我也没有找到涵盖此情况的示例。 - Hamilcar

在Gekko（或其他DAE软件包）中添加2阶或更高阶导数没有问题。只需为每个超过1阶的阶数定义一个新的微分状态即可。这是一个例子：https://apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential Gekko还可以解决更高（3+）指数的DAE，而许多其他DAE积分器仅支持指数-1（或某些指数-2 Hessenberg形式）：https://apmonitor.com/wiki/index.php/Apps/PendulumMotion - John Hedengren

Answer 2

如果代数运算失败，您可以选择进行约束的数值解法，例如在每个时间步长运行fsolve:

import sys
from numpy import linspace
from scipy.integrate import odeint
from scipy.optimize import fsolve

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

def F_r(a, v):
    return (((1 - a) / 3) ** 2 + (2 * (1 + a) / 3) ** 2) * v

def constraint(a, v):
    return (F_lon - F_r(a, v)) / mass - a

def integral(y, _):
    v = y[1]
    a, _, ier, mesg = fsolve(constraint, 0, args=[v, ], full_output=True)
    if ier != 1:
        print "I coudn't solve the algebraic constraint, error:\n\n", mesg
        sys.stdout.flush()
    return [v, a]

dydt = odeint(integral, y0, time)

显然，这会减慢您的时间积分速度。始终检查fsolve能否找到好的解，并刷新输出，以便您可以实时了解并停止模拟。

关于如何在以前的时间步骤中“缓存”变量的值，您可以利用默认参数仅在函数定义时计算这一事实，

from numpy import linspace
from scipy.integrate import odeint

#you can choose a better guess using fsolve instead of 0
def integral(y, _, F_l, M, cache=[0]):
    v, preva = y[1], cache[0]
    #use value for 'a' from the previous timestep
    F_r = (((1 - preva) / 3) ** 2 + (2 * (1 + preva) / 3) ** 2) * v 
    #calculate the new value
    a = (F_l - F_r) / M
    cache[0] = a
    return [v, a]

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

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

请注意，为了使技巧生效，cache参数必须是可变的，这就是为什么我使用一个列表的原因。如果您不熟悉默认参数值的工作方式，请参阅此链接。

请注意，这两个代码并不产生相同的结果，你应该非常小心地使用前一时间步的值，这涉及数值稳定性和准确性。然而第二个代码显然更快。

- gg349

---

嗨，谢谢flebol。对这两种解决方案进行基本的分析，结果如下：%timeit odeint(integral1, y0, time) 100次循环，3次中最好的时间为9.03毫秒 %timeit odeint(integral2, y0, time, args=(F_lon, mass)) 1000次循环，3次中最好的时间为972微秒 其中integral1是第一个示例。在这些条件下，这两个示例之间的差异相当小（对于1000个时间步长的数量级为10^-7），但是，改变一些值（例如'F_lon=1000; mass=100'）第二种方法会失败。因此，我可以接受时间惩罚。谢谢。 - Shaun_M

@Shaun_M，9.03毫秒和972微秒之间的差异是10的因数，第一个解决方案慢了10倍。 - gg349

抱歉表述不够清晰。我的意思是指解决方案之间的数值差异。我试图表达的观点是，在这种情况下，使用缓存积分可以接受精度，但稳定性不行。谢谢。 - Shaun_M

@Shaun_M 感谢您提供了一种有趣的解决 DAEs 的方法。您的解决方案是否可以适用于具有奇异“质量矩阵”的 ODE 问题，例如 M*du/dt=f(u)？ - Graham G