## Abstract

Lagrange interpolation estimates a function at a reference position *Χ* from known values of the function at distinct non-uniformly spaced points *x*_{1}, …, *x*_{n}. Here, the corresponding *n*-point finite-difference formulae are derived to estimate derivatives up to order *n*−1 at *Χ*. A recurrence relation is derived that permits the errors to be determined as a Taylor series to any accuracy. The error coefficient multiplying the *n*+*j*th derivative is a polynomial of order *j*+1 in the elementary symmetric functions for the displacements *x*_{1}−*Χ*,…,*x*_{n}−*Χ*. Appendices state finite-difference formulae to estimate the derivatives and the first four error terms for *n*=1,…,5.

## 1. Introduction

Computational engineering often requires the numerical solution of differential equations. A natural and direct way to construct finite-difference computational models is to replace the differential operators at some reference point *x*=*Χ* with discrete counterparts *D*_{d}, corresponding to derivatives of polynomial Lagrange interpolation from the function values at *n*>*d* distinct grid points *x*_{1},…,*x*_{n}. If the grid points are regularly spaced, then the *Χ*-dependence of the finite-difference operators *D*_{d} is explicitly known and tabulated (Abramowitz & Stegun 1965: eqn. 25.2.7 and eqns. 25.3.4 to 25.3.6, tables 25.1, 25.2). In applications the grid spacing might not be uniform (e.g. grid points to include sites where values are available or are sought). For non-uniform grids, Fornberg (1988, 1998) and Corless & Rokicki (1996) give neat computer algorithms that construct *D*_{d} for 0≤*d*<*n*. In §3 of the present paper, an explicit formula for *D*_{d} is derived in terms of elementary symmetric functions. Appendix A evaluates *D*_{d} in terms of the displacements *α*_{i}=*x*_{i}−*Χ* for the cases 0≤*d*<*n*, *n*=1,…,5.

In a term-by-term finite-difference model of a differential equation, the size of the errors is related to the worst of the errors that arise in replacing by *D*_{d}. For a computational scheme constructed in terms of *D*_{d}, it may be possible to make slight adjustments to the coefficients multiplying each of the *D*_{d} so that there is extra cancellation of the errors. Crandall (1955) performed such error cancellation with *n*=3 and a uniform grid for the diffusion equation at two levels in time. Mitchell & Griffiths (1980, ch. 2, table 1) demonstrate the leap in accuracy over the Crank–Nicolson scheme (1947). Smith (2000) extended the Crandall (1955) scheme to include grid non-uniformity via neat Taylor series for the *n*=3 errors in *D*_{0}, *D*_{1} and *D*_{2}. The motivation for the present paper is to derive error Taylor series for all *n*. Appendix B states the first four error terms for 0≤*d*<*n*, *n*=1,…,5. Computer algebra packages (e.g. Maple or Mathematica) make it straightforward to confirm the validity for *n*=1,…,5 of the neat error expressions.

Section 2 introduces elementary symmetric functions and states the main results, from which operators and errors can be constructed for any number of grid points. The subsequent four sections detail a direct derivation of the main results, involving generalized Vandermonde determinants and Schur functions (De Marchi 2001). Functions introduced by Shur in his 1901 thesis on groups of matrices are today called *S* or Schur functions (MacDonald 1995).

## 2. Elementary symmetric functions and main results

In this paper, *α* denotes the ordered set of displacements *α*_{i}=*x*_{i}−*Χ*. For the set *α*, the elementary symmetric functions are defined as the sum of all distinct permutations of order *i* over the set. An equivalent algebraic definition (Baker 1994; MacDonald 1995) is that for arbitrary *z*,(2.1)For indices *i*<0 or *i*>*n*, it is implicit that . The zero-order elementary symmetric function is . To minimize confusion with powers, the superscript indicating the set will usually be omitted. For example, with *n*=3,(2.2)The derivatives with respect to varied reference points are(2.3)If the chosen reference point *Χ* is the centroid, then there is the simplification(2.4)With *e*_{1}=0, equations (B 6*a*–*c*) correspond to Smith's (2000) eqns (3.3*a*) to (3.3*c*). If the reference point *Χ* coincides with any of the grid points, the simplification is(2.5)On a uniformly spaced grid with *Χ* chosen to be the centroid, *e*_{i}=0 for all odd *i*.

