Divided differences

In mathematics, divided differences is an algorithm, historically used for computing tables of logarithms and trigonometric functions. Charles Babbage's difference engine, an early mechanical calculator, was designed to use this algorithm in its operation.

Divided differences is a recursive division process. Given a sequence of data points $(x_{0},y_{0}),\ldots ,(x_{n},y_{n})$, the method calculates the coefficients of the interpolation polynomial of these points in the Newton form.

Definition

Given n + 1 data points

$(x_{0},y_{0}),\ldots ,(x_{n},y_{n})$

where the $x_{k}$ are assumed to be pairwise distinct, the forward divided differences are defined as:

$[y_{k}]:=y_{k},\qquad k\in \{0,\ldots ,n\}$

$[y_{k},\ldots ,y_{k+j}]:={\frac {[y_{k+1},\ldots ,y_{k+j}]-[y_{k},\ldots ,y_{k+j-1}]}{x_{k+j}-x_{k}}},\qquad k\in \{0,\ldots ,n-j\},\ j\in \{1,\ldots ,k\}.$

To make the recursive process of computation clearer, the divided differences can be put in tabular form, where the columns correspond to the value of j above, and each entry in the table is computed from the difference of the entries to its immediate lower left and to its immediate upper left, divided by a difference of corresponding x-values:

${\begin{matrix}x_{0}&y_{0}=[y_{0}]&&&\\&&[y_{0},y_{1}]&&\\x_{1}&y_{1}=[y_{1}]&&[y_{0},y_{1},y_{2}]&\\&&[y_{1},y_{2}]&&[y_{0},y_{1},y_{2},y_{3}]\\x_{2}&y_{2}=[y_{2}]&&[y_{1},y_{2},y_{3}]&\\&&[y_{2},y_{3}]&&\\x_{3}&y_{3}=[y_{3}]&&&\\\end{matrix}}$

Notation

Note that the divided difference $[y_{k},\ldots ,y_{k+j}]$ depends on the values $x_{k},\ldots ,x_{k+j}$ and $y_{k},\ldots ,y_{k+j}$, but the notation hides the dependency on the x-values. If the data points are given by a function ƒ,

$(x_{0},f(x_{0})),\ldots ,(x_{k},f(x_{n}))$

one sometimes writes

$f[x_{k},\ldots ,x_{k+j}]$

for the divided difference

$[f(x_{k}),\ldots ,f(x_{k+j})]$

Several other notations for the divided difference of the function ƒ on the nodes x₀, ..., x_n are also used, for instance:

$[x_{0},\ldots ,x_{n}]f,$

$[x_{0},\ldots ,x_{n};f],$

$D[x_{0},\ldots ,x_{n}]f$

Example

Divided differences for $k=0$ and the first few values of $j$:

${\begin{aligned}{\mathopen {[}}y_{0}]&=y_{0}\\{\mathopen {[}}y_{0},y_{1}]&={\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}\\{\mathopen {[}}y_{0},y_{1},y_{2}]&={\frac {{\mathopen {[}}y_{1},y_{2}]-{\mathopen {[}}y_{0},y_{1}]}{x_{2}-x_{0}}}={\frac {{\frac {y_{2}-y_{1}}{x_{2}-x_{1}}}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}}{x_{2}-x_{0}}}={\frac {y_{2}-y_{1}}{(x_{2}-x_{1})(x_{2}-x_{0})}}-{\frac {y_{1}-y_{0}}{(x_{1}-x_{0})(x_{2}-x_{0})}}\\{\mathopen {[}}y_{0},y_{1},y_{2},y_{3}]&={\frac {{\mathopen {[}}y_{1},y_{2},y_{3}]-{\mathopen {[}}y_{0},y_{1},y_{2}]}{x_{3}-x_{0}}}\end{aligned}}$

Properties

• Linearity

$(f+g)[x_{0},\dots ,x_{n}]=f[x_{0},\dots ,x_{n}]+g[x_{0},\dots ,x_{n}]$

