Gauss Laguerre integration
Consider
Equation 1
In the usual event that this integral decreases exponentially, it is possible to write it as a Gauss Laguerre polynomial. This enables us to re-write the integral as
Then a change of variable to 
makes this
Since the common forms of fun(r) decrease exponentially, it is reasonable to expand
Then
There is a very explicit sense in which this is a "best" fit to the function in this interval. The integral
Equation 2
And the L point sum of the appropriate weights times the value of f(x) and the zeroes of the L'th Laguerre polynomial is
The weight function appropriate for exponential integrals is
. The
following tables are from Stroud and Secrest[1]
|
|
|
|
X |
A |
|
0.07415878376 |
0.1767684749 |
|
0.3912686133 |
0.3004781436 |
|
0.9639573439 |
0.2675995470 |
|
1.796175582 |
0.1599133721 |
|
2.893651381 |
0.06824937997 |
|
4.264215539 |
0.02123930760 |
|
5.918141561 |
0.004841627351 |
|
7.868618915 |
0.0008049127473 |
|
10.13242371 |
0.9652472093E-4 |
|
12.73088146 |
0.8207305258E-5 |
|
15.69127833 |
0.4830566724E-6 |
|
19.04899320 |
0.1904991361E-7 |
|
22.85084976 |
0.4816684630E-9 |
|
27.16066932 |
0.7348258839E-11 |
|
32.06912225 |
0.6202275387E-13 |
|
37.71290580 |
0.2541430843E-15 |
|
44.31736279 |
0.4078861296E-18 |
|
52.31290245 |
0.1707750187E-21 |
|
62.80242315 |
0.6715064649E-26 |
|
|
|
|
X |
A |
|
0.07053988969 |
0.1687468018 |
|
0.3721268180 |
0.2912543620 |
|
0.9165821024 |
0.2666861028 |
|
1.707306531 |
0.1660024532 |
|
2.749199255 |
0.07482606466 |
|
4.048925313 |
0.02496441730 |
|
5.615174970 |
0.006202550844 |
|
7.459017453 |
0.001144962386 |
|
9.594392869 |
0.1557417730E-3 |
|
12.03880254 |
0.1540144086E-4 |
|
14.81429344 |
0.1086486366E-5 |
|
17.94889552 |
0.5330120909E-7 |
|
21.47878824 |
0.1757981179E-8 |
|
25.45170279 |
0.3725502402E-10 |
|
29.93255463 |
0.4767529251E-12 |
|
35.01343424 |
0.3372844243E-14 |
|
40.83305705 |
0.1155014339E-16 |
|
47.61999404 |
0.1539522140E-19 |
|
55.81079575 |
0.5286442725E-23 |
|
66.52441652 |
0.1656456612E-27 |
These tables in a form more appropriate for downloading are also in for\19pt.dat and for\20pts.dat. As coded in Fortran to check for my typos, see for\TGLAGU.FOR and for\GLAGU.FOR. There are two tables, on the grounds that if you can do it once, you can do it twice. Stroud goes up to N=68, and actually has 36 digit accuracy for A and x. If this seems profitable to you, contact me. The integral is
Note that if fun = 1 that the function at the last point is 7.7826x1028. This is about the size of the last weight. If the function is reasonably convergent by x=66, this naive approach may work.
It is not important to get the exponential right. One of the functions with the best convergence as a power series is the exponential function. To test this In the second part of my test for typos code above, I integrate
for various values of alpha.
|
Alpha |
Analytic |
20 pt |
20pt-19pt |
|
-1 |
1 |
1 |
10-10 |
|
-0.5 |
2 |
2 |
10-10 |
|
-0.2 |
5 |
4.999998 |
10-6 |
|
-0.1 |
10 |
9.9943 |
2x10-3 |
|
-0.05 |
20 |
19.507 |
0.105 |
|
-2 |
0.5 |
0.5 |
5x10-11 |
|
-4 |
0.25 |
0.25 |
5.7x10-9 |
|
-8 |
0.125 |
0.12495 |
2.9x10-5 |
|
-16 |
0.0625 |
0.05967 |
7x10-4 |
Details on how to attempt an inappropriate integration in
addition to consideration of 
à
Are in wsteve-GLintegration\wsteveGL.htm
In general there are 4 possibilities
a. I made a typo
I did check my code, but some of you are beginning to have doubts.
b. You made a typo
Test these by integrating 
=3,628,800. Note 66.5244165210 = 1.6975037x1018 -- barely tests the last number Also integrate
Since the last term to the 20th power = 2.88x1036 which will somewhat test the last term.
c. Both integrals give the same answer to 7-10 digits.
Home run! Use the 20 point value and assume that the error is on the order of the difference. I like statements like the results with 19 and 20 point Gauss quadrature agreed to 8 digits. Give Stroud and Secrest some credit.
d. The integrals differ in the third digit
Sorry! If a 38th degree polynomial does not represent the function, it probably cannot be done without some thought.
When this does not work as indicated by less than 4-digit agreement between the 19-point evaluation and the 20-point evaluation. Somewhat stronger methods are needed.
Variable change
The Gauss Laguerre points are not particularly magic. Any interval can be changed to the integral between 0 and 1. Let
Equation 3
So that
This makes the integral in equation 1 become
In general the only requirement on the function used to change variables is that it be positive definite. In general there is a bit of problem in finding r(y). In this case, however, equation 2 is
So that the integral becomes
This is ready for mid-point trapezoidal rule integration with a user-specified number of points.
I took the liberty of coding this in Fortran, so that I could draw some pictures. for\VARCHAN.FOR
H=1D0/M
AINT=0
DO I=1,M
YI=H*(I-.5D0)
RI=-LOG(1-YI)/ALPHA
FT=FUN(RI)/(1-YI)
AINT=AINT+FT
WRITE(2,*)YI,FT
WRITE(3,*)YI,RI
ENDDO
AINT=H*AINT/ALPHA
ANAL=1/BETA
The function passed to this routine is
FUNCTION FUN(R)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/PASS/BETA
FUN=EXP(-BETA*R)
Note that alpha and beta can be different. To no one's surprise
INPUT BETA, ALPHA, M
1 1 10
NUMERICAL INTEGRAL IS 1.000000000000000
ANALYTICAL INTEGRAL IS 1.000000000000000
Figure 1 Integrand versus y for alpha = beta = 1
Figure 2 R(y) versus y for alpha = beta =1
INPUT BETA, ALPHA, M
.1 1 10
NUMERICAL INTEGRAL IS 3.513015505058240
ANALYTICAL INTEGRAL IS 10.000000000000000
In the case where 
and the integration assumes e-r.
The ten-point value is off by almost a factor of 3.
Figure 3 Integrand for exp(-0.1 r) integrated with exp( - r).
Using 20 points helps only a little
.1 1 20
NUMERICAL INTEGRAL IS 3.947173835135680
ANALYTICAL INTEGRAL IS 10.000000000000000
.1 1 5000
NUMERICAL INTEGRAL IS 6.515255637935450
ANALYTICAL INTEGRAL IS 10.000000000000000
Recall that Gauss Laguerre for 20 points was nearly exact.
Figure 4 Integrand for trap rule with exp(-0.1 r) integrated using exp( -r)
This case will extrapolate as h rather than h2 since we are essentially integrating a singularity at y = 1.
Figure 5 Integrand of exp( - 0.1 r) using exp(-r) near y = 1.
This should immediately make you think of
the BLI assignment
On the other side integrating exp( -10 r) gives
Figure 6 Integration of exp(-10r) with exp(-r)
And the same works for the exp(-0.1 r) if I use alpha = 0.01.
.1 .01 128
NUMERICAL INTEGRAL IS 9.997711409803650
ANALYTICAL INTEGRAL IS 10.000000000000000
Figure 7 integrand versus y for exp(-0.1r) evaluated using exp(-0.01r)
The difference of course is the 1/(1-y) which makes singularities as r goes to infinity much worse than those as r goes to zero.
Two regions
Sometimes you just have to bite the bullet and admit that you have to look closely at the function in the range between 0 and R. Do not just stop at R. Write equation 1 as
The first integral is to be looked at in detail using mid-point trap rule or findfun or whatever. It is over a finite range and can be extrapolated to zero step size. It is I2 that needs to be examined here Change variables to y=r-R
Then change to x = alpha y
And of course using Gauss Laguerre quadrature
Two things differ from the last lecture. First there is no attempt to evaluate the entire integral in I2. Mostly it is present to ensure that we have gone out far enough with R.
This means that the 2 and 3 point Gauss Laguerre polynomials are appropriate.
|
|
|
|
X |
A |
|
0.5857864376 |
0.8535533906 |
|
3.414213562 |
0.1464466094 |
=
|
|
|
|
X |
A |
|
0.4157745568 |
0.7110930099 |
|
2.294280360 |
0.2785177336 |
|
6.289945083 |
0.01038925650 |
Second the spacing needed to determine the function "h" has been determined in evaluating the first region. This allows us to say
and thereby determine a reasonable guess for alpha. The quotation marks of course mean that the evaluation of region 1 may well have involved a much smaller h than is necessary at the end range of the function.
Make the Gauss-Laguerre integrator work for some simple functions.
a.
b. 