首页 文章

使用scipy.integrate.odeint时出现TypeError

提问于
浏览
3

我正试图用 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 回答

  • 2

    正如@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() (列表的总和是列表的串联,而不是总和列表) .

    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)
    

    enter image description here

相关问题