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?

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ First try: support.wolfram.com/kb/3442 . Esp turn off symbolic processing if speed is an issue. $\endgroup$
    – Searke
    Commented May 16, 2013 at 16:42
  • $\begingroup$ This is helpful. It takes 90 seconds now. Thanks!! Other ideas are welcomed. $\endgroup$
    – lagoa
    Commented May 16, 2013 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
    Commented May 17, 2013 at 2:24
  • $\begingroup$ They are important numerical values. I have edited the question with their values. Thanks for pointing it out. $\endgroup$
    – lagoa
    Commented May 17, 2013 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
    Commented May 17, 2013 at 14:07

1 Answer 1

Reset to default
3
$\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}  *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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