$(\lambda \cdot f)[x_{0},\dots ,x_{n}]=\lambda \cdot f[x_{0},\dots ,x_{n}]$

• Leibniz rule

$(f\cdot g)[x_{0},\dots ,x_{n}]=f[x_{0}]\cdot g[x_{0},\dots ,x_{n}]+f[x_{0},x_{1}]\cdot g[x_{1},\dots ,x_{n}]+\dots +f[x_{0},\dots ,x_{n}]\cdot g[x_{n}]=\sum _{r=0}^{n}f[x_{0},\ldots ,x_{r}]\cdot g[x_{r},\ldots ,x_{n}]$

• Divided differences are symmetric: If $\sigma :\{0,\dots ,n\}\to \{0,\dots ,n\}$ is a permutation then

$f[x_{0},\dots ,x_{n}]=f[x_{\sigma (0)},\dots ,x_{\sigma (n)}]$

• Polynomial interpolation in the Newton form: if $p$ is a polynomial function of degree $\leq n$, then

$p(x)=p[x_{0}]+p[x_{0},x_{1}](x-x_{0})+p[x_{0},x_{1},x_{2}](x-x_{0})(x-x_{1})+\cdots +p[x_{0},\ldots ,x_{n}](x-x_{0})(x-x_{1})\cdots (x-x_{n-1})$

• If $p$ is a polynomial function of degree $<n$, then

$p[x_{0},\dots ,x_{n}]=0.$

• Mean value theorem for divided differences: if $f$ is n times differentiable, then

$f[x_{0},\dots ,x_{n}]={\frac {f^{(n)}(\xi )}{n!}}$ for a number $\xi$ in the open interval determined by the smallest and largest of the $x_{k}$'s.

Matrix form

The divided difference scheme can be put into an upper triangular matrix:

$T_{f}(x_{0},\dots ,x_{n})={\begin{pmatrix}f[x_{0}]&f[x_{0},x_{1}]&f[x_{0},x_{1},x_{2}]&\ldots &f[x_{0},\dots ,x_{n}]\\0&f[x_{1}]&f[x_{1},x_{2}]&\ldots &f[x_{1},\dots ,x_{n}]\\0&0&f[x_{2}]&\ldots &f[x_{2},\dots ,x_{n}]\\\vdots &\vdots &&\ddots &\vdots \\0&0&0&\ldots &f[x_{n}]\end{pmatrix}}$.

Then it holds

• $T_{f+g}(x)=T_{f}(x)+T_{g}(x)$
• $T_{\lambda f}(x)=\lambda T_{f}(x)$ if $\lambda$ is a scalar
• $T_{f\cdot g}(x)=T_{f}(x)\cdot T_{g}(x)$

This follows from the Leibniz rule. It means that multiplication of such matrices is commutative. Summarised, the matrices of divided difference schemes with respect to the same set of nodes x form a commutative ring.

• Since $T_{f}(x)$ is a triangular matrix, its eigenvalues are obviously $f(x_{0}),\dots ,f(x_{n})$.
• Let $\delta _{\xi }$ be a Kronecker delta-like function, that is

$\delta _{\xi }(t)={\begin{cases}1&:t=\xi ,\\0&:{\mbox{else}}.\end{cases}}$

Obviously $f\cdot \delta _{\xi }=f(\xi )\cdot \delta _{\xi }$, thus $\delta _{\xi }$ is an eigenfunction of the pointwise function multiplication. That is $T_{\delta _{x_{i}}}(x)$ is somehow an "eigenmatrix" of $T_{f}(x)$: $T_{f}(x)\cdot T_{\delta _{x_{i}}}(x)=f(x_{i})\cdot T_{\delta _{x_{i}}}(x)$. However, all columns of $T_{\delta _{x_{i}}}(x)$ are multiples of each other, the matrix rank of $T_{\delta _{x_{i}}}(x)$ is 1. So you can compose the matrix of all eigenvectors of $T_{f}(x)$ from the $i$-th column of each $T_{\delta _{x_{i}}}(x)$. Denote the matrix of eigenvectors with $U(x)$. Example

