Scipy.integrate float error

Question

I'm trying to integrate this nasty integral that consists of couple of functions.

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

printing out the sig(R) gives me TypeError: a float is require.

I need to note that R is a radius of a object that is directly proportional to its mass. The masses are given by an structured array. I'm not sure if I am using the correct integral command to evaluate all the masses. Only the array has a set of values.

If you need anymore information, please tell, and i'll be happy to share whatever I can. Any advice will be appreciated.

Answer 1

1. Computing an integral with integrand.quad()

Read the documentation

• integrand.quad returns a length 2 tuple with the first value holding the estimated value of the integral. Make sure you just get the first element and not the whole tuple

• You don't need lambda, to pass the limits. lambda is used in place of an integrand function not for the limits.

• The args parameter is used for passing additional parameters to your integrand so you don't need to pass k since integrating over this parameter. You will need to pass w, as I will show next.

2. TypeError: a float is required

• You need function W(kR) to return a float for the integration to work. Let's loop through each value and calculate the integral.
• Remember you only need to pass arguments which the function will use. sig(R) does not use R and W(kR) does not use kR.
• Just compute and store the fourier transform once since it does not depend on the integral. There is no need to call W()each time.

• To solve the type error we will loop through the array from the fourier transform and compute the integral.

  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. Naming convention for variables and functions

Look at the answers here. Try to use meaningful variable and function names to minimize bugs and confusion. All local variables and functions should be lower case. Global variables that are constants should be upper case. Here is my attempt, but you will know the meaning of these variables and functions better.

#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 integrand(k, P)

Although integrand() is nested it can still search for function P in the global namespace so don't pass it.