Add gl node weight expansions

README.md

@@ -56,7 +56,7 @@ There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four
 
 4. 4th kind, weight function `w(x) = sqrt((1-x)/(1+x))`
 
-They are all have explicit simple formulas for the nodes and weights <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&pg=PA180#v=onepage&q&f=false">[4]</a>.
+They are all have explicit simple formulas for the nodes and weights <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&pg=PA180#v=onepage&q&f=false">[3]</a>.
 ## The algorithm for Gauss-Legendre
 Gauss quadrature for the weight function `w(x) = 1`.
 
@@ -71,28 +71,31 @@ 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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.)
+*  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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[2]</a>.)
 
 * 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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.
+The Gauss-Radau nodes and weights can be computed via the `(0,1)` Gauss-Jacobi nodes and weights <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[2]</a>.
 
 ## 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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.
+The Gauss-Lobatto nodes and weights can be computed via the `(1,1)` Gauss-Jacobi nodes and weights <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[2]</a>.
 
 ## The algorithm for Gauss-Laguerre
-Gauss quadrature for the weight function `w(x) = exp(-x)` on `[0,Inf)`
+Gauss quadrature for the weight function `w(x) = exp(-x)`, `w(x) = x.^alpha.*exp(-x)` or `w(x) = x.^alpha.*exp(-qm*x.^m)` on `[0,Inf)`
 
-* For `n<128`: Use the Golub-Welsch algorithm. 
+* For `n<128`: Use forward recurrence. 
 
-* For `method=GLR`: 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 <a href="http://epubs.siam.org/doi/pdf/10.1137/06067016X">[2]</a>.
+* For `method=GW`: Use the Golub-Welsh algorithm.
+
+* For `n>=128`: Use explicit asymptotic expansions [5] and this is O(sqrt(n)) when allowed to stop when the weights are below the smallest positive floating point number.
+
+* For `method=gen`: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [4], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight `w(x) = x^alpha exp(-q_m x^m)` and asymptotics of the orthogonal polynomials can also be called by the user.
 
-* For `n>=128`: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight `w(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 Gauss-Hermite
 Gauss quadrature for the weight function `w(x) = exp(-x^2)` on the real line.
@@ -118,15 +121,15 @@ The paper <a href="http://arxiv.org/abs/1410.5286">[7]</a> also derives an `O(n)
 
 ## References:
 [1] I. Bogaert, <a href="http://epubs.siam.org/doi/abs/10.1137/140954969">"Iteration-free computation of Gauss-Legendre quadrature nodes and weights"</a>, SIAM J. Sci. Comput., 36(3), A1008-A1026, 2014.
-
-[2] A. Glaser, X. Liu, and V. Rokhlin. <a href="http://epubs.siam.org/doi/pdf/10.1137/06067016X">"A fast algorithm for the calculation of the roots of special functions."</a> SIAM J. Sci. Comput., 29 (2007), 1420-1438.
 
-[3] N. Hale and A. Townsend, <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">"Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature
+[2] N. Hale and A. Townsend, <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">"Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature
        nodes and weights"</a>, SIAM J. Sci. Comput., 2012.
 
-[4] J. C. Mason and D. C. Handscomb, <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&dq=Mason%20and%20Handscomb&pg=PP1#v=onepage&q=Mason%20and%20Handscomb&f=false">"Chebyshev Polynomials"</a>, CRC Press, 2002.
+[3] J. C. Mason and D. C. Handscomb, <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&dq=Mason%20and%20Handscomb&pg=PP1#v=onepage&q=Mason%20and%20Handscomb&f=false">"Chebyshev Polynomials"</a>, CRC Press, 2002.
+
+[4] D. Huybrechs and P. Opsomer. Construction and implementation of asymptotic expansions for Laguerre-type orthogonal polynomials. IMA J. Numer. Anal., 2017, Accepted.
 
-[5] P. Opsomer, (in preparation).
+[5] D. Huybrechs and P. Opsomer. Arbitrary-order asymptotic expansions of generalized Gaussian quadrature rules, 2017, in preparation.
 
 [6] A. Townsend, <a href="http://math.mit.edu/~ajt/papers/QuadratureEssay.pdf"> The race for high order Gauss-Legendre quadrature</a>, in SIAM News, March 2015.  
 

src/FastGaussQuadrature.jl

@@ -7,6 +7,7 @@ using Compat, SpecialFunctions
 export gausslegendre
 export gausschebyshev
 export gausslaguerre
+export gaussfreud
 export gausshermite
 export gaussjacobi
 export gausslobatto
@@ -18,6 +19,7 @@ import SpecialFunctions: besselj, airyai, airyaiprime
 include("gausslegendre.jl")
 include("gausschebyshev.jl")
 include("gausslaguerre.jl")
+include("gaussfreud.jl")
 include("gausshermite.jl")
 include("gaussjacobi.jl")
 include("gausslobatto.jl")