Divided differences

The divided difference scheme can be put into an upper triangular matrix. Let ${\displaystyle 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}]\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\ldots &f[x_{n}]\end{pmatrix}}}$.

Then it holds

• ${\displaystyle T_{f+g}x=T_{f}x+T_{g}x}$
• ${\displaystyle 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 form a commutative ring.

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

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

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

${\displaystyle 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 ${\displaystyle T_{f}x}$ can be written as

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

${\displaystyle {\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 a polynomial function ${\displaystyle q}$ with ${\displaystyle q(\xi )=(\xi -x_{0})\cdots (\xi -x_{n})}$ this can be written as

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

Alternatively, we can allow counting backwards from the start of the sequence by defining ${\displaystyle x_{k}=x_{k+n+1}=x_{k-(n+1)}}$ whenever ${\displaystyle k<0}$ or ${\displaystyle n<k}$. This definition allows ${\displaystyle x_{-1}}$ to be interpreted as ${\displaystyle x_{n}}$, ${\displaystyle x_{-2}}$ to be interpreted as ${\displaystyle x_{n-1}}$, ${\displaystyle x_{-n}}$ to be interpreted as ${\displaystyle x_{0}}$, etc. The expanded form of the divided difference thus becomes

${\displaystyle f[x_{0},\dots ,x_{n}]=\sum _{j=0}^{n}{\frac {f(x_{j})}{\prod \limits _{k=j-n}^{j-1}(x_{j}-x_{k})}}+\sum _{j=0}^{n}{\frac {f(x_{j})}{\prod \limits _{k=j+1}^{j+n}(x_{j}-x_{k})}}}$

Yet another characterization utilizes limits:

${\displaystyle f[x_{0},\dots ,x_{n}]=\sum _{j=0}^{n}\lim _{x\rightarrow x_{j}}\left[{\frac {f(x_{j})(x-x_{j})}{\prod \limits _{k=0}^{n}(x-x_{k})}}\right]}$

You can represent partial fractions using the expanded form of divided differences. (This does not simplify computation, but is interesting in itself.) If ${\displaystyle p}$ and ${\displaystyle q}$ are polynomial functions, where ${\displaystyle \mathrm {deg} \ p<\mathrm {deg} \ q}$ and ${\displaystyle q}$ is given in terms of linear factors by ${\displaystyle q(\xi )=(\xi -x_{1})\cdot \dots \cdot (\xi -x_{n})}$, then it follows from partial fraction decomposition that

${\displaystyle {\frac {p(\xi )}{q(\xi )}}=\left(t\to {\frac {p(t)}{\xi -t}}\right)[x_{1},\dots ,x_{n}].}$

If limits of the divided differences are accepted, then this connection does also hold, if some of the ${\displaystyle x_{j}}$ coincide.

If ${\displaystyle f}$ is a polynomial function with arbitrary degree and it is decomposed by ${\displaystyle f(x)=p(x)+q(x)\cdot d(x)}$ using polynomial division of ${\displaystyle f}$ by ${\displaystyle q}$, then

${\displaystyle {\frac {p(\xi )}{q(\xi )}}=\left(t\to {\frac {f(t)}{\xi -t}}\right)[x_{1},\dots ,x_{n}].}$

The divided differences can be expressed as

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

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

This is called the Peano form of the divided differences and ${\displaystyle B_{n-1}}$ is called the Peano kernel for the divided differences, both named after Giuseppe Peano.

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 similar values. However we know, that difference quotients approximate the derivative and vice versa:

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

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

${\displaystyle 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 }$

${\displaystyle \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 ${\displaystyle y-x}$ by expanding the Taylor series at the center between ${\displaystyle x}$ and ${\displaystyle y}$:

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

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

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

${\displaystyle {\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 }$

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 ${\displaystyle f}$ to a divided difference ${\displaystyle 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: ${\displaystyle p_{n}(x)=x^{n}.}$

Regular Taylor series is a weighted sum of power functions: ${\displaystyle 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: ${\displaystyle 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 ${\displaystyle n}$ terms vanish, because we have a higher difference order than polynomial order, and in the following term the divided difference is one:

${\displaystyle {\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 ${\displaystyle {\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 ${\displaystyle f}$. The nice thing is, that there is a simpler way. It holds

${\displaystyle 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 ${\displaystyle p_{n}}$ by a division of formal power series. See how this reduces to the successive computation of powers when we compute ${\displaystyle p_{n}[h]}$ for several ${\displaystyle 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.

Given n data points

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

with

${\displaystyle x_{\nu }=x_{0}+\nu h,\ h>0,\ \nu =0,\ldots ,n-1}$

the divided differences can be calculated via forward differences defined as

${\displaystyle \Delta ^{(0)}y_{i}:=y_{i}}$

${\displaystyle \Delta ^{(k)}y_{i}:=\Delta ^{(k-1)}y_{i+1}-\Delta ^{(k-1)}y_{i},\ k\geq 1.}$

The relationship between divided differences and forward differences is[4]

${\displaystyle f[x_{0},x_{1},\ldots ,x_{k}]={\frac {1}{k!h^{k}}}\Delta ^{(k)}f(x_{0}).}$

${\displaystyle {\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}}}$