Numerical integration of a special integral

I want to numerically compute an integral of the following form $$\int_0^{\infty}\left[\int^1_0J_0(v\rho)\cos(\rho^2/2)\rho\,\mathrm{d}\rho \right]^2\,v\,\mathrm{d}v$$

where $$J_0(x)$$ is the zero-th order Bessel function, BesselJ[0,x] in Mathematica. As there is a square of integral in variable that needs to be done numerically before integrating another variable numerically, the function $$\textbf{NIntegrate}$$ seems do not work, anyone know how to deal with integrals of this form numerically using functions of Mathematica?

• Use iteratively NIntegrate[y*(Integrate[x*f[x,y],{x,0,1}])^2,{y,0,1}]. Integrate does numerical integration too, just first does analytical... Commented Nov 17, 2017 at 9:48
• @JoséAntonioDíazNavas Thank you, but this seems does not work. The description of the cause in Mathematica is similar as when you numerically integrate Cos(c x) for x in (0,1) and c is a constant. Commented Nov 17, 2017 at 10:14
• It could be helpful if you specify what is f[x,y]definition Commented Nov 17, 2017 at 10:17
• @JoséAntonioDíazNavas I have edited the question and $f$ is given in it. Commented Nov 17, 2017 at 10:24
• I have faced several times with this kind of integrals. Unfortunately, MMA does not help to solve that inside, I had to go over analytical techniques to solve it manually. Maybe someone can provide a solution I am not able to see... Commented Nov 17, 2017 at 10:31

Here's a way to compute the result:

(* pre-process integrand for Levin Rule - saves ~10% of time *)
li = NIntegrateLevinIntegrandReduce[BesselJ[0, v   r]  Cos[r^2/2] r, r];
int[v0_?NumericQ] := Block[{v = v0},
NIntegrate["ignored", {r, 0, 1},
"Amplitude" -> li@"Amplitude", "Kernel" -> li@"Kernel",
"DifferentialMatrix" -> First@li@"DifferentialMatrices",
Method -> {"GaussKronrodRule", "Points" -> 11}},
PrecisionGoal -> 8,
MaxRecursion -> 50]
];

NIntegrate[int[v]^2 v, {v, 0, Infinity},
MaxRecursion -> 20, PrecisionGoal -> 6,
Method -> {"DoubleExponential"}] // AbsoluteTiming
(* {307.02, 0.460368} *)


From Mariusz Iwaniuk's comment:

0.25 + Sin[1]/4
(* 0.460368 *)

• It seems that there is a bug with Sum. Using NSum I get the same result. Commented May 26 at 5:49

The problem you asked for can be converted into a single multi-dimensional integral. Considering that $$\left(\int_0^1 \rho \cos \left(\frac{\rho ^2}{2}\right) J_0(v \rho ) \, d\rho \right){}^2=\\\left(\int_0^1 \rho _1 \cos \left(\frac{\rho _1^2}{2}\right) J_0\left(v \rho _1\right) \, d\rho _1\right) \int_0^1 \rho _2 \cos \left(\frac{\rho _2^2}{2}\right) J_0\left(v \rho _2\right) \, d\rho _2=\\\int _0^1\int _0^1\rho _1 \rho _2 \cos \left(\frac{\rho _1^2}{2}\right) \cos \left(\frac{\rho _2^2}{2}\right) J_0\left(v \rho _1\right) J_0\left(v \rho _2\right)d\rho _2d\rho _1$$,

So your problem can be evaluated by

Integrate[
BesselJ[0, v  r1] Cos[r1^2/2] r1  BesselJ[0, v  r2] Cos[
r2^2/2] r2  v, {r1, 0, 1}, {r2, 0, 1}, {v, 0, \[Infinity]}]


Unfortunately, the above integral does not converge, just as the error report of the above code and the answer from Goofy. But now we can do the above integral when The integration interval of v is finite:

