View source for 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\).
<ref name="john">
http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html
John Burkardt.
GEN_LAGUERRE_RULE. Generalized Gauss-Laguerre Quadrature Rules.
2010.02.23.
</ref>.

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:
<poem><nomathjax><nowiki>
#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 );
 }
}
</nowiki></nomathjax></poem>

The code above produced the following output:
<poem><nomathjax><nowiki>
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
</nowiki></nomathjax></poem>

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==
<references/>

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

==Keywords==

[[LaguerreL]],

[[Category:English]]
[[Category:LaguerreL]]
[[Category:Numerical integration]]