In §3, the *n*-point finite-difference operator *D*_{d}, operating on a function *f*(*x*), is shown to be the weighted sum of the function values at the grid points(2.6)Extensive numerical tests confirm the agreement of this explicit equation (2.6) with results from the computational algorithms of Fornberg (1988, 1998) and Corless & Rokicki (1996). *D*_{0}[*f*] is an *n*-point Lagrange interpolation (Abramowitz & Stegun 1965: eqn 25.2.2) and *D*_{d}[*f*] is the *d* th derivative, with respect to *Χ* of the Lagrange interpolation (Fornberg 1988). A mathematical way of expressing the equivalence of the subscript *d* to the number of *Χ*-derivatives is the consistency relationship(2.7)In appendix A, the sign changes and the increasing factorial numerators between successive *D*_{0}, …, *D*_{n−1} can be explained from equations (2.3) and (2.7).

If the function *f*(*x*) is not a polynomial in *x* of degree, *n*−1, then an error will arise at degree *n* or beyond. For uniform spacing, series for differences in terms of derivatives are well known (Abramowitz & Stegun 1965: eqns 25.3.16–25.3.20). In §4 it is shown that the error terms from the weighted sum of a Taylor series about the reference point *Χ* can be written as a series involving Schur functions in the displacements,(2.8)After some technical preliminaries in §5, it is shown in §6 that the higher-order Schur functions can be calculated through the recurrence relation for *j*≥*n*:(2.9)Exact arithmetic avoids instability for large *j*. An interpretation of the left-hand side term in equation (2.8) gives the low-order error coefficients for 0≤*j*<*n*:(2.10)From these degree zero starting values in equation (2.10), at the *ℓ*th application the recurrence relation in equation (2.9) generates the *j*=*n*+*ℓ*−1 term, that has homogeneous degree *n*+*ℓ*−1−*d* in the displacements and is polynomial of order *ℓ* in *e*_{1}, …, *e*_{n} In particular, the leading four error terms presented in appendix B are respectively linear, quadratic, cubic and quartic in *e*_{1}, …, *e*_{n}.

For the errors, a consequence of the consistency relationship in equation (2.7) and(2.11)In appendix B, the sign changes and decreasing factorial denominators for the lowest-order error terms *f*^{(n)}(*Χ*) in *D*_{d} can be linked to equations (2.3) and (2.11).

## 3. Derivation of difference operators

Let the operator *D*_{d}[*f*] be the weighted sum of discrete values of a function *f*(*x*) over *n* distinct points so that(3.1)Taking the Taylor series of *f*(*x*_{i}) about the position *Χ* and writing *α*_{i}=*x*_{i}−*Χ*,(3.2)To avoid convergence considerations, the circle of convergence about *Χ* is assumed to include all the *x*_{i}. Let *D*_{d,m}[*f*] represent the truncated form of *D*_{d}[*f*] with the *j*-summation terminated at degree *m*−1. For finite-term truncations, the order of *i* and *j* summations can be exchanged:(3.3)There are *n* weights *w*_{i} to be selected. The truncated operator *D*_{d,n}[*f*] can be forced to represent the *d*th derivative operator for *d*<*n*:(3.4)With the standard notation for the Kronecker delta(3.5)then the unit column vector represents the derivative to be approximated. The system to be solved can thus be written in matrix form as(3.6)Cramer's rule states that any system *Aw*=*b* with non-zero det(*A*) has a general solution for each component *w*_{y} of *w*=(*w*_{1}, …, *w*_{n}),(3.7)where the unit row vector picks out the component *w*_{y} of the solution.

