# How to write a custom function of m-order difference? [duplicate]

In the text, the formula of calculating m-order difference with m + 1 nodes is as follows:

$$f\left[x_{k}, x_{k+1}\right]=\frac{f\left(x_{k+1}\right)-f\left(x_{k}\right)}{x_{k+1}-x_{k}}$$ $$f\left[x_{k}, x_{k+1}, \cdots, x_{k+m}\right]=\frac{f\left[x_{k+1}, \cdots, x_{k+m}\right]-f\left[x_{k}, \cdots, x_{k+m-1}\right]}{x_{k+m}-x_{k}}$$

The Newton interpolation formula with n+1 nodes can be obtained by using the above formula:

\begin{aligned} N_{n}(x)=& f\left(x_{0}\right)+f\left[x_{0}, x_{1}\right]\left(x-x_{0}\right)+f\left[x_{0}, x_{1}, x_{2}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots \\ &+f\left[x_{0}, x_{1}, \cdots, x_{n}\right]\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-1}\right) \end{aligned}

But the built-in function DifferenceQuotient can only calculate high-order difference for equal step h:

DifferenceQuotient[f[x], {x, 3, h}]


What can I do to write an ingenious m-order difference function to calculate Newton's interpolation polynomials?

Update: I tried to write a custom function Newton to solve this problem, I hope I can get other more ingenious methods $$\color{Gray} {\text{(武汉 岩石 数值分析-2007-2)}}$$:

ClearAll["*"]
f[{x_, y_}] := (f[y] - f[x])/(y - x)
f[x_List] := (f[Drop[x, {Length[x] - 1}]] -
f[Drop[x, {Length[x]}]])/(x[[Length[x]]] - x[[Length[x] - 1]])
ω[n_, s_] := Product[(s - xx[[i]]), {i, 1, n}]
xx = {0., 1., 2., 3.};
y = {2., 3., 0., -1.};
Evaluate[(f /@ xx)] = y;
Newton[n_?(# > 0 && IntegerQ[#] &), s_ : x] :=
f[xx[[1]]] +
Sum[f[Table[xx[[i]], {i, 1, m + 1}]]*ω[m, s], {m, 1, n}]

Newton[3, x] // Expand


Note: These definitions used are from page 99 of this book.

Or use the definition on page 5 of this book:

• What have you tried? – xzczd Sep 17 '20 at 3:45
• @xzczd Thank you. I have updated the question. – A little mouse on the pampas Sep 17 '20 at 5:02
• Your implementation is wrong, please double check it. – xzczd Sep 17 '20 at 5:46
• @xzczd Thank you for pointing out the mistake. I have added some necessary information. Equation (3.3) and equation (3.3') are equivalent. – A little mouse on the pampas Sep 17 '20 at 6:10
• I see, I've made a mistake in testing. – xzczd Sep 17 '20 at 6:42

## 1 Answer

newton[{X_, Y_}, x_] :=
Module[{f, asso = X -> Y // Thread // Association},
f@{a_, b___, c_} := f@{a, b, c} = (f@{b, c} - f@{a, b})/(c - a);
f@{a_} := asso@a;
Sum[f@# Times @@ (x - Most@#) &@X[[;; i]], {i, Length@X}]]

xx = {0., 1., 2., 3., 5};
y = {2., 3., 0., -1., -3};
newton[{xx, y}, x] // Expand
(* 2. + 6.5 x - 7.75 x^2 + 2.5 x^3 - 0.25 x^4 *)
`