From 9ecd26a43862541b1c341cb318472d56ad0b382b Mon Sep 17 00:00:00 2001
From: Michael Orlitzky
Date: Wed, 6 Feb 2013 14:47:19 -0500
Subject: [PATCH] Add the central_difference function and a test for it.
---
central_difference.m | 51 ++++++++++++++++++++++++++++++++
tests/central_difference_tests.m | 7 +++++
2 files changed, 58 insertions(+)
create mode 100644 central_difference.m
create mode 100644 tests/central_difference_tests.m
diff --git a/central_difference.m b/central_difference.m
new file mode 100644
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
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);
--
2.33.1