In this form, the system of equations (3.6), upon factoring out and cancelling factorials, has a solution for each component(3.8)The denominator in equation (3.8) is a Vandermonde determinant (De Marchi 2001), hereafter denoted by VDM(*α*), in terms of the ordered set *α*=(*α*_{1},…,*α*_{n}). It has value(3.9)The matrix in the numerator has zero last column except for the value *d*! at the position (*d*+1, *n*+1) and it has zero last row except for 1 at the position (*n*+1, *y*). As temporary notation within this section, let(3.10a)a length *n*−1 ordered set of displacements that excludes *α*_{y}, and let(3.10b)a length *n*−1 ordered set of integers excluding *d* at position *d*+1 that arise as powers of the displacements. By expansion down the last column then the last row, the numerator in (3.8) can be written(3.11)The denominator can be evaluated in a way that involves VDM(*β*(*y*))(3.12)This is non-zero because the grid points *x*_{j}, and therefore the displacements *α*_{j}, are distinct. Then, the quotient (3.8) takes the form(3.13)where *S*_{λ}(*β*) is a Schur function over *β*(*y*) with partition *λ* (Baker 1994; MacDonald 1995):(3.14)

Partitions can be calculated by taking the difference in the powers of the numerator and the denominator in equation (3.14), in reverse order (Baker 1994; MacDonald 1995). The powers in the numerator are *γ*=(0,…,*n*−1)\(*d*) and those in the denominator are (0,…,*n*−2), so that the partition *λ* is given by(3.15)For convenience, the notation *a*^{b} represents *b* occurrences of *a*, for example, (1^{4})=(1, 1, 1, 1). Trailing zeros in partitions are dropped as they are equivalent to multiplication of the Schur function by *e*_{0}=1. The conjugate of *λ* is obtained by transposing the diagram of *λ* to give *λ*′=(*n*−*d*−1) (Baker 1994; MacDonald 1995).

The Jacobi–Trudi identity for elementary symmetric functions states (MacDonald 1995) that for an arbitrary partition *λ* of length *ℓ*(3.16)In this particular case with *λ*′=(*n*−*d*−1), the Schur function has the simple form(3.17)This gives the explicit form of equation (3.13) as(3.18)The weighted summation in equation (3.1) over all *n* of the points gives the difference operator that approximates the *d*th derivative(3.19)Also,(3.20)where the definition in equation (2.1) of elementary symmetric functions and the result has been used, that is, *β*(*i*) is only of length *n*−1. Equating powers of *z* gives:(3.21)The temporary notation *β* can be replaced in equation (3.19) to give the result:(3.22)The displacement differences *α*_{i}−*α*_{j} can also be written as grid differences *x*_{i}−*x*_{j}. Thus, the denominators do not depend on *Χ*.

*D*_{0}[*f*](*Χ*) is a polynomial of degree *n*−1 in *Χ* and can be recognized as an *n*-point Lagrange interpolation of *f*(*Χ*) (Abramowitz & Stegun 1965: eqn 25.2.2). If a general function *f*(*Χ*) is replaced by *D*_{0}[*f*](*Χ*) then the grid-point values *f*(*x*_{i}) and operators *D*_{d}[*f*](*Χ*) are unchanged. That restriction to polynomials of degree *n*−1 permits *D*_{d,n} to be replaced by *D*_{d} in the derivative matching equation (3.4). The freedom to vary *Χ* implies that *D*_{d}[*f*](*Χ*) is the *d*th derivative with respect to *Χ* of *D*_{0}[*f*](*Χ*). Fornberg (1988) made that linkage the premise for an algorithm, rather than a consequence.

## 4. Derivation of error terms

