解隐式常微分方程（微分代数方程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]传回这个函数呢？

Answer 1

虽然这个问题有点老，但更新一下还是很有价值，可能对任何偶然看到这个问题的人都有帮助。目前在Python中，有不少可以解决隐式常微分方程（ODE）的包。

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)

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 参数必须是可变的，所以我用的是列表。如果你不太了解默认参数是怎么工作的，可以看看这个链接。

要注意，这两段代码的结果是不同的，使用上一个时间步的值时要非常小心，这关系到数值的稳定性和精确度。不过，第二段代码明显快得多。