odeint：无法根据规则“安全”将数组数据从 dtype('complex128') 转换为 dtype('float64')

Question

以下代码给出错误：Cannot cast array data from dtype('complex128') to dtype('float64') based on 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'

但是，如果我从 func 返回实值（而不是复数），例如：

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

然后odeint 工作没有任何错误。

谁能帮我确定/解决这个问题的根源吗？

提前致谢！

Answer 1

您正在更改数据类型并返回复数值，而 ODE 求解器需要实数值。

你也可以试着让你的输入变得复杂，这样至少这点不会产生矛盾

y0 = 0+0j

您对解决方案的期望究竟是什么？在任何合理的解释中，您都会得到零函数。

Answer 2

谢谢，路兹！基本上，我正在尝试解决以下（薛定谔）微分方程：

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 w.r.t 的一阶导数。时间，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)

这应该为 usol 提供 len(t) x 1 的形状

对于 k 中的每个值，我的最终解应该是 len(t) x n 的形状，因为 k 有 n 个元素。

我还在 Matlab 中使用 ode45 解决了这个问题。从 Python 获得的解决方案与我在 Matlab 中找到的解决方案大不相同。我知道 Matlab 解决方案是正确的。