At degree *n* and beyond, errors will arise. It is useful to be able to calculate the higher-order errors, for example, to extend high-order numerical schemes to non-uniform grids (Smith 2000). The general difference operator can be written(4.1)where(4.2)The derivation in §3 for the approximate derivatives ensures that with 0≤*j*<*n*(4.3)For *j*≥*n*, the expression of equation (3.8) for the weights *w*_{i} has the consequence(4.4)The denominator is VDM(*α*). From equation (3.5) and by expansion down the last column, the numerator is(4.5)where *Γ*=(0, 1, 2, …, *n*−1, *j*)\(*d*) (similar to *γ* but including *j* at position *n*).

Then(4.6)where the partition is(4.7)The conjugate partition is of length *j*−*n*+1. Inserting the above expression into equation (4.1) gives the explicit form for the general difference operator in terms of Schur functions as(4.8)With the initial for 0≤*j*<*n* defined as(4.9)equation (4.8) can also be written as(4.10)

## 5. Preliminary results

Before the recurrence relation of equation (2.9) is derived, some preliminary results are first obtained. As used earlier, the Jacobi–Trudi identity for the conjugate partition gives the Schur functions in terms of elementary symmetric functions(5.1)where denotes element *s* of the conjugate partition . The square matrix of size *j*−*n*+1, which gives the subscripts for the elementary symmetric functions in equation (5.1), is(5.2)By the definition of equation (2.1), *e*_{i}=0 when *i*>*n* so the highest subscript that yields a non-zero elementary symmetric function is given when the subscript *i*=*n*. The first element *n*−*d* of the conjugate partition gives the subscripts *n*−*d*−*s*+*t* on the first row. So, with *s*=1, the last non-zero elementary symmetric function *e*_{n} arises when *n*−*d*−1+*t*=*n*, that is, *t*=*d*+1. Since *j*−*n*+1≥*t*, then the first row consists of the elements *e*_{n−d}, …, *e*_{n} padded with zeros for *j*≥*n*+*d*; otherwise it consists of the elements *e*_{n−d}, …, *e*_{j−d}. Accordingly, the Schur function is considered over two intervals:(5.3)For convenience, the notation refers to the upper-triangular matrix *M*_{i} (of size *i*) with row *x* removed and the notation refers to *M*_{i} with row *x* and column *y* removed. The second row of equation (5.2), and hence the first row of *M*_{j−n+1}, has final element *e*_{n} when *j*−*n*=*n*, so that *j*=2*n*, giving(5.4)The values of this matrix are a direct consequence of (5.2). The matrix is upper-triangular since for *s*>*t*+1 in equation (5.2), that is, in the strictly lower triangular region of (5.4*a*), then the conjugate partition has elements 1+*s*−*t*<0 and *e*_{i}=0 for *i*<0. For other values of *j*≥*n*, equation (5.2) shows that the matrices *M*_{j−n+1} may be defined iteratively in terms of the above equation (5.4*a*). The last rows of equation (5.4*a*) and in the second case below, equation (5.4*b*), are chosen for compatibility in this iterative definition and as a result they preserve the upper-triangular nature of *M*_{j−n+1}.(5.4)In the first case, the final column goes up from 1 to *e*_{j−n}. In the second case, the zeros at the start of the final column are a consequence of *e*_{i}=0 when *i*>*n*.

With the first row and column removed, it is clear that for *i*>1(5.5)For brevity in the following derivations, this result is also assumed for the case *i*=1. From the iterative definition it is clear that(5.6)Since *M*_{i} is upper triangular with unit diagonal elements,(5.7)For 1≤*k*, *ℓ*≤*i* and *i*≥2(5.8)This result is due to the trailing 1s on the leading diagonal of *M*_{i}. The determinant can be expanded up the leading diagonal until the first of either row *k* or column *ℓ* is reached when the trailing 1s end and the expansion of the determinant stops.

