Trying to build a function to compute divided difference for arbitrary list of points

As the title says Im trying to build a function to compute divided difference for arbitrary list of points.

The divided difference of a set of points $$\{(x_1,y_1),\ldots, (x_n,y_n)\}$$ where $$x_j is defined recursively by

$$f[i,j]:=\frac{f[i+1,j]-f[i,j-1]}{x_j-x_i},\, f[i,i]:=y_i$$

for $$j\ge i$$.I tried to mimic this definition with the code

f[i_, j_, m_] := f[i, j, m] = (f[i + 1, j, m] - f[i, j - 1, m])/(
m[[j, 1]] - m[[i, 1]]) Boole[j > i] + m[[j, 2]] Boole[j == i]


where m is an arbitrary $$n\times 2$$ matrix. However I tried to see if this works calling f[1,1,{{1,2},{3,4}}]but it gives me a recursion limit error. How I can fix the above code? Thank you in advance.

• In recursive definitions, you need to specify a one or more stopping definitions. E.g. for the factorial function f[0] = 1; f[n_] := n f[n - 1] from the docs almost works, but better is something like f[0] = 1; f[n_Integer?Positive] := n f[n - 1] or f[0] = 1; f[n_Integer] /; n > 0 := n f[n - 1] – Michael E2 Jul 27 '19 at 12:11
• @MichaelE2 but these stops are already defined in the function, when you can something where $j=i$ – Masacroso Jul 27 '19 at 12:20
• FWIW, here's a procedural def. based on the algorithm in Atkinson, Intro. to Num. Anal.: dividedDifference[x_, f_] := Module[{d}, With[{n = Length[(x)]}, d = f; (* div. diff. computed in place *) Do[ d[[j]] = (d[[j]] - d[[j - 1]])/(x[[j]] - x[[j - k + 1]]), {k, 2, n}, {j, n, k, -1}]; d ]] -- x is a list of nodes and f a list of values at the nodes. For input m, the call should be dividedDifference @@ Transpose@m. – Michael E2 Jul 27 '19 at 12:23
• Masacroso, you show only one definition rule for f[i_, j_, m_], which calls f[i + 1, j, m] and f[i, j - 1, m], both of which in turn call f on new arguments ad infinitum. Unless there are definition rules that you are not showing, there are no stops to this process. The calls to f will occur even if Boole[] evaluates to zero, in case you think Boole[] might stop the process. It would be possible to put the stops in using an If[] statement. – Michael E2 Jul 27 '19 at 12:28
• @Michael's implementation is in fact a triangular recursion, and is really the most efficient method to build up the divided difference table. – J. M.'s discontentment Jul 28 '19 at 2:50

Well, instead of evaluating an individual entry as $$f_{ij}$$ once at a time, here I provide an implementation generating the whole upper-triangular $$f$$ matrix. So here we go:

Clear[fmat]
fmat[pts : {{_, _} ..}] := Module[{n = Length[pts], kernels, xs, ys, denomenators, entries},
{xs, ys} = pts\[Transpose];
kernels = NestList[Insert[#, 0, 2] &, {-1, 1}, n - 2];
denomenators = ListCorrelate[#, xs] & /@ kernels;
entries = FoldList[Differences[#]/#2 &, ys, denomenators];
SparseArray[MapIndexed[Band[Prepend[#2, 1]] -> # &, entries], n]
]


To check its validity, one can see a symbolic result and compare it with that on the Wiki page:

Then below calculation is trivial

Also, I have tested on lists of random reals with a larger shape, say $$1000\times 2$$. The efficiency of generating the $$1000\times 1000$$ $$f$$ matrix is acceptable, at least.

P.S. I learnt that the calculation of divided difference is related to Newton interpolation. So if your purpose is to do so interpolations, you can directly use built-in functions like Interpolation, InterpolatingPolynomial, etc.

ClearAll[dividedDif]
dividedDif = NestList[#/(#[[1, 1]] + #[[-1, -1]] &[Denominator /@ (List @@ #)]) & /@
Differences[#] &, Divide @@ (Differences /@ #), Length[#[[1]]] - 2][[All, 1 ]] &;

Prepend[dividedDif @ #, #[[1, 1]]] & @ {{r, s, t, u, v}, {a, b, c, d, e}}


• Nice result, but I think you meant to write "dividedDif" on the last line with the Prepend. When I do that, I get the answer you show. Oddly, if you use the values {{1,2},{3,4}}, you don't get what was in the previous answer although I might have input that incorrectly. – Mark R Jul 27 '19 at 19:59
• Thank you @MarkR. Fixed it. – kglr Jul 27 '19 at 20:00