central_difference.m

   1 function coefficients = central_difference(xs, x)
   2   ##
   3   ## The first order central difference at x1 is,
   4   ##
   5   ##   f'(x1) = (f(x2) - f(x0))/2
   6   ##
   7   ## where the index x1 is of course arbitrary but x2, x0 are adjacent
   8   ## to x1. The coefficients we seek are the coefficients of f(xj) for
   9   ## j = 1,...,N-2, where N is the length of ``xs``. We omit the first
  10   ## and last coefficient because at x0 and xN, the previous/next
  11   ## value is not available.
  12   ##
  13   ## This should probably take an 'order' parameter as well; see
  14   ## forward_euler().
  15   ##
  16   ## INPUT:
  17   ##
  18   ##   * ``xs`` - The vector of x-coordinates.
  19   ##
  20   ##   * ``x`` - The point `x` at which you'd like to evaluate the
  21   ##   derivative of the specified `integer_order`. This should be an
  22   ##   element of `xs`.
  23   ##
  24   ## OUTPUT:
  25   ##
  26   ##   * ``coefficients`` - The vector of coefficients, in order, of
  27   ##   f(x0), f(x1), ..., f(xn).
  28   ##
  29
  30   if (length(xs) < 3)
  31     ## We need at least one point other than the first and last.
  32     coefficients = NA;
  33     return;
  34   end
  35
  36   x_idx = find(xs == x);
  37
  38   if (x_idx == 1 || x_idx == length(xs))
  39     ## You asked for the difference at the first or last element, which
  40     ## we can't do.
  41     coefficients = NA;
  42     return;
  43   end
  44
  45   ## Start with a vector of zeros.
  46   coefficients = zeros(1, length(xs));
  47
  48   ## And fill in the two values that we know.
  49   coefficients(x_idx - 1) = -1/2;
  50   coefficients(x_idx + 1) = 1/2;
  51 end