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

如果代数操作失败，可以对约束进行数值求解，例如在每个时间步运行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是否找到了一个好的解决方案，并刷新输出，这样您就可以实现它，并停止模拟。在

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

^{pr2}$

请注意，为了使技巧生效，cache参数必须是可变的，这就是我使用list的原因。如果您不熟悉默认参数的工作方式，请参阅this链接。在

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