$U(x_{0},x_{1},x_{2},x_{3})={\begin{pmatrix}1&{\frac {1}{(x_{1}-x_{0})}}&{\frac {1}{(x_{2}-x_{0})\cdot (x_{2}-x_{1})}}&{\frac {1}{(x_{3}-x_{0})\cdot (x_{3}-x_{1})\cdot (x_{3}-x_{2})}}\\0&1&{\frac {1}{(x_{2}-x_{1})}}&{\frac {1}{(x_{3}-x_{1})\cdot (x_{3}-x_{2})}}\\0&0&1&{\frac {1}{(x_{3}-x_{2})}}\\0&0&0&1\end{pmatrix}}$

The diagonalization of $T_{f}(x)$ can be written as

$U(x)\cdot \operatorname {diag} (f(x_{0}),\dots ,f(x_{n}))=T_{f}(x)\cdot U(x)$.

Polynomials and power series

The matrix

$J={\begin{pmatrix}x_{0}&1&0&0&\cdots &0\\0&x_{1}&1&0&\cdots &0\\0&0&x_{2}&1&&0\\\vdots &\vdots &&\ddots &\ddots &\\0&0&0&0&\;\ddots &1\\0&0&0&0&&x_{n}\end{pmatrix}}$

contains the divided difference scheme for the identity function with respect to the nodes $x_{0},\dots ,x_{n}$, thus $J^{m}$ contains the divided differences for the power function with exponent $m$. Consequently, you can obtain the divided differences for a polynomial function $p$ by applying $p$ to the matrix $J$: If

$p(\xi )=a_{0}+a_{1}\cdot \xi +\dots +a_{m}\cdot \xi ^{m}$

and

$p(J)=a_{0}+a_{1}\cdot J+\dots +a_{m}\cdot J^{m}$

then

$T_{p}(x)=p(J).$

This is known as Opitz' formula.

Now consider increasing the degree of $p$ to infinity, i.e. turn the Taylor polynomial into a Taylor series. Let $f$ be a function which corresponds to a power series. You can compute the divided difference scheme for $f$ by applying the corresponding matrix series to $J$: If

$f(\xi )=\sum _{k=0}^{\infty }a_{k}\xi ^{k}$

and

$f(I)=\sum _{k=0}^{\infty }a_{k}J^{k}$

then

$T_{f}(x)=f(J).$

Alternative characterizations

Expanded form

${\begin{aligned}f[x_{0}]&=f(x_{0})\\f[x_{0},x_{1}]&={\frac {f(x_{0})}{(x_{0}-x_{1})}}+{\frac {f(x_{1})}{(x_{1}-x_{0})}}\\f[x_{0},x_{1},x_{2}]&={\frac {f(x_{0})}{(x_{0}-x_{1})\cdot (x_{0}-x_{2})}}+{\frac {f(x_{1})}{(x_{1}-x_{0})\cdot (x_{1}-x_{2})}}+{\frac {f(x_{2})}{(x_{2}-x_{0})\cdot (x_{2}-x_{1})}}\\f[x_{0},x_{1},x_{2},x_{3}]&={\frac {f(x_{0})}{(x_{0}-x_{1})\cdot (x_{0}-x_{2})\cdot (x_{0}-x_{3})}}+{\frac {f(x_{1})}{(x_{1}-x_{0})\cdot (x_{1}-x_{2})\cdot (x_{1}-x_{3})}}+{\frac {f(x_{2})}{(x_{2}-x_{0})\cdot (x_{2}-x_{1})\cdot (x_{2}-x_{3})}}+\\&\quad \quad {\frac {f(x_{3})}{(x_{3}-x_{0})\cdot (x_{3}-x_{1})\cdot (x_{3}-x_{2})}}\\f[x_{0},\dots ,x_{n}]&=\sum _{j=0}^{n}{\frac {f(x_{j})}{\prod _{k\in \{0,\dots ,n\}\setminus \{j\}}(x_{j}-x_{k})}}\end{aligned}}$

