Add the central_difference function and a test for it.

diff --git a/central_difference.m b/central_difference.m
new file mode 100644 (file)
index 0000000..ad58fd4
--- /dev/null
+++ b/central_difference.m
@@ -0,0 +1,51 @@
+function coefficients = central_difference(xs, x)
+  ##
+  ## The first order central difference at x1 is,
+  ##
+  ##   f'(x1) = (f(x2) - f(x0))/2
+  ##
+  ## where the index x1 is of course arbitrary but x2, x0 are adjacent
+  ## to x1. The coefficients we seek are the coefficients of f(xj) for
+  ## j = 1,...,N-2, where N is the length of ``xs``. We omit the first
+  ## and last coefficient because at x0 and xN, the previous/next
+  ## value is not available.
+  ##
+  ## This should probably take an 'order' parameter as well; see
+  ## forward_euler().
+  ##
+  ## INPUT:
+  ##
+  ##   * ``xs`` - The vector of x-coordinates.
+  ##
+  ##   * ``x`` - The point `x` at which you'd like to evaluate the
+  ##   derivative of the specified `integer_order`. This should be an
+  ##   element of `xs`.
+  ##
+  ## OUTPUT:
+  ##
+  ##   * ``coefficients`` - The vector of coefficients, in order, of
+  ##   f(x0), f(x1), ..., f(xn).
+  ##
+
+  if (length(xs) < 3)
+    ## We need at least one point other than the first and last.
+    coefficients = NA;
+    return;
+  end
+
+  x_idx = find(xs == x);
+
+  if (x_idx == 1 || x_idx == length(xs))
+    ## You asked for the difference at the first or last element, which
+    ## we can't do.
+    coefficients = NA;
+    return;
+  end
+
+  ## Start with a vector of zeros.
+  coefficients = zeros(1, length(xs));
+
+  ## And fill in the two values that we know.
+  coefficients(x_idx - 1) = -1/2;
+  coefficients(x_idx + 1) = 1/2;
+end

diff --git a/tests/central_difference_tests.m b/tests/central_difference_tests.m
new file mode 100644 (file)
index 0000000..017a7ad
--- /dev/null
+++ b/tests/central_difference_tests.m
@@ -0,0 +1,7 @@
+xs = [0, 0.5, 1];
+expected = [-0.5, 0, 0.5];
+actual = central_difference(xs, 0.5);
+
+unit_test_equals("Three point central difference works", ...
+                expected, ...
+                actual);