Move homework #1's octave code into its directory.

Add the homework #1 directory to the divided_difference() load path.

diff --git a/bisect.m b/bisect.m
deleted file mode 100644 (file)
index 3570e57..0000000
--- a/bisect.m
+++ /dev/null
@@ -1,68 +0,0 @@
-function [root,iterations] = bisect(f, a, b, epsilon)
-  ## Find a root of the function `f` on the closed interval [a,b].
-  ##
-  ## It is assumed that there are an odd number of roots within [a,b].
-  ## The intermediate value theorem is used to determine whether or not
-  ## there is a root on a given interval. So, for example, cos(-2) and
-  ## cos(2) are both less than zero; the I.V.T. would not conclude that
-  ## cos(x) has a root on -2<=x<=2.
-  ##
-  ## If `f` has more than one root on [a,b], the smallest root will be
-  ## returned. This is an implementation detail.
-  ##
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``f`` - The function whose root we seek.
-  ##
-  ##   * ``a`` - The "left" endpoint of the interval in which we search.
-  ##
-  ##   * ``b`` - The "right" endpoint of the interval in which we search.
-  ##
-  ##   * ``epsilon`` - How close we should be to the actual root before we
-  ##     halt the search and return the current approximation.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``root`` - The root that we found.
-  ##
-  ##   * ``iterations`` - The number of bisections that we performed
-  ##     during the search.
-  ##
-
-  ## Store these so we don't have to recompute them.
-  fa = f(a);
-  fb = f(b);
-
-  ## We first check for roots at a,b before bothering to bisect the
-  ## interval.  
-  if (fa == 0)
-    root = a;
-    iterations = 0;
-    return;
-  end
-
-  if (fb == 0)
-    root = b;
-    iterations = 0;
-    return;
-  end
-
-  ## Bisect the interval.
-  c = (a+b)/2;
-  iterations = 1;
-  
-  ## If we're within the prescribed tolerance, we're done.
-  if (b-c < epsilon)
-    root = c;
-    return;
-  end
-
-  if (has_root(fa,f(c)))
-    [root,subiterations] = bisect(f,a,c,epsilon);
-    iterations = iterations + subiterations;
-  else
-    [root,subiterations] = bisect(f,c,b,epsilon);
-    iterations = iterations + subiterations;
-  end
-end

diff --git a/divided_difference.m b/divided_difference.m
index a9039b10c042caa773c4b162f0d6ee94e95734fe..785b66f27e812ff0d908bdbf317b05077e15d694 100644 (file)
--- a/divided_difference.m
+++ b/divided_difference.m
@@ -15,7 +15,8 @@ function dd = divided_difference(f, xs)
   ## OUTPUTS:
   ##
   ##   * ``dd`` - The divided difference f[xs(1), xs(2),...]
   ## OUTPUTS:
   ##
   ##   * ``dd`` - The divided difference f[xs(1), xs(2),...]

-  ## 
+  ##
+  addpath('../homework1/src');

   order = length(xs) - 1;
   
   if (order < 0)
   order = length(xs) - 1;
   
   if (order < 0)


@@ -23,10 +24,10 @@ function dd = divided_difference(f, xs)
     dd = NA;
   elseif (order == 0)
     ## Our base case.
     dd = NA;
   elseif (order == 0)
     ## Our base case.

-    dd = f(xs(1))
+    dd = f(xs(1));

   else
     ## Order >= 1.
   else
     ## Order >= 1.

-    cs = divided_difference_coefficients(xs)
+    cs = divided_difference_coefficients(xs);

     dd = dot(cs, f(xs));
   end
 end
     dd = dot(cs, f(xs));
   end
 end

diff --git a/divided_difference_coefficients.m b/divided_difference_coefficients.m
deleted file mode 100644 (file)
index 53d5799..0000000
--- a/divided_difference_coefficients.m
+++ /dev/null
@@ -1,26 +0,0 @@
-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)) == dd. Used to solve linear systems.
-  ##
-
-  coefficients = [];
-    
-  for xj = xs
-    this_coeff = 1;
-    for xi = xs
-      if (xi != xj)
-       ## Append (xj - xi) to the vector of coefficients.
-       this_coeff = this_coeff * (1 / (xj - xi));
-      end
-    end
-    coefficients(end+1) = this_coeff;
-  end
-end

