Difference between revisions of "Gauss-Laguerre quadrature"

From TORI
Jump to navigation Jump to search
 
m (Text replacement - "\$([^\$]+)\$" to "\\(\1\\)")
 
Line 1: Line 1:
 
[[Gauss-Laguerre quadrature]] is approximation of integral of exponentially decaying function with the finite sum defined as follows:
 
[[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$
+
\(\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,
+
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)$
+
\(\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)
+
\(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\) 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]],
+
\(x_n=\mathrm{LaguerreLzero}[M,n]\) is \(n\)th zero of \(M\)th polynomial [[LaguerreL]],
   
$\displaystyle
+
\(\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 [[Gauss-Laguerre quadrature]].
 
The scheme described in the preamble seems to be simplest case of the [[Gauss-Laguerre 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
+
\(\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 $\alpha$ and $b$ are constants.
+
where \(\alpha\) and \(b\) are constants.
   
The formula has the same form, the approximation $S$ is written as
+
The formula has the same form, the approximation \(S\) is written as
   
$\displaystyle S=\sum_{n=1}^M \, w_n \, f(x_n)$
+
\(\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 $x_n$ and \weights $w_n$ for $0<n\le M$,
+
coordiantes \(x_n\) and \weights \(w_n\) for \(0<n\le M\),
at given $\alpha$, $b$ and $M$.
+
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 $M$
+
number of points in the rule, denoted above with \(M\)
   
Order of the corresponding lagiuerr, denoted above with $\alpha$
+
Order of the corresponding lagiuerr, denoted above with \(\alpha\)
   
$A$ is the left endpoint (typically 0). (above it is explicitly set to zero)
+
\(A\) is the left endpoint (typically 0). (above it is explicitly set to zero)
   
$B$ is the exponential scale factor, typically 1.
+
\(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 $\alpha\!=\!0$ ===
+
===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 $x$ and $w$ are printed to verify the correct reading;
+
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$
+
\(\displaystyle_o^\infty x^m \exp(-x) \mathrm d x\)
   
for $m$ from 1 to 17.
+
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.
+
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",
+
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 $m\!=\!16$, the precision drops down, as it is supposed to do.
+
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

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\). [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 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

  1. 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.

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

Keywords

LaguerreL,