By expansion up the leading diagonal in equations (5.4*a*,*b*), when 1≤*k*, *ℓ*≤*j*−*n* and *j*−*n*≥2,(5.9)In the first case the matrix can be reduced to size by equation (5.8). Since the row removed *j*−*n*−*k*+1 is less than the column removed *ℓ*, it can be seen by considering equation (5.4) that the last row is all zero, giving the zero determinant. In the second case, when the row removed *j*−*n*−*k*+1 is greater than or equal to the column removed *ℓ*, then the matrix can be reduced to size by equation (5.8).

Expanding the determinant along the first row in equation (5.3) gives(5.10)By expansion of the determinant up the final column in equations (5.4*a*,*b*), when *ℓ*≤*j*−*n* (i.e. not removing the final column),(5.11)Strictly speaking, with *j*=*n*+1, so for compatibility with the first case above it is assumed that .

## 6. Construction of the recurrence relation

The results of §5 form the building blocks used in deriving the recurrence relation. In accordance with the intervals over which these results are valid, is considered for (a) low-order error terms *n*≤*j*≤*n*+*d*, (b) moderate-order error terms *n*+*d*<*j*≤2*n* and (c) high-order error terms *j*≥2*n*. The initial values of are defined on the interval 0≤*j*<*n* as in equation (4.9)(6.1)It is left to show that with these initial values the Schur functions can be calculated for all *j*≥*n* through the recurrence relation(6.2)

### (a) Low-order error terms: *n*≤*j*≤*n*+*d*

For the interval *n*≤*j*≤*n*+*d*, the first case in equation (5.10) gives(6.3)where the summation is split up as(6.4)The last case is with *ℓ*=*j*−*n*+1 and equations (5.6) and (5.7) have been used to simplify the determinant. When *j*=*n*, it is clear from the summation in equation (6.3) that *σ*_{1}(*j*)+*σ*_{2}(*j*)=0, since these terms do not arise. For the remaining *j*>*n*, the first case of equation (5.11) is used to give(6.5)The order of summation is exchanged to give(6.6)The inner summation of the sum *σ*_{1}(*n*)+*σ*_{2}(*j*)=0 is split such that(6.7)When *j*=*n*+1 then *σ*_{2}(*j*)=0 as *k* only takes one value in the outer summation, hence *ℓ* takes all the values in the inner summation. For *j*>*n*+1 the remaining part of the split is given by(6.8)Using the second case of equation (5.9), since from the inner summation *j*−*n*−*k*+1≥*ℓ*,(6.9)The outer summation implies that *n*≤*j*−*k*. Since *k*≥1 and *j*≤*n*+*d*, for this interval, *j*−*k*≤*n*+*d*−1. Together these inequalities imply that *n*≤*j*−*k*≤*n*+*d*, so the first case of equation (5.10) may be inserted with *j* replaced by *j*−*k* to give(6.10)Using the first case of equation (5.9), since from the inner summation *j*−*n*−*k*+1<*ℓ*,(6.11)The initial conditions of equation (6.1) are used to rewrite *σ*_{3}(*j*) as(6.12)As *j*≥*n* for this interval and *n*≥*k* from the upper limit, *j*−*k*≥0. Combining this with the lower limit gives 0≤*j*−*k*<*n*. From the initial conditions of equation (6.1), the only non-zero initial value for arises when *j*−*k*=*d* so that as required.

Finally, from equation (6.3), the recurrence relation over the interval *n*≤*j*≤*n*+*d* is(6.13)

### (b) Moderate-order error terms: *n*+*d*<*j*≤2*n*

