Add diffusion_matrix_sparse() and its tests.

[octave.git] / forward_euler.m

diff --git a/forward_euler.m b/forward_euler.m
index ec678a15c741ee29bb7aec50929e6330959a7adc..ff9c5eb2a498bc2be65d366fa9ba497d568354b7 100644 (file)
--- a/forward_euler.m
+++ b/forward_euler.m
@@ -5,13 +5,45 @@ function coefficients = forward_euler(integer_order, xs, x)
   ##
   ##   xs = [x0,x1,x2,x3,x4]
   ##
-  ##   f'(x=x1) ~= [f(x2)-f(x1)]/(x2-x1)
+  ##   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.
+  ## 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.
+    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
 
@@ -37,7 +69,10 @@ function coefficients = forward_euler(integer_order, xs, x)
   trailing_zeros = zeros(1, trailing_zero_count);
 
   targets = xs(first_nonzero_idx : last_nonzero_idx);
-  cs = divided_difference_coefficients(targets);
 
-  coefficients = horzcat(leading_zeros, cs, trailing_zeros);  
+  # 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