Algorithm for computing polynomial coefficients
In mathematics, divided differences is an algorithm, historically used for computing tables of logarithms and trigonometric functions.[citation needed] Charles Babbage's difference engine, an early mechanical calculator, was designed to use this algorithm in its operation.[1]
Divided differences is a recursive division process. Given a sequence of data points
, the method calculates the coefficients of the interpolation polynomial of these points in the Newton form.
Definition
Given n + 1 data points
![{\displaystyle (x_{0},y_{0}),\ldots ,(x_{n},y_{n})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2b31b2f2436513ff671f30a751d6c5f6189b7ca5)
where the
![{\displaystyle x_{k}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6d2b88c64c76a03611549fb9b4cf4ed060b56002)
are assumed to be pairwise distinct, the
forward divided differences are defined as:
![{\displaystyle {\begin{aligned}{\mathopen {[}}y_{k}]&:=y_{k},&&k\in \{0,\ldots ,n\}\\{\mathopen {[}}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}}},&&k\in \{0,\ldots ,n-j\},\ j\in \{1,\ldots ,n\}.\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2a88275cd95e7358b1eebecb1d8bddb800676026)
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:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d60ef92b00038923edcedaecccbadd25373b07c1)
Notation
Note that the divided difference
depends on the values
and
, but the notation hides the dependency on the x-values. If the data points are given by a function f,
![{\displaystyle (x_{0},y_{0}),\ldots ,(x_{k},y_{n})=(x_{0},f(x_{0})),\ldots ,(x_{n},f(x_{n}))}](https://wikimedia.org/api/rest_v1/media/math/render/svg/71ef4071f5c9fb8e0a87e97c5bea8586af94a9e1)
one sometimes writes the divided difference in the notation
![{\displaystyle f[x_{k},\ldots ,x_{k+j}]\ {\stackrel {\text{def}}{=}}\ [f(x_{k}),\ldots ,f(x_{k+j})]=[y_{k},\ldots ,y_{k+j}].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d4cab5dd7621f97867cb52b2176197b00a80dafc)
Other notations for the divided difference of the function
ƒ on the nodes
x0, ...,
xn are:
![{\displaystyle f[x_{k},\ldots ,x_{k+j}]={\mathopen {[}}x_{0},\ldots ,x_{n}]f={\mathopen {[}}x_{0},\ldots ,x_{n};f]=D[x_{0},\ldots ,x_{n}]f.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8c20017aea8a25909a2c547f84cacbce55615c49)
Example
Divided differences for
and the first few values of
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c7c5a021041c0405a351383f2e34f0de88e4e58b)
Thus, the table corresponding to these terms upto two columns has the following form:
![{\displaystyle {\begin{matrix}x_{0}&y_{0}&&\\&&{y_{1}-y_{0} \over x_{1}-x_{0}}&\\x_{1}&y_{1}&&{{y_{2}-y_{1} \over x_{2}-x_{1}}-{y_{1}-y_{0} \over x_{1}-x_{0}} \over x_{2}-x_{0}}\\&&{y_{2}-y_{1} \over x_{2}-x_{1}}&\\x_{2}&y_{2}&&\vdots \\&&\vdots &\\\vdots &&&\vdots \\&&\vdots &\\x_{n}&y_{n}&&\\\end{matrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a41df0780e824759efadd4de5fcd73e8b95f07e5)
Properties
- Linearity
![{\displaystyle {\begin{aligned}(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}]\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b3b56f7d22b2b850570cca476f8668b0f19a0618)
- Leibniz rule
![{\displaystyle (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}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7eb7cdd2db0caf752548b0e82bc8cd744a30cb18)
- Divided differences are symmetric: If
is a permutation then ![{\displaystyle f[x_{0},\dots ,x_{n}]=f[x_{\sigma (0)},\dots ,x_{\sigma (n)}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/02ce355ce667f1d9c4825733b42075bb74702610)
- Polynomial interpolation in the Newton form: if
is a polynomial function of degree
, and
is the divided difference, then ![{\displaystyle P_{n-1}(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})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/993d474f737e93d7bc0249c6c6458d5395f3e44b)
- If
is a polynomial function of degree
, then ![{\displaystyle p[x_{0},\dots ,x_{n}]=0.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2c7e5d88eb9d0afc0358a14c9b337b291064de44)
- Mean value theorem for divided differences: if
is n times differentiable, then ![{\displaystyle f[x_{0},\dots ,x_{n}]={\frac {f^{(n)}(\xi )}{n!}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b21b9d4d2d31605f2a2049acd94166b8ca3ac5b8)
for a number
in the open interval determined by the smallest and largest of the
's.
Matrix form
The divided difference scheme can be put into an upper triangular matrix:
![{\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}]\\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}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9779253726b0fdf9b90a533cd7cb1561be1d7b34)
Then it holds
![{\displaystyle T_{f+g}(x)=T_{f}(x)+T_{g}(x)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f97de77053d4d4b3e39a50443baef6a874b56117)
if
is a scalar
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
is a triangular matrix, its eigenvalues are obviously
. - Let
be a Kronecker delta-like function, that is ![{\displaystyle \delta _{\xi }(t)={\begin{cases}1&:t=\xi ,\\0&:{\mbox{else}}.\end{cases}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/75d8040f34849330d6ae47b0622bad1ffcda24d0)
Obviously
, thus
is an eigenfunction of the pointwise function multiplication. That is
is somehow an "eigenmatrix" of
:
. However, all columns of
are multiples of each other, the matrix rank of
is 1. So you can compose the matrix of all eigenvectors of
from the
-th column of each
. Denote the matrix of eigenvectors with
. 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})(x_{2}-x_{1})}}&{\frac {1}{(x_{3}-x_{0})(x_{3}-x_{1})(x_{3}-x_{2})}}\\0&1&{\frac {1}{(x_{2}-x_{1})}}&{\frac {1}{(x_{3}-x_{1})(x_{3}-x_{2})}}\\0&0&1&{\frac {1}{(x_{3}-x_{2})}}\\0&0&0&1\end{pmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/84b0b9f077aa604fba6e3e4be6024a0b52d80f48)
The diagonalization of
can be written as ![{\displaystyle U(x)\cdot \operatorname {diag} (f(x_{0}),\dots ,f(x_{n}))=T_{f}(x)\cdot U(x).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dad127dd2098dbe6df407819b77ffb5aca059594)
Polynomials and power series
The matrix
![{\displaystyle 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9476776f7a5a5df0f14e860a315b8810c15805fb)
contains the divided difference scheme for the
identity function with respect to the nodes
![{\displaystyle x_{0},\dots ,x_{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/db64119863b01d965d94edaf33a384f1462f7adf)
, thus
![{\displaystyle J^{m}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/99b7529b2d9d206785626cf1811aba4d70ac9e5a)
contains the divided differences for the
power function with
exponent ![{\displaystyle m}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0a07d98bb302f3856cbabc47b2b9016692e3f7bc)
. Consequently, you can obtain the divided differences for a
polynomial function ![{\displaystyle p}](https://wikimedia.org/api/rest_v1/media/math/render/svg/81eac1e205430d1f40810df36a0edffdc367af36)
by applying
![{\displaystyle p}](https://wikimedia.org/api/rest_v1/media/math/render/svg/81eac1e205430d1f40810df36a0edffdc367af36)
to the matrix
![{\displaystyle J}](https://wikimedia.org/api/rest_v1/media/math/render/svg/359e4f407b49910e02c27c2f52e87a36cd74c053)
: If
![{\displaystyle p(\xi )=a_{0}+a_{1}\cdot \xi +\dots +a_{m}\cdot \xi ^{m}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/831ed179760178d377bd7c7425bb21bc4ae0cc65)
and
![{\displaystyle p(J)=a_{0}+a_{1}\cdot J+\dots +a_{m}\cdot J^{m}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bc74d2c0b78f1d1c09abc7a282a17f88956d568c)
then
![{\displaystyle T_{p}(x)=p(J).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a0a37f33db5e904c7e7b23b0b2e8126ed159263e)
This is known as
Opitz' formula.
[2][3] Now consider increasing the degree of
to infinity, i.e. turn the Taylor polynomial into a Taylor series. Let
be a function which corresponds to a power series. You can compute the divided difference scheme for
by applying the corresponding matrix series to
: If
![{\displaystyle f(\xi )=\sum _{k=0}^{\infty }a_{k}\xi ^{k}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f04768f93c8aded7909131ede0083eb99ffb6843)
and
![{\displaystyle f(J)=\sum _{k=0}^{\infty }a_{k}J^{k}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1c2e79817d4f446ea10f10b231973919b57aa81a)
then
![{\displaystyle T_{f}(x)=f(J).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/78e638feeebc5d9fceffc758bad3ce9af30868f7)
Alternative characterizations
Expanded form
![{\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})}}+\\&\quad \quad {\frac {f(x_{2})}{(x_{2}-x_{0})\cdot (x_{2}-x_{1})\cdot (x_{2}-x_{3})}}+{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/14ec0849d6411a16b8e4d30ab46fa6465bc76928)
With the help of the polynomial function
this can be written as
![{\displaystyle f[x_{0},\dots ,x_{n}]=\sum _{j=0}^{n}{\frac {f(x_{j})}{\omega '(x_{j})}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/23345b9df9f4974173732a145c0555de91569917)
Peano form
If
and
, the divided differences can be expressed as[4]
![{\displaystyle f[x_{0},\ldots ,x_{n}]={\frac {1}{(n-1)!}}\int _{x_{0}}^{x_{n}}f^{(n)}(t)\;B_{n-1}(t)\,dt}](https://wikimedia.org/api/rest_v1/media/math/render/svg/71b2827e1bbfd4857051b4ea15edbaa5f49538b6)
where
![{\displaystyle f^{(n)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9dfb1963ccde0e87eb3838f51dc19041e2ff3816)
is the
![{\displaystyle n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a601995d55609f2d9f5e233e36fbe9ea26011b3b)
-th
derivative of the function
![{\displaystyle f}](https://wikimedia.org/api/rest_v1/media/math/render/svg/132e57acb643253e7810ee9702d9581f159a1c61)
and
![{\displaystyle B_{n-1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/35a38fdd982d39522072e328364b4ccec132f6dc)
is a certain
B-spline of degree
![{\displaystyle n-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fbd0b0f32b28f51962943ee9ede4fb34198a2521)
for the data points
![{\displaystyle x_{0},\dots ,x_{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/db64119863b01d965d94edaf33a384f1462f7adf)
, given by the formula
![{\displaystyle B_{n-1}(t)=\sum _{k=0}^{n}{\frac {(\max(0,x_{k}-t))^{n-1}}{\omega '(x_{k})}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f962c0435f9c41159ed2fa3a3bbc3d839a398422)
This is a consequence of the Peano kernel theorem; it is called the Peano form of the divided differences and
is the Peano kernel for the divided differences, all named after Giuseppe Peano.
Forward and backward differences
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
![{\displaystyle (x_{0},y_{0}),\ldots ,(x_{n},y_{n})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2b31b2f2436513ff671f30a751d6c5f6189b7ca5)
with
![{\displaystyle x_{k}=x_{0}+kh,\ {\text{ for }}\ k=0,\ldots ,n{\text{ and fixed }}h>0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a4dfd0c650c116ee1a216757cb6580575a72cb58)
the forward differences are defined as
![{\displaystyle {\begin{aligned}\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.\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c687ca3b8d960b102e901fa0f747b6261894b733)
whereas the backward differences are defined as:
![{\displaystyle {\begin{aligned}\nabla ^{(0)}y_{k}&:=y_{k},\qquad k=0,\ldots ,n\\\nabla ^{(j)}y_{k}&:=\nabla ^{(j-1)}y_{k}-\nabla ^{(j-1)}y_{k-1},\qquad k=0,\ldots ,n-j,\ j=1,\dots ,n.\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/57e142f7e965cfde8e4f08fa670dae0d2c42b573)
Thus the forward difference table is written as:
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aa01535bab7531c27c6b1f6648a55883f2493925)
whereas the backwards difference table is written as:
![{\displaystyle {\begin{matrix}y_{0}&&&\\&\nabla y_{1}&&\\y_{1}&&\nabla ^{2}y_{2}&\\&\nabla y_{2}&&\nabla ^{3}y_{3}\\y_{2}&&\nabla ^{2}y_{3}&\\&\nabla y_{3}&&\\y_{3}&&&\\\end{matrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/be0329cd07f4a6613e1fa9e90c1c0e537101f9bf)
The relationship between divided differences and forward differences is[5]
![{\displaystyle [y_{j},y_{j+1},\ldots ,y_{j+k}]={\frac {1}{k!h^{k}}}\Delta ^{(k)}y_{j},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/324e8711d27bb1c041f2bd268e916e7a65b2f68a)
whereas for backward differences:
[citation needed]![{\displaystyle [{y}_{j},y_{j-1},\ldots ,{y}_{j-k}]={\frac {1}{k!h^{k}}}\nabla ^{(k)}y_{j}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6bb191ad42a2ed1475b80c1e0f3e8b65deede658)
See also
References
- ^ Isaacson, Walter (2014). The Innovators. Simon & Schuster. p. 20. ISBN 978-1-4767-0869-0.
- ^ de Boor, Carl, Divided Differences, Surv. Approx. Theory 1 (2005), 46–69, [1]
- ^ Opitz, G. Steigungsmatrizen, Z. Angew. Math. Mech. (1964), 44, T52–T54
- ^ Skof, Fulvia (2011-04-30). Giuseppe Peano between Mathematics and Logic: Proceeding of the International Conference in honour of Giuseppe Peano on the 150th anniversary of his birth and the centennial of the Formulario Mathematico Torino (Italy) October 2-3, 2008. Springer Science & Business Media. p. 40. ISBN 978-88-470-1836-5.
- ^ Burden, Richard L.; Faires, J. Douglas (2011). Numerical Analysis (9th ed.). Cengage Learning. p. 129. ISBN 9780538733519.
- 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.