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'(x=x1) ~= [f(x2)-f(x1)]/(x2-x1)
   9   ##
  10   ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids
  11   ## the solution of linear systems.
  12   ##
  13   if (integer_order == 0)
  14     df = x;
  15     return;
  16   end
  17
  18   if (even(integer_order))
  19     offset_b = integer_order / 2;
  20     offset_f = offset_b;
  21   else
  22     ## When the order is odd, we need one more "forward" point than we
  23     ## do "backward" points.
  24     offset_b = (integer_order - 1) / 2;
  25     offset_f = offset_b + 1;
  26   end
  27
  28   ## Zero out the coefficients for terms that won't appear. We compute
  29   ## where `x` is, and we just computed how far back/forward we need to
  30   ## look from `x`, so we just need to make the rest zeros.
  31   x_idx = find(xs == x);
  32   first_nonzero_idx = x_idx - offset_b;
  33   last_nonzero_idx = x_idx + offset_f;
  34   leading_zero_count = first_nonzero_idx - 1;
  35   leading_zeros = zeros(1, leading_zero_count);
  36   trailing_zero_count = length(xs) - last_nonzero_idx;
  37   trailing_zeros = zeros(1, trailing_zero_count);
  38
  39   targets = xs(first_nonzero_idx : last_nonzero_idx);
  40   cs = divided_difference_coefficients(targets);
  41
  42   coefficients = horzcat(leading_zeros, cs, trailing_zeros);
  43 end