Gauss-Laguerre quadrature

Gauss-Laguerre quadrature is approximation of integral of exponentially decaying function with the finite sum defined as follows:

\(\displaystyle \int_0^\infty f(x)\, \exp(-x)\, \mathrm d x \approx S\)

where \(f(x)\, \exp(-x)\) is integrand, \(f\) is supposed to be smooth (and preferably holomorphic) function,

\(\displaystyle S=\sum_{n=1}^M \, w_n \, f(x_n)\)

\(M\) is integer number (and should be for the precise approximation at non–polynomial \(f\), number \(M\) should be large)

\(x_n\) and \(w_n\) for integer \(n\) are coordinates of nodes and corresponding weights, characterised in the following:

\(x_n=\mathrm{LaguerreLzero}[M,n]\) is \(n\)th zero of \(M\)th polynomial LaguerreL,

\(\displaystyle w_n = \frac{x_n}{(M\!+\!1)\, |\mathrm{LaguerreL}(M\!+\!1,x_n)|\,} \)

Generalized Gauss Laguerre

The scheme described in the preamble seems to be simplest case of the Gauss-Laguerre quadrature. It can be generalised. For the application in atomic physics, the following integral is approximated:

\(\displaystyle \int_0^\infty x^\alpha\, \exp(-b x)\, f(x)\, \mathrm d x \approx S\)

where \(\alpha\) and \(b\) are constants.

The formula has the same form, the approximation \(S\) is written as

\(\displaystyle S=\sum_{n=1}^M \, w_n \, f(x_n)\)

John Burkardt

2010.02.23, John Burkardt suggests the C++ routine to evaluate coordiantes \(x_n\) and \weights \(w_n\) for \(0<n\le M\), at given \(\alpha\), \(b\) and \(M\). .

The code is loaded as http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.cpp and can be downloaded and distributed under the the GNU LGPL license (the use and free distribution is allowed, but the distributors should not modify the code).

The code reducers the following parameters:

number of points in the rule, denoted above with \(M\)

Order of the corresponding lagiuerr, denoted above with \(\alpha\)

\(A\) is the left endpoint (typically 0). (above it is explicitly set to zero)

\(B\) is the exponential scale factor, typically 1.

Here is the example of the use:

./gen_laguerre_rule 8 0 0 1 8n0a0a1

gerenates

file 8n0a0a1e_x.txt with 8 values of coordinates and

file 8n0a0a1_r.txt with corresponding weights.

However, the poetic users can give the output file some more meaningful name, than just "8n0a0a1".

If the user forgets to specify the parameters, then, the program by John Burkardt asks the user to specify them, and gets them from the command line.

Test with \(\alpha\!=\!0\)

Assuming, that gen_laguerre_rule is already compiled, the primitive test can be performed with C++ code below:

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
int main(){ int m,n,M=8; DB r,s,f,x[8],w[8]; FILE *o;
system("./gen_laguerre_rule 8 0 0 1 8n0a0a1");
 o=fopen("8n0a0a1_x.txt","r"); DO(n,8) { fscanf(o,"%lf", &r); x[n]=r; printf("%2d %20.16lf\n",n,r);} fclose(o);
o=fopen("8n0a0a1_w.txt","r"); DO(n,8) { fscanf(o,"%lf", &r); w[n]=r; printf("%2d %20.16lf\n",n,r);} fclose(o);
f=1.;
for(m=1;m<18;m++)
        { f*=m; s=0.;
        DO(n,8){ r=x[n]; s+=w[n]*pow(r,m); } printf("%2d %20.16le %20.16lf\n", m, s, s/f );
 }
}

The code above produced the following output:

03 August 2016 08:04:24 PM

