Lessons from last time

The 3-d integral of interest is

This should not be evaluated using the variable mu, but should be evaluated using theta so that

Mid-point trap rule applied to the phi integration is Gauss quadrature for a periodic function. Mid-point trap rule applied to the theta integration is Gauss Tchebyshev quadrature for a function defined over a fixed region.

Gauss quadrature means that if the function is expanded in the appropriate set of orthogonal polynomials as

The integral using N points is exact for P_j with j <= 2N

For the phi integration the P's are not polynomials but are rather cosines and sines. For the theta integration the polynomials are Tchebyshev polynomials. To be specific the phi integration is

The theta integration is

The theta dependence of M

The circle swept out as a line, depends on the value of . For r along the z axis, the same point is sampled over and over again as phi changes. It is thus quite reasonable to make M a function of phi. That is to define

This approximately doubles the effective number of phi evaluations in the important region. Combined these yield

Equation 1

The r dependence of N and M_max.

The third integral in the integration over volume is

Equation 2

It differs from the first two in that it is of infinite range and that each value of f(r) involves M x N operations. since it is usually a an integral over phi followed by an integration over theta of 20 to 100 points.

Figure 1 Region with only a few centers.

In the case of a few atoms clustered near the origin as shown in figure 1, the number of points in the angular integration can be left roughly constant as r is increased, since the information needed to accurately estimate these integrals does not increase as r increases. Thus the values of N and M_max in equation 1 can be left constant as a function of r. For large r . In this equation the value of alpha is "essentially" determined by the electronic states of the electrons making up the solid regions.

Figure 2 extended region filling all of space.

In the case of a region filling all of space as in figure 2, the values of N and M_max must increase with increasing r. The mid point trap rule integrations in N and M are "essentially" exact for about 3 to 5 points across each region. If the natural units of Angstroms or Bohr radii have been chosen, the atoms are on the order of 2-5 apart. Thus N and M_max need to be approximately .

The value for large r is given by as before, but now alpha is the Debye length or simply a very small integrating factor inserted to make the integral convergent. Note that for r = 100 a_B. That N = M_max= 300 so that approximately 45000 evaluations are needed for each value of fun(r). For r = 1000 a_B , 4.5 x 10⁶ evaluations are needed. While even this last number is reasonable with today's computers, there is a premium on using few r values and rapidly estimating the overall behavior of fun(r).

Gauss Laguerre integration

Since the common forms of fun(r) decrease exponentially, it is reasonable to expand

Then

There is a very explicit sense, to be detailed in a future lecture, in which this is a "best" fit to the function in this interval. The integral

Equation 3

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. In C writing the code will be good for you. 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.7826x10²⁸. 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

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.52441652¹⁰ = 1.6975037x10¹⁸ -- barely tests the last number Also integrate

Since the last term to the 20th power = 2.88x10³⁶ 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.

Some thought

The first thing to do is to write r and f(r) from the integrator, to plot these and see where the points are relative to the value of the function. Then note that a change of variables to makes .