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(), and gausshermite(). This package is heavily influenced by Chebfun.

An introduction to Gauss quadrature can be found here.

• The fastest Julia code for Gauss quadrature nodes and weights (without tabulation).

• Change the perception that Gauss quadrature rules are expensive to compute.

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

tic(), gausschebyshev( 100000 ); toc()
elapsed time: 0.007636825 seconds

tic(), gausslegendre( 100000 ); toc() 
elapsed time: 0.017749388 seconds

tic(), gaussjacobi( 100000, .9, -.1 ); toc() 
elapsed time: 4.670444327 seconds

tic(), gaussradau( 100000 ); toc() 
elapsed time: 4.51240011 seconds

tic(), gausslobatto( 100000 ); toc() 
elapsed time: 3.989099163 seconds

tic(), gausshermite( 100000 ); toc()
elapsed time: 0.482970027 seconds

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

tic(), gausslegendre( 1000000001 ); toc()
elapsed time: 2526.19562756 seconds

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

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

1. 1st kind, weight function w(x) = 1/sqrt(1-x^2)

2. 2nd kind, weight function w(x) = sqrt(1-x^2)

3. 3rd kind, weight function w(x) = sqrt((1+x)/(1-x))

4. 4th kind, weight function w(x) = sqrt((1-x)/(1+x))

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

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].

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 [2].)

Warning: a and b need to be relatively small (-1<a,b<10).

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 [2].

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 [2].

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

• 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 [4].

The paper [4] 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.

[1] I. Bogaert, "Iteration-free computation of Gauss-Legendre quadrature nodes and weights", SIAM J. Sci. Comput., 36(3), A1008-A1026, 2014.

[2] N. Hale and A. Townsend, "Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights", SIAM J. Sci. Comput., 2012.

[3] J. C. Mason and D. C. Handscomb, "Chebyshev Polynomials", CRC Press, 2002.

[4] A. Townsend, T. Trogdon, and S. Olver, "Fast computation of Gauss quadrature nodes and weights on the whole real line", submitted, 2014.