X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=forward_euler.m;fp=forward_euler.m;h=7d7e12de4b8aeb06b29e0f48b236103ff91c50b6;hp=0000000000000000000000000000000000000000;hb=876e38fb99680748dfdf334ba450f633566d9b6a;hpb=62d652799ded51169bda744d8728e1d33582fa5f

diff --git a/forward_euler.m b/forward_euler.m
new file mode 100644
index 0000000..7d7e12d
--- /dev/null
+++ b/forward_euler.m
@@ -0,0 +1,78 @@
+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.
+  ##
+  ##
+  ## INPUTS:
+  ##
+  ##   * ``integer_order`` - The order of the derivative which we're
+  ##     approximating.
+  ##
+  ##   * ``xs`` - The vector of x-coordinates.
+  ##
+  ##   * ``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`.
+  ##
+  ##
+  ## 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.
+    df = NA;
+    return;
+  end
+
+  if (integer_order == 0)
+    df = x;
+    return;
+  end
+
+  if (length(xs) < 2)
+    ## You can't approximate a derivative of order greater than zero
+    ## with zero or one points!
+    df = NA
+    return;
+  end
+
+  if (even(integer_order))
+    offset_b = integer_order / 2;
+    offset_f = offset_b;
+  else
+    ## When the order is odd, we need one more "forward" point than we
+    ## do "backward" points.
+    offset_b = (integer_order - 1) / 2;
+    offset_f = offset_b + 1;
+  end
+
+  ## 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);
+
+  targets = xs(first_nonzero_idx : last_nonzero_idx);
+
+  # The multiplier comes from the Taylor expansion.
+  multiplier = factorial(integer_order);
+  cs = divided_difference_coefficients(targets) * multiplier;
+
+  coefficients = horzcat(leading_zeros, cs, trailing_zeros);
+end