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: 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 ); } }
 * 1) include
 * 2) include
 * 3) include
 * 4) define DB double
 * 5) define DO(x,y) for(x=0;x<y;x++)

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.

Keywords
LaguerreL,