divided_difference.m

   1 #!/usr/bin/octave --silent
   2
   3 function dd = divided_difference(f, xs)
   4   ## Compute divided difference of `f` at points `xs`. The argument `xs`
   5   ## is assumed to be a vector containing at least one element. If it
   6   ## contains n elements, the (n-1)st divided difference will be
   7   ## calculated.
   8   order = length(xs);
   9
  10   if (order < 1)
  11     ## Can't do anything here. Return nothing.
  12     dd = NA;
  13   elseif (order == 1)
  14     ## Our base case.
  15     dd = f(xs(1))
  16   else
  17     ## Order > 1, recurse.
  18
  19     ## f[x0,...,x_n-1]
  20     f0 = divided_difference(f, xs(1:end-1));
  21     ## f[x1,...,x_n]
  22     f1 = divided_difference(f, xs(2:end));
  23
  24     # http://mathworld.wolfram.com/DividedDifference.html
  25     dd = (f0 - f1)/(xs(1) - xs(end))
  26   end
  27 end