With the help of the polynomial function $\omega (\xi )=(\xi -x_{0})\cdots (\xi -x_{n})$ this can be written as

$f[x_{0},\dots ,x_{n}]=\sum _{j=0}^{n}{\frac {f(x_{j})}{\omega '(x_{j})}}.$

Peano form

If $x_{0}<x_{1}<\cdots <x_{n}$ and $n\geq 1$, the divided differences can be expressed as

$f[x_{0},\ldots ,x_{n}]={\frac {1}{(n-1)!}}\int _{x_{0}}^{x_{n}}f^{(n)}(t)\;B_{n-1}(t)\,dt$

where $f^{(n)}$ is the $n$-th derivative of the function $f$ and $B_{n-1}$ is a certain B-spline of degree $n-1$ for the data points $x_{0},\dots ,x_{n}$, given by the formula

$B_{n-1}(t)=\sum _{k=0}^{n}{\frac {(\max(0,x_{k}-t))^{n-1}}{\omega '(x_{k})}}$

This is a consequence of the Peano kernel theorem; it is called the Peano form of the divided differences and $B_{n-1}$ is the Peano kernel for the divided differences, all named after Giuseppe Peano.

Taylor form

First order

If nodes are cumulated, then the numerical computation of the divided differences is inaccurate, because you divide almost two zeros, each of which with a high relative error due to differences of approximations to similar values. However we know, that difference quotients approximate the derivative and vice versa:

${\frac {f(y)-f(x)}{y-x}}\approx f'(x)$ for $x\approx y$

This approximation can be turned into an identity whenever Taylor's theorem applies.

$f(y)=f(x)+f'(x)\cdot (y-x)+f''(x)\cdot {\frac {(y-x)^{2}}{2!}}+f'''(x)\cdot {\frac {(y-x)^{3}}{3!}}+\dots$

$\Rightarrow {\frac {f(y)-f(x)}{y-x}}=f'(x)+f''(x)\cdot {\frac {y-x}{2!}}+f'''(x)\cdot {\frac {(y-x)^{2}}{3!}}+\dots$

You can eliminate the odd powers of $y-x$ by expanding the Taylor series at the center between $x$ and $y$:

$x=m-h,y=m+h$, that is $m={\frac {x+y}{2}},h={\frac {y-x}{2}}$

$f(m+h)=f(m)+f'(m)\cdot h+f''(m)\cdot {\frac {h^{2}}{2!}}+f'''(m)\cdot {\frac {h^{3}}{3!}}+\dots$

$f(m-h)=f(m)-f'(m)\cdot h+f''(m)\cdot {\frac {h^{2}}{2!}}-f'''(m)\cdot {\frac {h^{3}}{3!}}+\dots$

${\frac {f(y)-f(x)}{y-x}}={\frac {f(m+h)-f(m-h)}{2\cdot h}}=f'(m)+f'''(m)\cdot {\frac {h^{2}}{3!}}+\dots$

Higher order

The Taylor series or any other representation with function series can in principle be used to approximate divided differences. Taylor series are infinite sums of power functions. The mapping from a function $f$ to a divided difference $f[x_{0},\dots ,x_{n}]$ is a linear functional. We can as well apply this functional to the function summands.

Express power notation with an ordinary function: $p_{n}(x)=x^{n}.$

Regular Taylor series is a weighted sum of power functions: $f=f(0)\cdot p_{0}+f'(0)\cdot p_{1}+{\frac {f''(0)}{2!}}\cdot p_{2}+{\frac {f'''(0)}{3!}}\cdot p_{3}+\dots$

Taylor series for divided differences: $f[x_{0},\dots ,x_{n}]=f(0)\cdot p_{0}[x_{0},\dots ,x_{n}]+f'(0)\cdot p_{1}[x_{0},\dots ,x_{n}]+{\frac {f''(0)}{2!}}\cdot p_{2}[x_{0},\dots ,x_{n}]+{\frac {f'''(0)}{3!}}\cdot p_{3}[x_{0},\dots ,x_{n}]+\dots$

