# How to calculate and plot this integral? [closed]

Any help to integrate this equation:

$$\begin{eqnarray} p(k,q)&=&-16 H^2 q^2\int_{q_a}^q \int_{q_a}^q \frac{1}{\sqrt{H^2-k_1{^2}}}\sinh(\sqrt{H^2-k_1{^2}}(q-b))\frac{1}{H b} \Big[m^2 y(b)^2 \nonumber\\ &-&\frac{1}{H b} y^{\prime 2}(b)\Big]+\frac{1}{\sqrt{H^2-k_2{^2}}}\sinh(\sqrt{H^2-k_2{^2}}(q-f)) \frac{1}{H b} \Big[m^2 y^2(f) \nonumber\\ &-&\frac{1}{H f} Y^{\prime 2} (f) \Big] d b ~d f \end{eqnarray}$$

where:

$$y(b) = c_1 b^{\frac{1}{2}+\sqrt{1+\frac{m^2}{H^2}}}+ c_2 b^{\frac{1}{2}-\sqrt{1+\frac{m^2}{H^2}}},$$

$$y(f) = c_1 f^{\frac{1}{2}+\sqrt{1+\frac{m^2}{H^2}}}+ c_2 f^{\frac{1}{2}-\sqrt{1+\frac{m^2}{H^2}}}.$$

Here is my naive trial:

I let $$c_1=c_2=1$$. Also let $$k_1=k_2=k$$. Let $$m=H=100$$ GeV. Also take $$q_a=-100$$.

y[b_] := c1*b^(1/2 + Sqrt[1 + m^2/H^2]) + c2*b^(1/2 - Sqrt[1 + m^2/H^2]);

y[f_] := c1*f^(1/2 + Sqrt[1 + m^2/H^2]) + c2*f^(1/2 - Sqrt[1 + m^2/H^2]);

eq[k_,b_,f_,q_] :=
-16 H^2 q^2 ( (1/Sqrt[H^2 - k1^2])*Sinh[Sqrt[H^2 - k1^2]*(q - b)]*(1/(H*b))*
(m^2*y[b]^2 - (1/(H*b))*D[y[b], b]^2) + (1/Sqrt[H^2 - k2^2])*
Sinh[Sqrt[H^2 - k2^2]*(q - f)]*(1/(H*b))*(m^2*y[f]^2 -  (1/(H*f))*D[y[f], f]^2) )


Then the integration should be taken in the limit where $$q\to 0$$

i1[k_?NumericQ, f_?NumericQ, q_?NumericQ] := NIntegrate[eq[k, b, f, q], {b, -100, q}]

i2[k_?NumericQ, q_?NumericQ] :=
NIntegrate[i1[k, f, q], {f, -100, q}]

eqq[k_] = Limit[i2[k, q], q -> 0]

Plot[eqq[k], {k, 0, 1}]


But this code doesn't work!

In the equation's reference 2309.11694 $$H=H^\star=10^{14}=m=-q_a$$ GeV, and the plot of $$p(k,q)$$ as a function of $$k/H$$ is given by:

• In your code functions y1,y2 are not defined. Oct 22, 2023 at 14:04
• @AlexTrounev. You mean D[y[b], b] ? Ok, edit it. There are only y[b] and y[f]. Oct 22, 2023 at 14:29
• Upper-case Y in expression is not defined. Maybe a typo.
– josh
Oct 22, 2023 at 15:06
• Additional context from Physics.SE: Calculating the power spectrum of the gravitational waves Oct 22, 2023 at 15:54
• You're calling NIntegrate on an integral with symbolic parameters. You need to specify your parameter values c1 = c2 = 1; m = H = 100; , and set k1 and k2 to k. I am also not sure if Limit will work with inexact input
– ydd
Oct 22, 2023 at 16:28