5
$\begingroup$
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]

plot

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

$\endgroup$
  • 2
    $\begingroup$ First try: support.wolfram.com/kb/3442 . Esp turn off symbolic processing if speed is an issue. $\endgroup$ – Searke May 16 '13 at 16:42
  • $\begingroup$ This is helpful. It takes 90 seconds now. Thanks!! Other ideas are welcomed. $\endgroup$ – lagoa May 16 '13 at 17:34
  • 3
    $\begingroup$ 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}] $\endgroup$ – Michael E2 May 17 '13 at 2:24
  • $\begingroup$ They are important numerical values. I have edited the question with their values. Thanks for pointing it out. $\endgroup$ – lagoa May 17 '13 at 12:19
  • 1
    $\begingroup$ Have you seen the relevant section for "LevinRule" ? It seems to be spot on for oscillatory integrals with small oscillation around a kernel. $\endgroup$ – gpap May 17 '13 at 14:07
2
$\begingroup$

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}  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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