5
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 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, 2019 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, 2019 at 12:20
  • 2
    $\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, 2019 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, 2019 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.'s lack of A.I. ♦
    Jul 28, 2019 at 2:50

2 Answers 2

Reset to default
7
$\begingroup$

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:

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.

Improve this answer
$\endgroup$
4
$\begingroup$ 
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}}

enter image description here

Improve this answer
$\endgroup$
2
  • 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, 2019 at 19:59
  • $\begingroup$ Thank you @MarkR. Fixed it. $\endgroup$
    – kglr
    Jul 27, 2019 at 20:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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