MA332 lecture notes, Fall 1998
Week 3, Thursday

Reading: F&B Chapter 4

Romberg integration = trapezoid rule (with a trick to reuse data)
		    + Richardson's extrapolation


Richardson's extrap says that if you have a series of estimates that
are converging on and answer, and you know _how_ it's converging,
you might be able to guess where it's going.

For example:  x0 = 1;
	      x1 = 2;

I don't know where this series is going, but I do know that
x1's error is expected to be half of x0's error.  What do I guess?

     x1 + 2 * (x1 - x0) = 3

More general example:

M = exact answer
N(h) = estimate based on interval size h

Assume that we use a numerical integration technique with the
following behavior:

M = N(h) + K1 h + K2 h^2 + ...

For small h, the dominant term of the error is K1 h, so we say
that the error is O(h).

If we halve the size of the interval, we get

M = N(h/2) + K1 h/2 + K2 h^2/4 + ...

If the O(h) term is dominant, then the error has been approximately
halved.  Thus, if we take

	 x0 = N(h)
	 x1 = N(h/2)

Then the extrapolation is

	 x1 + 2 (x1-x0), or what the book calls N2 (h).

Since the O(h) terms cancel, we get

M = N2 (h) - K2/2 h^2 - K3 h^3

Again, for small h the h^2 term dominates the error, so we say that
the error is O(h^2).


1) We can repeat this process to cancel the O(h^2) terms, and so on.

   x0 = M + K h^2
   x1 = M + K H^2 / 4

   M = x1 + (x1 - x0) / 3

   In general

   M = x1 + (x1 - x0) / (2^j - 1)

   where j>=1 is the order of the error term you are cancelling,
   and assuming that you are halving h each time.

2) The limit is the point where the subtraction of adjacent estimates
   starts losing accuracy.


How do we apply this to integration?

1) Make successive estimates using the trapezoid rule, taking
   advantage of the fact that to go from n intervals to 2n intervals
   you only have to evaluate n new points (not true for Midpoint
   or Simpson).

2) Find an estimate of the error terms that is a power of h multipled
   by a constant (that does not vary as we change the number (size) of
   subintervals).

3) Use Richardson extrapolation.

From page 119 we know that the error term for Composite Trapezoid
Rule is

     (b-a) h^2 / 12 * f''(xi),  where xi is in the interval (a,b)

Well, that's no good because there is no reason to expect xi to
be the same for different values of h.

On page 130, they claim (unproved) that if f(x) and all its
derivatives are continuous, then the error term can be written
as a summation of

   Ki * h^2i     for i = 1 .. infinity

In other words, the dominant term of the error is  K1 * h^2, which
is exactly the form we need for Richardson's extrapolation.

Furthermore, since the error terms only have even powers, once we
cancel the O(h^2) term, the next dominant term is O(h^4)!!

Example:

Use Maple's trapezoid function to generate successive
estimates of the integral of sqrt (1+ x^2), the arrange the
Romberg estimates into a triangle.

1) For the most part the estimates converge as expected, although
there is one example of overshooting the mark.

Break


Gaussian Quadrature
-------------------

Choose any polynomial of order 3.  There are 4 parameters.

Evaluate it at sqrt(3)/3 and -sqrt(3)/3 and add 'em up.

That sum is the exact integral of the polynomial over (-1,1)!

Proof:

Assume that f(x) = a0 + a1*x + a2*x^2 + a3*x^3

Then int(f(x)) = a0 * int(1) + a1 * int(x) + a2 * int(x^2) + a3 * int(x^3)

If we assume that int (f(x)) = c1 f(x1) + c2 f(x2)
= a0 (c1+c2) + a1 (c1 x1+c2 x2) + a2 (c1 x1^2 + c2 x2^2)...

Since we don't know the values of the a's, we have to match up
all the terms:

int(1) = c1 + c2 = 2
int(x) = c1 x1 + c2 x2 = 0
int(x^2) = c1 x1^2 + c2 x2^2 = 2/3
int(x^3) = c1 x1^3 + c2 x2^3 = 0

Solve that nonlinear 4x4 system!

c1 = c2 = 0    x1 = sqrt(3)/3   x2 = -sqrt(3)/3

So, how can this possibly work?

Using Maple, we can see that a1 and a3 have no effect on the integral
because odd functions integrate to zero on (-1,1).  Thus there
are really only two parameters (not four), and we only need 2 points
to determine them.

1) What about intervals that don't happen to be (-1,1)?

Just use a change of variables:  x = (2x'-a-b) / (b-a)

x' = a  ->  x = -1
x' = b  ->  x = 1

x1' = ((b-a) x1 + a + b) / 2
x2' = ((b-a) x2 + a + b) / 2

dx' = (b-a)/2 dx

int(f(x), a..b) = (b-a)/2 (f(x1') + f(x2'))

2) Can we do the same thing for higher-order polynomials?

You bet, see the table on page 137.

3) Why should we expect a method that is exact for 3rd order polynomials
   to be accurate in general.

   Taylor's theorem?  Smoothness?  Because we hope so?