# Is there any way to make this numerical integral evaluate faster?

k = 100;
\[Delta]pz = 0;
ks = 1.18;
int[x0_?NumericQ, X_?NumericQ] := NIntegrate[ 2 Sqrt[\[Pi]] E^(-x^2/2) E^(-(x - X)^2/2) r/Sqrt[(x - x0)^2 + r^2]E^(-ks Sqrt[r^2 + (x - x0 )^2])BesselJ[0, Sqrt[(X + 2*k)^2 + \[Delta]pz^2] r], {x, -10, 10}, {r, 0, 10}, PrecisionGoal -> 3, AccuracyGoal -> 8, MinRecursion -> 1]


The evaluation of int[0,0] takes around 6 seconds. I would like to considerably reduce this time. Here is the plot of the integrand

(*plot*)
Plot3D[Evaluate[2 N[Sqrt[\[Pi]], 30] E^(-x^2/2) E^(-(x - X)^2/2) r/Sqrt[(x - x0)^2 + r^2]E^(-ks Sqrt[r^2 + (x - x0 )^2])BesselJ[0, Sqrt[(X + 2*k)^2 + \[Delta]pz^2] r] /. {x0 -> 0, X -> 0}], {x, -5, 5}, {r, 0, 10}, PlotRange -> All]


which has a series of oscillating peaks. What is the best method to integrate this kind of function?

• First try: support.wolfram.com/kb/3442 . Esp turn off symbolic processing if speed is an issue. May 16, 2013 at 16:42
• This is helpful. It takes 90 seconds now. Thanks!! Other ideas are welcomed. May 16, 2013 at 17:34
• What are k, ks and \[Delta]pz? For instance, this takes less than a second to evaluate: NIntegrate[2 Sqrt[\[Pi]] E^(-x^2/2) E^(-(x - X)^2/2) r/Sqrt[(x - x0)^2 + r^2] E^(-ks Sqrt[r^2 + (x - x0)^2]) BesselJ[0, Sqrt[(X + 2*k)^2 + \[Delta]pz^2] r] /. {k -> 1, ks -> 1, \[Delta]pz -> 1, x0 -> 0, X -> 0}, {x, -10, 10}, {r, 0, 10}] May 17, 2013 at 2:24
• They are important numerical values. I have edited the question with their values. Thanks for pointing it out. May 17, 2013 at 12:19
• Have you seen the relevant section for "LevinRule" ? It seems to be spot on for oscillatory integrals with small oscillation around a kernel.
– gpap
May 17, 2013 at 14:07

The Levin rule by itself doesn't help much, but if we increase the number of collocation points some, we get a little increase in speed. If we integrate from -Infinity to Infinity, we get some more (in case -10 to 10 is meant to be a finite approximation of an infinite interval, a common numeric trick that is usually not advisable in Mathematica).

k = 100;
δpz = 0;
ks = 118/100;
integrand =
2 Sqrt[π] E^(-x^2/2) E^(-(x - X)^2/2) r/
Sqrt[(x - x0)^2 + r^2] E^(-ks Sqrt[r^2 + (x - x0)^2]) BesselJ[0,
Sqrt[(X + 2*k)^2 + δpz^2] r];

ClearAll[int];
int[x00_?NumericQ, X0_?NumericQ, opts___] :=
Block[{x0 = x00, X = X0},
NIntegrate[integrand, {x, -10, 10}, {r, 0, 10}, opts,
PrecisionGoal -> 3, AccuracyGoal -> 8, MinRecursion -> 1]
];

int[0, 0] // AbsoluteTiming
(*  {2.32213, 0.00017724}  *)

int[0, 0, Method -> {"LevinRule"}] // AbsoluteTiming
(*  {2.45073, 0.00017724}  *)

int[0, 0, Method -> {"LevinRule", "Points" -> 11}] // AbsoluteTiming
(*  {1.61575, 0.00017724}  *)

Block[{x0 = 0, X = 0},
NIntegrate[integrand,
{x, -Infinity , Infinity}, {r, 0, 10},
PrecisionGoal -> 3, AccuracyGoal -> 8,
Method -> {"LevinRule", "Points" -> 11}]
] // AbsoluteTiming
(*  {1.11198, 0.000177237}  *)