Numerical integration of a special integral

Question

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?

Answer 1

Here's a way to compute the result:

(* pre-process integrand for Levin Rule - saves ~10% of time *)
li = NIntegrate`LevinIntegrandReduce[BesselJ[0, v   r]  Cos[r^2/2] r, r];
int[v0_?NumericQ] := Block[{v = v0},
   NIntegrate["ignored", {r, 0, 1}, 
    Method -> {"LevinRule", "AdditiveTerm" -> li@"AdditiveTerm", 
      "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 *)

Answer 2

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]}] //
   TraditionalForm // TeXForm

$$\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.