forward_euler.m

   1 function coefficients = forward_euler(integer_order, xs, x)
   2   ##
   3   ## Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
   4   ## Take for example a first order approximation, with,
   5   ##
   6   ##   xs = [x0,x1,x2,x3,x4]
   7   ##
   8   ##   f'(x1) ~= [f(x2)-f(x1)]/(x2-x1)
   9   ##
  10   ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids the
  11   ## solution of linear systems.
  12   ##
  13   ##
  14   ## INPUTS:
  15   ##
  16   ##   * ``integer_order`` - The order of the derivative which we're
  17   ##     approximating.
  18   ##
  19   ##   * ``xs`` - The vector of x-coordinates.
  20   ##
  21   ##   * ``x`` - The point `x` at which you'd like to evaluate the
  22   ##     derivative of the specified `integer_order`. This should be an
  23   ##     element of `xs`.
  24   ##
  25   ##
  26   ## OUTPUTS:
  27   ##
  28   ##   * ``coefficients`` - The vector of coefficients, in order, of
  29   ##     f(x0), f(x1), ..., f(xn).
  30   ##
  31
  32   if (integer_order < 0)
  33     ## You have made a grave mistake.
  34     coefficients = NA;
  35     return;
  36   end
  37
  38   if (integer_order == 0)
  39     coefficients = x;
  40     return;
  41   end
  42
  43   if (length(xs) < 2)
  44     ## You can't approximate a derivative of order greater than zero
  45     ## with zero or one points!
  46     coefficients = NA
  47     return;
  48   end
  49
  50   if (even(integer_order))
  51     offset_b = integer_order / 2;
  52     offset_f = offset_b;
  53   else
  54     ## When the order is odd, we need one more "forward" point than we
  55     ## do "backward" points.
  56     offset_b = (integer_order - 1) / 2;
  57     offset_f = offset_b + 1;
  58   end
  59
  60   ## Zero out the coefficients for terms that won't appear. We compute
  61   ## where `x` is, and we just computed how far back/forward we need to
  62   ## look from `x`, so we just need to make the rest zeros.
  63   x_idx = find(xs == x);
  64   first_nonzero_idx = x_idx - offset_b;
  65   last_nonzero_idx = x_idx + offset_f;
  66   leading_zero_count = first_nonzero_idx - 1;
  67   leading_zeros = zeros(1, leading_zero_count);
  68   trailing_zero_count = length(xs) - last_nonzero_idx;
  69   trailing_zeros = zeros(1, trailing_zero_count);
  70
  71   targets = xs(first_nonzero_idx : last_nonzero_idx);
  72
  73   # The multiplier comes from the Taylor expansion.
  74   multiplier = factorial(integer_order);
  75   cs = divided_difference_coefficients(targets) * multiplier;
  76
  77   coefficients = horzcat(leading_zeros, cs, trailing_zeros);
  78 end