quadrature_least_squares, a FORTRAN90 code which computes weights for "sub-interpolatory" quadrature rules.
A large class of quadrature rules may be computed by specifying a set of N abscissas, or sample points, X(1:N), determining the Lagrange interpolation basis functions L(1:N), and then setting a weight vector W by
W(i) = I(L(i))
after which, the integral of any function f(x) is estimated by
I(f) \approx Q(f) = sum ( 1 <= i <= N ) W(i) * f(X(i))
We call this an interpolatory rule because the function f(x) has first been interpolated by
f(x) \approx sum ( 1 <= i <= N ) L(i) * f(X(i))
after which, we apply the integration operator:
I(f) \approx I(sum ( 1 <= i <= N ) L(i) * f(X(i)))
= sum ( 1 <= i <= N ) I(L(i)) * f(X(i))
= sum ( 1 <= i <= N ) W(i) * f(X(i)).
For badly chosen sets of X, or high values of N, or unruly functions f(x), interpolation may be a bad way to approximate the function. An alternative is to seek a polynomial interpolant of degree D < N-1, and then integrate that. We might call this a "sub-interpolatory" rule.
As it turns out, a natural way to seek such a rule is to write out the N by D+1 Vandermonde matrix and use a least squares solver. Even though the N by N Vandermonde matrix is ill-conditioned for Gauss elimination, a least squares approach can produce usable solutions from the N by D+1 matrix.
The outline of this procedure was devised by Professor Mac Hyman of Tulane University.
The computer code and data files described and made available on this web page are distributed under the MIT license
quadrature_least_squares is available in a C version and a C++ version and a Fortran90 version and a MATLAB version and an Octave version.
clenshaw_curtis_rule, a FORTRAN90 code which defines a multiple dimension Clenshaw Curtis quadrature rule.
QR_SOLVE, a FORTRAN90 code which computes the least squares solution of a rectangular linear system A*x=b.
QUADMOM, a FORTRAN90 code which computes a Gaussian quadrature rule for a weight function rho(x) based on the Golub-Welsch procedure that only requires knowledge of the moments of rho(x).
QUADRATURE_GOLUB_WELSCH, a FORTRAN90 code which computes the points and weights of a Gaussian quadrature rule using the Golub-Welsch procedure, assuming that the points have been specified.
QUADRULE, a FORTRAN90 code which defines quadrature rules for approximating an integral over a 1D domain.
QUADRULE_FAST, a FORTRAN90 code which defines efficient versions of a few 1D quadrature rules.
TEST_INT, a FORTRAN90 code which defines test integrands for 1D quadrature rules.
TRUNCATED_NORMAL_RULE, a FORTRAN90 code which computes a quadrature rule for a normal probability density function (PDF), also called a Gaussian distribution, that has been truncated to [A,+oo), (-oo,B] or [A,B].
VANDERMONDE, a FORTRAN90 code which carries out certain operations associated with the Vandermonde matrix.