GEN_LAGUERRE_RULE
  C++ version

  Compute a generalized Gauss-Laguerre rule for approximating
    Integral ( a <= x < +oo ) |x-a|^ALPHA exp(-B*x(x-a)) f(x) dx
  of order ORDER.

  The user specifies ORDER, ALPHA, A, B, and FILENAME.

  ORDER is the number of points in the rule:
  ALPHA is the exponent of |X|:
  A is the left endpoint (typically 0).
  B is the exponential scale factor, typically 1:
  FILENAME is used to generate 3 files:
  * filename_w.txt - the weight file
  * filename_x.txt - the abscissa file.
  * filename_r.txt - the region file.

  ORDER = 8
  ALPHA = 0
  A = 0
  B = 1
  FILENAME "8n0a0a1".

  Creating quadrature files.

  Root file name is "8n0a0a1".

  Weight file will be "8n0a0a1_w.txt".
  Abscissa file will be "8n0a0a1_x.txt".
  Region file will be "8n0a0a1_r.txt".

GEN_LAGUERRE_RULE:
  Normal end of execution.

03 August 2016 08:04:24 PM
 0 0.1702796323051008
 1 0.9037017767993794
 2 2.2510866298661290
 3 4.2667001702876570
 4 7.0459054023934637
 5 10.7585160101809905
 6 15.7406786412780004
 7 22.8631317368892688
 0 0.3691885893416376
 1 0.4187867808143426
 2 0.1757949866371721
 3 0.0333434922612156
 4 0.0027945362352257
 5 0.0000907650877336
 6 0.0000008485746716
 7 0.0000000010480012
 1 9.9999999999999956e-01 0.9999999999999996
 2 1.9999999999999987e+00 0.9999999999999993
 3 5.9999999999999938e+00 0.9999999999999990
 4 2.3999999999999972e+01 0.9999999999999988
 5 1.1999999999999987e+02 0.9999999999999989
 6 7.1999999999999909e+02 0.9999999999999988
 7 5.0399999999999945e+03 0.9999999999999989
 8 4.0319999999999949e+04 0.9999999999999988
 9 3.6287999999999930e+05 0.9999999999999981
10 3.6287999999999907e+06 0.9999999999999974
11 3.9916799999999873e+07 0.9999999999999968
12 4.7900159999999815e+08 0.9999999999999961
13 6.2270207999999723e+09 0.9999999999999956
14 8.7178291199999603e+10 0.9999999999999954
15 1.3076743679999939e+12 0.9999999999999953
16 2.0921164185599906e+13 0.9999222999222954
17 3.5545170124799856e+14 0.9993372640431424

First, the output by John Burkardt describes, what does the routine do; then, the tables of \(x\) and \(w\) are printed to verify the correct reading; then, the integral of power function is evaluated numerically,

\(\displaystyle_o^\infty x^m \exp(-x) \mathrm d x\)

for \(m\) from 1 to 17. As the result is expected to be close to \(m!\), ration of the first column to \(m!\) is printed in the last column.

while \(m\!<\! 2M\!-\!1\), the Gauss quadrature integrates "exactly", id est, with precision close to the maximal precision achievable with double variables; since \(m\!=\!16\), the precision drops down, as it is supposed to do.

It is not recommended to look inside the code: First, the license does not allow to modify it; Second, there seem to be misprint in the comment, it contains line

Compute a generalized Gauss-Laguerre rule for approximating
Integral ( a <= x < +oo ) |x-a|^ALPHA exp(-B*x(x-a)) f(x) dx
of order ORDER.

that seems to be wrong. (Extra "x" appears in the regiment of the exponential.) However, the misprint in the comment is not suppose to affect the calculus.

Application

Gauss-Laguerre quadrature seems to be efficient in the numerical calculations with approximation of atomic wave functions in spherical coordinates.

Also, the Gauss-Laguerre quadrature seems to be efficient for calculation of the Green function of the propagation go light in cylindrical coordinates.

References

https://en.wikipedia.org/wiki/Gauss–Laguerre_quadrature

Keywords

LaguerreL,