FastGaussQuadrature

A Julia package to compute n-point Gauss quadrature nodes and weights to 16-digit accuracy and in O(n) time. So far the package includes gausschebyshev(), gausslegendre(), gaussjacobi(), gaussradau(), gausslobatto(), gausslaguerre(), and gausshermite(). This package is heavily influenced by Chebfun.

An introduction to Gauss quadrature can be found here. For a quirky account on the history of computing Gauss-Legendre quadrature, see [4].

Our Aims

The fastest Julia code for Gauss quadrature nodes and weights (without tabulation).
Change the perception that Gauss quadrature rules are expensive to compute.

Examples

Here we compute 100000 nodes and weights of the Gauss rules. Try a million or ten million.

@time gausschebyshev( 100000 );
0.002681 seconds (9 allocations: 1.526 MB, 228.45% gc time)

@time gausslegendre( 100000 ); 
0.007110 seconds (17 allocations: 2.671 MB)

@time gaussjacobi( 100000, .9, -.1 );
1.782347 seconds (20.84 k allocations: 1.611 GB, 22.89% gc time)

@time gaussradau( 100000 );
1.849520 seconds (741.84 k allocations: 1.625 GB, 22.59% gc time)

@time gausslobatto( 100000 );
1.905083 seconds (819.73 k allocations: 1.626 GB, 23.45% gc time)

@time gausslaguerre( 100000 )
.891567 seconds (115.19 M allocations: 3.540 GB, 3.05% gc time)

@time gausshermite( 100000 );
0.249756 seconds (201.22 k allocations: 131.643 MB, 4.92% gc time)

The paper [1] computed a billion Gauss-Legendre nodes. So here we will do a billion + 1. This is (probably) a world record:

@time gausslegendre( 1000000001 );
131.392154 seconds (17 allocations: 26.077 GB, 1.17% gc time)

(The nodes near the endpoints coalesce in 16-digits of precision.)

The algorithm for Gauss-Chebyshev

There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions:

1st kind, weight function w(x) = 1/sqrt(1-x^2)
2nd kind, weight function w(x) = sqrt(1-x^2)
3rd kind, weight function w(x) = sqrt((1+x)/(1-x))
4th kind, weight function w(x) = sqrt((1-x)/(1+x))

They are all have explicit simple formulas for the nodes and weights [4].

The algorithm for Gauss-Legendre

Gauss quadrature for the weight function w(x) = 1.

For n<=5: Use an analytic expression.
For n<=60: Use Newton's method to solve Pn(x)=0. Evaluate Pn and Pn' by 3-term recurrence. Weights are related to Pn'.
For n>60: Use asymptotic expansions for the Legendre nodes and weights [1].

The algorithm for Gauss-Jacobi

Gauss quadrature for the weight functions w(x) = (1-x)^a(1+x)^b, a,b>-1.

For n<=100: Use Newton's method to solve Pn(x)=0. Evaluate Pn and Pn' by three-term recurrence.
For n>100: Use Newton's method to solve Pn(x)=0. Evaluate Pn and Pn' by an asymptotic expansion (in the interior of [-1,1]) and the three-term recurrence O(n^-2) close to the endpoints. (This is a small modification to the algorithm described in [3].)
For max(a,b)>5: Use the Golub-Welsch algorithm requiring O(n^2) operations.

The algorithm for Gauss-Radau

Gauss quadrature for the weight function w(x)=1, except the endpoint -1 is included as a quadrature node.

The Gauss-Radau nodes and weights can be computed via the (0,1) Gauss-Jacobi nodes and weights [3].

The algorithm for Gauss-Lobatto

Gauss quadrature for the weight function w(x)=1, except the endpoints -1 and 1 are included as nodes.

The Gauss-Lobatto nodes and weights can be computed via the (1,1) Gauss-Jacobi nodes and weights [3].

The algorithm for Gauss-Laguerre

Gauss quadrature for the weight function w(x) = exp(-x) on [0,Inf)

For n<=128: Use the Golub-Welsch algorithm.
For n>128: Use the Glaser-Lui-Rohklin algorithm. Evaluate Ln and Ln' by using Taylor series expansions near roots generated by solving the second-order differential equation that Ln satisfies, see [2].

The algorithm for Gauss-Hermite

Gauss quadrature for the weight function w(x) = exp(-x^2) on the real line.

For n<200: Use Newton's method to solve Hn(x)=0. Evaluate Hn and Hn' by three-term recurrence.
For n>=200: Use Newton's method to solve Hn(x)=0. Evaluate Hn and Hn' by a uniform asymptotic expansion, see [6]. * The paper [6] also derives an O(n) algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form exp(-V(x)), where V(x) is a real polynomial.

Example usage

@time nodes, weights = gausslegendre( 100000 );
0.007890 seconds (19 allocations: 2.671 MB)

# integrates f(x) = x^2 from -1 to 1
@time dot( weights, nodes.^2 )
0.004264 seconds (7 allocations: 781.484 KB)
0.666666666666666