odeint() 中的 Sympy 表达式给出推导错误

Question

我正在尝试自动化几个过程，例如使用 Sympy 从拉格朗日量生成 ODE，并使用 Numpy 和 Scipy 对它们进行数值积分。 完整代码在最后。 作为使用 solve() 生成 ODE 的结果，我得到了一个包含 Sympy 表达式的字典，如下所示：

{Derivative(lambda1(t), t): (y(t) + 1)/(x(t)*y(t)),
 Derivative(z(t), t): x(t),
 Derivative(x(t), t): y(t)*z(t),
 Derivative(y(t), t): -x(t)*z(t)
 }

然后据此我想将Scipy的微分方程组与odeint()积分。为此，我需要从 def Field(Q,t): 中的字典中提取表达式（例如 lambdify），以引入为 odeint(Field,Q_0,t_array)。这是我运行遇到困难的地方：

我第一次尝试

def Equ2(nQ,t,Q,Field):
x1,y1,z1,lamb1 = nQ
dQ =[]
for f in Q:
    dQ.append(lambdify(Q, Field[f.diff(t)],'numpy' )(x1,y1,z1,lamb1))
return dQ[0:len(nQ)]

但它不能转到 odeint()，因为它需要带有 2 个参数的字段，我试图将它传递到 odeint() 的可选 arga=() 中，给了我一个（长）错误：

ValueError                   Traceback (most recent call last)
<ipython-input-20-63f086b8a252> in Equ2(nQ, t, Q, Field)
      20     dQ =[]
      21     for f in Q:
 ---> 22         dQ.append(lambdify(Q, Field[f.diff(t)],'numpy' )(x1,y1,z1,lamb1))
      23     return dQ[0:len(Q)-1]

[...]
ValueError: 
Can't calculate 1st derivative wrt 14.0430379424125.

所以我尝试了基本相同但没有循环，

def Equ1(nQ,t):
x1,y1,z1,lamb1 = nQ
dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy' )(x1,y1,z1,lamb1)
dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy' )(x1,y1,z1,lamb1)
dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy' )(x1,y1,z1,lamb1)
dlam = lambdify((x,y,z,lam[0]), field[lam[0].diff(t)],'numpy' )(x1,y1,z1,lamb1)
return [dx,dy,dz]

并且有（我认为）同样的问题：

ValueError                                Traceback (most recent call last)
<ipython-input-20-63f086b8a252> in Equ1(nQ, t)
      9 def Equ1(nQ,t):
     10     x1,y1,z1,lamb1 = nQ
---> 11     dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy' )(x1,y1,z1,lamb1)
     12     dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy' )(x1,y1,z1,lamb1)
     13     dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy' )(x1,y1,z1,lamb1)

[...]
ValueError: 
Can't calculate 1st derivative wrt 17.6326726993661.

如果我简单地尝试：

def Equ0(nQ,t):
x,y,z,lamb = nQ
dx = y*z
dy = -x*z
dz = x
dlam = (y+1.)/(x*y)
return [dx,dy,dz]

集成工作正常。另外，如果我调用 EquX() 函数时使用与 odeint() 内部类似的参数，它们就可以正常工作。

完整代码

from sympy import *
from sympy.physics.mechanics import dynamicsymbols
from numpy import linspace, sin, cos
from scipy.integrate import odeint

t = Symbol('t')
x = Function('x')(t)
y = Function('y')(t)
z = Function('z')(t)

lam = dynamicsymbols('lambda1:{0}'.format(5))
f = x.diff(t)- y*z
eq = Matrix([x.diff(t) - lam[0].diff(t)*y*x*z+z,
         y.diff(t) +x*z,
         z.diff(t)-x
        ])

field = solve(list(eq)+[f],[x.diff(t),y.diff(t),z.diff(t),lam[0].diff(t)])


def Equ0(nQ,t):
    x,y,z,lamb = nQ
    dx = y*z
    dy = -x*z
    dz = x
    dlam = (y+1.)/(x*y)
    return [dx,dy,dz]

