author Michael Orlitzky Sat, 15 Sep 2012 03:51:30 +0000 (23:51 -0400) committer Michael Orlitzky Sat, 15 Sep 2012 03:51:30 +0000 (23:51 -0400)
 forward_euler.m patch | blob | history

@@ -1,26 +1,20 @@
-function df = forward_euler(integer_order, h, f, x)
+function coefficients = forward_euler(integer_order, xs, x)
##
-  ## Use the forward Euler method to compute the derivative of `f` at
-  ## a point `x`.
+  ## Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
+  ## Take for example a first order approximation, with,
##
-  ## INPUTS:
+  ##   xs = [x0,x1,x2,x3,x4]
##
-  ##   * ``integer_order`` - The order of the derivative.
+  ##   f'(x=x1) ~= [f(x2)-f(x1)]/(x2-x1)
##
-  ##   * ``h`` - The step size.
+  ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids
+  ## the solution of linear systems.
##
-  ##   * ``f`` - The function whose derivative we're computing.
-  ##
-  ##   * ``x`` - The point at which to compute the derivative.
-  ##
-
if (integer_order == 0)
df = x;
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;
@@ -31,15 +25,19 @@ function df = forward_euler(integer_order, h, f, x)
offset_f = offset_b + 1;
end

-  backward_xs = [x-(offset_b*h) : h : x];
-
-  ## 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)];
+  ## 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;