divided_difference.m

   1 function dd = divided_difference(f, xs)
   2   ## Compute divided difference of `f` at points `xs`. The argument `xs`
   3   ## is assumed to be a vector containing at least one element. If it
   4   ## contains n elements, the (n-1)st divided difference will be
   5   ## calculated.
   6   ##
   7   ## INPUTS:
   8   ##
   9   ##   * ``f`` - The function whose divided differences we want.
  10   ##
  11   ##   * ``xs`` - A vector containing x-coordinates. The length of `xs`
  12   ##   determines the order of the divided difference.
  13   ##
  14   ##
  15   ## OUTPUTS:
  16   ##
  17   ##   * ``dd`` - The divided difference f[xs(1), xs(2),...]
  18   ##
  19   order = length(xs) - 1;
  20
  21   if (order < 0)
  22     ## Can't do anything here. Return nothing.
  23     dd = NA;
  24   elseif (order == 0)
  25     ## Our base case.
  26     dd = f(xs(1))
  27   else
  28     ## Order >= 1, recurse.
  29
  30     ## f[x0,...,x_n-1]
  31     f0 = divided_difference(f, xs(1:end-1));
  32     ## f[x1,...,x_n]
  33     f1 = divided_difference(f, xs(2:end));
  34
  35     # http://mathworld.wolfram.com/DividedDifference.html
  36     dd = (f0 - f1)/(xs(1) - xs(end));
  37   end
  38 end