0
$\begingroup$

I intended to use Mathematica to perform the following numerical integration arising in the context of cosmology. I wish to evaluate the function $\sigma(R)$ at some specific points and plot those points.

$$\sigma(R) = \frac{1}{2\pi^2R^6}\int_{0}^{\infty}\frac{dk}{k^4}P(k) (-3kR\cos(kR)+3\sin(kR))^2$$

$$P(k)=2\pi^2\delta_H^2\frac{k}{H_0^4}T(k)^2$$

The function $T(k)$ is defined via

$$T(x) = \frac{\log[1+0.171x]}{0.171x}[1+0.284x+(1.18x)^2+(0.399x)^3+(0.490x)^4]^{-0.25}$$

$$x = \frac{k}{k_{eq}}$$

where $\delta_H,H_0,k_{eq}$ are some constant cosmological parameters

Apparently, my naive Mathematica code led to the error message,

NIntegrate::inumr: The integrand has evaluated to non numerical values for all sampling points in the region with boundaries {{∞, 0.}}

Might I ask for some help in correcting my code?

H0 := 0.5*100/(3*10^5) (*Hubble rate in unit of Mpc-1, h=0.5 is used*)
keq := 0.0731*0.5^2 
T[k_] := 
  (Log[1 + 0.171 (k/keq)]/(0.171*(k/keq)))*(1 + 0.284 (k/keq) + (1.18*(k/keq))^2 + 
  (0.399*k/keq)^3 + (0.490 k/keq)^4)^(-0.25)
P[k_] := 2*Pi^2*DeltaH^2*(k*H0^(-4))*T[k]^2
f[R_] := 
  (1/(2*Pi^2*R^6))*
     NIntegrate[(P[k]/k^4)*(3*(-k*R*Cos[k*R] + Sin[k*R]))^2, {k, 0, Infinity}]

I attach a screen capture detailing my code, definition of the variables and the error message.

enter image description here

$\endgroup$

closed as off-topic by Michael E2, m_goldberg, MarcoB, Henrik Schumacher, José Antonio Díaz Navas Apr 30 '18 at 9:27

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Michael E2, m_goldberg, MarcoB, Henrik Schumacher, José Antonio Díaz Navas
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 2
    $\begingroup$ Some tips: 1) please post the code rather than an image. 2) Use = instead of := for constants (e.g., H0). 3) Prevent evaluation by using T[k_?NumericQ]:= instead of T[k_]:= (and the same for P and f). 4) Those squares $\square$ in the error messages probably mean you mistyped something (e.g., Square does not square a number, but represents the geometrical figure instead; use ^2 for squaring). 5) Use, whenever possible, rationalised units (e.g., $G=c=\hbar=1$ or something like that) so as to have numbers of order one. This makes numerical routines much more reliable and stable $\endgroup$ – AccidentalFourierTransform Apr 29 '18 at 22:57
  • $\begingroup$ Thanks for your suggestions! Unfortunately, the issues still remain unsolved.. $\endgroup$ – Heng Fai Chang Apr 29 '18 at 23:22
  • $\begingroup$ Thanks for the comments. Yes, I do mean x^2 by Square [..]. After making the corrections, the issue still remains unsolved. $\endgroup$ – Heng Fai Chang Apr 29 '18 at 23:53
  • $\begingroup$ I get that DeltaH is undefined (that is, it has no numeric value and leads to the error). $\endgroup$ – Michael E2 Apr 30 '18 at 0:09
  • 2
    $\begingroup$ This kind of error is easy to avoid if you just try out your functions before using them. Simply typing P[1] would have revealed your trouble. $\endgroup$ – John Doty Apr 30 '18 at 0:37
1
$\begingroup$

Your pasted code did not include the definition of DeltaH

H0 = 0.5*100/(3*10^5) // 
  Rationalize;(*Hubble rate in unit of Mpc-1, h=0.5 is used*)

keq = 0.0731*0.5^2 // Rationalize;

DeltaH = 19*10^-6;

T[k_] = (Log[1 + 0.171 (k/keq)]/(0.171*(k/keq)))*(1 + 
       0.284 (k/keq) + (1.18*(k/keq))^2 + (0.399*
          k/keq)^3 + (0.490 k/keq)^4)^(-0.25) // Rationalize[#, 0] &;

P[k_] = 2*Pi^2*DeltaH^2*(k*H0^(-4))*T[k]^2;

Since f is defined using a numeric technique, its argument should be restricted to numeric values.

f[R_?NumericQ] := (1/(2*Pi^2*R^6))*
  NIntegrate[(P[k]/k^4)*(3*(-k*R*Cos[k*R] + Sin[k*R]))^2, {k, 0, Infinity}]

f[16]

(* 1.72279 *)

Plot[f[R], {R, 2, 20}]

enter image description here

Or a LogPlot

LogPlot[f[R], {R, .01, 20}]

enter image description here

$\endgroup$
  • $\begingroup$ Let me check. Many thanks for your suggestions! $\endgroup$ – Heng Fai Chang Apr 30 '18 at 0:22

Not the answer you're looking for? Browse other questions tagged or ask your own question.