Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Is it possible to make NIntegrate faster?

• To: mathgroup at smc.vnet.net
• Subject: [mg105212] Re: [mg105163] Is it possible to make NIntegrate faster?
• From: Leo Alekseyev <dnquark at gmail.com>
• Date: Tue, 24 Nov 2009 05:50:47 -0500 (EST)
• References: <200911231150.GAA21785@smc.vnet.net>

Good point, DrMajorBob; I seem to have underestimated the telepathic
powers of  Mathgroup.  As you will see below, the problem with my
NIntegrate was that Mathematica tried to do SymbolicProcessing where
it didn't need to.

Andrew Moylan writes:
> You integrand contains a Bessel function over an infinite range. Specialized
> methods exist for numerical integrals of that sort, and NIntegrate automatically
> selects one.  However, since your integral decays very fast, the oscillatory
> nature of the integrand is not very important. Therefore the default methods are
> quite adequate (and faster). You can choose to use only default methods with
> this option: Method -> {Automatic, SymbolicProcessing -> 0} Of course if, for
> some other parameters, your integral doesn't decay as fast and the rapidly
> oscillatory nature becomes important, then you may still wish to use the
> specialized method for Bessel integrands.

Using this suggestion, and playing with AccuracyGoal settings, I was
able to get a factor of 30-40 performance improvement.  This is still
about 5 times slower compared to an alternative system (but fast
enough so that I no longer feel compelled to rewrite my Mathematica
code).  I will include the integrands below, and if off the top of
your head you can suggest any other optimizations, please do.
Personally, I'm out of ideas -- everything is wrapped in N[] and no
SetDelayed is used anywhere.

Here is the short version, to illustrate the sort of integral I'm working with:

ZDipoleFieldTIntegrandRSC[xi_][x_, y_, z_, h_, k0_, e1_, e2par_,
 e2perp_] :=
 Block[{r = Sqrt[x^2 + y^2],  m1 = 1, m2 = 1,
   qz1, qz2, rr, tt, ex1 = e1, ex2 = e2par,
  nintOpts}, qz1 = Sqrt[e1 m1 k0^2 - (xi/r)^2];
 qz2 = Sqrt[e2par (m2 - (xi/r)^2/e2perp)];
 tt = (2 ex2 qz1)/(ex2 qz1 + ex1 qz2);
 I/(2 Pi)*
  tt ((xi/r)^2*((xi/r)*BesselJ[0, xi] {0, 0, 1/e2perp} +
      qz2 {x/r/ex2, y/r/ex2, 0}*I*BesselJ[1, xi]) Exp[-I qz2 z] Exp[
      I qz1 h]/(2 r qz1))]

Module[{xi, nintOpts = {MaxRecursion -> 400, AccuracyGoal -> 5}},
 With[{
  z = -0.2, x = 0.1, y = 0, k0 = 1, h = 0.5, e1 = 1, e2par = =
1,
  e2perp = 1
  },
 NIntegrate[
  ZDipoleFieldTIntegrandRSC[xi][x, y, z, h, k0, e1, e2par,
   e2perp], {xi, N[Sqrt[x^2 + y^2] Sqrt[e1] k0], Infinity},
  Evaluate[Sequence @@ nintOpts]]]
 ]


And here's the complete version (involves sums of a couple of terms
like the one above and a condition that changes the integrand based on
whether or not z>0)


ZDipoleFieldTotalIntegrand2[kkk_, rsc_][x_, y_, z_, h_, k0_, e1_,
  e2par_, e2perp_] :=
 Module[{r, intlim = Infinity, sgn, m1 = 1, m2 = 1, qz1, qz2,
   ex1 = e1, ex2 = e2par, nintOpts, qz1fn, qz2fn, rr, tt, k, corr},
  r = Sqrt[x^2 + y^2]; If[r == 0, r = eps];
  sgn = If[z < h, 1., -1.];
  If[rsc == False, k = kkk; corr = 1., k = kkk/r; corr = 1./r];
  qz1 = Sqrt[e1 m1 k0^2 - k^2] // N;
  qz2 = Sqrt[e2par (m2 k0^2 - k^2/e2perp)] // N;
  rr = -((-ex2 qz1 + ex1 qz2)/(ex2 qz1 + ex1 qz2)) // N;
  tt = (2 ex2 qz1)/(ex2 qz1 + ex1 qz2) // N;
  corr*If[z >= 0,
    N[I/(2 Pi) 1/
        e1 (k^2*(k*BesselJ[0, k r] {0, 0, 1} +
           sgn*qz1 {x/r, y/r, 0}*I*BesselJ[1, k r]) Exp[
           I qz1 Abs[z - h]]/(2 qz1))] +
     N[I/(2 Pi e1)*
       rr (k^2*(k*BesselJ[0, k r] {0, 0, 1} -
           qz1 {x/r, y/r, 0}*I*BesselJ[1, k r]) Exp[
           I qz1 (z + h)]/(2 qz1))],
    N[I/(2 Pi)*
      tt (k^2*(k*BesselJ[0, k r] {0, 0, 1/e2perp} +
          qz2 {x/r/ex2, y/r/ex2, 0}*I*
           BesselJ[1, k r]) Exp[-I qz2 z] Exp[I qz1 h]/(2 qz1))]]]

