From: Michael Orlitzky Date: Sat, 15 Sep 2012 03:51:30 +0000 (-0400) Subject: Rewrite the forward_euler() function to return a vector of coefficients rather than... X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=ba154db78b02163be343b235438e3add73cba526;p=octave.git Rewrite the forward_euler() function to return a vector of coefficients rather than assuming that we will really have a function to evaluate. --- diff --git a/forward_euler.m b/forward_euler.m index c6ef217..ec678a1 100644 --- a/forward_euler.m +++ b/forward_euler.m @@ -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; + leading_zeros = zeros(1, leading_zero_count); + trailing_zero_count = length(xs) - last_nonzero_idx; + trailing_zeros = zeros(1, trailing_zero_count); - xs = horzcat(backward_xs, forward_xs); + targets = xs(first_nonzero_idx : last_nonzero_idx); + cs = divided_difference_coefficients(targets); - ## 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