def Equ1(nQ,t):
    x1,y1,z1,lamb1 = nQ
    dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy' )(x1,y1,z1,lamb1)
    dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy' )(x1,y1,z1,lamb1)
    dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy' )(x1,y1,z1,lamb1)
    dlam = lambdify((x,y,z,lam[0]), field[lam[0].diff(t)],'numpy' )(x1,y1,z1,lamb1)
    return [dx,dy,dz]


def Equ2(nQ,t,Q,Field):
    x1,y1,z1,lamb1 = nQ
    dQ =[]
    for f in Q:
        dQ.append(lambdify(Q, Field[f.diff(t)],'numpy' )(x1,y1,z1,lamb1))
    return dQ[0:len(Q)-1]

q = [x,y,z,lam[0]]
nq = [1,2,3,4]
time=linspace(0,10,10)

### This line works just fine:
print Equ0(nq,t), Equ1(nq,t), Equ2(nq,t,q,field) #They give the same output 

sol0 = odeint(Equ0,nq,time)
sol1 = odeint(Equ1,nq,time) #Errors here
sol2 = odeint(Equ2,nq,time,args=(q,field)) #And here

最后是完整的错误：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
 <ipython-input-20-63f086b8a252> in Equ1(nQ, t)
  9 def Equ1(nQ,t):
 10     x1,y1,z1,lamb1 = nQ
 ---> 11     dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy' )(x1,y1,z1,lamb1)
 12     dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy' )(x1,y1,z1,lamb1)
 13     dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy' )(x1,y1,z1,lamb1)

/usr/local/lib/python2.7/dist-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions)
2864         new_symbols = list(map(sympify, symbols))  # e.g. x, 2, y, z---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-63f086b8a252> in Equ1(nQ, t)
  9 def Equ1(nQ,t):
 10     x1,y1,z1,lamb1 = nQ
 ---> 11     dx = lambdify((x,y,z,lam[0]), field[x.diff(t)],'numpy' )(x1,y1,z1,lamb1)
 12     dy = lambdify((x,y,z,lam[0]), field[y.diff(t)],'numpy' )(x1,y1,z1,lamb1)
 13     dz = lambdify((x,y,z,lam[0]), field[z.diff(t)],'numpy' )(x1,y1,z1,lamb1)

  /usr/local/lib/python2.7/dist-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions)
  2864         new_symbols = list(map(sympify, symbols))  # e.g. x, 2, y, z
2865         assumptions.setdefault("evaluate", True)
---> 2866         return Derivative(self, *new_symbols, **assumptions)
2867 
2868     ###########################################################################

/usr/local/lib/python2.7/dist-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions)
1068                 ordinal = 'st' if last_digit == 1 else 'nd' if last_digit == 2 else 'rd' if last_digit == 3 else 'th'
1069                 raise ValueError(filldedent('''
---> 1070                 Can\'t calculate %s%s derivative wrt %s.''' % (count, ordinal, v)))
1071 
1072             if all_zero and not count == 0:

ValueError: 
Can't calculate 1st derivative wrt 0.0.

2865         assumptions.setdefault("evaluate", True)
---> 2866         return Derivative(self, *new_symbols, **assumptions)
2867 
2868     ###########################################################################

/usr/local/lib/python2.7/dist-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions)
1068                 ordinal = 'st' if last_digit == 1 else 'nd' if last_digit == 2 else 'rd' if last_digit == 3 else 'th'
1069                 raise ValueError(filldedent('''
---> 1070                 Can\'t calculate %s%s derivative wrt %s.''' % (count, ordinal, v)))
1071 
1072             if all_zero and not count == 0:

ValueError: 
Can't calculate 1st derivative wrt 0.0.

TL;DR odeint() 内部出现一些推导错误，我无法使用自定义函数在 odeint() 外部重现。

Answer 1

如果您想使用 t 作为符号，您应该避免在函数声明中将 t 声明为浮点数。尝试用另一个名称 s 或 tt 或 ...

替换浮点数 t

odeint() 中的 Sympy 表达式给出推导错误

Sympy expression in odeint() gives derivation error

python

scipy

sympy

ode

derivative