我正试图用 scipy.integrate.odeint
求解一组耦合微分方程 . 但是,当我尝试运行该程序时,我收到以下错误:
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe' odepack.error: Result from function call is not a proper array of floats.
这是我使用的代码:
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计算的,因此可能是 sypmpy
表达式?如果是这样,我能做些什么吗?我尝试使用lambdify,但没有成功 .
1 回答
正如@Warren Weckesser所说,并且正如您所怀疑的那样,您需要首先对表达式进行lambdify,以便各种偏导数
dV/dvar[j]
返回浮点值 .更一般地说,你的
afleiden
函数的一个问题是它评估V
的分析导数,而不计算这些表达式的值 . 我假设v,a,c
是你的问题的参数,deacribing潜在的功能V(x,y,z)
. 我也假设你的o.d.e.是dX/dt = dV/dX(x,y,z)
,其中
X=[x,y,z]
是您的变量列表 .如果是这种情况,那么你有3个差异 . 方程式,而不是
6
,如func()
(列表的总和是列表的串联,而不是总和列表) .与你的潜力
V
的表达式一致,原点是一个驱逐者,并且点y(t)
在模拟过程中趋于无穷大 . 如果在表达式的开头添加减号,则原点变为吸引子,解决方案收敛于0
: