Introduction

                Last time we considered

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

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 2

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 findfun work of a few weeks ago.

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.