In a slow and (in my case, uninformed) effort to add some numerical quadrature capabilities to PanAxiom, I’ve been looking at QUADPACK and in particular the GSL implementation of it, and the GSLL wrapper of GSL with Lisp. (Note that QUADPACK as a Fortran 90 library is available through netlib.)
My interest has been taken up with Gauss-Kronrod quadrature, which up until recently was no more than a name to me – I knew nothing of its theory or its practice.
As you know, Gaussian quadrature takes the form
where the weights and the nodes (or abscissae) are chosen so that the rule is exact for all polynomials of degree . This means that Gaussian rules will in general be far more accurate that Newton-Cotes rules with the same number of nodes, since a Newton-Cotes rule with nodes will be accurate only for polynomials up to degree . (Although Newton-Cotes rules with odd numbers of nodes will in fact be accurate for polynomials up to degree ). The formula above can be easily transformed to an arbitrary interval on integration by a linear transformation.
In fact, more generally Gaussian quadrature is defined with a weight functions , so that we are approximating the integral
but we are interested here in the weight function . For the Gaussian formula above, it turns out that the nodes are the roots of the Legendre polynomials of the first type. These are a sequence of polynomials which satisfy several useful properties:
- They are orthogonal over the interval with respect to the weight function , in that if then
- The roots of are all real, distinct, and lie in the region .
- The Legendre polynomial can be computed as being the unique polynomial (up to a multiplicative constant) which satisfies
The weights for Gaussian integration also can be computed in terms of the Legendre polynomials; if is a root of , then the corresponding weight is
There’s a problem with Gaussian integration and that is error calculation. When you attempt to calculate an integral numerically, you need some way of estimating the error. The most efficient way (which is easy with Newton-Cotes rules), is to compare the results of two computations: the nodes of one being a subset of the nodes of the other. So we might compute our integral at nodes , (with ) and then compute the integral again at nodes which include a new set of nodes interlaced with the first set:
This means that instead of having to compute the function first at nodes and again at nodes, we only have to compute the function at the new nodes .
But the nodes of one Gaussian rule will not be subsets of nodes of a higher order rule. In 1964, the Russian mathematician Alexander Kronrod had the idea of extending an -node Gaussian rule to nodes, and choosing the new nodes so that the rules would be exact for all polynomials of degree up to .
It turns out, although Kronrod didn’t know this himself, that the new nodes were the roots of another set of polynomials called Stieltjes polynomials, denoted and which can be uniquely determined, up to a multiplicative constant, by
Stieltjes discussed these in a letter to Charles Hermite in 1894, and conjectured that the roots were all real, in the interval , and interlaced with the zeros of the Legendre polynomials. These conjectures were all proved by Gabriel Szegö in 1931.
The computation of Gauss-Kronrod nodes and weights, and extensions to other weight functions is an active area of research, especially in regards to their asymptotic properties. But in fact we can throw all of the theory away and compute nodes and weights easily by using a modern computer algebra system and the above definitions.
For example, consider the Gaussian seven point rule, and the extra eight nodes required to form a Kronrod extension. First we need to compute the nodes and weights, using the Legendre polynomials. This can be done easily in PanAxiom:
(1) -> L7(x) == legendreP(7,x)
(2) -> D7 == D(L7,x)
(3) -> n7 := solve(L7,1.0e-10)
(4) -> w7:=[2/((1-x^2)*(D7(x))^2) for x in n7]
For ease of display we are just showing 10 decimal places, but we can easily compute, and show, more. You can check these values against the ones here, for example.
Now for the Kronrod extension, and the first thing is to determine the Stieltjes polynonial . We start by forming a generic polynomial with unknown coefficients:
(5) -> P:=reduce(+,[a[i]*x^i for i in 0..8])
and then create the list of equations relating to the integral definition above
(6) -> II:List POLY FRAC INT:=[integrate(x^k*legendreP(7,x)*P,x=-1..1)::EXPR FRAC INT for k in 0..7]
In order to obtain a single polynomial, we shall set the coefficient of to 1:
sol:=solve(concat(II,a-1),[a[i] for i in 0..8])
Now we can turn this into a polynomial:
E:=reduce(+,[x^(i-1)*rhs(sol.i) for i in 1..9])
and solve it:
These are the Kronrod nodes which, as you see, interlace with the initial Gaussian nodes obtain a roots of the Legendre polynomial. So, put ’em all together:
To determine the weights, we want an integral using all the nodes to be exact for all powers of up to and including . So first a tiny function to provide the value of that integral over :
int(i)==if rem(i,2)=0 then 2/(i+1) else 0
Now, the equations for the weights:
wteqs:=[reduce(+,[w[i]*(k15.i)^m for i in 1..15])-int(m) for m in 0..14]
and their solution:
solve(wteqs,[w[i] for i in 1..15]).1
Again, these can be checked at the Gauss-Kronrod wikipedia page.
The use of a computer algebra system means that the above computations can be done to arbitrary precision – as long as we ensure that the working precision is greater than the displayed precision, to account for rounding errors. If for example in PanAxiom we set digits(50) and outputGeneral(25) and run through the above computations again, we find the same results as those computed by Pavel Holoborodko.
(Note: Pavel is the author of the QuickLaTeX extension for WordPress, which allows complete LaTeX code in a post – I am using it here and in all my posts. As far as I’m concerned, QuickLaTeX is an extension worth shifting from wordpress.com to one’s own wordpress.org. No other extension comes close for its handling of mathematics.)
I suspect that this somewhat bone-headed method will fall apart, or simply get drowned by rounding errors, for large . The now more standard approach is based on some work by Dirk Laurie in the late 1990’s. The Gaussian nodes can be computed first of all as the eigenvalues of a particular tridiagonal matrix; Laurie showed how to extend this result so as to be able to compute the Kronrod nodes. Versions of his method are now standard for developing Gauss-Kronrod rules.
Note also that the whole business of trying to interlace nodes with Gaussian rules is ameliorated by using Clenshaw-Curtis quadrature where the nodes of any orders are automatically nested in the nodes of higher orders. And Clenshaw-Curtis quadrature seems to be roughly as accurate as Gaussian quadrature.