Difference between revisions of "GaussLaguerre quadrature"
m (Text replacement  "\$([^\$]+)\$" to "\\(\1\\)") 

Line 1:  Line 1:  
[[GaussLaguerre quadrature]] is approximation of integral of exponentially decaying function with the finite sum defined as follows: 
[[GaussLaguerre 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 
+  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)\,} 
w_n = \frac{x_n}{(M\!+\!1)\, \mathrm{LaguerreL}(M\!+\!1,x_n)\,} 

+  \) 

−  $ 

==Generalized Gauss Laguerre== 
==Generalized Gauss Laguerre== 

The scheme described in the preamble seems to be simplest case of the [[GaussLaguerre quadrature]]. 
The scheme described in the preamble seems to be simplest case of the [[GaussLaguerre quadrature]]. 

Line 21:  Line 21:  
For the application in atomic physics, the following integral is approximated: 
For the application in atomic physics, the following integral is approximated: 

−  +  \(\displaystyle 

\int_0^\infty 
\int_0^\infty 

−  x^\alpha\, \exp(b x)\, f(x)\, \mathrm d x \approx S 
+  x^\alpha\, \exp(b x)\, f(x)\, \mathrm d x \approx S\) 
−  where 
+  where \(\alpha\) and \(b\) are constants. 
−  The formula has the same form, the approximation 
+  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 == 
== John Burkardt == 

2010.02.23, John Burkardt suggests the [[C++]] routine to evaluate 
2010.02.23, John Burkardt suggests the [[C++]] routine to evaluate 

−  coordiantes 
+  coordiantes \(x_n\) and \weights \(w_n\) for \(0<n\le M\), 
−  at given 
+  at given \(\alpha\), \(b\) and \(M\). 
<ref name="john"> 
<ref name="john"> 

http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html 
http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html 

Line 50:  Line 50:  
The code reducers the following parameters: 
The code reducers the following parameters: 

−  number of points in the rule, denoted above with 
+  number of points in the rule, denoted above with \(M\) 
−  Order of the corresponding lagiuerr, denoted above with 
+  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: 
Here is the example of the use: 

Line 73:  Line 73:  
John Burkardt asks the user to specify them, and gets them from the command line. 
John Burkardt asks the user to specify them, and gets them from the command line. 

−  ===Test with 
+  ===Test with \(\alpha\!=\!0\) === 
Assuming, that gen_laguerre_rule is already compiled, 
Assuming, that gen_laguerre_rule is already compiled, 

Line 171:  Line 171:  
First, the output by John Burkardt describes, what does the routine do; 
First, the output by John Burkardt describes, what does the routine do; 

−  then, the tables of 
+  then, the tables of \(x\) and \(w\) are printed to verify the correct reading; 
then, the integral of power function is evaluated numerically, 
then, the integral of power function is evaluated numerically, 

−  +  \(\displaystyle_o^\infty x^m \exp(x) \mathrm d x\) 

−  for 
+  for \(m\) from 1 to 17. 
−  As the result is expected to be close to 
+  As the result is expected to be close to \(m!\), ration of the first column to \(m!\) is printed in the last column. 
−  while 
+  while \(m\!<\! 2M\!\!1\), the Gauss quadrature integrates "exactly", 
id est, with precision close to the maximal precision achievable with double variables; 
id est, with precision close to the maximal precision achievable with double variables; 

−  since 
+  since \(m\!=\!16\), the precision drops down, as it is supposed to do. 
It is not recommended to look inside the code: 
It is not recommended to look inside the code: 
Latest revision as of 18:44, 30 July 2019
GaussLaguerre 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)\,} \)
Contents
Generalized Gauss Laguerre
The scheme described in the preamble seems to be simplest case of the GaussLaguerre 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\). ^{[1]}.
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 GaussLaguerre rule for approximating
Integral ( a <= x < +oo ) xa^ALPHA exp(B*x(xa)) 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.9999999999999956e01 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 GaussLaguerre rule for approximating Integral ( a <= x < +oo ) xa^ALPHA exp(B*x(xa)) 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
GaussLaguerre quadrature seems to be efficient in the numerical calculations with approximation of atomic wave functions in spherical coordinates.
Also, the GaussLaguerre quadrature seems to be efficient for calculation of the Green function of the propagation go light in cylindrical coordinates.
References
 ↑ http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html John Burkardt. GEN_LAGUERRE_RULE. Generalized GaussLaguerre Quadrature Rules. 2010.02.23.
https://en.wikipedia.org/wiki/Gauss–Laguerre_quadrature