求解隐式ODE（微分代数方程DAE）

求解隐式ODE（微分代数方程DAE）

如果代数运算失败，则可以

fsolve

求出约束的数值解，例如在每个时间步上运行：

import sysfrom numpy import linspacefrom scipy.integrate import odeintfrom scipy.optimize import fsolvey0 = [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) * vdef constraint(a, v):    return (F_lon - F_r(a, v)) / mass - adef 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:nn", mesg        sys.stdout.flush()    return [v, a]dydt = odeint(integral, y0, time)

显然，这会减慢您的时间整合。始终检查

fsolve

找到一个好的解决方案，并刷新输出，以便您可以在发生时立即意识到并停止模拟。

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

from numpy import linspacefrom scipy.integrate import odeint#you can choose a better guess using fsolve instead of 0def 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

参数必须是可变的，这就是为什么我使用列表。如果您不熟悉默认参数的工作方式，请参见此链接。

请注意，这两个代码不会产生相同的结果，因此在数值上的稳定性和精度上，您应该非常小心地使用上一个时间步长的值。第二个显然要快得多。