使用scipy.integrate.odeint时出现TypeError

Question

我正在尝试用 scipy.integrate.odeint 来解决一组耦合的微分方程。但是当我运行程序时，出现了以下错误：

TypeError: 无法将数组数据从 dtype('O') 转换为 dtype('float64')，根据规则 'safe' odepack.error: 函数调用的结果不是一个正确的浮点数组。

这是我使用的代码：

V = "v**2/2*log(1+x**2 + (y/a)**2 + (z/c)**2)" 
var = ['x','y','z']

def afleiden(func, var):
    f = sympify(func)
    partAfg = [f.diff(var[i]) for i in range(len(var))]
    return partAfg



init=[0.3,0.2,0.9,0.2,0.6,0.7] 

def func(rv, t, pot, var):
    return rv[3:6] + afleiden(pot,var) 

# rv is a list with 6 elements of witch the last 3 are part of the diff equations

t = np.arange(0,10,0.01)

y = odeint(func, init, t, args=(V, var,))

这可能是因为 afleiden 中的方程是用 Sympy 计算的，所以可能是 sympy 表达式？如果是这样，我有什么办法解决这个问题吗？我尝试过使用 lambdify，但没有成功。

Answer 1

正如@Warren Weckesser所说的，以及你所怀疑的，你需要先将表达式转换成可以计算的形式，这样各种偏导数dV/dvar[j]才能返回一个浮点数值。

更一般来说，你的afleiden函数有一个问题，就是它计算了V的解析导数，但并没有计算这些表达式的具体值。我假设v,a,c是你问题中的参数，用来描述一个势能函数V(x,y,z)。我还假设你的常微分方程是

dX/dt = dV/dX(x,y,z)，

其中X=[x,y,z]是你的变量列表。

如果是这样的话，那么你有3个微分方程，而不是你在func()中提到的6个（列表的和是将列表连接起来，而不是求和后的列表）。

import numpy as np
from sympy import lambdify, sympify
from scipy.integrate import odeint

var = ['x', 'y', 'z']
V = sympify("v**2/2*log(1+x**2 + (y/a)**2 + (z/c)**2)")
dVdvar_analytical = [V.diff(var[i]) for i in range(len(var))]
dVdvar = [lambdify(('x', 'y', 'z', 'v', 'a', 'c'), df) for df in dVdvar_analytical]

def afleiden(variables, _, params, dVdvar):
    x, y, z = variables
    v, a, c = params
    return [dVdvarj(x, y, z, v, a, c) for dVdvarj in dVdvar ]

variables0, params = [0.3, 0.2, 0.9], [0.2, 0.6, 0.7]

t = np.arange(0, 10, .1)
y = odeint(afleiden, variables0, t, args=(params, dVdvar))
plot(t, y)

根据你对势能V的表达，原点是一个排斥点，而在模拟中，点y(t)在很长时间后会趋向于无穷大。如果你在表达式前面加上一个负号，原点就变成了一个吸引点，解会收敛到0：

#example with minus sign
V = sympify("-v**2/2*log(1+x**2 + (y/a)**2 + (z/c)**2)")
t = np.arange(0, 100, .1)