Recursive difference formulas for derivatives can be obtained by canceling
the truncation error at each order of numerical approximation. This method is
called Richardson extrapolation. It can be used only if the data values
are equally spaced with step size h. We consider first the simplest cases
in this hierarchy and then construct the general algorithm.
Denote the three-point central difference approximation as D1(h). It has the truncation error of order O(h2). By canceling the O(h2) error, we define a new central difference:
D2(h) = ( 4D1(h) - D1(2h)) / 3
The new approximation D2(h) is in fact the five-point central difference:
D2(h) = ( -I2 + 8 I1 - 8 I-1 + I-2 ) / 12 h
The five-point central difference is more accurate than the three-point difference: it has the truncation error of order O(h4)!
Denote the two-point forward difference approximation as D1(h). It has the truncation error of order O(h). By canceling the O(h) error, we define a new forward difference:
D2(h) = 2D1(h) - D1(2h)
The new approximation D2(h) is in fact the three-point forward difference:
D2(h) = ( -I2 + 4 I1 - 3 I0) / 2 h
The three-point forward difference is more accurate than the two-point difference: it has the truncation error of order O(h2)!
Construct the following mapping from Dk(h) to Dk+1(h):
Dk+1 = Dk(h) + (Dk(h) - Dk(2h)) / (4k - 1)
In order to compute the central difference approximation up to the mth order, one needs to compute the central difference approximations of lower order with larger step sizes: 2h, 4h, ..., up to 2m-1h. These approximations can be arranged in a table of recursive derivatives (similar to the table of Newton's divided differences):
| step size | three-point difference |
five-point difference |
nine-point difference |
17-point difference |
|---|---|---|---|---|
| h | D1(h) | D2(h) | D3(h) | D4(h) |
| 2h | D1(2h) | D2(2h) | D3(2h) | |
| 4h | D1(4h) | D1(4h) | ||
| 8h | D1(8h) |
The first row of this table gives a sequence of higher-order central difference approximations for the first derivative I'(t0), where t = t0 is the central point surrounded by the points t = t0 - 2k-1h and t = t0 + 2k-1h, for k = 1,2,...,m. The approximation Dk(h) has the truncation error O(h2k). If h is small, the truncation error rapidly decreases with larger k. However, a higher-order central difference does not lead to the most optimal approximation for the derivative due to the round-off error. There exists usually an optimal order k of the central difference approximation that gives the minimum of the total of the truncation and round-off errors.
Similar hierarchies can be constructed for central difference approximations of second-order and higher-order derivatives. The starting approximation for these hierarchies should be the central difference formula of order of O(h2)), e.g. for the second derivative:
D1(h) = ( I1 - 2 I0 + I-1 ) / h2
Hierarchies for forward and backward differences can be constructed from the algorithm above by replacing the numerical factor 4k by 2k.
