求解隐式ODE(微分代数方程DAE) [英] Solve an implicit ODE (differential algebraic equation DAE)

问题描述

我正在尝试使用来自scipy的odeint解决二阶ODE.我遇到的问题是该函数隐式地耦合到二阶项，如简化的代码段所示(请忽略示例的假装物理特性):

I'm trying to solve a second order ODE using odeint from scipy. The issue I'm having is the function is implicitly coupled to the second order term, as seen in the simplified snippet (please ignore the pretend physics of the example):

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的求值之间存在很多逻辑，代数运算失败.

in this case I realise it is possible to algebraically solve for the implicit variable, however in my actual scenario there is a lot of logic between F_r and the evaluation of a and algebraic manipulation fails.

我相信可以使用MATLAB的 ode15i 来解决DAE.功能，但我尽一切可能避免这种情况.

I believe the DAE could be solved using MATLAB's ode15i function, but I'm trying to avoid that scenario if at all possible.

我的问题是-有办法解决python(最好是scipy)中的隐式ODE函数(DAE)吗?并有更好的方法解决上述问题吗?

My question is - is there a way to solve implicit ODE functions (DAE) in python( scipy preferably)? And is there a better way to pose the problem above to do so?

作为最后的选择，从上一个时间步长开始通过a是可以接受的.在每个时间步之后如何将dydt[1]传递回函数?

As a last resort, it may be acceptable to pass a from the previous time-step. How could I pass dydt[1] back into the function after each time-step?

推荐答案

如果代数运算失败，则可以寻求约束的数值解决方案，例如在每个时间步长运行fsolve:

if algebraic manipulation fails, you can go for a numerical solution of your constraint, running for example fsolve at each timestep:

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找到了一个好的解决方案，并刷新输出，以便您可以在发生时立即意识到并停止模拟.

Clearly this will slow down your time integration. Always check that fsolve finds a good solution, and flush the output so that you can realize it as it happens and stop the simulation.

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

About how to "cache" the value of a variable at a previous timestep, you can exploit the fact that default arguments are calculated only at the function definition,

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参数必须是可变的，这就是我使用列表的原因.请参阅此链接不熟悉默认参数的工作原理.

Notice that in order for the trick to work the cache parameter must be mutable, and that's why I use a list. See this link if you are not familiar with how default arguments work.

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

Notice that the two codes DO NOT produce the same result, and you should be very careful using the value at the previous timestep, both for numerical stability and precision. The second is clearly much faster though.

这篇关于求解隐式ODE(微分代数方程DAE)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！