Objective

In this lecture, I intend to demonstrate that Gauss Tchebyshev Integration the "best" rule for the second mu integral of simple three-d integration. It is "best" because the weights are equal and because the values of x_i are easily found and varied. It is better than mid-point trapezoidal rule for integration regions in which the derivatives of the end points differ. It superiority to mid-point trap is primarily due to the fact that the points are more closely spaced at the ends than in the mid-region. It shares with mid-point trap the fact that the integrand never needs to be evaluated at either end.

An expansion in orthogonal polynomials

Last time we found that a periodic function,

in the interval (-1,1) has an integral given by . We also found that the error in a mid-point trap rule estimate of this function is . This says that the integration error is approximately equal to 2c_N which is the 2N+1 'th coefficient in a Fourier expansion of the function. For N on the order of 100 or so this allows fun to be rather well represented by this expansion. The coefficients can be found by utilizing the orthogonality over this interval to arrive at

The Fourier expansion is limited to periodic functions, but any function can be expanded as

If the P's are orthogonal in the sense that

Then

Equation 1

The function w(x) is at our disposal. It can be made exponential so that the interval can go from 0 to infinity in the case of Laguerre polynomials. It can be made a exp(-alpha x*x) so that the interval runs from -infinity to plus infinity in the case of Hermite polynomials. It can be made 1 in the case of Legendre polynomials.

Tchebyshev Polynomials

The Tchebyshev polynomial has the property that it is always less than or equal to 1.

so that

This means that the in the expansion.

the error in approximating the function is approximately equal the last term dropped. A good rule of thumb is to assume that it is less than the last term held. There is a sub industry called Chebyshev economization in which a polynomial expansion is made about x_i which is more accurate than needed near x_i in order to maximize its range. It is then re-expanded as a the above. The coefficients larger than those needed for the desired accuarcy are dropped. Then the series is rewritten as a lower order polynomial.

Evaluating the integrals for the coefficients (and finding w(x))

The orthogonality condition is

Change variables to

So that the condition becomes

Thus define

to find

This yields

Note that the expansion coefficients in the series are given by equation 1 as

which seems to be a somewhat nasty integral, especially if one wants to use an end point method to evaluate it. As above change variables to so that

This simplifies the integral for the coefficients to

Gauss Quadrature

The function of interest is

The sum in principle includes an infinite number of terms. The integral is

and the obgect is to find A’s and x_i’s such that the error appears in as high a c_j as possible. For the first N terms

as long as the x_i are at distinct locations it is a simple matter of solving the N simultaneous equations for the N values of A_i. Then the N+1’st equation is simply

solve this one by letting x_i be the N zeroes of P_N(x), which for a positive definite w(x) can be shown to all lie between 0 and 1. Now look at the

N+2 nd equation

and note that by adding just the correct amounts of the first N equations, this can be factored into

Which is also solved by this same set of x_i’s. Continuing with this the last 2N+1 th equation however becomes

which for the set of x_i’s at the zeroes of P_N(x) is as wrong as it can be. This is the same as found for periodic functions in the last lecture. The conclusion is that

where the x_i's are the zeroes of P_N(x), and the A_i's in general need to be solved for. Codes for finding these are in Press and should be used on a second pass.

	Gauss Legendre
X=+-value	A, N=20 same for +- values
0.9931285992	0.01761400714
0.9639719273	0.04060142980
0.9122344283	0.06267204833
0.8391169718	0.08327674158
0.7463319065	0.1019301198
0.6360536807	0.1181945320
0.5108670020	0.1316886384
0.3737060887	0.1420961093
0.2277858511	0.1491729865
0.07652652113	0.1527533871

The 21 point one counts the central point only once

	Gauss Legendre
X=+-value	A, N=20 same for +- values
0.9937521706	0.01601722826
0.9672268386	0.03695378977
0.9200993342	0.05713442542
0.8533633646	0.07610011363
0.7684399635	0.09344442346
0.6671388041	0.1087972992
0.5516188359	0.1218314161
0.4243421202	0.1322689386
0.2880213168	0.1398873948
0.1455618542	0.1445244040
0.00000	0.1460811336

The reason that there are two tables, one with 20 and one with 21 terms, is that integral should always be evaluated with both. When the values are the same, this says that the 41 and 43 terms in the expansion of the function are both effectively zero and that the answer is correct. When the answers are different, neither should be trusted. Note that the weights at the ends are much smaller than those in the middle and that the points at the ends are much closer spaced than the ones in the middle.

Gauss Tchebyshev Quadrature

The desired integral is

The zeroes of the N'th Tchebyshev polynomial occur at the points where

This requires

Recall of course that we need to have x_i between -1 and 1

The A's turn out to be constant = as can easily be seen by letting fun = w so that

Let so that

In the more general case

Equation 2

Wait a minute.

Equation 2 is just a mid-point trap rule evaluation of

Equation 3

This can be arrived at beginning with

Equation 4

By changing variables to

h=3.1415927/50 open(1,file='test.pts') do i=1,50 y=h*(i-.5) x=cos(y) f=x**4 write(1,*)x,f enddo close(1)

The above analysis shows that a mid-point trap rule of Equation 3 is much more accurate, than a mid point trap rule of equation 4. I used the simple code to the left to generate the data plotted in figure 1.

Figure 1 The points calculated in a Chebyshec evaluation of fun(x) = x⁴.

The point is that sin(y) is the weight, and cos(y) for y spaced as a mid-point trap rule is the set of y_i's needed for a Gauss Tchebyshev integration. Equation 3 is more accurate than equation 4, because the end points which have data on only one side are interpolated with a much closer spacing.

Assignment

Let fun=sin((p/2)x*x)

beg=0

end=

That is Calculate the integral to 6+ digit accuracy. [Abramowitz and Stegun give 0.567822]

This is a maximum of the integral. This is a Fresnel integral which arises in diffraction problems. It is also a chirp which arises when two stars collide.

To get the most out of this assignment, use mid-point trap rule and Gauss Tchebyshev with a varynig number of points.

Try to extrapotate to an infinite number of points. Also use 20 and 21 point Gauss quadrature. -- be sure to check for type's both yours and mine.