RFC: quadgk integration routine

[stevengj]

This patch adds a new routine quadgk (inspired by the Matlab quadgk function) to Julia:

I,E = quadgk(f, a,b,c...; reltol=eps*100, abstol=0, maxevals=10^7, order=7)

which performs adaptive 1d Gauss-Kronrod integration (addressing issue #1235), returning an estimated integral I and error E.

For example, quadgk(cos, 0,1) integrates cos from 0 to 1, returning sin(1). Complex-valued functions and complex endpoints (contour integration) are also supported, as is arbitrary-precision arithmetic. e.g. quadgk(cos, 0, BigFloat(1), order=40) integrates cos using BigFloat arithmetic, with a higher-order rule to get exponential accuracy.

The code is all written from scratch by me, with the most complicated routine being the kronrod function to compute Kronrod points and weights; I implemented this from a pseudo-code description in a well-known 1997 paper by Laurie.

[stevengj]

Update: quadgk now also handles infinite and semi-infinite real intervals (via the usual coordinate-transformation trick).

[StefanKarpinski]

Phenomenal. Since this is new functionality, I'm in favor of just merging it. The only question to me is if we should have this in base or a package.

[ViralBShah]

@stevengj How come you wrote the Cubature routines as a package? Would it best to have one place for all the numerical integration routines? I do not have a view one way or another, but just wanted to know. People are used to having some quadrature routines in the base library due to Matlab's inclusion of various quad routines in their base library.

[kmsquire]

@ViralBShah, the Cubature.jl package is GPL, presumably because the C library it depends on is GPL.

(Since @stevengj wrote both, he may have the ability to change the license on both, but perhaps he doesn't want to, or can't because the C library includes code copyrighted by others.)

[stevengj]

I don't have the ability to change the license of the Cubature package, since it uses some code from other GPL projects.

In any case, it seemed good to have a native Julia implementation of at least 1d quadrature, so that it can handle arbitrary numeric types, including BigFloat; that will never happen in a C library like the Cubature package. And it seems like at least 1d integration should be bundled with Julia proper, as that is such a basic task.

(Integration is such an important task, and there are so many different algorithms for different types of problems, that it seems both inevitable and desirable that we eventually have multiple external packages as well. e.g. it would be nice to have a package wrapping the (GPL) Cuba library, too. We already have the GSL package to wrap the GSL routines, too.)

[stevengj]

Regarding Matlab's various quad* routines, it seems like quadgk is enough since, unlike Matlab, the order is an adjustable parameter. e.g. if you use order=1 then it is not too different from adaptive Simpson (for low-accuracy integration of nonsmooth integrands). I don't see much point in having both quadl and quadgk; the two seem more-or-less redundant to me, with quadgk being nicer since it doesn't evaluate the integrand at the endpoints (which makes it easier to handle singular integrands).

[ViralBShah]

Thanks for the explanation. Please merge when ready.

[coveralls]

Changes Unknown when pulling 949421f on stevengj:quadgk into * on JuliaLang:master*.

[IainNZ]

Hah what the heck, why/how did Coveralls just run

[Keno]

Clicking on the link, that doesn't appear to be our repository? API key collision?

[staticfloat]

Hmmm, yeah, that's pretty bad. We should let somebody know about that.