python - Scipy.integrate浮点错误

Question

我正在尝试整合这个由几个函数组成的令人讨厌的积分。

import matplotlib.pyplot as plt
import numpy as np
import os
import scipy.integrate as integrate

path = '/users/username/Desktop/untitled folder/python files/document/'

os.chdir( path )
data = np.load('msii_phasespace.npy',mmap_mode='r')
# data.size: 167197
# data.shape: (167197,)
# data.dtype: dtype([('x', '<f4'), ('y', '<f4'), ('z', '<f4'),
  # ('velx', '<f4'), ('vely', '<f4'), ('velz', '<f4'), ('m200', '<f4')])

## Constant
rho_m = 3.3e-14
M = data['x'] Mass of dark matter haloes
R = ((3*M)/(rho_m*4*(3.14)))**(1.0/3.0) # Km // Radius of sphere 
k = 0.001 # Mpc h^-1 // Wave Dispersion relation 
delt_c = 1.686 
h = 0.73 # km s^-1 Mpc^-1 
e = 2.718281 # Eulers number
T_CMB = 2.725   
Omega_m = 0.27 
kR = k*R


def T(k): 
    q = k/((Omega_m)*h**2)*((T_CMB)/27)**2
    L = np.log(e+1.84*q)
    C = 14.4 + 325/(1+60.5*q**1.11)
    return L/(L+C*q**2)

def P(k): 
    A = 0.75 
    n = 0.95 
    return A*k**n*T(k)**2

def W(kR): # Fourier Transfrom in the top hat function
    return 3*(np.sin(k*R)-k*R*np.cos(k*R))/(k*R)**3

def sig(R): 
    def integrand(k,P,W):
        return k**2*P(k)*W(kR)**2
I1 = integrate.quad(integrand, lambda k: 0.0, lambda k: np.Inf, args=(k,))
return ((1.0/(2.0*np.pi**2)) * I1)**0.5

打印出 sig(R) 会给出 TypeError: a float is require。

我需要注意的是，R 是物体的半径，与其质量成正比。质量由结构化数组给出。我不确定我是否使用正确的积分命令来评估所有质量。只有数组有一组值。

如果您需要更多信息，请告知，我很乐意分享我能分享的一切。任何建议将不胜感激。

Answer 1

<强>1。使用 integrand.quad()

计算积分

阅读documentation

• integrand.quad 返回一个长度为 2 的元组，其中第一个值保存积分的估计值。确保您只获取第一个元素而不是整个元组

• 您不需要 lambda 来通过限制。 lambda 用于代替被积函数而不是限制。

• args 参数用于将其他参数传递给被积函数，因此您无需传递 k，因为在此参数上进行积分。您需要传递 w，正如我接下来将展示的那样。

<强>2。类型错误:需要 float

• 您需要函数 W(kR) 返回一个浮点值，积分才能正常工作。让我们循环遍历每个值并计算积分。
• 请记住，您只需要传递函数将使用的参数。 sig(R) 不使用 R 并且 W(kR) 不使用 kR。
• 只需计算并存储一次傅里叶变换，因为它不依赖于积分。无需每次都调用W()。

• 为了解决类型错误，我们将循环遍历傅里叶变换的数组并计算积分。

  fourier_transform = 3*(np.sin(k*R)-k*R*np.cos(k*R))/(k*R)**3
  def sig():
      I_one = np.zeros(len(fourier_transform))
      def integrand(k, w):
          return k**2*p_func(k)*w
      for i, w in enumerate(fourier_transform):
          I_one[i], err = integrate.quad(integrand, 0.0, np.inf, args = w)
          print I_one[i]
      return ((1.0/(2.0*np.pi**2)) * I_one)**0.5
  

<强>3。变量和函数的命名约定

查看答案 here 。尝试使用有意义的变量和函数名称，以尽量减少错误和困惑。所有局部变量和函数都应该小写。作为常量的全局变量应该是大写的。这是我的尝试，但是你会更好地了解这些变量和函数的含义。

#global constants should be capitalized
RHO_M = 3.3e-14
H = 0.73 # km s^-1 Mpc^-1 
EULERS_NUM = 2.718281 # Eulers number
T_CMB = 2.725   
OMEGA_M = 0.27 

#nonconstant variables should be lower case 
mass = data['x'] #Mass of dark matter haloes
radius = ((3*mass)/(RHO_M*4*(3.14)))**(1.0/3.0) # Km // Radius of sphere 
#not sure if this is a constant???
mpc_h = 0.001 # Mpc h^-1 // Wave Dispersion relation 
mpc_radius = mpc_h*radius

#this stays constant for all integrations so let's just compute it once
fourier_transform = 3*(np.sin(mpc_h*radius)*mpc_h*radius*np.cos(mpc_h*radius))/(mpc_h*radius)**3

def t_calc(k): 
    q = k/((OMEGA_M)*H**2)*((T_CMB)/27)**2
#I didn't change this because lower case L creates confusion
    L = np.log(EULERS_NUM+1.84*q)
    c = 14.4 + 325/(1+60.5*q**1.11)
    return L/(L+c*q**2)

def p_calc(k): 
    a = 0.75 
    n = 0.95 
    return a*k**n*t_calc(k)**2

def sig():
    I_one = np.zeros(len(fourier_transform))
    def integrand(k, w):
        return k**2*p_func(k)*w
    for i, w in enumerate(fourier_transform):
        I_one[i], err = integrate.quad(integrand, 0.0, np.inf, args = w)
        print I_one[i]
    return ((1.0/(2.0*np.pi**2)) * I_one)**0.5

<强>5。 def 被积函数(k, P)

虽然 integrand() 是嵌套的，但它仍然可以在全局命名空间中搜索函数 P，因此不要传递它。