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.35.1