FastGaussQuadrature.jl
A Julia package to compute n
point Gauss quadrature nodes and weights to 16digit 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 GaussLegendre quadrature, see [6].
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 GaussLegendre nodes. So here we will do a billion + 1.
@time gausslegendre( 1000000001 );
131.392154 seconds (17 allocations: 26.077 GB, 1.17% gc time)
(The nodes near the endpoints coalesce in 16digits of precision.)
The algorithm for GaussChebyshev
There are four kinds of GaussChebyshev quadrature rules, corresponding to four weight functions:

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

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

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

4th kind, weight function
w(x) = sqrt((1x)/(1+x))
They are all have explicit simple formulas for the nodes and weights [4].
The algorithm for GaussLegendre
Gauss quadrature for the weight function w(x) = 1
.

For
n<=5
: Use an analytic expression. 
For
n<=60
: Use Newton's method to solvePn(x)=0
. EvaluatePn
andPn'
by 3term recurrence. Weights are related toPn'
. 
For
n>60
: Use asymptotic expansions for the Legendre nodes and weights [1].
The algorithm for GaussJacobi
Gauss quadrature for the weight functions w(x) = (1x)^a(1+x)^b
, a,b>1
.

For
n<=100
: Use Newton's method to solvePn(x)=0
. EvaluatePn
andPn'
by threeterm recurrence. 
For
n>100
: Use Newton's method to solvePn(x)=0
. EvaluatePn
andPn'
by an asymptotic expansion (in the interior of[1,1]
) and the threeterm recurrenceO(n^2)
close to the endpoints. (This is a small modification to the algorithm described in [3].) 
For
max(a,b)>5
: Use the GolubWelsch algorithm requiringO(n^2)
operations.
The algorithm for GaussRadau
Gauss quadrature for the weight function w(x)=1
, except the endpoint 1
is included as a quadrature node.
The GaussRadau nodes and weights can be computed via the (0,1)
GaussJacobi nodes and weights [3].
The algorithm for GaussLobatto
Gauss quadrature for the weight function w(x)=1
, except the endpoints 1
and 1
are included as nodes.
The GaussLobatto nodes and weights can be computed via the (1,1)
GaussJacobi nodes and weights [3].
The algorithm for GaussLaguerre
Gauss quadrature for the weight function w(x) = exp(x)
on [0,Inf)

For
n<128
: Use the GolubWelsch algorithm. 
For
method=GLR
: Use the GlaserLuiRohklin algorithm. EvaluateLn
andLn'
by using Taylor series expansions near roots generated by solving the secondorder differential equation thatLn
satisfies, see [2]. 
For
n>=128
: Use a Newton procedure on RiemannHilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weightw(x) = x^alpha exp(q_m x^m)
and this is O(sqrt(n)) when allowed to stop when the weights are below the smallest positive floating point number.
The algorithm for GaussHermite
Gauss quadrature for the weight function w(x) = exp(x^2)
on the real line.

For
n<200
: Use Newton's method to solveHn(x)=0
. EvaluateHn
andHn'
by threeterm recurrence. 
For
n>=200
: Use Newton's method to solveHn(x)=0
. EvaluateHn
andHn'
by a uniform asymptotic expansion, see [7].
The paper [7] also derives an O(n)
algorithm for generalized GaussHermite 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
References:
[1] I. Bogaert, "Iterationfree computation of GaussLegendre quadrature nodes and weights", SIAM J. Sci. Comput., 36(3), A1008A1026, 2014.
[2] A. Glaser, X. Liu, and V. Rokhlin. "A fast algorithm for the calculation of the roots of special functions." SIAM J. Sci. Comput., 29 (2007), 14201438.
[3] N. Hale and A. Townsend, "Fast and accurate computation of GaussLegendre and GaussJacobi quadrature nodes and weights", SIAM J. Sci. Comput., 2012.
[4] J. C. Mason and D. C. Handscomb, "Chebyshev Polynomials", CRC Press, 2002.
[5] P. Opsomer, (in preparation).
[6] A. Townsend, The race for high order GaussLegendre quadrature, in SIAM News, March 2015.
[7] A. Townsend, T. Trogdon, and S. Olver, "Fast computation of Gauss quadrature nodes and weights on the whole real line", to appear in IMA Numer. Anal., 2014.
[8] M. Vanlessen, "Strong asymptotics of LaguerreType orthogonal polynomials and applications in Random Matrix Theory", Constr. Approx., 25:125175, 2007.