TypeError when using scipy.integrate.odeint

Question

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.

Answer 1

as @Warren Weckesser says and as you suspected, you need to lambdify the expressions first, so that the various partial derivatives dV/dvar[j] return a floating point value.

More in general, one problem of your afleiden function is that it evaluates the analytical derivatives of V, without calculating the value of those expressions. I assume that v,a,c are parameters of your problems, deacribing a potential function V(x,y,z). I also assume that your o.d.e. is

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

where X=[x,y,z] is your list of variables.

If this is the case, then you have 3 diff. equations, and not 6 as in your func() (a sum of lists is the concatenation of the lists, not the list of the sums).

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)

Consistently with your expression for the potential V, the origin is a repellor and the point y(t) tends to infinity in the simulation for large times. if you add a minus sign at the begin of the expression, the origin becomes an attractor and the solution converges to 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)