Proposal: gaussian quadrature

[hsnyder]

Hi everyone.

Would there be interest in adding methods to compute Gaussian quadrature points and weights (plus possibly Gauss-Lobatto points and weights)? I'd be happy to offer the implementation if there is interest.

[arjenmarkus]

Sofar the library contains the trapezoid and Simpson rules only. More sophisticated rules are certainly welcome, I'd say. Op di 22 dec. 2020 om 06:13 schreef Harris Snyder <notifications@github.com

…

: Hi everyone. Would there be interest in adding methods to compute Gaussian quadrature points and weights (plus possibly Gauss-Lobatto points and weights)? I'd be happy to offer the implementation if there is interest. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#277>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN6YR7REA2IL4CCH7ACJGLSWATIJANCNFSM4VFEBWXA> .

[ivan-pi]

I'd be certainly interested in these. In a previous reply #112 (comment) I suggested these might be more suitable for a separate package. On second thought, since many Fortran users are physicists and engineers, at least for the classic Gaussian rules it might be beneficial to provide them here. If nothing else, the feedback from other stdlib developers and page visitors can help design a tidy API.

I have some very old Gaussian quadrature routines dating back to 1966 available here: https://github.com/ivan-pi/stroud_quad. If I am not mistaken they are accurate up to 10-11 digits for most cases.

[jvdp1]

Thank you for your proposition.
I agree wth @arjenmarkus and @ivan-pi: it would be good to have them in stdlib ;)

[hsnyder]

Okay, sounds like the idea has some support. I'll wait another couple of days (holidays coming up anyway) and if there aren't any objections raised we can discuss an API

[certik]

@hsnyder do not wait. We want this in stdlib, as we agreed the scope of stdlib can be roughly as scipy, and thus it definitely includes quadrature rules. Go ahead and start proposing the API. See also my comment with an implementation accurate to full double precision: #112 (comment)

[certik]

I think this is a duplicate of #112. But I'll keep it open, we can close both #277 and #112 once this is implemented in stdlib.

[ivan-pi]

The NAG and IMSL mathematical libraries both offer Gaussian quadrature weights:

• NAG
• IMSL

[hsnyder]

Thanks for the input everyone. Lots of options re: implementation it seems (depending on licenses)

Regarding API, I would strongly prefer returning the points and weights directly and having the user use sum(f(x)*w) rather than having the user pass in a function and having the library perform the integration.

One idea regarding the API: something virtually identical to what @certik has here: https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/quadrature.f90
i.e.
function gauss_pts(n) and function gauss_weights(n)

but since you'll always want both, I could see the argument for a subroutine
subroutine gauss_xw(n,x,w), where x and w would be intent(out)

We should provide analagous functions for lobatto points and weights, as well, and perhaps others. Thoughts?

[jvdp1]

I am not a typical user for such functions. However, I suggest to use a similar API to trapz and trapz_weights

[certik]

@hsnyder yes, we can lump them together to return points and weights at once. The same for lobatto. Indeed, returning them to the user is the way to go.

[hsnyder]

@jvdp1 I think that since the use case for gaussian quadrature is a bit different (usually you're integrating a function rather than discrete samples), it makes more sense to return the points and weights instead.

@certik sounds good, I like your tabulated implementation from hfsolver. I propose to use that as the basis for an implementation, and fall back on another implementation in the n > 64 case. From what various people have posted in this thread and in #112, it seems like there are several options. Is there a favourite among them? Quadpack perhaps?

[hsnyder]

I've started working on an implementation... Should we provide separate implementations for getting the points and weights in real32, or support double precision only and let the user convert if they care? If we do want to offer real32 versions, is it important that the whole implementation be real32, or can the computation be done in 64 and then down-converted?

I'd be surprised if there's a lot of use of 32 bit platforms anymore, but figured I'd ask.

[arjenmarkus]

Wrt 32 bits platforms: probably not that much indeed, but that is a different matter than 32 bits versus 64 bits reals ;).

For the sort of problems I work on professionally, 32 bits (or, I prefer to talk of single precision, the two are not entirely equivalent), offer in general quite enough precision. Single-precision means 5 or 6 decimals. If you consider such quantities as the water depth (to borrow from my field), that would typically mean better than 1 mm precision (assuming water depths up to 1 km). That is far more precise than any measurements in that field. Using single-precision means a more economic use of memory and disk space. That is certainly not to say that double-precision does not enter my field of expertise - I merely mean that both types of real precision are still useful.

That said, I would say that for the purposes of this part of the standard library it suffices to have double-precision points and weights and have a convenient interface to convert to single-precision.

[hsnyder]

Hi @arjenmarkus

Good points - it was lazy of me to conflate 32 bit reals with 32 bit platforms. As per your suggestion, I'll provide interfaces for both kinds of reals, but the back-end computation of the nodes and weights will be done in 64 bit.