The proofs over the remaining intervals are much the same with differing summation indices. For the interval *n*+*d*<*j*≤2*n*, the second case in equation (5.10) and the first case in equation (5.11) give(6.14)On exchanging the order of summation(6.15)where the notation *σ*_{1}(*j*), *σ*_{2}(*j*) and *σ*_{3}(*j*) is reused to again denote a split in the summation. The first part of the split is(6.16)When *j*−*n*−*d*=*j*−*n*, that is, *d*=0, the split does not arise, hence, *σ*_{2}(*j*)+*σ*_{3}(*j*)=0. The remaining terms for *d*>0 are split in the inner summation to give(6.17)and(6.18)The last case of equation (5.10) with *j* replaced by *j*−*k* gives(6.19)since from the outer summation *j*−*k*≥*n*+*d*. The second case of equation (5.9) is used on the inner summation since the limits give *j*−*n*−*k*+1≥*ℓ*, so that(6.20)Then(6.21)where the first case of equation (5.10) has been used with *j* replaced by *j*−*k*, since from the outer summation *k*≤*j*−*n* and *k*≥*j*−*n*−*d*+1 so that *n*≤*j*−*k*≤*n*+*d*−1. The remaining part of the summation for *d*≥1 is(6.22)by the first case in equation (5.9) as, from the inner summation, *ℓ*≥*j*−*n*−*k*+2 so that *j*−*n*−*k*+1<*ℓ*. When *j*<2*n* the initial conditions of equation (6.1) give(6.23)This result is elementary, since *j*>*n*+*d* for this interval and from the outer summation *n*≥*k* giving *j*−*k*>*d*, similarly *n*>*j*−*k*, and so . Then,(6.24)With equation (6.23) used as required to extend the upper limit of the summation, the recurrence relation for *n*+*d*<*j*≤2*n* is(6.25)

### (c) High-order error terms: *j*≥2*n*

For the interval *j*≥2*n*, the second cases in equations (5.10) and (5.11) give(6.26)Exchanging the order of summation and splitting the summations into three parts with further reuse of the *σ*_{1}(*j*), *σ*_{2}(*j*) and *σ*_{3}(*j*) notation,(6.27)

The first part of the split summation is(6.28)where the second case of equation (5.9) has been used since from the summation limits and the second case of equation (5.10) has been used since from the outer summation *j*−*k*≥*n*+*d*. For *j*≥2*n*+*d*, *σ*_{2}(*j*)+*σ*_{3}(*j*)=0 since *σ*_{2}(*j*) and *σ*_{3}(*j*) do not arise in this case. For *j*<2*n*+*d*(6.29)where the second case of equation (5.9) has been used since from the inner summation . Then(6.30)where the first case of equation (5.10) has been used since in this interval *j*≥2*n*, from the upper limit *n*≥*k* so that *n*≤*j*−*k* and, from the lower limit, , which combined give . The last split of the summation is(6.31)where the first case of equation (5.9) has been used since from the inner summation . Finally, the recurrence relation for *j*≥2*n* is(6.32)completing the proof of the recurrence relation of equation (6.2) with initial conditions of equation (6.1).

## 7. Concluding remarks

Explicit one-dimensional difference operators *D*_{d} have been derived that mimic derivative operators at a reference point *Χ* for any number *n* of distinct points *x*_{1}, …, *x*_{n} over an irregular grid and for any derivative *d*<*n*. Along with these, a recurrence relation has been derived that allows calculation of Taylor series for the errors. The *n*+*j*th derivative error terms are polynomials of order *j*+1 in the elementary symmetric functions for the displacements *x*_{1}−*Χ*,…, *x*_{n}−*Χ*.

The Taylor series for the errors makes it elementary to obtain the error from a linear sum of *D*_{d} terms, for example, when selecting coefficients in a finite-difference scheme to mimic a differential equation. At all accuracy levels, the error coefficients involve polynomials in the *n* non-constant elementary symmetric functions *e*_{1},…,*e*_{n} for the set of displacements. The difference operators *D*_{d}, together with the elementary symmetric functions, are a natural combination of tools with which to extend high-order numerical schemes from uniform to non-uniform grids.

## Acknowledgements

M.K.B. would like to thank the Natural Environment Research Council for Ph.D. funding.

## Footnotes

- Received July 9, 2004.
- Accepted October 11, 2004.

- © 2005 The Royal Society