diff --git a/even.m b/even.m
deleted file mode 100644 (file)
index f4f3767..0000000
--- a/even.m
+++ /dev/null
@@ -1,14 +0,0 @@
-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
deleted file mode 100644 (file)
index ec678a1..0000000
--- a/forward_euler.m
+++ /dev/null
@@ -1,43 +0,0 @@
-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'(x=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.
-  ##
-  if (integer_order == 0)
-    df = x;
-    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);
-  cs = divided_difference_coefficients(targets);
-
-  coefficients = horzcat(leading_zeros, cs, trailing_zeros);  
-end

diff --git a/has_root.m b/has_root.m
deleted file mode 100644 (file)
index ffd9666..0000000
--- a/has_root.m
+++ /dev/null
@@ -1,34 +0,0 @@
-function has_root = has_root(fa, fb)
-  ## Use the intermediate value theorem to determine whether or not some
-  ## function has an odd number of roots on an interval. If the function
-  ## in question has an even number of roots, the result will be
-  ## incorrect.
-  ##
-  ## Call the function whose roots we're concerned with 'f'. The two
-  ## parameters `fa` and `fb` should correspond to f(a) and f(b).
-  ##
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``fa`` - The value of `f` at one end of the interval.
-  ##
-  ##   * ``fb`` - The value of `f` at the other end of the interval.
-  ##
-  ## OUTPUTS:
-  ##
-  ##  * ``has_root`` - True if we can use the I.V.T. to conclude that
-  ##    there is a root on [a,b], false otherwise.
-  ##
-  
-  ## If either f(a) or f(b) is zero, the product of their signs will be
-  ## zero and either a or b is a root. If the product of their signs is
-  ## negative, then f(a) and f(b) are non-zero and have opposite sign,
-  ## so there must be a root on (a,b). The only case we don't want is
-  ## when f(a) and f(b) have the same sign; in this case, the product of
-  ## their signs would be one.
-  if (sign(fa) * sign(fb) != 1)
-    has_root = true;
-  else
-    has_root = false;
-  end
-end

diff --git a/odd.m b/odd.m
deleted file mode 100644 (file)
index b0bc626..0000000
--- a/odd.m
+++ /dev/null
@@ -1,15 +0,0 @@
-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
deleted file mode 100644 (file)
index 96ded5e..0000000
--- a/partition.m
+++ /dev/null
@@ -1,28 +0,0 @@
-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

diff --git a/solve_poisson.m b/solve_poisson.m
deleted file mode 100644 (file)
index 946d1d8..0000000
--- a/solve_poisson.m
+++ /dev/null
@@ -1,61 +0,0 @@
-function u = solve_poisson(integerN, f)
-  ##
-  ## Numerically solve the poisson equation,
-  ##
-  ##   -u_xx(x) = f(x)
-  ##
-  ## in one dimension, subject to the boundary conditions,
-  ##
-  ##   u(0) = 0
-  ##   u(1) = 1
-  ##
-  ## over the interval [0,1]. It is assumed that the function `f` is at
-  ## least continuous on [0,1].
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``integerN`` - An integer representing the number of subintervals
-  ##     we should use to approximate `u`.
-  ##
-  ##   * ``f`` - The function on the right hand side of the poisson
-  ##     equation.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * u - The function `u` that approximately solves -u_xx(x) = f(x).
-  ##
-
-  [xs,h] = partition(integerN, 0, 1);
-
-  ## We cannot evaluate u_xx at the endpoints because our
-  ## differentiation algorithm relies on the points directly to the left
-  ## and right of `x`.
-  differentiable_points = xs(2:end-1)
-
-  ## Prepend and append a zero row to f(xs).
-  ## Oh, and negate the whole thing (see the formula).
-  F = - cat(2, 0, f(differentiable_points), 0)
-
-  ## These are the coefficient vectors for the u(x0) and u(xn)
-  ## constraints. There should be N zeros and a single 1.
-  the_rest_zeros = zeros(1, integerN);
-  u_x0_coeffs = cat(2, 1, the_rest_zeros);
-  u_x4_coeffs = cat(2, the_rest_zeros, 1);
-
-  ## Start with the u(x0) row.
-  A = u_x0_coeffs;
-
-  for x = differentiable_points
-    ## Append each row obtained from the forward Euler method to A.
-    u_row = forward_euler(2, xs, x);
-    A = cat(1, A, u_row);
-  end
-
-  ## Finally, append the last row for x4.
-  A = cat(1, A, u_x4_coeffs);
-
-  ## Solve the system. This gives us the value of u(x) at xs.
-  U = A \ F';
-
-  plot(xs,U)
-end