Re: Is it possible to make NIntegrate faster?
Leo Alekseyev
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 <btr...@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 <dnq...@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
>>
>
>
> --
> DrMaj...@yahoo.com
>
DrMajorBob
Leo Alekseyev
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 <btr...@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
>> z = -0.2, x = 0.1, y = 0, k0 = 1, h = 0.5, e1 = 1, e2par = 1,
> --
> DrMaj...@yahoo.com
>
DrMajorBob
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_] =
DrMajorBob
it's not surprising to find cases where it isn't.
The dissected version might be more readable in some situations, too.
Bobby
On Wed, 25 Nov 2009 01:31:29 -0600, Leo Alekseyev <dnq...@gmail.com>
wrote:
> 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 <btr...@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_] >