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
>
- References:
- Is it possible to make NIntegrate faster?
- From: Leo Alekseyev <dnquark@gmail.com>
- Is it possible to make NIntegrate faster?