GLQ: A Gauss-Laguerre Quadrature solver

Lightweight C++ library for calculating Gauss-Laguerre quadrature abscissas and weights.

This package is based on the algorithms described in:

• W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, Third Edition, Cambridge University Press, 2007.

Gauss-Laguerre Quadrature

The Gauss-Laguerre approximates the following improper integral:

$$\int_{0}^{+\infty}x^\alpha e^{-x} f(x)dx \approx \sum_{i=1}^N w_i f(x_i)$$

Where $x_i$ is the root of the Laguerre polynomial defined as:

$$L_0^\alpha(x) = 1$$ $$L_1^\alpha(x) = 1 + \alpha - x$$ $$L_{n}^\alpha(x) = \frac{\left(2n - 1 + \alpha - x\right)L_{n-1}^\alpha(x) - (n - 1 + \alpha)L_{n-2}^\alpha(x)}{n}$$

For $n \geq 2$.

and $w_i$ the corresponding weight from:

$$w_i = \frac{\Gamma(n+\alpha+1)x_i}{n!\left(n+1\right)^2\left[L_{n+1}^{\alpha}(x_i)\right]^2}$$

Root finding method

The Laguerre polynomial roots are determined with the Newton-Raphson method through:

$$x_{k+1} = x_k - \frac{L_n^\alpha(x_k)}{L_n^{'\alpha}(x_k)}$$

Where $L_n^{'\alpha}(x_k)$ is the polynomial derivative through:

$$L_n^{'\alpha}(x_k) = \frac{nL_n^{\alpha}(x_k) - (n + \alpha)L_{n-2}^{\alpha}(x_k)}{x_k}$$ Once the approximation changes less than the defined NEWTON_EPSILON, the root is saved as $x_i$.

Initial guess

The Newton-Raphson method requires an initial guess defined as the following:

For $i=0$:

$$x'_i = \frac{(1 + \alpha)(3 + 0.92\alpha)}{1 + 2.4N + 1.8\alpha}$$

For $i=1$:

$$x'_i = x_{i-1} + \frac{15 + 6.25\alpha}{1 + 0.9\alpha + 2.5N}$$

For $i\geq2$:

$$x'_i = x_{i-1} + \frac{x_{i-1} - x_{i-2}}{1 + 0.3\alpha}\left(\frac{1 + 2.55(i-1)}{1.9(i-1)}+\frac{1.26(i-1)}{1 + 3.5(i-1)}\right)$$

Usage

#include "glq.hpp"
#include <iostream>

int main() {
  int num_nodes = 16;
  glq::Quadrature q = glq::quadrature(num_nodes);

  std::cout << "Abscissas:\n";
  for (double x : q.abscissas) {
    std::cout << x << " ";
  }
  std::cout << "\nWeights:\n";
  for (double w : q.weights) {
    std::cout << w << " ";
  }
  std::cout << std::endl;
}