[SciPy-User] scipy integration question

[Ralf Gommers]

On Wed, Jan 12, 2011 at 5:21 AM, Robert Thompson <rthompsonj at gmail.com>wrote:

>  I'm just learning python today so please forgive me. I come from a C
> background (although I wouldn't call myself good) so I thought migrating
> some programs over to python would be a good start. I'm trying to do an
> integration involving multiple variables and it's giving me completely
> different answers from my C program and I have no idea why.
>
> here's a snippet from my c code that I'm trying to emulate:
> for(j=0; j<kbins; j++)
> {
>   k[j] = pow(10.,log10(kmin)+j*dk);
>   q = k[j]/(Om*h*h);
>   W = 3.*(sin(k[j]*R)-k[j]*R*cos(k[j]*R))/pow(k[j]*R,3.)
>   T = bunch of random stuff, it's a function of q alone though;
>   P[j] = k[j]*T*T;
>   ptemp = P[j];
>   ktemp = k[j];
>   sigma2[j] = qromb(k_integrate,0,kmax);
> }
>
> double k_integrate(double ktemp,double ptemp, double W)
> {
>   integrand_k = ktemp*ktemp*ktemp*ptemp*W*W;
>   return integrand_k;
> }
>
> And here is my corresponding python code:
> from scipy.integrate import quad
> from pylab import *
>
> Odm = 0.26
> Ob  = 0.04
> Om  = Odm+Ob
> Ok  = 0.0
> Ol  = 0.7
> O0  = Ol+Om
> h   = 0.7
>
> rho0 = 1.462e11 #Msun/Mpc^3
> z = 0
> kbins = 10000
> kmin  = 0.000001
> kmax  = 10000
> dk = (log10(kmax)-log10(kmin))/kbins
>
> Ho = h*100*1.02272e-12
> c  = 3.067e-7
>
> def k_integrand(k,W,P):
>     return k**3*W**2*P
>
> M = 14.4963
> R = ((3*10**(M))/(4*pi*rho0))**(1./3.)
>
> powerspec, wavenumber, sigma2 = [],[],[]
>
> for j in range(0,kbins):
>     k = 10**(log10(kmin)+j*dk)
>     wavenumber.append(k)
>
>     q = k/(Om*h**2)
>     W = 3.*(sin(k*R)-k*R*cos(k*R))/(k*R)**3
>     T = log(1.+2.34*q)/(2.34*q) *
> (1.+3.89*q+(16.1*q)**2.+(5.46*q)**3.+(6.71*q)**4.)**(-1./4.)
>     P = k*T**2
>     powerspec.append(P)
>
>     sigmaresult, err = quad(k_integrand,0,kmax,args=(W,P))
>     sigma2.append(sigmaresult)
>
>
> If I stick in k=1 the C program returns log10(sigma2)=16.076 but the python
> script returns log10(sigma2)=6.777. Any ideas on what is causing this huge
> discrepancy?  I checked q,W,T,&P and they are in agreement so I know it's
> something to do with the integration.  Thanks in advance for your time!
>
> I'm not sure what the for loop and the lists in your Python code are for,
but I suspect you want something like this:

from scipy.integrate import quad
from numpy import *

Odm = 0.26
Ob  = 0.04
Om  = Odm + Ob
h   = 0.7
rho0 = 1.462e11 #Msun/Mpc^3

M = 14.4963
R = ((3*10**(M))/(4*pi*rho0))**(1./3.)

def k_integrand(k):
    q = k/(Om*h**2)
    W = 3.*(sin(k*R)-k*R*cos(k*R))/(k*R)**3
    T = log(1.+2.34*q)/(2.34*q) *
(1.+3.89*q+(16.1*q)**2.+(5.46*q)**3.+(6.71*q)**4.)**(-1./4.)
    P = k*T**2
    return k**3*W**2*P

kmax = 1e4
sigmaresult, err = quad(k_integrand, 0, kmax)


This gives a very small number because your k_integrand is <1e-11 over most
of the range. If you replace it with a known function (like integrate f(x) =
x from 0 to kmax) you'll see if you get the right answer or not.

Cheers,
Ralf
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110115/7e7bb0de/attachment.html>