0

'numpy.float64' object cannot be interpreted as an integer运行以下代码时遇到此错误消息：

from scipy.integrate import nquad
from numpy import sqrt, sin, cos, pi, arcsin, maximum, log, log10, exp
from matplotlib.pyplot import xscale, yscale, plot, show, savefig

#define constants
msun=1500
m1=30*msun
m2=30*msun
M=m1+m2
eta=m1*m2/M**2
mSMBH=4.2*10**6*msun
v=5/10**4
pc=3.09*10**16
r0=1.005*10**21              #nfw parameter#
rho0=1.27/10**49             #nfw parameter#
pmtpy=3*10**8*365*24*3600    #rate from per-meter to per-yr
rhoc=1.35*10**(-53)
cm=2.16

#define functions
def g(x):
    return log(1+x)-x/(1+x)
def y(x):
    return 1.43*10**12*1500/x
def sig(x):
    return 13.216*y(x)**0.41/(1+1.102*y(x)**0.2+6.22*y(x)**0.333)
def c(x):
    return 2.881*(1+(sig(x)/1.257)**1.022)*exp(0.06/sig(x)**2)
def r0(x):
    return (3*x/(800*c(x)**3*rhoc*pi))**(1/3)
def rho0(x):
    return 200*rhoc*c(x)**3/3/g(c(x))
def sigma(x):
    return sqrt(4*pi*rho0(x)*r0(x)**2*g(cm)/cm)
def mSMBH(x):
    return 1500*10**(8.15+4.38*log10(1500*sigma(x)))
def rcap(x):
    return 8*mSMBH(x)

def Rsp(x):
    return 0.122*sqrt(mSMBH(x)/rho0(x)/r0(x))
def rhosp(r,x):
    return 3*rho0(x)*r0(x)*Rsp(x)**(4/3)*(1-rcap(x)/r)**3/r**(7/3)

#x,y,z are to be integrated, and m is the x-axis parameter to be searched through using while loop
def v(z,m):
    return sqrt(mSMBH(m)/z)
def vrel(y,z,m):
    return 2*v(z,m)*sin(y/2)
def vcm(y,z,m):
    return 2*v(z,m)*cos(y/2)
def bmin(y,z,m):
    return 2*M/vrel(y,z,m)
def bmax(y,z,m):
    return maximum(bmin(y,z,m), M*(340*eta/(3*cos(1)*vrel(y,z,m)**9))**(1/7))
def a(x,y,z,m):
    if bmax(y,z,m)>bmin(y,z,m):
        return M/(vrel(y,z,m)**2*((bmax(y,z,m)/x)**7-1))
    else:
        return 0
def rt(x,y,z,m):
    return a(x,y,z,m)*(mSMBH(m)/M)**(1/3)
def rd(x,y,z,m):
    return maximum(rcap(m),rt(x,y,z,m))
def rz(x,y,z,m):
    return rd(x,y,z,m)*(sqrt(1+8*mSMBH(m)/(rd(x,y,z,m)*vcm(y,z,m)**2))-1)/2
def thetal(x,y,z,m):
    if rd(x,y,z,m)<z and rz(x,y,z,m)<z:
        return arcsin(rd(x,y,z,m)*sqrt(2*(vcm(y,z,m)**2/2-mSMBH(m)/z+mSMBH(m)/rd(x,y,z,m)))/(vcm(y,z,m)*z))
    else:
        return pi/2

#prepare looping tools    
m=10**9*msun
M=[]           #x value array
Rate=[]        #result estimation array
RateUp=[]      #result upper bound array
RateDown=[]    #result lower bound array
options={'limit':200}

while m<10**14*msun:
    def zbound():
        return [rcap(m),Rsp(m)]
    def ybound(z_foo):
        return [10**(-10),pi-10**(-10)]
    def xbound(y,z):
        return [bmin(y,z,m)-10**(-6),bmax(y,z,m)-10**(-6)]
    def f(x,y,z):
        return pmtpy*8*pi**2*v(z,m)/m1**2*rhosp(z,m)**2*z**2*sin(y/2)**2*cos(y/2)*x*cos(thetal(x,y,z,m))
    ans=nquad(f, [xbound,ybound,zbound],opts=[options,options,options])   
    M.append(m)
    Rate.append(ans[0])
    RateUp.append(ans[0]+ans[1])
    RateDown.append(ans[0]-ans[1])
    m+=10**15

xscale('log')
yscale('log')
plot(M,Rate,M,RateUp,M,RateDown)
show()

基本上，我正在尝试评估取决于参数 m 的三重积分。该错误具体发生在

    File "C:/Users/user/.spyder-py3/temp.py", line 83, in xbound
        return [bmin(y,z,m)-10**(-6),bmax(y,z,m)-10**(-6)]
    File "C:/Users/user/.spyder-py3/temp.py", line 52, in bmax
        return maximum(bmin(y,z,m), M*(340*eta/(3*cos(1)*vrel(y,z,m)**9))**(1/7))

TypeError: 'numpy.float64' object cannot be interpreted as an integer

但是我真的不知道为什么在我的代码中的“bmax”中的任何地方都需要一个整数。知道可能是什么问题吗？非常感谢您的帮助！

4

1 回答 1

0

问题的出现是因为你有M=[]，这使它成为一个list. 这完全覆盖了M = m1 + m2您在顶部所做的。然后在第 55 行，你有M * BLAH，其中M不是一个 numpy 数组（它是一个列表），因此该*运算符被解释为“重复列表 BLAH 次”。这就是问题所在。

PS：更改它以np.array(M) * (340*eta/(3*cos(1.)*vrel(y,z,m)**9))**(1/7)解决此问题，但会产生其他问题。例如，如果您使用的是 python2，那么(1/7) = 0，您可能想要(1./7.)，等等。np.arange您预先了解群众 - 因此使用创建完整M数组然后处理它可能会更快。

PS：由于您有与天文学相关的计算，请查看astropy软件包，它可以完成您正在寻找的一些事情。

于 2018-05-30T10:56:38.963 回答