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

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

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

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


python scipy differential-equations odeint ode45
person Shoaibur Rahman    schedule 13.04.2020    source источник
comment
В конце концов, этот вопрос касается той же PDE Шредингера, что и в stackoverflow.com/questions/29803342/, только вместо использования RK4 здесь используется scipy решатель. Функции производных в обоих вариантах остаются неизменными, можно подумать о том, чтобы сделать вариант разделения оператора (vhat) более локальным, чтобы избежать фазовых ошибок при уменьшении больших чисел по модулю 2*pi.   -  person Lutz Lehmann    schedule 17.04.2020

Ответы (2)


arrow_upward
0
arrow_downward

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

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

y0 = 0+0j

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

person Lutz Lehmann    schedule 13.04.2020

arrow_upward
0
arrow_downward

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

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 - первая производная от u относительно время, k - волновое число. Ниже мой исправленный / измененный код.

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)

Это должно дать усолу форму len (t) x 1.

Ожидается, что для каждого значения в k мое окончательное решение будет иметь форму len (t) x n, поскольку k имеет n элементов.

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

person Shoaibur Rahman    schedule 16.04.2020
comment
Взгляните на stackoverflow.com/questions/ 29803342 / и stackoverflow.com/questions/29617089/ для обработки аналогичного уравнения. Вы уверены в форме уравнения, вы пропустили множитель i в первом уравнении ?? - person Lutz Lehmann; 17.04.2020
comment
Какую версию Python вы используете? Разница в том, что когда n/2 является целочисленным делением, а в более новой версии 3 - делением с плавающей запятой. - person Lutz Lehmann; 17.04.2020
comment
@LutzLehmann Большое спасибо! Да, я пропустил множитель i во втором члене первого уравнения. Связанный код полезен. - person Shoaibur Rahman; 17.04.2020
comment
Я использую Python 3. - person Shoaibur Rahman; 17.04.2020
comment
y0, t0 = ut[0], t[0] не имеет смысла, вы хотите, чтобы y0 был полным вектором состояния ut. Вы не инициализируете usol = [ut] в этом коде. Цикл времени должен быть for i in t[1:]: usol.append( utsol.integrate(i) ). - person Lutz Lehmann; 17.04.2020
comment
Спасибо! Это сработало! Как я могу принять ваш ответ? Я здесь новенький, поэтому могу пропустить этот вариант. - person Shoaibur Rahman; 18.04.2020