Objective

In this lecture, I intend to demonstrate that mid-point trap rule is the "best" rule for the innermost phi integral of simple three-d integration. In addition it is also the "best" method for functions that are known to start and end with zero derivatives.

Mid point trap rule -- The simplest numerical integration method

The integral from beg to end can be written as

Expand the function about the values x_i

So that

In each of the N integrals let t=x-x_i

When integrated from -h/2 to h/2, all of the odd terms integrate to zero and the rest are only a little harder so that by expanding about the mid point of each region an integral in general can be written as
Note that the end points are never evaluated. Simplifying a bit

The first term is the mid-point trapezoidal rule estimate of the integral.

The code accurate to order h² for mid point trap rule is

h=(end-beg)/N this lower case combination of

ai=0 Fortran,basic, and C is not a new language

for I=1 to N BOB but merely my attempt to indicate

x=i x h - h/2 code that you may implement

ai=ai+fun(x) in any language

enddo

ai=ai x h

Write this estimate for the integral of f''

Then notice that the integral on the left is easy to evaluate

Solve for the sum on f''

Then substitute into the equation for the original integral to find

And indeed the last term can be treated in the same way to yield

Equation 1

Do not fear! We almost never evaluate the derivatives. What we do is Richardson extrapolation. In an easy version of this note that the integral is quadratic in h. We plot the value of the trap rule integral estimate versus h².

For small enough values of h² the integral estimate is linear as a function of h², for larger values, it is quadratic. The quadratic values can be extrapolated to zero using Lagrange interpolation.

With

Even easier, the last two points can be extrapolated to zero using

Press

The integration techniques here are discussed in Press in Chapter 4, section 4.3 in particular uses POLINT, which is very similar to our uneven Lagrange method for extrapolating. Section 4.4 is similar to our mid-point methods. Press relies a bit too heavily on a single difference to estimate the errors in my opinion and his stopping when satisfied is to my mind always suspect. Press also goes to too high an order in this code. Ease of implementation does not justify high orders. Error estimates tend to oscillate. Also see section on periodic functions below.

A third estimate

A totally incorrect error estimate can result from an oscillating function that gets hit only at the peaks. This can be reduced dramatically by adding one more estimator corresponding to the line between points 1 and 3 in the picture.

Each integral estimate point has an error of order h². Each two-point extrapolation has an error of order h⁴.

The three-point extrapolation has an error of order h⁶.

Assume that estimates are made with 100, 200, and 400 points. The integral estimate makes an error of order C x (0.01)² = C x 10^-4. The two-point estimator uses the second derivative, which is large by a factor of 2, but in an expansion that divides by a factor of 2. Thus the two point estimator makes an error of order C x 10^-8, while the three point estimator makes an error of order C x 10^-12. Press's method goes to even higher orders, but at this point, the entire method may have problems. A mid-point trap rule for complex functions using N, 3N, 9N and 27N points which enables the same points to be re-used complete with Lagrange extrapolation to the integral estimate is in for\mptrap.for

End point trap rule

The use of the end points complicates the evaluation in two ways. First the end values need to be divided by 2, secondly they need to be calculated and this may involve analytically expanding functions such as (1-cos(x))/x². There is an advantage; placing an equal number of points in the center of each region can halve the spacing h. With mid-point trap rule it is necessary to place twice as many points in a more complex scheme in order to take advantage of the first set of points to cut the size by a factor of 3.

The code accurate to order h² for end point trap rule is

h=(end-beg)/N

ai=(fun(beg)+fun(end))/2

for I=1 to N-1

x=beg+i x h

ai=ai+fun(x)

enddo

ai=ai ´ h

Assignment

Write a code to evaluate an integral of our familiar Lorentzian -- Note that the analytical integral is
In G/R p.60

Periodic Functions

Frequently the integral of interest is of the form

For today the integral of interest is

The integrand has the property that

The first error estimate in equation 1 is

So is the second error estimate

As are all higher order estimates.

Trap rule error in periodic functions

A periodic function can always be expanded as

Note that c_j is complex

Note that

Since the first exponential is always 1. The integral from 0 to P of fun is

That is the actual value of the integral is given by the c₀ term alone. All other terms in the expansion integrate to zero.

The mid point trap rule approximation to this integral is

note interchange of summation orders
where h=P/N and x_k=h/2+k h so that

Let z=exp(2pij/N),so the last term can be written as a sum of z^k and we can use the familiar relation for z¹1

So that

which separating out the j=0, and j=mN terms yields

And we note that exp(2pij)=1 for all j so that the final term is always zero.

The error in the trap rule integration of fun is just given by c_N the N th coefficient in a complex expansion of fun. This implies that the first 2N coefficients, counting the real and imaginary parts of c_j as two coefficients, are exactly accounted for in a simple midpoint trap rule of a periodic function.

Note that the terms that enter are not uniform. If for some accidental reason c₁₀ = 0, 10 point integration will seem exact up to c₂₀. Changing N from 10 to 20 will give a false accuracy reading. The correct change is from 10 to 11 or 13, Do not simply double the points. --- Press pay attention.

Periodic extension

Note that any function which is zero at 0 and 0 at P where P is a large number can be treated as though it were periodic. In particular the two-body correlation function for a classical system is defined to be

Changing variables in the r2 integration can make this into

Then finally

For a Lennard Jones potential the wave function contains a exp(-1/r⁵). For a Coulomb potential it is a less steep, but still the two-body correlation function

Contains a short range part that goes as exp(-b/r⁵)

And a long-range part that eventually becomes 1+exp(-ar). Thus the derivatives at both ends of the integration region are zero making the mid point trap rule extremely accurate for expectation values which can be reduced to integrals such as