Add the divided_difference Octave function.

[octave.git] / divided_difference.m

diff --git a/divided_difference.m b/divided_difference.m
new file mode 100755 (executable)
index 0000000..139fa47
--- /dev/null
+++ b/divided_difference.m
@@ -0,0 +1,27 @@
+#!/usr/bin/octave --silent
+
+function dd = divided_difference(f, xs)
+  ## Compute divided difference of `f` at points `xs`. The argument `xs`
+  ## is assumed to be a vector containing at least one element. If it
+  ## contains n elements, the (n-1)st divided difference will be
+  ## calculated.
+  order = length(xs);
+  
+  if (order < 1)
+    ## Can't do anything here. Return nothing.
+    dd = NA;
+  elseif (order == 1)
+    ## Our base case.
+    dd = f(xs(1))
+  else
+    ## Order > 1, recurse.
+
+    ## f[x0,...,x_n-1]
+    f0 = divided_difference(f, xs(1:end-1));
+    ## f[x1,...,x_n]
+    f1 = divided_difference(f, xs(2:end));
+
+    # http://mathworld.wolfram.com/DividedDifference.html
+    dd = (f0 - f1)/(xs(1) - xs(end))
+  end
+end