Stieltjes polynomials and integration

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

    \[ \int_{-1}^1f(x)\,dx\approx \sum_{k=1}^{n}w_kf(x_k) \]

where the weights w_k and the nodes (or abscissae) x_k are chosen so that the rule is exact for all polynomials of degree 2n-1.  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 n nodes will be accurate only for polynomials up to degree n.  (Although Newton-Cotes rules with odd numbers n of nodes will in fact be accurate for polynomials up to degree n+1).   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 w(x), so that we are approximating the integral

    \[ \int_{-1}^1w(x)f(x)\,dx \]

but we are interested here in the weight function w(x)=1.  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 P_0(x), P_1(x), P_2(x),\ldots which satisfy several useful properties:

  1. They are orthogonal over the interval [-1,1] with respect to the weight function w(x)=1, in that if m\ne n then

        \[ \int^1_{-1}P_n(x)P_m(x)\,dx=0. \]

  2. The roots of P_n(x) are all real, distinct, and lie in the region [-1,1].
  3. The Legendre polynomial P_n(x) can be computed as being the unique polynomial (up to a multiplicative constant) which satisfies

        \[ \int^1_{-1}P_n(x)x^k\,dx=0\mbox{\quad for $k=0,1,2,\ldots n-1$} \]

In fact they satisfy many more properties; see the wikipedia page or the Wolfram mathworld page.

The weights for Gaussian integration also can be computed in terms of the Legendre polynomials; if x_i is a root of P_n(x), then the corresponding weight is

    \[ w_i=\frac{2}{(1-x_i^2)[P'_n(x_i)]^2}. \]

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 a_1,a_2,a_3,\ldots a_n, (with a_i<a_{i+1}) and then compute the integral again at 2n+1 nodes which include a new set of n+1 nodes interlaced with the first set:

    \[ b_0<a_1<b_1<a_2<b_2<\cdots <b_{n-1}<a_n<b_n. \]

This means that instead of having to compute the function first at n nodes and again at 2n+1 nodes, we only have to compute the function at the n+1 new nodes b_i.

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 n-node Gaussian rule to 2n+1 nodes, and choosing the new nodes so that the rules would be exact for all polynomials of degree up to 2n+1.

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 E_n(x) and which can be uniquely determined, up to a multiplicative constant, by

    \[ \int^1_{-1}P_n(x)E_{n+1}(x)x^k\,dx = 0\mbox{\quad $k=0,1,2,\ldots n$}. \]

Stieltjes discussed these in a letter to Charles Hermite in 1894, and conjectured that the roots were all real, in the interval [-1,1], 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 w(x) 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)

[- 0.4058451514, - 0.7415311856, - 0.9491079123, 0.0, 0.9491079123, 0.7415311856, 0.4058451514]

(4) -> w7:=[2/((1-x^2)*(D7(x))^2) for x in n7]

[0.3818300505, 0.2797053915, 0.1294849662, 0.4179591837, 0.1294849662, 0.2797053915, 0.3818300505]

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 E_8(x).  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 x^8 to 1:

sol:=solve(concat(II,a[8]-1),[a[i] for i in 0..8])

    \[ \left[a_0 =\frac{52932681}{4854324041},a_1=0,a_2=\frac{202548}{653429},a_3=0,a_4=\frac{7794}{5491},a4=0,a_6=\frac{36}{17},a_7=0,a_8=1\right] \]

Now we can turn this into a polynomial:

E:=reduce(+,[x^(i-1)*rhs(sol.i) for i in 1..9])

and solve it:

 k8:=map(rhs,solve(E,1.0e-15))

[- 0.207784955, - 0.5860872355, - 0.8648644234, - 0.9914553711, 0.9914553711, 0.8648644234, 0.5860872355, 0.207784955]

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:

k15:=concat(n7,map(rhs,k8))

To determine the weights, we want an integral using all the nodes to be exact for all powers of x up to and including x^{15}. So first a tiny function to provide the value of that integral over [-1,1]:

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

[w_1 = 0.0229353220_1, w_2 = 0.0630920926_3, w_3 = 0.1047900103, w_4 = 0.1406532597, w_5 = 0.1690047266, w_6 = 0.1903505781, w_7 = 0.2044329401, w_8 = 0.2094821411, w_9 = 0.2044329401, w_{10} = 0.1903505781, w_{11} = 0.1690047266, w_{12} = 0.1406532597, w_{13} = 0.1047900103, w_{14} = 0.0630920926_3, w_{15} = 0.0229353220_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 n. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *