-function df = forward_euler(integer_order, h, f, x)
+function coefficients = forward_euler(integer_order, xs, x)
+ ##
+ ## Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
+ ## Take for example a first order approximation, with,
+ ##
+ ## xs = [x0,x1,x2,x3,x4]
+ ##
+ ## f'(x1) ~= [f(x2)-f(x1)]/(x2-x1)
+ ##
+ ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids the
+ ## solution of linear systems.
##
- ## Use the forward Euler method to compute the derivative of `f` at
- ## a point `x`.
##
## INPUTS:
##
- ## * ``integer_order`` - The order of the derivative.
+ ## * ``integer_order`` - The order of the derivative which we're
+ ## approximating.
##
- ## * ``h`` - The step size.
+ ## * ``xs`` - The vector of x-coordinates.
##
- ## * ``f`` - The function whose derivative we're computing.
+ ## * ``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`.
##
- ## * ``x`` - The point at which to compute the derivative.
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``coefficients`` - The vector of coefficients, in order, of
+ ## f(x0), f(x1), ..., f(xn).
##
+ if (integer_order < 0)
+ ## You have made a grave mistake.
+ coefficients = NA;
+ return;
+ end
+
if (integer_order == 0)
- df = x;
+ coefficients = x;
+ return;
+ end
+
+ if (length(xs) < 2)
+ ## You can't approximate a derivative of order greater than zero
+ ## with zero or one points!
+ coefficients = NA
return;
end
-
- ## We need a few points around `x` to compute the derivative at `x`.
- ## The number of points depends on the order.
+
if (even(integer_order))
offset_b = integer_order / 2;
offset_f = offset_b;
offset_f = offset_b + 1;
end
- backward_xs = [x-(offset_b*h) : h : x];
+ ## Zero out the coefficients for terms that won't appear. We compute
+ ## where `x` is, and we just computed how far back/forward we need to
+ ## look from `x`, so we just need to make the rest zeros.
+ x_idx = find(xs == x);
+ first_nonzero_idx = x_idx - offset_b;
+ last_nonzero_idx = x_idx + offset_f;
+ leading_zero_count = first_nonzero_idx - 1;
+ leading_zeros = zeros(1, leading_zero_count);
+ trailing_zero_count = length(xs) - last_nonzero_idx;
+ trailing_zeros = zeros(1, trailing_zero_count);
- ## We'll always have at least one forward point, so start this vector
- ## from (x + h) and include `x` itself in the backward points.
- forward_xs = [x+h : h : x+(offset_f*h)];
+ targets = xs(first_nonzero_idx : last_nonzero_idx);
- xs = horzcat(backward_xs, forward_xs);
+ # The multiplier comes from the Taylor expansion.
+ multiplier = factorial(integer_order);
+ cs = divided_difference_coefficients(targets) * multiplier;
- ## Now that we have all of the points that we need in xs, we can use
- ## the divided_difference function to compute the derivative.
- df = divided_difference(f, xs);
+ coefficients = horzcat(leading_zeros, cs, trailing_zeros);
end