--- /dev/null
+function coefficients = divided_difference_coefficients(xs)
+ ## Compute divided difference coefficients of `f` at points `xs`.
+ ##
+ ## INPUTS:
+ ##
+ ## * ``xs`` - A vector containing x-coordinates.
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``coefficients`` - The vector of coefficients such that
+ ## dot(coefficients, f(xs)) == f[xs]. Used to solve linear systems.
+ ##
+
+ coefficients = [];
+
+ for xj = xs
+ this_coeff = 1;
+ for xi = xs
+ if (xi != xj)
+ this_coeff = this_coeff * (1 / (xj - xi));
+ end
+ end
+ coefficients(end+1) = this_coeff;
+ end
+end
--- /dev/null
+function even = even(integer_n)
+ ## Returns true if its argument is even; false otherwise.
+ ##
+ ## INPUTS:
+ ##
+ ## * ``integer_n`` - The integer whose parity you're determining.
+ ##
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``even`` - True if `integer_n` is even, false otherwise.
+ ##
+ even = rem(integer_n, 2) == 0;
+end
--- /dev/null
+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
--- /dev/null
+function odd = odd(integer_n)
+ ## Returns true if its argument is odd; false otherwise.
+ ##
+ ## INPUTS:
+ ##
+ ## * ``integer_n`` - The integer whose parity you're determining.
+ ##
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``odd`` - True if `integer_n` is odd, false otherwise.
+ ##
+
+ odd = !even(integer_n);
+end
--- /dev/null
+function [p,delta] = partition(integerN, a, b)
+ ## Partition the interval [a,b] into integerN subintervals. We do not
+ ## requite that a<b.
+ ##
+ ## INPUTS:
+ ##
+ ## * integerN - The number of subintervals.
+ ##
+ ## * a - The "left" endpoint of the interval to partition.
+ ##
+ ## * b - The "right" endpoint of the interval to partition.
+ ##
+ ##
+ ## OUTPUTS:
+ ##
+ ## * p - The resulting partition, as a vector of length integerN+1.
+ ##
+ ## * delta - The distance between x_i and x_{i+1} in the partition.
+ ##
+ ##
+
+ ## We don't use abs() here because `b` might be less than `a`. In that
+ ## case, we want delta negative so that when we add it to `a`, we move
+ ## towards `b`.
+ delta = (b - a)/integerN;
+
+ p = [a : delta : b];
+end
unit_test_equals("sin[0, pi, 2*pi] == 0", ...
0, ...
divided_difference(@sin, [0,pi,2*pi]));
+
+unit_test_equals("zero order divided_difference_coefficients", ...
+ [1], ...
+ divided_difference_coefficients([0]));
+
+unit_test_equals("first order divided_difference_coefficients", ...
+ [-1, 1] / pi, ...
+ divided_difference_coefficients([0, pi]));
+
+unit_test_equals("second order divided_difference_coefficients", ...
+ [1, -2, 1] / (2*pi^2), ...
+ divided_difference_coefficients([0, pi, 2*pi]));
+
+
+unit_test_equals("1 is odd", ...
+ true, ...
+ odd(1));
+
+unit_test_equals("1 is not even", ...
+ false, ...
+ even(1));
+
+unit_test_equals("2 is not odd", ...
+ false, ...
+ odd(2));
+
+unit_test_equals("2 is even", ...
+ true, ...
+ even(2));
projects / octave.git / commitdiff