New LevinIntegrate package for highly oscillatory integration
Andrew Moylan
My LevinIntegrate package, for numerical integration of highly oscillatory
functions, is now available at
http://andrew.j.moylan.googlepages.com/levinintegrate.
LevinIntegrate is an automatic integrator based on Levin's method for highly
oscillatory integration. It outperforms NIntegrate (and other automatic
integrators) for a variety of irregularly oscillatory functions:
LevinIntegrate[Exp[I*Cosh[x]], {x, 0, 10}] // Timing
>> {0.02, -0.138718 + 1.20194 * I}
Compare with:
NIntegrate[Exp[I*Cosh[x]], {x, 0, 10}, MaxRecursion -> 50] // Timing
>> (warnings)
>> {4.546, -0.138718 + 1.20194 * I}
For more information or to download the package, see
http://andrew.j.moylan.googlepages.com/levinintegrate.
This is an early release of LevinIntegrate; it is still under development.
If you have any comments, suggestions, or problems, please do let me know of
them!
Andrew Moylan
antononcube
plug-in mechanism. It would be best if the package is made in such a
way that it provides a new integration strategy. (I.e. the plug-in is
its ultimate goal.)
Below is an example how the hook up can be done just using
LevinIntegrate. Note that with that new strategy, LevinOscillatory,
piecewise functions are properly integrated. E.g. this works
In[21]:= NIntegrate[If[x < 5, Sqrt[x], Exp[I*Cosh[x]]], {x, 0, 10},
Method -> LevinOscillatory] // Timing
Out[21]= {0.021987,-0.138718+1.20194 I}
but LevinIntegrate gives up:
In[23]:= LevinIntegrate[If[x < 5, Sqrt[x], Exp[I*Cosh[x]]], {x, 0,
10}] // Timing
During evaluation of In[23]:= LevinIntegrate::unknownoscillator: \
LevinIntegrate only operates on integrands containing oscillatory \
factors matching E^_|BesselJ[_Integer,_]|Sin[_]|Cos[_]. If your \
integrand can still be integrated using a Levin-type method, use \
GeneralisedLevinIntegrate.
Out[23]= {0.001462,Null}
Plug-in code:
Needs["Moylan`LevinIntegrate`"]
Clear[LevinOscillatory]
LevinOscillatory /:
NIntegrate`InitializeIntegrationStrategy[LevinOscillatory, nfs_,
ranges_, strOpts_, allOpts_] :=
Block[{},
LevinOscillatory[{First /@ ranges, strOpts}]
];
LevinOscillatory[{vars_, strOpts_}]["Algorithm"[regions_, opts___]] :=
Module[{integrands, ranges, res},
integrands = (#@"Integrand")@"FunctionExpression" & /@ regions;
ranges =
First@Outer[Prepend[#1, #2] &, #@"Boundaries", vars, 1] & /@
regions;
res = MapThread[
LevinIntegrate[#1, Sequence @@ #2,
Sequence @@ strOpts] &, {integrands, ranges}];
If[FreeQ[res, Null],
Total@res,
Total@
MapThread[
If[NumberQ[#1], #1,
NIntegrate[#2, Sequence @@ #3 // Evaluate,
DeleteCases[opts, Method -> _] // Evaluate]] &, {res,
integrands, ranges}]
]
];
Anton Antonov,
Wolfram Research, Inc.
On Oct 10, 3:51 am, "Andrew Moylan" <andrew.j.moy...@gmail.com> wrote:
> Hi all,
>
> My LevinIntegrate package, for numerical integration of highly oscillatory
> functions, is now available athttp://andrew.j.moylan.googlepages.com/levinintegrate.
>
> LevinIntegrate is an automatic integrator based on Levin's method for highly
> oscillatory integration. It outperformsNIntegrate(and other automatic
> integrators) for a variety of irregularly oscillatory functions:
>
> LevinIntegrate[Exp[I*Cosh[x]], {x, 0, 10}] // Timing
>
> >> {0.02, -0.138718 + 1.20194 * I}
>
> Compare with:
>
> NIntegrate[Exp[I*Cosh[x]], {x, 0, 10}, MaxRecursion -> 50] // Timing
>
> >> (warnings)
> >> {4.546, -0.138718 + 1.20194 * I}
>
> For more information or to download the package, seehttp://andrew.j.moylan.googlepages.com/levinintegrate.
Andrew Moylan
Yeah I agree: the correct place for LevinIntegrate is as a rule or
strategy in NIntegrate. Gaining access to NIntegrate's pre-processing
(for piecewise functions, a la your example (but shouldn't Out[21] be
7.46604+0.00482644*I?)) is a good example of why that will be
valuable.
Thanks for your example demonstrating how to set up an integration
strategy that calls LevinIntegrate. It makes me think that setting up
LevinIntegrate as a rule or strategy ("LevinOscillatory"---I like it)
won't be all that tricky. Is there documentation available to help
authors of integration strategies?
Related matters:
1. Do there exist sensible rules for *automatically* selecting the
LevinOscillatory strategy during pre-processing (a la the rules that
NIntegrate pre-processing presently uses to select the existing
Oscillatory Clenshaw-Curtis rules when the frequency of the Sin/Cos/
Exp[I*] in question is large enough)?
1a. Is it difficult (for a package author) to extend NIntegrate's
pre-processing / rule selection procedure?
The answer to 1. is unclear to me. The whole problem with irregularly
oscillatory integrals is that it's no longer easy to judge "how"
oscillatory they are over a given range of the integration parameter:
Exp[I*g[x]] <-- hard to decide whether this is sufficiently
oscillatory to warrant applying LevinIntegrate rather than just
GaussKronrod, especially since g may be g[x_?NumericQ]:=bla.
Andrew
Robert Dodier
> My LevinIntegrate package, for numerical integration of highly oscillatory
> functions, is now available athttp://andrew.j.moylan.googlepages.com/levinintegrate.
Andrew, thanks for writing this package.
I'll recommend that you include some kind of license so that
other people know the limits of their rights, if any, wrt your
package.
With no explicit license statement, the rights of others are
limited to what is called "fair use" in the US; dunno what the
default might be in other parts of the world. You've already done
us all a big favor simply by creating the package; may I suggest
you take it one step further and tell us what we can do with it.
Sorry if you've already covered this ground; I looked in the
package and on your web page without finding a license statement.
best
Robert Dodier
P_ter
I would like to have an explanation of this so called "so that other people know the limits of their rights". It is filling in for other people and they did not even say anything yet.
The benevolent attitude to people exists everywhere. Also the people who see that a human is a dangerous animal.
I am thankful for those initiatives which just give. And I do not like systems nor arguments which are negative to these initiatives. I cite as an example for what I see as benevolent, the 100$ PC from Nicholas Negroponte. Many people have the same attitude on a smaller scale.
with friendly greetings,
P_ter