odeint: невозможно преобразовать данные массива из dtype ('complex128') в dtype ('float64') в соответствии с правилом 'safe' - PullRequest
Новая
       127

odeint: невозможно преобразовать данные массива из dtype ('complex128') в dtype ('float64') в соответствии с правилом 'safe'

0 голосов
/ 14 апреля 2020

Следующий код выдает ошибку: Невозможно преобразовать данные массива из dtype ('complex128') в dtype ('float64') согласно правилу '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'

Однако, если Я возвращаю реальные значения (вместо комплексных) из func, например: 

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

Тогда odeint работал без ошибок.

Может кто-нибудь помочь мне определить источник / решить эту проблему?

Заранее спасибо!

Ответы [ 2 ]

0 голосов
/ 16 апреля 2020

Спасибо, Лутц! По сути, я пытаюсь решить следующее дифференциальное уравнение (Шредингера): 

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

Принимая преобразование Фурье с обеих сторон, мы получим: 

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

, где u_t - первая производная из вас по времени, к является волновым числом. Ниже приведен мой исправленный / измененный код. 

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)

Это должно предоставить usol форму len (t) x 1

Для каждого значения в k мое окончательное решение должно быть форма len (t) xn, так как k имеет n элементов.

Я также решил эту проблему, используя ode45 в Matlab. Решение, полученное из Python, сильно отличается от того, что я нашел в Matlab. Я знаю, что решение Matlab верное.

0 голосов
/ 14 апреля 2020

Вы изменяете тип данных и возвращаете комплексные значения, когда решатель ODE ожидает реальные значения.

Вы также можете попытаться сделать свой ввод сложным, чтобы по крайней мере эта точка не вызывала противоречий 

y0 = 0+0j

Что именно вы ожидаете в качестве решения? В любой разумной интерпретации вы получите нулевую функцию.

...