Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

We used Richardson extrapolation to improve our approximations of integrals. This section shows how we can use Richardson extrapolation to improve the approximations found using the centred divided-difference formula.

Background

Useful background for this topic includes:

5. Interpolation
13.1 Centred Divided-Difference Formulae

References

Bradie, Section 6.3, Richardson Extrapolation, p.447.
Mathews, Section 6.1, Richardson's Extrapolation, p.330.
Chapra, Section 23.2, Richardson Extrapolation, p.635.

Theory

One significant problem with approximating a derivative is subtractive cancellation. In the formula (f(x + h) − f(x − h))/(2 h), for very small values of h, the two function evaluations f(x + h) and f(x − h) will be approximately equal, and therefore subtractive cancellation will occur. Therefore, it would be undesirable, and dangerous, to use ever smaller values of h.

Instead of using a smaller value of h, suppose we are attempting to approximate an exact value e with an approximation a(h). In this case, e is the derivative f⁽¹⁾(x) and the approximation is a(h) = (f(x + h) − f(x − h))/(2 h). Suppose now that the error of the approximation is defined by a Taylor series of the form:

e = a(h) + K hⁿ + o(hⁿ)

Consider now the approximation using h/2:

e = a(h/2) + K (h/2)ⁿ + o((h/2)ⁿ)
= a(h/2) + K/2ⁿ hⁿ + o(hⁿ)

Multiplying this second expression by 2ⁿ and subtracting off the first equation yeilds

2ⁿe − e = 2ⁿa(h/2) − a(h) + K/2ⁿ hⁿ − K hⁿ + o(hⁿ)

Note that the hⁿ term cancels and we are left with

(2ⁿ − 1)e = 2ⁿa(h/2) − a(h) + o(hⁿ)

If we look at the full Taylor series for the centred divided-difference formula, we note that the error terms are of the form K_nhⁿ:

We can write this as:

where K₁ = −1/6 f⁽³⁾(x)h², etc.

This is not the case if you use the backward divided-difference formula, which is of the form:

To convince yourself of this, you may look at this Maple worksheet which demonstrates this fact.

For the centred divided-difference formula, this is identical to the pattern for the composite trapezoidal rule, and therefore, we can use Richardson extrapolation to get a better answer. Note, for integration, applying Richardson extrapolation to the results of the composite trapezoidal rule is termed Romberg integration. There is no special name for Richardson extrapolation as applied to derivative formula.

HOWTO

Problem

Approximate the derivative of a univariate function f(x) at a point x₀. We will assume that we can calculate f(x) for arbitrary values of x.

Assumptions

We will assume that the function is sufficiently differentiable.

Tools

We will use interpolation and extrapolation.

Initial Values

We will start with an initial value h.

Process

For i = 0, 1, 2, let R_{i, 0} be the centred divided- difference formula with a step of h/2ⁱ. Then, for j = 1, 2, ..., i, calculate

Halting Conditions

There are two conditions which may cause the iteration process to halt:

We halt if |R_{i, i} − R_{i − 1, i − 1}| < ε_step, or
If we have iterated some maximum number of times, say N, and have not met Condition 1, we halt and indicate that a solution was not found.

If we halt due to Condition 1, we state that R_{i, i} is our approximation to the derivative.

If we halt due to Condition 2, we may state that a solution may not exist.

Examples

Example 1

Approximate the derivative of the function f(x) = e^-x sin(x) at the point x = 1.0 using the centred divided-difference formula and Richardson extrapolation starting with h = 0.5 and continuing using ε_step = 0.0001.

We calculate the derivative using h = 0.5 and h = 0.25:

-0.06821507210
-0.1001894112

Neither approximation is very close to the correct answer of -0.1107937653, however, one step of Richardson extrapolation yields:

-0.06821507210
-0.1001894112 -0.1108475242

Next, we approximate the derivative using h = 0.125, that is, -0.1081453352, and we fill in the third row of the table:

-0.06821507210
-0.1001894112 -0.1108475242
-0.1081453352 -0.1107973099 -0.1107939623

The difference |-0.1108475242 − (-0.1107939623)| = 0.000053 < ε_step, and thus -0.1107939623 is our approximation.

The relative error of our approximation is 1.8 × 10^-5.

Example 2

Repeat Exercise 1, but use the backward divided-difference formula.

Again, using h = 0.5 and h = 0.25 and applying one step of Richardson extrapolation, we get:

-0.2344655259
-0.1369349794 -0.1044247972

We cannot continue because the error for the backward-difference formula is 1/4 f⁽⁴⁾(ξ)h³, and therefore we cannot resort to using more than one step of Richardson extrapolation. However, we can still generate a table using only one step:

-0.2344655259
-0.1369349794 -0.1044247972
-0.1167038388 -0.1099601253
-0.1121938272 -0.1106904900
-0.1111341968 -0.1107809867

The difference between the last two, |R_4,1 − R_3,1| = 0.000090 < ε_step, so we are finished and we use -0.1107809867 as our approximation.

The relative error of our last answer is 1.2 × 10^-4.

Exercise 3

Use Richardson extrapolation to help approximate the 2nd derivative of the function given in Exercise 1, using the same values.

Using the formula for the second derivative, we calculate:

-0.4230489884 -0.4039640736 -0.3976024353 -0.3991434368 -0.3975365579 -0.3975321661

The first column is (f(x+h) - 2 f(x) + f(x-h))/h² for h = 0.5, 0.25, and 0.125. The next two columns are Richard extrapolations and the difference between the last two values |R_2,2 − R_1,1| = 0.000070 < ε_step.

Engineering

To be completed.

Error

The error is discussed in the theory section: the terms R_{i, j} have an error which reduces according to O(h^{2 + 2j}).

Questions

Question 1

Use h = 0.5 and ε_step = 0.00001 to approximate the derivative of f(x) = tan(x) at x = 1 using the centred divided- difference formula and Richardson extrapolation.

Answer: R_5,5 = 3.425518831

Question 2

Repeat Question 1 but for the function f(x) = sin(x)/x.

Answer: R_3,3 = -.3011686798

Question 3

Does it significantly hurt our approximation to use Richardson extrapolation in a situation where it does not apply?

Answer: recall that each value R_i,j approximates the integral, so any weighted average of these will still approximate the same integral.

Matlab

Matlab may be used as follows to find the derivative using Richardson extrapolation. Note that matrices start at (1,1), so we must make an adjustment. We assume a function y = f(x) is defined.

x = 1;
h = 0.5;
eps_step = 0.00001;
R(1, 1) = (f(x + h) - f(x - h))/(2*h);
for i=1:100
   h = h/2;

   R(i + 1, 1) = (f(x + h) - f(x - h))/(2*h);

   for j=1:i
      R(i + 1, j + 1) = (4^j*R(i + 1, j) - R(i, j))/(4^j - 1);
   end

   if ( abs( R(i + 1, i + 1) - R(i, i) ) < eps_step )
      break;
   elseif ( i == 100 )
      error( 'Richardson extrapolation failed to converge' );
   end
end

Maple

Maple may be used as follows to find the derivative using Richardson extrapolation:

f := x -> exp(-x) * sin(x) - 0.5 * exp(-x) * cos(x);
x := 1;
h := 0.5;
eps_step := 0.00001;
R[0, 0] := (f(x + h) - f(x - h))/(2*h);
for i from 1 to 100 do
   h := h/2;

   R[i, 0] := (f(x + h) - f(x - h))/(2*h);

   for j from 1 to i do
      R[i, j] := (4^j*R[i, j - 1] - R[i - 1, j - 1])/(4^j - 1);
   end do;

   if abs( R[i,i] - R[i - 1, i - 1] ) < eps_step then
      break;
   elif i = 100 then
      error "Richardson extrapolation failed to converge";
   end if;
end do;

Topic 12.3: Richardson Extrapolation