emsr/quadrature

This is a C++ quadrature library reengineered from GSL and with new things added such as double exponential methods.

quadrature

Note: I originally envisioned providing this library as an extension to the gnu C++ standard library. I've given up on that idea. The standard C++ library is difficult enough to keep up with as it is. I've changed the namespace from __gnu_cxx to emsr (if nlohmann can do it so can I :-)). It stands for Excellent Math System - Reborn. Really!

This is a C++ quadrature library reengineered from GSL and with new things added. This was hived off from tr29124_test.

This library began as a way to test the orthogonal polynomials that were added to the Gnu implementation of the C++ standard library as part of the implementation of TR1 Special Math Functions. This library was dropped by the original author some years back. Then and now, the foundation of this library is the Gnu Scientific Library (GSL) quadrature package which itself is based on the venerable QUADPACK library.

This library is not a simple translation - that would not add real value.

The goals and ideal of this implementation are:

• Type genericity - allow different number systems to be used with the library including multi-precision ones. Currently, float, double, long double, and __float128 have been tested.
• Use C++ containers and algorithms throughout. The adaptive algorithms in GSL use what is esentially a priority queue to decide which chunk to work on next. C++ obviously has a lot to offer here.
• Add new quadrature algorithms. The GSL was focused on what I think is the middle of the spectrum: Gauss-Kronrod, Fejer and Clenshaw-Curtis, etc. This library adds simple recursive midpoint, trapezoid, and Simpson rules at the "low-end" and double-exponential, sinh-tanh rules at the "high-end".
• Support contour integration in the complex plane and for vector spaces.

Currently, this is a one-dimensional library. Cubature is a future goal. Monte-Carlo is another subject left out for now. The standard C++ library was developed in part to support this use-case and it is a worthy thing to have in a C++ library. After we get these done.

quadrature's Issues

Make sure that a == b returns zero quickly, and that b < a (b -> a) gives -(a -> b).

If a == b return the appropriate type of zero quickly.
If the number system is ordered then we can do the re-ordering.

type deduction fails for std::complex return types for qagX integration methods

Greetings,

I was testing out some of the integration routines for functions of complex type and real argument. The qagX_integrate algorithms fail to type deduce because the error handling functions cannot distinguish argument type return type. The same is true in the squad_integrate, however the fix there is simple, since one can take std::isinf(iv.fx[i]) || std::isnan(iv.fx[i]) to
std::isinf(std::abs(iv.fx[i])) || std::isnan(std::abs(iv.fx[i])) in all occurrences of the inf/nan test since the standard seems to only accept floating point types. The fixes for the qagX routines are a bit more involved I believe, unless I am missing something.

Cheers, Andrew

The qcheb algorithm is a mystery and is inflexible.

There is actually a good bit online about this. it's one of those ancient algos that gets passed around.

c***purpose  this routine computes the chebyshev series expansion
c            of degrees 12 and 24 of a function using a
c            fast fourier transform method
c            f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)),
c            f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)),
c            where t(k,x) is the chebyshev polynomial of degree k.
c        parameters
c          on entry
c           x      - real
c                    vector of dimension 11 containing the
c                    values cos(k*pi/24), k = 1, ..., 11
c
c           fval   - real
c                    vector of dimension 25 containing the
c                    function values at the points
c                    (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24,
c                    where (a,b) is the approximation interval.
c                    fval(1) and fval(25) are divided by two
c                    (these values are destroyed at output).
c
c          on return
c           cheb12 - real
c                    vector of dimension 13 containing the
c                    chebyshev coefficients for degree 12
c
c           cheb24 - real
c                    vector of dimension 25 containing the
c                    chebyshev coefficients for degree 24

So the polynomial expansion is in terms of T_n(x) and the roots are those of another Chebyshev polynomial U_{n-1}(x).

Let's have the ability to have expansions in T_{2m}(x) (and T_m(x)) at roots of U_{2m-1}(x).

In fact, I want a general Clenshaw thing for arbutrary polynomials with a recurrence relation.

Add Lebedev quadrature for 3D spherical (2D surface) integration.

See https://en.wikipedia.org/wiki/Lebedev_quadrature
There are implementations all over in all languages.
See gsl.
Lebedev, V. I. (1976). "Quadratures on a sphere". Zh. Vȳchisl. Mat. Mat. Fiz. 16 (2): 293–306. doi:10.1016/0041-5553(76)90100-2

Perhaps the tanh_sinh integral could use an open midpoint rule?

Perhaps the tanh_sinh integral could use an open midpoint rule?
This would only be useful if the trapezoid rule hits the endpoints.
OTOH, de rules don't work on bad endpoints anyway.

Oscillatory integrals and tests...

A paper: E. Sermutlu, H.T. Eyyubo~
glu / Applied Mathematics and Computation 189 (2007) 452–461
has analgorithm and a battery of 25 integrals.
Check out both the algorithm and the battery of tests for possible inclusion.

Make canned double exp integral rules.

Maybe max of 128. Lower power of 2 rules could skip in that single list. The weights will differ.

Try to decouple this from upstream projects like specfun (tr28124) and polynomial.

Maybe I can have polynomial as a submodule.
The upper incomplete gamma __gnu_cxx::tgamma is just cxx_math.

The double-exponential methods need error estimates.

They are recursive methods IIRC so it should be relatively easy.

__gnu_cxx::factorial<_Tp> not defined in the --gnu_cxx namespace

The factorial function is not defined in the __gnu_cxx namespace . i checked in jacobi.h and integrator.h and then searched the entire codebase to see where its defined but i could not find it anywhere

No member named 'factorial' in namespace '__gnu_cxx'

is it possible to add it