Objective

In this lecture, mid-point trap rule will be shown to be the "best" rule for periodic functions. The actual error for this rule is derived in the section #Trap rule error in periodic functions. This is actually a bit of a simplification of the integration method discussed in BLI_integration.htm.

Mid point trap rule -- The simplest numerical integration method

The integral from beg to end can be written as
[1]

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
Text Box: The code accurate to order h2 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
for I=1 to N language BOB but merely my
x=i x h - h/2 attempt to indicate code that
ai=ai+fun(x) you may implement in any
enddo language
ai=ai x h

Note that the end points are never evaluated. Simplifying a bit

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

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

Extrapolate to zero in h². 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^2, the integral estimate is linear as a function of h² so that the values recorded from the simple trap rule sum as Int₁, Int₂ and Int₃ can be written as

Where A is the value of the integral and B is a slope that we don’t even care about.

Richardson extrapolation. is simply solving the two equations. Multiply the first by h₂² and the second by h₁² and subtract to find

In particular if h₁ = 2h₂ , this becomes which looks surprisingly like Simpson’s rule. We however have computers and can evaluate this for (h₃,h₁), (h₂,h₁) and (h₃,h₂). The average is a good estimate of the value, while the sum of the squares of the differences is a fair estimate of the error.

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. for/mptrap.for

End point trap rule See also BLI_integration.htm

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 in the box to the left[2].

Periodic Functions

A periodic function is of the form

If the function is integrated from 0 to T, the first error estimate in 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 and that

Since the first exponential is always 1. The integral from 0 to T of f 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.

Let so that the mid point trap rule approximation to this integral is

note interchange of summation orders. Let , so that

The last term can be written as a sum of z^k and we can use the familiar relation for z¹1

So that

Note that so that the the numerator in is always zero. The denominator is also zero for j=mN for which . In this case the last sum in is N so that

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.

Assignment

Evaluate Using various integration points ranging from N=m/2 to N=4m. Plot the integration error as a function of 1/N². If the method “breaks” plot the details where it breaks.