ZDipoleFieldTotal2[x_?NumericQ, y_?NumericQ, z_?NumericQ, h_?NumericQ,
   k0_?NumericQ, e1_?NumericQ, e2par_?NumericQ, e2perp_?NumericQ,
  opts___?OptionQ] :=
 Module[{r = Sqrt[x^2 + y^2], intlim = Infinity,
   sgn = If[z < h, 1, -1], k, nintOpts = Flatten[{opts}]},
  NIntegrate[
    ZDipoleFieldTotalIntegrand2[k, False][x, y, z, h, k0, e1, e2par,
     e2perp], {k, 0, Sqrt[e1] k0}, Evaluate[Sequence @@ nintOpts]] +
   NIntegrate[
    ZDipoleFieldTotalIntegrand2[k, True][x, y, z, h, k0, e1, e2par,
     e2perp], {k, N[r Sqrt[e1] k0], Infinity},
    Evaluate[Sequence @@ nintOpts]]]

With[{
  z = -0.2,
  x = 0.5,
  xvec = Linspace[.5, 1, 11],
  y = 0, k0 = 1, h = 0.5, e1 = 1, e2par = 1, e2perp = 1},
 ZDipoleFieldTotal2[x, y, z, h, k0, e1, e2par, e2perp]
 ]:


On Mon, Nov 23, 2009 at 10:19 PM, DrMajorBob <btreat1 at austin.rr.com> wrote:
> You told us everything but what we need to know: the integrand.
>
> It may be defined with SetDelayed, where Set could be much faster. There may
> be other issues.
>
> But in a vacuum... who knows?
>
> Bobby
>
> On Mon, 23 Nov 2009 05:50:40 -0600, Leo Alekseyev <dnquark at gmail.com> wrote:
>
>> Dear Mathgroup,
>>
>> Recently I have been using NIntegrate fairly extensively.  I am
>> dealing with an oscillatory integral that has a singularity.
>> NIntegrate is able to treat it reasonably well -- the only default I
>> had to change was increasing MaxRecursion.  However, it is slow.  10
>> points of my integrand take about 40 seconds to evaluate.  After I
>> ported my code to another system, this same integral took about a second
>> using
>> the Gauss-Kronrod method (quadgk in the other system).  Furthermore, by
>> increasing the absolute and relative tolerance values, I could improve
>> the speed without losing too much precision, so currently the
>> integrals evaluate in 0.4 seconds.
>>
>> I have been playing with various NIntegrate parameters to try to
>> improve the speed, to no effect.  My integrands are straightforward
>> (although long) algebraic expressions involving a few Bessel functions
>> and exponentials, wrapped inside a Module; all subexpressions use N[]
>> so that nothing should be symbolic...  Ideally I hoped to find some
>> sort of a speed/accuracy tradeoff, but that hasn't happened.
>>
>> I read the numerical integration tutorial in the docs, but am finding
>> it hard to figure out how to improve the efficiency of my integration.
>>  I would expect Mathematica to get to at least within an order of
>> magnitude of the other system using the same integration strategy.  The
>> current
>> performance isn't satisfactory -- but neither is the solution of
>> porting perfectly good Mathematica code to the other system...
>>
>> I would much appreciate any suggestions.
>>
>> Thanks,
>> --Leo
>>
>
>
> --
> DrMajorBob at yahoo.com
>

• References:
  • Is it possible to make NIntegrate faster?
    • From: Leo Alekseyev <dnquark@gmail.com>