Add the homework #1 directory to the divided_difference() load path.
+++ /dev/null
-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
## OUTPUTS:
##
## * ``dd`` - The divided difference f[xs(1), xs(2),...]
- ##
+ ##
+ addpath('../homework1/src');
order = length(xs) - 1;
if (order < 0)
dd = NA;
elseif (order == 0)
## Our base case.
- dd = f(xs(1))
+ dd = f(xs(1));
else
## Order >= 1.
- cs = divided_difference_coefficients(xs)
+ cs = divided_difference_coefficients(xs);
dd = dot(cs, f(xs));
end
end
+++ /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)) == 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
+++ /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'(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
+++ /dev/null
-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
+++ /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
+++ /dev/null
-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
projects / octave.git / commitdiff