SciPy 无法将复数转换为浮点数

Question

我遇到了一个问题，我在下面对此进行了简化。

所以，我有一个函数（矩阵），它返回矩阵的特征值。

第二个函数 (deriv) 求特征值关于 delta 的导数。 SciPy自己的导函数很慢，所以我用了复差分法。

然后我对 r 和 theta 取被积函数的二重积分。在对 theta 进行一维积分之后，我求解了 r 的微分方程，而不是采用双积分，这使得计算速度更快。

最后，积分作品。但是，如果我尝试找到满足等式的特定 delta，则会出现错误：无法将复数转换为浮点数，我不明白虚数的来源。

我是 Python 世界的新手，任何帮助将不胜感激。谢谢。

from scipy.integrate import quad
from numpy import linalg as LA
import numpy as np
from pylab import *

from scipy.integrate import odeint
from scipy.integrate import ode
from scipy.optimize import fsolve


#This routine is for the matrix and returns the eigenvalues
def matrix(r, theta, delta1, delta2):

    mat = np.array([[r**2-1,r*np.cos(theta),0,delta1],r*np.cos(theta),r**2 - 1, -delta1,0],[0,-delta2,-r**2+1,-r*np.cos(theta)],[delta2,0,-r*np.cos(theta),-r**2+1]])
    return np.sort(LA.eigvals(mat))[2:4]

#This routine takes the derivatives of eigenvalues with respect to the parameter delta. I set delta1 = delta2.
def deriv(r, theta, delta):

    h = 0.00000000001

    return np.imag(matrix(r, theta, delta, delta+h*1j))/h

#This is the integrand that we need to integrate over r and theta
def integrand(theta,r, beta, delta):
    ee = matrix(r, theta, delta, delta)
    dd = deriv(r, theta, delta)
    return (np.tanh(0.5*beta*ee[0])*dd[0]+np.tanh(0.5*beta*ee[1])*dd[1])*r*r*np.sin(theta)

#Just integrate over theta
def polarint(y,r,beta,delta):
    return quad(integrand,0,np.pi,args = (r,beta,delta))[0]

#Instead of integrating directly over r, solve the differential equation.
#Lambda is the range of r.
def givesclen(delta, beta, lam):
    y0 = 0
    t_out = np.linspace(0, lam, 2500);
    rt = odeint(polarint,y0,t_out,args=(beta,delta))
    temp = (rt[-1]/delta)/(4*np.pi**2)-t_out[-1]/(2*np.pi**2)
    return temp[0]

#The goal is to find the delta; given sl, lam, beta
#Such that the result of the integration is equal to sl
def equationgap(delta, beta, lam,sl):
    return givesclen(delta, beta, lam)*4*np.pi - sl

#Test if the equationgap works, the result should be close to zero!
print equationgap(.5,40,500,.1744)

#Now use the fsolve function should find the delta to be .5!
#beta = 40
#lam = 500
#sl = 0.174
#fsolve(equationgap,.6,args = (beta, lam, sl))

编辑：

错误信息是：

Traceback (most recent call last):
File "test.py", line 38, in polarint
return quad(integrand,0,np.pi,args = (r,beta,delta))[0]
File "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 281, in quad
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
File "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 345, in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
File "test.py", line 30, in integrand
dd = deriv(r, theta, delta)
File "test.py", line 22, in deriv
return np.imag(matrix(r, theta, delta, delta+h*1j))/h
File "test.py", line 14, in matrix
mat = np.array([[r**2-1,r*np.cos(theta),0,delta1],[r*np.cos(theta),r**2 - 1, -delta1,0],[0,-delta2,-r**2+1,-r*np.cos(theta)],[delta2,0,-r*np.cos(theta),-r**2+1]])
TypeError: can't convert complex to float

Answer 1

我可以通过以下方式生成您的错误消息：

x = np.ones((3,))
x[1] = 1 + 1j

x 是一个 dtype=float 数组。尝试插入复杂值会产生错误。

如果在评估中产生错误：

mat = np.array([[r**2-1,r*np.cos(theta),0,delta1],
                [r*np.cos(theta),r**2 - 1, -delta1,0],
                [0,-delta2,-r**2+1,-r*np.cos(theta)],
                [delta2,0,-r*np.cos(theta),-r**2+1]])

问题是 - r、theta、delta1、delta2 的哪些值会产生此错误。 1st 3 是实标量，最后一个复数，我得到一个 (4,4) 复数数组。我怀疑其中一个变量是我们没有说明的。

---

我可以将pvs 的最小示例缩短为：

np.array([1, np.array([1j])])   # error
np.array([[1],np.array([1j])])  # ok

切换术语的顺序，结果不一样：

np.array([np.array([1j]),1])
# array([array([ 0.+1.j]), 1], dtype=object)

显然，错误是尝试从标量和数组的混合构建更大的数组的结果，其中一个是复杂的。错误消息不明确，可能是意外问题的结果。它可能值得一个错误报告。

解决方案是将所有0s 转换为列表[0]，或者确保delta2 值是标量，而不是数组。

Answer 2

这是失败的最小代码

from numpy import array
q = [[-1.0, 0.0, 0, array([ 0.6])], [0.0, -1.0, array([-0.6]), 0], [0, array([-0.6 -1.00000000e-11j]), 1.0, -0.0], [array([ 0.6 +1.00000000e-11j]), 0, -0.0, 1.0]]
array(q)

这有点奇怪。但是，请注意，其中有一个元素与普通数字混合的数组。这可以通过在给array()之前打印出你给它的对象来发现。

要修复它，请在您提供给 fsolve 的函数中将 delta 更改为 delta[0]：

def equationgap(delta, beta, lam,sl):
    return givesclen(delta[0], beta, lam)*4*np.pi - sl

因为其余代码期望 delta 是单个数字，而不是数组。 fsolve 将为要优化的函数提供一个数字数组，即使只有一个数字。