MathGroup Archive 2009

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

Search the Archive

Re: Is it possible to make NIntegrate faster?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg105230] Re: [mg105163] Is it possible to make NIntegrate faster?
  • From: Leo Alekseyev <dnquark at gmail.com>
  • Date: Wed, 25 Nov 2009 02:31:29 -0500 (EST)
  • References: <200911231150.GAA21785@smc.vnet.net>

Apologies for missing some definitions in the code; they are
eps=3.*^-16;
xvec = Table[x, {x, 0.5, 1, 0.5/10}] (* not linspace *)

Thanks DrMajorBob for the suggestion, but implementing the integrands
as you said (creating branches with a Set symbolic result) and thus
removing every If[] statement inside the body of the function results
in absolutely no measurable performance gain...  This might actually
be a good thing -- this means that Mathematica is working at the
optimal level without the need for extensive "manual labor" in terms
of carefully defining the function.

--Leo

On Tue, Nov 24, 2009 at 1:42 AM, DrMajorBob <btreat1 at austin.rr.com> wrote:
> The complete version leaves "eps" and "Linspace" undefined, so I still can't
> test it much. Perhaps "eps" should be a parameter in
> ZDipoleFieldTotalIntegrand2?
>
> I'd also define branches of the function separately, remove N from the body
> of it, and Set a symbolic result for each branch (perhaps after Simplify).
> That could greatly reduce redundancy in applying the same logic on every
> call. In the first set of brackets defining the function, I'd place the
> variables that determine a logical branch, so that there's no logic on the
> right hand side.
>
> For instance, for rsc == True, x == y == 0, z < 0 (with r taking the place
> of "eps"):
>
> ZDipoleFieldTotalIntegrand2[True, 0 | 0., 0 | 0., z_?Negative, r_] >  Function[{k, h, k0, e1, ex2, e2perp},
>  Evaluate@Simplify@Module[{qz1, qz2, rr, tt},
>     qz1 = Sqrt[e1 k0^2 - k^2];
>     qz2 = Sqrt[ex2 (k0^2 - k^2/e2perp)];
>     rr = -((-ex2 qz1 + e1 qz2)/(ex2 qz1 + e1 qz2));
>     tt = (2 ex2 qz1)/(ex2 qz1 + e1 qz2);
>     (I/(2 Pi)*
>       tt (k^2*(k*BesselJ[0, k r] {0, 0, 1/e2perp}) Exp[-I qz2 z] Exp[
>           I qz1 h]/(2 qz1)))]
>  ]
>
> I eliminated some extraneous variables and algebraic terms in the Module,
> but it wouldn't matter if you didn't bother, since the logic and
> substitutions only execute once... rather than every time the case is
> encountered.
>
> N can be applied in the function that calls this one, but NIntegrate will
> pass Real variables anyway, so there's probably no need.
>
> Along these lines, I'm betting you could reach timings a lot closer to those
> of the alternative system.
>
> Bobby
>
> On Mon, 23 Nov 2009 23:18:54 -0600, Leo Alekseyev <dnquark at gmail.com> wrote:
>
>> 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
>>>
>
>
> --
> DrMajorBob at yahoo.com
>


  • Prev by Date: Re: Kernel crashes in ReplaceAll - with or without Maximize failure
  • Next by Date: Re: problem about calling Mathematica in Fortran
  • Previous by thread: Re: Is it possible to make NIntegrate faster?
  • Next by thread: Re: Re: Is it possible to make NIntegrate