Abstract
Methods for the computation of classical Gaussian quadrature rules are described which are effective both for small and large degree. These methods are reliable because the iterative computation of the nodes has guaranteed convergence, and they are fast due to their fourth-order convergence and its asymptotic exactness for an appropriate selection of the variables. For Gauss–Hermite and Gauss–Laguerre quadratures, local Taylor series can be used for computing efficiently the orthogonal polynomials involved, with exact initial values for the Hermite case and first values computed with a continued fraction for the Laguerre case. The resulting algorithms have almost unrestricted validity with respect to the parameters. Full relative precision is reached for the Hermite nodes, without any accuracy loss and for any degree, and a mild accuracy loss occurs for the Hermite and Laguerre weights as well as for the Laguerre nodes. These fast methods are exclusively based on convergent processes, which, together with the high order of convergence of the underlying iterative method, makes them particularly useful for high accuracy computations. We show examples of very high accuracy computations (of up to 1000 digits of accuracy).
This is a preview of subscription content, log in via an institution to check access.
Similar content being viewed by others
Notes
Notice that we have changed the notation for the nodes with respect to the previous section and the index runs differently since we are considering the positive nodes.
Of course, we could use other moments, as for instance \(\displaystyle \int _{-\infty }^{\infty } e^{-x^2}dx=\sqrt{\pi }= \delta _{1,k}w_0 f(0)+2\displaystyle \sum _{j=1}^{\lfloor n/2 \rfloor } w_j f(\alpha _j)\), \(k=n-2\lfloor n/2\rfloor \).
This is a consequence of the absolute error relation (see [24, Eq. (2.13)]) \(x^{(k+1)}-\alpha \approx \frac{\displaystyle {1}}{\displaystyle {12}}A'(\alpha )(x^{(k)}-\alpha )^4\), with \(\alpha \) the root that is computed, together with the fact that for the Hermite equation \(A(x)=2n+1-x^2\), which implies that \(1-\frac{\displaystyle {x^{(k+1)}}}{\displaystyle {\alpha }} \approx \frac{\displaystyle {1}}{\displaystyle {6}}(x^{(k)}-\alpha )^4\approx \frac{\displaystyle {1}}{\displaystyle {6}}(x^{(k)}-x^{(k+1)})^4\)
A fair comparison of efficiency between different methods should be always made by using implementations in a same platform and programming language (as is the case of the codes in this paper and [14]). This means that the different codes (if available) should be translated to a same language and carefully optimized. It would be certainly interesting to coherently benchmark the different available approaches, but this is outside the scope of the present paper.
Note that all the z-values generated by \(T_{-1}\) are an increasing sequence and therefore the stopping rule is safe.
We have \(z_e-z_{l}=(\alpha ^2-1/4)^{1/4}(1+\mathcal{O}(n^{-1/2}))\) and the maximum value of A(z) is reached at \(z=z_e\), where \(A(z_e)=2(L-\sqrt{\alpha ^2-1/4})\). Therefore, the distance between consecutive nodes in the z variable can be bounded by \(\varDelta z=z_{i+1}-z_{i}>\pi /\sqrt{A(z_e)}\); with this we estimated that the number of zeros smaller than \(z_e\) is \((z_e-z_l)\sqrt{A(z_e)}/\pi \sim \frac{\displaystyle {2(\alpha ^2-1/4)^{1/4}}}{\displaystyle {\pi }}\sqrt{n}\), which is an upper bound.
Or, if \(h(z_e)\) is very large, we can instead take \({\bar{y}}(z_e)=1\) and \(\dot{{\bar{y}}}(z_e)=1/h(z_e)\) to prevent overflows.
References
Bogaert, I.: Iteration-free computation of Gauss–Legendre quadrature nodes and weights. SIAM J. Sci. Comput. 36(3), A1008–A1026 (2014). https://doi.org/10.1137/140954969
Bogaert, I., Michiels, B., Fostier, J.: O(1) computation of Legendre polynomials and Gauss–Legendre nodes and weights for parallel computing. SIAM J. Sci. Comput. 34(3), C83–C101 (2012). https://doi.org/10.1137/110855442
Bremer, J.: On the numerical calculation of the roots of special functions satisfying second order ordinary differential equations. SIAM J. Sci. Comput. 39(1), A55–A82 (2017). https://doi.org/10.1137/16M1057139
Cash, J.: A note on the numerical solution of linear recurrence relations. Numer. Math. 34, 371–386 (1980). https://doi.org/10.1007/BF01403675
Davis, P., Rabinowitz, P.: Abscissas and weights for Gaussian quadratures of high order. J. Res. Nat. Bur. Stand. 56, 35–37 (1956)
Deaño, A., Huybrechs, D., Opsomer, P.: Construction and implementation of asymptotic expansions for Jacobi-type orthogonal polynomials. Adv. Comput. Math. 42(4), 791–822 (2016). https://doi.org/10.1007/s10444-015-9442-z
Deaño, A., Gil, A., Segura, J.: New inequalities from classical Sturm theorems. J. Approx. Theory 131(2), 208–230 (2004). https://doi.org/10.1016/j.jat.2004.09.006
Deaño, A., Segura, J.: Global Sturm inequalities for the real zeros of the solutions of the Gauss hypergeometric differential equation. J. Approx. Theory 148(1), 92–110 (2007). https://doi.org/10.1016/j.jat.2007.02.005
Dimitrov, D.K., Nikolov, G.P.: Sharp bounds for the extreme zeros of classical orthogonal polynomials. J. Approx. Theory 162(10), 1793–1804 (2010). https://doi.org/10.1016/j.jat.2009.11.006
Gautschi, W.: Gauss–Radau formulae for Jacobi and Laguerre weight functions. Math. Comput. Simul. 54(4–5), 403–412 (2000). https://doi.org/10.1016/S0378-4754(00)00179-8. 1999 International Symposium on Computational Sciences, to honor John R. Rice (West Lafayette, IN)
Gautschi, W.: Generalized Gauss–Radau and Gauss–Lobatto formulae. BIT 44(4), 711–720 (2004). https://doi.org/10.1007/s10543-004-3812-0
Gil, A., Segura, J., Temme, N.M.: Numerical Methods for Special Functions. Society for Industrial and Applied Mathematics (SIAM), Philadelphia (2007). https://doi.org/10.1137/1.9780898717822
Gil, A., Segura, J., Temme, N.M.: Recent software developments for special functions in the Santander–Amsterdam project. Sci. Comput. Program. 90A, 42–54 (2014)
Gil, A., Segura, J., Temme, N.M.: Asymptotic approximations to the nodes and weights of Gauss–Hermite and Gauss–Laguerre quadratures. Stud. Appl. Math. 140(3), 298–332 (2018). https://doi.org/10.1111/sapm.12201
Gil, A., Segura, J., Temme, N.M.: Noniterative computation of Gauss–Jacobi quadrature. SIAM J. Sci. Comput. 41(1), A668–A693 (2019). https://doi.org/10.1137/18M1179006
Glaser, A., Liu, X., Rokhlin, V.: A fast algorithm for the calculation of the roots of special functions. SIAM J. Sci. Comput. 29(4), 1420–1438 (2007). https://doi.org/10.1137/06067016X
Golub, G.H., Welsch, J.H.: Calculation of Gauss quadrature rules. Math. Comput. 23(1969), 221–230. addendum, ibid. 23(106, loose microfiche suppl), A1–A10 (1969)
Hale, N., Townsend, A.: Fast and accurate computation of Gauss–Legendre and Gauss–Jacobi quadrature nodes and weights. SIAM J. Sci. Comput. 35(2), A652–A674 (2013). https://doi.org/10.1137/120889873
Huybrechs, D., Opsomer, P.: Construction and implementation of asymptotic expansions for Laguerre-type orthogonal polynomials. IMA J. Numer. Anal. 38(3), 1085–1118 (2018). https://doi.org/10.1093/imanum/drx030
Johansson, F., Mezzarobba, M.: Fast and rigorous arbitrary-precision computation of Gauss–Legendre quadrature nodes and weights. arXiv preprint arXiv:1802.03948 (2018)
Kreuser, P.: Über das Verhalten der Integrale homogener linearer Differenzengleichungen im Unendlichen. Dissertation. Tübingen, 48 S (1914)
Lowan, A.N., Davids, N., Levenson, A.: Table of the zeros of the Legendre polynomials of order 1–16 and the weight coefficients for Gauss’ mechanical quadrature formula. Bull. Am. Math. Soc. 48, 739–743 (1942)
Olver, F.W.J.: Asymptotics and Special Functions. AKP Classics. A K Peters, Ltd., Wellesley (1997). (Reprint of the 1974 original [Academic Press, New York; MR0435697 (55 #8655)])
Segura, J.: Reliable computation of the zeros of solutions of second order linear ODEs using a fourth order method. SIAM J. Numer. Anal. 48(2), 452–469 (2010). https://doi.org/10.1137/090747762
Segura, J., Temme, N.M.: Numerically satisfactory solutions of Kummer recurrence relations. Numer. Math. 111(1), 109–119 (2008). https://doi.org/10.1007/s00211-008-0175-5
Swarztrauber, P.N.: On computing the points and weights for Gauss–Legendre quadrature. SIAM J. Sci. Comput. 24(3), 945–954 (2002). https://doi.org/10.1137/S1064827500379690. (electronic)
Townsend, A., Trogdon, T., Olver, S.: Fast computation of Gauss quadrature nodes and weights on the whole real line. IMA J. Numer. Anal. 36(1), 337–358 (2016)
Wilf, H.S.: Mathematics for the Physical Sciences. Dover Publications, Inc., New York (1978). (Reprinting of the 1962 original, Dover Books in Advanced Mathematics)
Yakimiw, E.: Accurate computation of weights in classical Gauss–Christoffel quadrature rules. J. Comput. Phys. 129(2), 406–430 (1996). https://doi.org/10.1006/jcph.1996.0258
Acknowledgements
The authors thank the anonymous referees for their constructive comments and suggestions. NMT thanks CWI for scientific support.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work was supported by Ministerio de Ciencia, Innovación y Universidades, projects MTM2015-67142-P (MINECO/FEDER, UE) and PGC2018-098279-B-I00 (MCIU/AEI/FEDER,UE).
Rights and permissions
About this article
Cite this article
Gil, A., Segura, J. & Temme, N.M. Fast, reliable and unrestricted iterative computation of Gauss–Hermite and Gauss–Laguerre quadratures. Numer. Math. 143, 649–682 (2019). https://doi.org/10.1007/s00211-019-01066-2
Received:
Revised:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-019-01066-2