Natural Integration range
Consider
let x=(a+b)/2+t(b-a)/2; dx=dt(b-a)/2 so that
Where
Orthogonal polynomials
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
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
an 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. For now leave it equal to
1.
Gauss Quadrature
The function of interest is
The
sum in principle includes an infinite number of terms. The integral is
Note
that the w is absorbed into the A. The
object is to find A’s and xi’s such that the error appears in as
high a cj as possible. For
the first N terms
As
long as the xi are at distinct locations it
is a simple matter of solving the N simultaneous equations for the N values of
Ai. Then the N+1’st equation
is simply
solve
this one by letting xi be the N zeroes of PN(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
This is also solved by this same set of xi’s. Continuing with this the last 2N+1 th equation however becomes
which for the set of xi’s at the zeroes of
PN(x) is as wrong as it can
be. The integral is exact for 2N terms,
then wrong for the 2N+1 th term. This is the same as found for periodic
functions in the last lecture. The
conclusion is that
where the xi's are the zeroes of PN(x), and the Ai's
in general need to be solved for. Codes
for finding these are in Press[1]
and should be used on a second pass.
|
xi |
AI Gauss Legendre |
|
|
±0.93246951420315202781 |
0.17132449237917034504 |
|
|
±0.66120938646626451366 |
0.36076157304813860756 |
|
|
±0.23861918608319690863 |
0.46791393457269104738 |
|
Tables with x, A values for N from 2 to 10 are in GaussLeg.htm
|
|
Gauss Legendre[2] |
|
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.
Assignment 1
Typos are the biggest problem
with these types of tables. Integrate a
constant between -1 and 1 to find that it gives 2, and then integrate t2
to find that it gives 2/3. Then integrate
tN for N = 10,11,12,13 and so using the 6
term Gauss Legendre formula given above.
Repeat with the 20 and 21 term form.
Bli Integration 2 The Lagrange Polynomial
This is POLYL in for\lagrange1.for
Where
By virtue of the fact that one term is always left out of
the product that covers N values in which x always appears once, it is easy to
see that Lm is a polynomial of degree N-1 and that sums of the Lm values
times constants will also be of degree N-1.
The value of Lk(xm¹k) where xm¹k is
any of the data points other than the k’th one will have a zero in the product
for the j=m term 
so that Lk(xm¹k)=0. The value of Lk(xm=k)
is quite different. Each term in the
product is of the form 
and by construction
the dangerous term with xj=xk has been left out so that Lk(xm=k)=1.
The Lagrange interpolation implicitly
uses the function values to construct the N-1 derivatives in the region with an
error equal to the N the derivative anywhere in the region and then uses these
to form the polynomial approximating function.
This means that the error term in the Lagrange polynomial approximation
for an internal value of x is 
where x is any point between the first and last point and w is
on the order of the size of the region.
For values of x outside the region w becomes on the order of the distance
between x and the most distant point in the interpolation and x is any point in that same region.
Figure 1 shows point used by PolyL and definition of mbeg,
nbeg
The data shown is as returned by BLI. The x axis contains values of xdat, while the y axis contains values of fdat. The 1, 2, …, 5 refers to the value of j in . In the lagrange code these x values are referred to as J+MBEG. The integration of interest is from x=nbeg to x=nbeg+1.
Where
The value nbeg is returned by locate, when it is called as CALL LOCATE(X,XDAT,NDAT,NBEG).
SUBROUTINE UELAG(NL,X,MBEG,NSKIP,ALAG,XDAT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ALAG(*),XDAT(*)
DO M=1,NL
IF(M+MBEG.EQ.NSKIP)THEN
ç NSKIP is for error
analysis.
ALAG(M)=0
ELSE
ALAG(M)=1
DO J=1,NL
IF(J.NE.M.AND.J+MBEG.NE.NSKIP)THEN
ALAG(M)=ALAG(M)*(X-XDAT(J+MBEG))/(XDAT(M+MBEG)-
2
XDAT(J+MBEG))
ENDIF
ENDDO
ENDIF
ENDDO
RETURN
END
FUNCTION POLYL(NL,X)
C NL
is the number of data points in the interpolation
C X
is the coordinate of the output value
C XDAT(NBEG) is the data point just below x
C XDAT(NSKIP) is not included in the interpolation
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NDAT=55,NLAG=12)
DIMENSION ALAG(NLAG),DLAG(NLAG),DDLAG(NLAG)
DIMENSION XDAT(NDAT),FDAT(NDAT)
SAVE NC,XDAT,FDAT
DATA NC/0/
IF(NC.EQ.0)THEN
OPEN(1,FILE='H2.TXT')
DO I=1,NDAT
READ(1,*)XDAT(I),FDAT(I)
ENDDO
CLOSE(1)
ENDIF
IF(NL.GT.NLAG)THEN
PRINT*,' EXCEEDED ALAG DIMENSION IN
POLYL'
READ(*,*)ITEST
STOP
ENDIF
CALL LOCATE(X,XDAT,NDAT,NBEG)
MBEG=MAX0(0,NBEG-NL/2)
MBEG=MIN0(NDAT-NL,MBEG)
NSKIP=0
CALL UELAG(NL,X,MBEG,NSKIP,ALAG,XDAT)
POLYL=0
DO I=1,NL
POLYL=POLYL+ALAG(I)*FDAT(I+MBEG)
ENDDO
RETURN
END
Assignment 2
Consider again the BLI integrals between 0 and 1
- ftest = x4
Integral should be 0.886226925455
Integral 0 to 1 should be 1.64493406685
error integral.
Use the BLI to find
the best set of points, xdat, and fdat.
Then use Gauss Quadrature to find the function Gdat
Continue with Locate
and Poly to write a function such that
using
Lagrange interpolation.
What is the order of
the error? How would you recommend
estimating the error?