odeint: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

0
votes

The following code gives the error: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

import numpy as np
from numpy.fft import fft
from scipy.integrate import odeint

t = np.linspace(0,9,10)

def func(y, t):
    k = 0
    dydt = fft(y)
    return dydt

y0 = 0
y = odeint(func, y0, t)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-4885da912033> in <module>
     10 
     11 y0 = 0
---> 12 y = odeint(func, y0, t)

~\AppData\Local\Continuum\anaconda3\envs\udacityDL\lib\site-packages\scipy\integrate\odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    243                              full_output, rtol, atol, tcrit, h0, hmax, hmin,
    244                              ixpr, mxstep, mxhnil, mxordn, mxords,
--> 245                              int(bool(tfirst)))
    246     if output[-1] < 0:
    247         warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

TypeError: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

However, if I return real values (instead of complex) from func like:

def func(y, t):
    k = 0
    dydt = fft(y)
    return np.abs(dydt)

Then the odeint worked without any error.

Could anyone please help me identifying the source of/solving this problem?

Thanks in advance!

2
pythonscipydifferential-equationsodeintode45
In the end, this question is about the same Schrödinger PDE as in stackoverflow.com/questions/29803342/…, only instead of the use of RK4 there here a scipy solver is used. The derivatives functions in both variants stay the same, one might think about making the operator splitting variant (vhat) more local to avoid the phase errors when reducing large numbers modulo 2*pi. - Lutz Lehmann

2 Answers

0
votes

You are changing the data type and return complex values, where the ODE solver expects real values.

You could also try to make your input complex, so that at least this point does not produce contradictions

y0 = 0+0j

What exactly do you expect as solution? In any reasonable interpretation you will get the zero function.

0
votes

Thanks, Lutz! Basically, I am trying to solve the following (Schrodinger) differential equation:

u_t = i/2 u_xx + |u|^2 u

Taking Fourier transform on both sides, we get:

fft(u_t) = -i/2 k^2 fft(u) + i fft(|u|^2 u)

where, u_t is the first derivative of u w.r.t. time, k is the wave number. Below is my corrected/modified code.

import numpy as np
from numpy.fft import fft
from scipy.integrate import odeint

def f(t, ut, k):
    u = ifft(ut)
    return -(1j/2) * k**2 * ut + 1j * fft( np.abs(u)**2 * u )



# Simulation
L = 30
n = 512
x2 = np.linspace(-L/2, L/2, n+1)
x = x2[0:n]

t = np.linspace(0, 2*np.pi, 41)
dt = t[1]-t[0]

k = 2*np.pi/L * np.concatenate( [np.arange(0, n/2), np.arange(-n/2,0,1)] )
print(x.shape, t.shape, k.shape)

# Initial solutions
u = 1 / np.cosh(x)
ut = fft(u)

# Solution for first value of k
y0, t0 = ut[0], t[0]
utsol = ode(f).set_integrator('zvode', method='bdf')
utsol.set_initial_value(y0, t0).set_f_params(k[0])
for i in t:
    usol.append( utsol.integrate(i+dt) ) # no idea if using i+dt is correct!
usol = np.array(usol)

This should provide usol with a shape of len(t) x 1

For each value in k, my final solution is expected to be a shape of len(t) x n, since k has n elements.

I also solved this using ode45 in Matlab. The solution obtained from Python is vastly different than what I found in Matlab. I know Matlab solution is correct.