3

I'm trying to solve a set of coupled differential equations using scipy.integrate.odeint. However when I try to run the program I get the following error:

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.

Here is the code I use:

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,))

Could it be because the equations from afleiden are calculated using Sympy and thus probably sypmpy expressions? If so, is there anything i can do about it? I tried using lambdify but that didn't work out.

4

1 回答 1

2

正如@Warren Weckesser 所说，正如您所怀疑的那样，您需要首先对表达式进行lambdify，以便各种偏导数dV/dvar[j]返回一个浮点值。

更一般地说，您的afleiden函数的一个问题是它评估 的解析导数V，而不计算这些表达式的值。我假设这v,a,c是你的问题的参数，描述了一个潜在的功能V(x,y,z)。我还假设你的颂歌是

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

X=[x,y,z]你的变量列表在哪里。

如果是这种情况，那么您有 3 个差异。方程式，而不是6您的func()（列表的总和是列表的串联，而不是总和的列表）。

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)

在此处输入图像描述

于 2014-04-30T21:59:08.493 回答