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. coefficients = NA; return; end if (integer_order == 0) 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 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