We know that the first $n$ terms vanish, because we have a higher difference order than polynomial order, and in the following term the divided difference is one:

${\begin{array}{llcl}\forall j<n&p_{j}[x_{0},\dots ,x_{n}]&=&0\\&p_{n}[x_{0},\dots ,x_{n}]&=&1\\&p_{n+1}[x_{0},\dots ,x_{n}]&=&x_{0}+\dots +x_{n}\\&p_{n+m}[x_{0},\dots ,x_{n}]&=&\sum _{a\in \{0,\dots ,n\}^{m}{\text{ with }}a_{1}\leq a_{2}\leq \dots \leq a_{m}}\prod _{j\in a}x_{j}.\\\end{array}}$

It follows that the Taylor series for the divided difference essentially starts with ${\frac {f^{(n)}(0)}{n!}}$ which is also a simple approximation of the divided difference, according to the mean value theorem for divided differences.

If we would have to compute the divided differences for the power functions in the usual way, we would encounter the same numerical problems that we had when computing the divided difference of $f$. The nice thing is, that there is a simpler way. It holds

$t^{n}=(1-x_{0}\cdot t)\dots \cdot (1-x_{n}\cdot t)\cdot (p_{0}[x_{0},\dots ,x_{n}]+p_{1}[x_{0},\dots ,x_{n}]\cdot t+p_{2}[x_{0},\dots ,x_{n}]\cdot t^{2}+\dots ).$

Consequently, we can compute the divided differences of $p_{n}$ by a division of formal power series. See how this reduces to the successive computation of powers when we compute $p_{n}[h]$ for several $n$.

If you need to compute a whole divided difference scheme with respect to a Taylor series, see the section about divided differences of power series.

Forward differences

Further information: Finite difference

When the data points are equidistantly distributed we get the special case called forward differences. They are easier to calculate than the more general divided differences.

Given n+1 data points

$(x_{0},y_{0}),\ldots ,(x_{n},y_{n})$

with

$x_{k}=x_{0}+kh,\ {\text{ for }}\ k=0,\ldots ,n{\text{ and fixed }}h>0$

the forward differences are defined as

$\Delta ^{(0)}y_{k}:=y_{k},\qquad k=0,\ldots ,n$

$\Delta ^{(j)}y_{k}:=\Delta ^{(j-1)}y_{k+1}-\Delta ^{(j-1)}y_{k},\qquad k=0,\ldots ,n-j,\ j=1,\dots ,n.$

${\begin{matrix}y_{0}&&&\\&\Delta y_{0}&&\\y_{1}&&\Delta ^{2}y_{0}&\\&\Delta y_{1}&&\Delta ^{3}y_{0}\\y_{2}&&\Delta ^{2}y_{1}&\\&\Delta y_{2}&&\\y_{3}&&&\\\end{matrix}}$

The relationship between divided differences and forward differences is

$[y_{0},y_{1},\ldots ,y_{n}]={\frac {1}{n!h^{n}}}\Delta ^{(n)}y_{0}.$

See also

• Difference quotient
• Neville's algorithm
• Polynomial interpolation
• Mean value theorem for divided differences
• Nörlund–Rice integral
• Pascal's triangle

References

• Louis Melville Milne-Thomson (2000) [1933]. The Calculus of Finite Differences. American Mathematical Soc. Chapter 1: Divided Differences. ISBN 978-0-8218-2107-7.
• Myron B. Allen; Eli L. Isaacson (1998). Numerical Analysis for Applied Science. John Wiley & Sons. Appendix A. ISBN 978-1-118-03027-1.
• Ron Goldman (2002). Pyramid Algorithms: A Dynamic Programming Approach to Curves and Surfaces for Geometric Modeling. Morgan Kaufmann. Chapter 4:Newton Interpolation and Difference Triangles. ISBN 978-0-08-051547-2.

External links

• An implementation in Haskell.