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<x_{j+1}$ 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.

  • 1
    $\begingroup$ 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] $\endgroup$ – Michael E2 Jul 27 at 12:11
  • $\begingroup$ @MichaelE2 but these stops are already defined in the function, when you can something where $j=i$ $\endgroup$ – Masacroso Jul 27 at 12:20
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Jul 27 at 12:23
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Jul 27 at 12:28
  • 1
    $\begingroup$ @Michael's implementation is in fact a triangular recursion, and is really the most efficient method to build up the divided difference table. $\endgroup$ – J. M. will be back soon Jul 28 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:

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:

enter image description here

Then below calculation is trivial

enter image description here

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.

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}}

enter image description here

  • 1
    $\begingroup$ 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. $\endgroup$ – Mark R Jul 27 at 19:59
  • $\begingroup$ Thank you @MarkR. Fixed it. $\endgroup$ – kglr Jul 27 at 20:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.