f[v0_] :=
NIntegrate[
BesselJ[0, v  r1] Cos[r1^2/2] r1  BesselJ[0, v  r2] Cos[
r2^2/2] r2  v, {r1, 0, 1}, {r2, 0, 1}, {v, 0, v0}]

f[10]


Edit: Another way to find the value of the above integral.

First, mathematica cannot evaluate the principal value of the above integral directly:

Integrate[
BesselJ[0, v  r1] Cos[r1^2/2] r1  BesselJ[0, v  r2] Cos[
r2^2/2] r2  v, {r1, 0, 1}, {r2, 0, 1}, {v, 0, \[Infinity]},
PrincipalValue -> True]


returns nothing.

Then I try to calculate $$\int_0^1 r \cos \left(\frac{r^2}{2}\right) J_0(r v) \, dr$$ using Taylor series:

Since

$$\cos(r^2/2)=\underset{k=0}{\overset{\infty }{\sum }}\frac{\left(-\frac{1}{4}\right)^k \left(r^2\right)^{2 k}}{(2 k)!}$$,

using

int1=Assuming[k >= 0 && k \[Element] Integers && v > 0,
Integrate[
BesselJ[0, v  r] ((-(1/4))^k (r^2)^(2 k))/(2 k)! r, {r, 0, 1}]]

(*(-1)^k 2^(-1 - 2 k)
HypergeometricPFQRegularized[{1 + 2 k}, {1, 2 + 2 k}, -(v^2/4)]*)


we know that $$\int_0^1 r \cos \left(\frac{r^2}{2}\right) J_0(r v) \, dr=\sum _{k=0}^{\infty } (-1)^k 2^{-2 k-1} \, _1\tilde{F}_2\left(2 k+1;1,2 k+2;-\frac{v^2}{4}\right)$$

and so

Sum[int1*(int1 /. k -> j), {k, 0, \[Infinity]}, {j, 0, \[Infinity]}] //


$$\left(\int_0^1 r \cos \left(\frac{r^2}{2}\right) J_0(r v) \, dr\right){}^2=\sum _{k=0}^{\infty } \sum _{j=0}^{\infty } (-1)^{j+k} 2^{-2 j-2 k-2} \, _1\tilde{F}_2\left(2 j+1;1,2 j+2;-\frac{v^2}{4}\right) \, _1\tilde{F}_2\left(2 k+1;1,2 k+2;-\frac{v^2}{4}\right)$$

Then we do the final integration. By

sumpart =
Assuming[k >= 0 && j >= 0,
Integrate[int1*(int1 /. k -> j)*v, {v, 0, \[Infinity]}]]
(*((-1)^(j + k) 2^(
1 - 2 (1 + j + k)))/((1 + 2 j + 2 k) Gamma[1 + 2 j] Gamma[1 + 2 k])*)

Sum[sumpart, {k, 0, \[Infinity]}, {j, 0, \[Infinity]}]

(*Sin[1]/4*)


We obtain a simple value Sin[1]/4. It should be the principal value of the question.

Update: bugs with Sum

Instead of using Sum, NSum gives a different value:

NSum[sumpart, {j, 0, \[Infinity]}, {k, 0, \[Infinity]}]
(*0.460368 - 1.99638*10^-45 I*)

Sum[sumpart, {j, 0, \[Infinity]}, {k, 0, \[Infinity]}] // N
(*0.210368*)



Now NSum returns the value Sin[1]/4+1/4 match the comment by @MariuszIwaniuk and the answer by Goofy.

• It should be: 1/4 + Sin[1]/4` ? Commented May 25 at 13:52
• @MariuszIwaniuk How did you get that result? Commented May 25 at 14:28
• There was a typo in my code. The inside integral is $\sim 1/(4v^2)$ as $v \rightarrow \infty$, and therefore it is not divergent. The divergence you see in your answer is from changing the order of integration, I think. Commented May 25 at 19:08