CAMB python convergence power spectrum code

Click For Summary
SUMMARY

The discussion centers on the Python CAMB code for computing the convergence power spectrum, specifically the integration method used in the code. The integration is performed using a dot product of two arrays instead of the traditional scipy.integrate.quad function. The first array, Delta-Chi, and the second, the integrand, are combined in a for-loop to achieve numerical integration. The results show that this method yields accurate results comparable to those obtained from scipy's quad function, particularly with N set to 1000 or higher.

PREREQUISITES
  • Familiarity with Python programming
  • Understanding of numerical integration techniques
  • Knowledge of the CAMB code and its applications in cosmology
  • Experience with NumPy and SciPy libraries
NEXT STEPS
  • Explore the CAMB documentation for advanced features and examples
  • Learn about numerical integration methods in Python, focusing on the trapezoidal rule and Simpson's rule
  • Investigate the use of NumPy's dot product for numerical computations
  • Study the implications of convergence power spectrum in cosmological simulations
USEFUL FOR

Researchers in cosmology, data scientists working with numerical simulations, and Python developers interested in advanced numerical methods.

sunrah
Messages
191
Reaction score
22
I'm trying to understand this python CAMB code: http://camb.readthedocs.io/en/latest/CAMBdemo.html
Scroll down to In[29] and In[30] to see it.

It's an integration over chi (comoving distance), yet scipy.integrate.quad is not called. It seems that the fun stuff happens in the last for-loop in In[30], where they use numpy.dot to take dot product of two arrays. What's going on here? The first array is Delta-Chi (see In[29]) and the second is simply the integrand. So how is this integration? Thanks
 
Technology news on Phys.org
I'v been messing around with this and found something interesting:

Python:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

xmax = 10
N = 1000
xx = np.linspace(0,xmax,N)

dx = (xx[2:] - xx[:-2])/2
xx = xx[1:-1]
Nx = len(xx)
yy1 = np.arange(0,Nx, dtype=np.float64)
yy2 = []

f = lambda z: z**3

for ix, x in enumerate(xx):
    foo = np.linspace(min(xx),x,N-2)**3
    yy1[ix] = np.dot(dx, foo)
    yy2.append(quad(f,0,x)[0])

plt.plot(xx,yy1,label='dot')
plt.plot(xx,yy2,label='quad')

If you run this script you'll find pretty good agreement between the integration via quad and integration using the dot product of dx and integrand. The value of N controls accuracy, with N>1000 not producing much noticeable difference. But why does this work? What is this method called thanks
 

Similar threads

  • · Replies 6 ·
Replies
6
Views
6K
  • · Replies 6 ·
Replies
6
Views
4K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 23 ·
Replies
23
Views
4K
Replies
5
Views
47K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 2 ·
Replies
2
Views
4K
  • · Replies 6 ·
Replies
6
Views
5K