我正在尝试使用来自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]
传递回函数中? 最佳答案
如果代数运算失败,则可以寻求约束的数值解决方案,例如在每个时间步上运行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
参数必须是可变的,这就是我使用列表的原因。如果您不熟悉默认参数的工作原理,请参见this链接。请注意,这两个代码不会产生相同的结果,因此在数值上的稳定性和精度上,您应该非常小心地使用上一个时间步长中的值。第二个显然要快得多。
关于python - 求解隐式ODE(微分代数方程DAE),我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/23578596/