Re: Is it possible to make NIntegrate faster?
- To: mathgroup at smc.vnet.net
- Subject: [mg105218] Re: [mg105163] Is it possible to make NIntegrate faster?
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Wed, 25 Nov 2009 02:29:15 -0500 (EST)
- References: <200911231150.GAA21785@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
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?