From 876e38fb99680748dfdf334ba450f633566d9b6a Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 15 Oct 2012 14:47:03 -0400 Subject: [PATCH] Move several functions out of the homework1/src directory and into the top-level octave directory. --- divided_difference_coefficients.m | 25 ++++++++++ even.m | 14 ++++++ forward_euler.m | 78 +++++++++++++++++++++++++++++++ odd.m | 15 ++++++ partition.m | 28 +++++++++++ run-tests.m | 29 ++++++++++++ 6 files changed, 189 insertions(+) create mode 100644 divided_difference_coefficients.m create mode 100644 even.m create mode 100644 forward_euler.m create mode 100644 odd.m create mode 100644 partition.m diff --git a/divided_difference_coefficients.m b/divided_difference_coefficients.m new file mode 100644 index 0000000..a40f559 --- /dev/null +++ b/divided_difference_coefficients.m @@ -0,0 +1,25 @@ +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 diff --git a/even.m b/even.m new file mode 100644 index 0000000..f4f3767 --- /dev/null +++ b/even.m @@ -0,0 +1,14 @@ +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 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 diff --git a/odd.m b/odd.m new file mode 100644 index 0000000..b0bc626 --- /dev/null +++ b/odd.m @@ -0,0 +1,15 @@ +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 diff --git a/partition.m b/partition.m new file mode 100644 index 0000000..96ded5e --- /dev/null +++ b/partition.m @@ -0,0 +1,28 @@ +function [p,delta] = partition(integerN, a, b) + ## Partition the interval [a,b] into integerN subintervals. We do not + ## requite that a