author Michael Orlitzky Fri, 22 Mar 2013 17:39:31 +0000 (13:39 -0400) committer Michael Orlitzky Fri, 22 Mar 2013 17:39:31 +0000 (13:39 -0400)
56 files changed:

@@ -1,82 +1,82 @@
function S = advection_matrix(integerN, x0, xN)
-  ##
-  ## The numerical solution of the advection-diffusion equation,
-  ##
-  ##   -d*u''(x) + v*u'(x) + r*u = f(x)
-  ##
-  ## in one dimension, subject to the boundary conditions,
-  ##
-  ##   u(x0)  = u(xN)
-  ##
-  ##   u'(x0) = u'(xN)
-  ##
-  ## over the interval [x0,xN] gives rise to a linear system:
-  ##
-  ##   AU = h^2*F
-  ##
-  ## where h = 1/n, and A is given by,
-  ##
-  ##   A = d*K + v*h*S + r*h^2*I.
-  ##
-  ## We will call the matrix S the "advection matrix," and it will be
-  ## understood that the first row (corresponding to j=0) is to be
-  ## omitted; since we have assumed that when j=0, u(xj) = u(x0) =
-  ## u(xN) and likewise for u'. ignored (i.e, added later).
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``integerN`` - An integer representing the number of
-  ##   subintervals we should use to approximate `u`. Must be greater
-  ##   than or equal to 2, since we have at least two values for u(x0)
-  ##   and u(xN).
-  ##
-  ##   * ``f`` - The function on the right hand side of the poisson
-  ##   equation.
-  ##
-  ##   * ``x0`` - The initial point.
-  ##
-  ##   * ``xN`` - The terminal point.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``S`` - The NxN matrix of coefficients for the vector [u(x1),
-  ##   ..., u(xN)].
-  ##
-  ## EXAMPLES:
-  ##
-  ## For integerN=4, x0=0, and x1=1, we will have four subintervals:
-  ##
-  ##   [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]
-  ##
-  ## The first row of the matrix 'S' should compute the "derivative"
-  ## at x1=0.25. By the finite difference formula, this is,
-  ##
-  ##   u'(x1) = (u(x2) - u(x0))/2
-  ##
-  ##          = (u(x2) - u(x4))/2
-  ##
-  ## Therefore, the first row of 'S' should look like,
-  ##
-  ##   2*S1 = [0, 1, 0, -1]
-  ##
-  ## and of course we would have F1 =  on the right-hand side.
-  ## Likewise, the last row of S should correspond to,
-  ##
-  ##   u'(x4) = (u(x5) - u(x3))/2
-  ##
-  ##          = (u(x1) - u(x3))/2
-  ##
-  ## So the last row of S will be,
-  ##
-  ##   2*S4 = [1, 0, -1, 0]
-  ##
-  ## Each row 'i' in between will have [-1, 0, 1] beginning at column
-  ## (i-1). So finally,
-  ##
-  ##   2*S = [0, 1, 0, -1]
-  ##        [-1, 0, 1, 0]
-  ##         [0, -1, 0, 1]
-  ##         [1, 0, -1, 0]
+  %
+  % The numerical solution of the advection-diffusion equation,
+  %
+  %   -d*u''(x) + v*u'(x) + r*u = f(x)
+  %
+  % in one dimension, subject to the boundary conditions,
+  %
+  %   u(x0)  = u(xN)
+  %
+  %   u'(x0) = u'(xN)
+  %
+  % over the interval [x0,xN] gives rise to a linear system:
+  %
+  %   AU = h^2*F
+  %
+  % where h = 1/n, and A is given by,
+  %
+  %   A = d*K + v*h*S + r*h^2*I.
+  %
+  % We will call the matrix S the "advection matrix," and it will be
+  % understood that the first row (corresponding to j=0) is to be
+  % omitted; since we have assumed that when j=0, u(xj) = u(x0) =
+  % u(xN) and likewise for u'. ignored (i.e, added later).
+  %
+  % INPUTS:
+  %
+  %   * ``integerN`` - An integer representing the number of
+  %   subintervals we should use to approximate `u`. Must be greater
+  %   than or equal to 2, since we have at least two values for u(x0)
+  %   and u(xN).
+  %
+  %   * ``f`` - The function on the right hand side of the poisson
+  %   equation.
+  %
+  %   * ``x0`` - The initial point.
+  %
+  %   * ``xN`` - The terminal point.
+  %
+  % OUTPUTS:
+  %
+  %   * ``S`` - The NxN matrix of coefficients for the vector [u(x1),
+  %   ..., u(xN)].
+  %
+  % EXAMPLES:
+  %
+  % For integerN=4, x0=0, and x1=1, we will have four subintervals:
+  %
+  %   [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]
+  %
+  % The first row of the matrix 'S' should compute the "derivative"
+  % at x1=0.25. By the finite difference formula, this is,
+  %
+  %   u'(x1) = (u(x2) - u(x0))/2
+  %
+  %          = (u(x2) - u(x4))/2
+  %
+  % Therefore, the first row of 'S' should look like,
+  %
+  %   2*S1 = [0, 1, 0, -1]
+  %
+  % and of course we would have F1 =  on the right-hand side.
+  % Likewise, the last row of S should correspond to,
+  %
+  %   u'(x4) = (u(x5) - u(x3))/2
+  %
+  %          = (u(x1) - u(x3))/2
+  %
+  % So the last row of S will be,
+  %
+  %   2*S4 = [1, 0, -1, 0]
+  %
+  % Each row 'i' in between will have [-1, 0, 1] beginning at column
+  % (i-1). So finally,
+  %
+  %   2*S = [0, 1, 0, -1]
+          [-1, 0, 1, 0]
+  %         [0, -1, 0, 1]
+  %         [1, 0, -1, 0]

if (integerN < 2)
S = NA;
@@ -85,28 +85,28 @@ function S = advection_matrix(integerN, x0, xN)

[xs,h] = partition(integerN, x0, xN);

-  ## We cannot evaluate u_xx at the endpoints because our
-  ## differentiation algorithm relies on the points directly to the
-  ## left and right of `x`. Since we're starting at j=1 anyway, we cut
-  ## off two from the beginning.
+  % We cannot evaluate u_xx at the endpoints because our
+  % differentiation algorithm relies on the points directly to the
+  % left and right of `x`. Since we're starting at j=1 anyway, we cut
+  % off two from the beginning.
differentiable_points = xs(3:end-1);

-  ## These are the coefficient vectors for the u(x0) and u(xn)
-  ## constraints. There should be N zeros and a single 1.
+  % 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 - 3);
u_x0_coeffs = cat(2, the_rest_zeros, [0.5, 0, -0.5]);
u_xN_coeffs = cat(2, [0.5, 0, -0.5], the_rest_zeros);

-  ## Start with the u(x0) row.
+  % Start with the u(x0) row.
S = u_x0_coeffs;

for x = differentiable_points
-    ## Append each row obtained from the forward Euler method to S.
-    ## Chop off x0 first.
+    % Append each row obtained from the forward Euler method to S.
+    % Chop off x0 first.
u_row = central_difference(xs(2:end), x);
S = cat(1, S, u_row);
end

-  ## Finally, append the last row for xN.
+  % Finally, append the last row for xN.
S = cat(1, S, u_xN_coeffs);
end
index a89536cbd11859212464d840f6e369e418418d61..dea3e77fb0531d019e901fcde58f9317300a51d2 100644 (file)
@@ -1,29 +1,29 @@
-  ##
-  ## Sparse version of the advection_matrix function. See
-  ##
+  %
+  % Sparse version of the advection_matrix function. See
+  %

if (integerN < 2)
S = NA;
return
end

-  ## The ones directly above the diagonal.
+  % The ones directly above the diagonal.
top = [ [zeros(integerN-1, 1), speye(integerN-1)]; ...
zeros(1, integerN)];

-  ## The negative ones directly below the diagonal.
+  % The negative ones directly below the diagonal.
bottom = [ [zeros(1, integerN-1); ...
-speye(integerN-1)    ], zeros(integerN, 1)];

-  ## Combine the top and bottom.
+  % Combine the top and bottom.
S = top + bottom;

-  ## Fill in the entries in the corner.
+  % Fill in the entries in the corner.
S(1, integerN) = -1;
S(integerN, 1) = 1;

-  ## And divide the whole thing by 2.
+  % And divide the whole thing by 2.
S = (1/2)*S;
end
@@ -1,24 +1,24 @@
function ip = c_inner_product(w, a, b, v1, v2)
-  ##
-  ## The usual inner product defined on the space of continuous
-  ## functions over the interval [a,b].
-  ##
-  ## INPUT:
-  ##
-  ##   * ``w`` -- The weight function.
-  ##
-  ##   * ``a`` -- The left endpoint of the interval.
-  ##
-  ##   * ``b`` -- The right endpoint of the interval.
-  ##
-  ##   * ``v1`` -- The first vector.
-  ##
-  ##   * ``v2`` -- The second vector.
-  ##
-  ## OUTPUT:
-  ##
-  ## The inner product <v1, v2>.
-  ##
+  %
+  % The usual inner product defined on the space of continuous
+  % functions over the interval [a,b].
+  %
+  % INPUT:
+  %
+  %   * ``w`` -- The weight function.
+  %
+  %   * ``a`` -- The left endpoint of the interval.
+  %
+  %   * ``b`` -- The right endpoint of the interval.
+  %
+  %   * ``v1`` -- The first vector.
+  %
+  %   * ``v2`` -- The second vector.
+  %
+  % OUTPUT:
+  %
+  % The inner product <v1, v2>.
+  %
integrand = @(x) w(x)*v1(x)*v2(x);
end
index 865ff8b472a566b0dfd3f5237939b575206352ea..5a4fdc16419f0acb02b07357f3fc9b5b6010eb85 100644 (file)
--- a/c_norm.m
+++ b/c_norm.m
@@ -1,20 +1,20 @@
function c_norm = c_norm(w, a, b, v)
-  ##
-  ## The norm on C[a,b] induced by c_inner_product.
-  ##
-  ## INPUT:
-  ##
-  ##   * ``w`` -- The weight function.
-  ##
-  ##   * ``a`` -- The left endpoint of the interval.
-  ##
-  ##   * ``b`` -- The right endpoint of the interval.
-  ##
-  ##   * ``v`` -- The vector.
-  ##
-  ## OUTPUT:
-  ##
-  ## The norm of `v`; that is, the inner product sqrt(<v, v>).
-  ##
+  %
+  % The norm on C[a,b] induced by c_inner_product.
+  %
+  % INPUT:
+  %
+  %   * ``w`` -- The weight function.
+  %
+  %   * ``a`` -- The left endpoint of the interval.
+  %
+  %   * ``b`` -- The right endpoint of the interval.
+  %
+  %   * ``v`` -- The vector.
+  %
+  % OUTPUT:
+  %
+  % The norm of `v`; that is, the inner product sqrt(<v, v>).
+  %
c_norm = sqrt(c_inner_product(w, a, b, v, v));
end
@@ -1,34 +1,34 @@
function coefficients = central_difference(xs, x)
-  ##
-  ## The first order central difference at x1 is,
-  ##
-  ##   f'(x1) = (f(x2) - f(x0))/2
-  ##
-  ## where the index x1 is of course arbitrary but x2, x0 are adjacent
-  ## to x1. The coefficients we seek are the coefficients of f(xj) for
-  ## j = 1,...,N-2, where N is the length of ``xs``. We omit the first
-  ## and last coefficient because at x0 and xN, the previous/next
-  ## value is not available.
-  ##
-  ## This should probably take an 'order' parameter as well; see
-  ## forward_euler().
-  ##
-  ## INPUT:
-  ##
-  ##   * ``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`.
-  ##
-  ## OUTPUT:
-  ##
-  ##   * ``coefficients`` - The vector of coefficients, in order, of
-  ##   f(x0), f(x1), ..., f(xn).
-  ##
+  %
+  % The first order central difference at x1 is,
+  %
+  %   f'(x1) = (f(x2) - f(x0))/2
+  %
+  % where the index x1 is of course arbitrary but x2, x0 are adjacent
+  % to x1. The coefficients we seek are the coefficients of f(xj) for
+  % j = 1,...,N-2, where N is the length of ``xs``. We omit the first
+  % and last coefficient because at x0 and xN, the previous/next
+  % value is not available.
+  %
+  % This should probably take an 'order' parameter as well; see
+  % forward_euler().
+  %
+  % INPUT:
+  %
+  %   * ``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`.
+  %
+  % OUTPUT:
+  %
+  %   * ``coefficients`` - The vector of coefficients, in order, of
+  %   f(x0), f(x1), ..., f(xn).
+  %

if (length(xs) < 3)
-    ## We need at least one point other than the first and last.
+    % We need at least one point other than the first and last.
coefficients = NA;
return;
end
@@ -36,16 +36,16 @@ function coefficients = central_difference(xs, x)
x_idx = find(xs == x);

if (x_idx == 1 || x_idx == length(xs))
-    ## You asked for the difference at the first or last element, which
-    ## we can't do.
+    % You asked for the difference at the first or last element, which
+    % we can't do.
coefficients = NA;
return;
end

-  ## Start with a vector of zeros.
+  % Start with a vector of zeros.
coefficients = zeros(1, length(xs));

-  ## And fill in the two values that we know.
+  % And fill in the two values that we know.
coefficients(x_idx - 1) = -1/2;
coefficients(x_idx + 1) = 1/2;
end
index dbe4e2a53a0c4af68a0374f7782dd9c5b7baab85..600363c904bec4eefb8073c2bbc94f285ded8d51 100644 (file)
@@ -1,26 +1,26 @@
function K = diffusion_matrix_sparse(integerN)
-  ##
-  ## A sparse representation of the matrix K in the advection-diffusion
-  ## equation. See advection_matrix.m for details.
-  ##
+  %
+  % A sparse representation of the matrix K in the advection-diffusion
+  % equation. See advection_matrix.m for details.
+  %

if (integerN < 2)
K = NA;
return
end

-  ## The negative ones directly above the diagonal.
+  % The negative ones directly above the diagonal.
top = [ [zeros(integerN-1, 1), -speye(integerN-1)]; ...
zeros(1, integerN)];

-  ## The negative ones directly below the diagonal.
+  % The negative ones directly below the diagonal.
bottom = [ [zeros(1, integerN-1); ...
-speye(integerN-1)    ], zeros(integerN, 1)];

-  ## Combine the top and bottom.
+  % Combine the top and bottom.
K = top + bottom + 2*speye(integerN);

-  ## Fill in the entries in the corner.
+  % Fill in the entries in the corner.
K(1, integerN) = -1;
K(integerN, 1) = -1;
end
index 4705f0a3d8d1ff3cf699cca2dd21c04d30a3c2fa..261abeded4ccf9eabf468f817d1b9fa9eafe0994 100644 (file)
@@ -1,32 +1,32 @@
function dd = divided_difference(f, xs)
-  ## Compute divided difference of `f` at points `xs`. The argument `xs`
-  ## is assumed to be a vector containing at least one element. If it
-  ## contains n elements, the (n-1)st divided difference will be
-  ## calculated.
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``f`` - The function whose divided differences we want.
-  ##
-  ##   * ``xs`` - A vector containing x-coordinates. The length of `xs`
-  ##   determines the order of the divided difference.
-  ##
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``dd`` - The divided difference f[xs(1), xs(2),...]
-  ##
+  % Compute divided difference of `f` at points `xs`. The argument `xs`
+  % is assumed to be a vector containing at least one element. If it
+  % contains n elements, the (n-1)st divided difference will be
+  % calculated.
+  %
+  % INPUTS:
+  %
+  %   * ``f`` - The function whose divided differences we want.
+  %
+  %   * ``xs`` - A vector containing x-coordinates. The length of `xs`
+  %   determines the order of the divided difference.
+  %
+  %
+  % OUTPUTS:
+  %
+  %   * ``dd`` - The divided difference f[xs(1), xs(2),...]
+  %

order = length(xs) - 1;

if (order < 0)
-    ## Can't do anything here. Return nothing.
+    % Can't do anything here. Return nothing.
dd = NA;
elseif (order == 0)
-    ## Our base case.
+    % Our base case.
dd = f(xs(1));
else
-    ## Order >= 1.
+    % Order >= 1.
cs = divided_difference_coefficients(xs);
dd = dot(cs, f(xs));
end
@@ -1,15 +1,15 @@
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.
-  ##
+  % 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 = [];

index e54f8185f49a1452846411ee654888f01f3cd1c4..fd668421debc745c92fb15db086537e2ff79cbef 100644 (file)
@@ -1,26 +1,26 @@
function envelope = envelope(A)
-  ## Compute the envelope of the matrix ``A``. The envelope of a matrix
-  ## is defined as the set of indices,
-  ##
-  ##   E = { (i,j) : i < j, A(k,j) != 0 for some k <= i }
-  ##
+  % Compute the envelope of the matrix ``A``. The envelope of a matrix
+  % is defined as the set of indices,
+  %
+  %   E = { (i,j) : i < j, A(k,j) != 0 for some k <= i }
+  %
if (!issymmetric(A) && !is_upper_triangular(A))
-    ## The envelope of a matrix is only defined for U-T or symmetric
-    ## matrices.
+    % The envelope of a matrix is only defined for U-T or symmetric
+    % matrices.
envelope = {NA};
return;
end

-  ## Start with an empty result, and append to it as we find
-  ## satisfactory indices.
+  % Start with an empty result, and append to it as we find
+  % satisfactory indices.
envelope = {};

for j = [ 1 : columns(A) ]
-    ## Everything below the first non-zero element in a column will be
-    ## part of the envelope. Since we're moving from top to bottom, we
-    ## can simply set a flag indicating that we've found the first
-    ## non-zero element. Thereafter, everything we encounter should be
-    ## added to the envelope.
+    % Everything below the first non-zero element in a column will be
+    % part of the envelope. Since we're moving from top to bottom, we
+    % can simply set a flag indicating that we've found the first
+    % non-zero element. Thereafter, everything we encounter should be
+    % added to the envelope.
found_nonzero = false;

for i = [ 1 : j-1 ]
index eb3b61c27729331e7b76828625f3874a8c44f04f..a70a1ac74db3c361b11c6d5e9bf3dea16af92dd4 100644 (file)
@@ -1,21 +1,21 @@
function [fixed_point, iterations] = fixed_point_method(g, epsilon, x0)
-  ## Find a fixed_point of the function `g` with initial guess x0.
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``g`` - The function to iterate.
-  ##
-  ##   * ``epsilon`` - We stop when two successive iterations are within
-  ##     epsilon of each other, taken under the infinity norm. halt the
-  ##     search and return the current approximation.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``fixed_point`` - The fixed point that we found.
-  ##
-  ##   * ``iterations`` - The number of iterations that we performed
-  ##     during the search.
-  ##
+  % Find a fixed_point of the function `g` with initial guess x0.
+  %
+  % INPUTS:
+  %
+  %   * ``g`` - The function to iterate.
+  %
+  %   * ``epsilon`` - We stop when two successive iterations are within
+  %     epsilon of each other, taken under the infinity norm. halt the
+  %     search and return the current approximation.
+  %
+  % OUTPUTS:
+  %
+  %   * ``fixed_point`` - The fixed point that we found.
+  %
+  %   * ``iterations`` - The number of iterations that we performed
+  %     during the search.
+  %

iterations = 0;
prev = x0;
index ff9c5eb2a498bc2be65d366fa9ba497d568354b7..6984cfc8dc9c14963bc1f8f1eb83b67f36869419 100644 (file)
@@ -1,36 +1,36 @@
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).
-  ##
+  %
+  % 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.
+    % You have made a grave mistake.
coefficients = NA;
return;
end
@@ -41,8 +41,8 @@ function coefficients = forward_euler(integer_order, xs, x)
end

if (length(xs) < 2)
-    ## You can't approximate a derivative of order greater than zero
-    ## with zero or one points!
+    % You can't approximate a derivative of order greater than zero
+    % with zero or one points!
coefficients = NA
return;
end
@@ -51,15 +51,15 @@ function coefficients = forward_euler(integer_order, xs, x)
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.
+    % 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.
+  % 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;
@@ -70,7 +70,7 @@ function coefficients = forward_euler(integer_order, xs, x)

targets = xs(first_nonzero_idx : last_nonzero_idx);

-  # The multiplier comes from the Taylor expansion.
+  % The multiplier comes from the Taylor expansion.
multiplier = factorial(integer_order);
cs = divided_difference_coefficients(targets) * multiplier;

index 3e1061d48340b1b3db9c4c4c45c1ced56528b82e..2c1fb22b85d174a512749908e8c1a0e45e4fcd2a 100644 (file)
@@ -1,20 +1,20 @@
function y = forward_euler1(x0, y0, f, h)
-  ## Compute one iteration of the forward Euler method.
-  ##
-  ## INPUT:
-  ##
-  ##   * ``x0`` - The initial x-coordinate.
-  ##
-  ##   * ``y0`` - The initial value y(x0).
-  ##
-  ##   * ``f`` - The function y' = f(x,y) of two variables.
-  ##
-  ##   * ``h`` - The step size.
-  ##
-  ## OUTPUT:
-  ##
-  ## The approximate value of y(x) at x = x0+h.
-  ##
+  % Compute one iteration of the forward Euler method.
+  %
+  % INPUT:
+  %
+  %   * ``x0`` - The initial x-coordinate.
+  %
+  %   * ``y0`` - The initial value y(x0).
+  %
+  %   * ``f`` - The function y' = f(x,y) of two variables.
+  %
+  %   * ``h`` - The step size.
+  %
+  % OUTPUT:
+  %
+  % The approximate value of y(x) at x = x0+h.
+  %
y_prime = f(x0,y0);
y = y0 +  h * y_prime;
end
index b78991d8b07683600389e8d4622900a986d9a317..a5042e1a7cde417ae6b9a9fb7e6c51f4d2495a6a 100644 (file)
@@ -1,4 +1,6 @@
function isUT = is_upper_triangular(A)
-  ## Returns true if ``A`` is upper triangular, false otherwise.
+  %
+  % Returns true if ``A`` is upper triangular, false otherwise.
+  %
isUT = isequal(A, triu(A));
end
index 0481d47be73d3af8c660b9a195b3fd4a375d16e2..36f67617a453124d1b3992aa7d6f45298333dc98 100644 (file)
@@ -1,25 +1,25 @@
function P = legendre_p(n)
-  ## Return the `n`th Legendre polynomial.
-  ##
-  ## INPUT:
-  ##
-  ##   * ``n`` - The index of the polynomial that we want.
-  ##
-  ## OUTPUT:
-  ##
-  ##   * ``P`` - A polynomial function of one argument.
-  ##
+  % Return the `n`th Legendre polynomial.
+  %
+  % INPUT:
+  %
+  %   * ``n`` - The index of the polynomial that we want.
+  %
+  % OUTPUT:
+  %
+  %   * ``P`` - A polynomial function of one argument.
+  %
if (n < 0)
-    ## Can't do anything here. Return nothing.
+    % Can't do anything here. Return nothing.
P = NA;
elseif (n == 0)
-    ## One of our base cases.
+    % One of our base cases.
P = @(x) 1;
elseif (n == 1)
-    ## The second base case.
+    % The second base case.
P = @(x) x;
else
-    ## Not one of the base cases, so use the recursive formula.
+    % Not one of the base cases, so use the recursive formula.
prev = legendre_p(n-1);
prev_prev = legendre_p(n-2);
P = @(x) (1/n).*( (2*n - 1).*x.*prev(x) - (n-1).*prev_prev(x) );
index 2c9d2ab08916107fe11ca42ba863f32e4341117f..3bb30be664003fcab28ca19259d7b8f4821f045b 100644 (file)
@@ -1,24 +1,24 @@
function P_tilde = legendre_p_tilde(n, a, b)
-  ## Return the `n`th Legendre polynomial scaled to the interval [a,b].
-  ##
-  ## INPUT:
-  ##
-  ##   * ``n`` - The index of the polynomial that we want.
-  ##
-  ##   * ``a`` - The left endpoint of the interval.
-  ##
-  ##   * ``b`` - The right endpoint of the interval.
-  ##
-  ## OUTPUT:
-  ##
-  ##   * ``P_tilde`` - A polynomial function of one argument.
-  ##
+  % Return the `n`th Legendre polynomial scaled to the interval [a,b].
+  %
+  % INPUT:
+  %
+  %   * ``n`` - The index of the polynomial that we want.
+  %
+  %   * ``a`` - The left endpoint of the interval.
+  %
+  %   * ``b`` - The right endpoint of the interval.
+  %
+  % OUTPUT:
+  %
+  %   * ``P_tilde`` - A polynomial function of one argument.
+  %
if (n < 0)
-    ## Can't do anything here. Return nothing.
+    % Can't do anything here. Return nothing.
P = NA;
else
-    ## Compute the Legendre polynomial over [-1,1] and mangle it to fit
-    ## the interval [a,b].
+    % Compute the Legendre polynomial over [-1,1] and mangle it to fit
+    % the interval [a,b].
P = legendre_p(n);
P_tilde = @(x) P( (2/(b-a)).*x + 1 - (2*b)/(b-a) );
end
index f4f37674e3d380c22695e4070fbc7f9caea46cd4..c251e2432b74d7314fb9a242020d0b35d32c5396 100644 (file)
@@ -1,14 +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.
-  ##
+  % 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
index b0bc6267be4db4e29615e2d6915353a018b4e4c6..c754a12a8a39e6b83c266e61831c6b1074402ee4 100644 (file)
@@ -1,15 +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.
-  ##
-
+  % 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
index 8bdd090945b49b1e99f5c7d9790f73ac218909da..2d33c1a771b56329243603a94456a8bbbd147664 100644 (file)
@@ -1,25 +1,25 @@
function [root, iterations] = newtons_method(f, f_prime, epsilon, x0)
-  ## Find a root of the function `f` with initial guess x0.
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``f`` - The function whose root we seek. Must return a column
-  ##     vector or scalar.
-  ##
-  ##   * ``f_prime`` - The derivative or Jacobian of ``f``.
-  ##
-  ##   * ``epsilon`` - We stop when the value of `f` becomes less
-  ##     than epsilon.
-  ##
-  ##   * ``x0`` - An initial guess. Either 1d or a column vector.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``root`` - The root that we found.
-  ##
-  ##   * ``iterations`` - The number of iterations that we performed
-  ##     during the search.
-  ##
+  % Find a root of the function `f` with initial guess x0.
+  %
+  % INPUTS:
+  %
+  %   * ``f`` - The function whose root we seek. Must return a column
+  %     vector or scalar.
+  %
+  %   * ``f_prime`` - The derivative or Jacobian of ``f``.
+  %
+  %   * ``epsilon`` - We stop when the value of `f` becomes less
+  %     than epsilon.
+  %
+  %   * ``x0`` - An initial guess. Either 1d or a column vector.
+  %
+  % OUTPUTS:
+  %
+  %   * ``root`` - The root that we found.
+  %
+  %   * ``iterations`` - The number of iterations that we performed
+  %     during the search.
+  %

if (columns(x0) > 1)
root = NA;
@@ -31,8 +31,8 @@ function [root, iterations] = newtons_method(f, f_prime, epsilon, x0)
f_val = f(x0);

while (norm(f_val, Inf) > epsilon)
-    ## This uses the modified algorithm (2.11.5) in Atkinson.
-    ## Should work in 1d, too.
+    % This uses the modified algorithm (2.11.5) in Atkinson.
+    % Should work in 1d, too.
delta = f_prime(xn) \ f_val;
xn = xn - delta;
f_val = f(xn);
index 60f7cd491431b94cf183cb7caa375ebff19e929f..87062f3e2d0bdbb10804851c0a814a68e8cb6894 100644 (file)
@@ -1,40 +1,40 @@
function alpha = step_length_cgm(r, A, p)
-  ##
-  ## Compute the step length for the conjugate gradient method (CGM).
-  ## The CGM attempts to solve,
-  ##
-  ##   Ax = b
-  ##
-  ## or equivalently,
-  ##
-  ##   min[phi(x) = (1/2)<Ax,x> - <b,x>]
-  ##
-  ## where ``A`` is positive-definite. In the process, we need to
-  ## compute a number of search directions ``p`` and optimal step
-  ## lengths ``alpha``; i.e.,
-  ##
-  ##   x_{k+1} = x_{k} + alpha_{k}*p_{k}
-  ##
-  ## This function computes alpha_{k} in the formula above.
-  ##
-  ## INPUT:
-  ##
-  ##   - ``r`` -- The residual, Ax - b, at the current step.
-  ##
-  ##   - ``A`` -- The matrix ``A`` in the formulation above.
-  ##
-  ##   - ``p`` -- The current search direction.
-  ##
-  ## OUTPUT:
-  ##
-  ##   - ``alpha`` -- The minimizer of ``f(x) = x + alpha*p`` along ``p`.
-  ##
-  ## NOTES:
-  ##
-  ## All vectors are assumed to be *column* vectors.
-  ##
+  %
+  % Compute the step length for the conjugate gradient method (CGM).
+  % The CGM attempts to solve,
+  %
+  %   Ax = b
+  %
+  % or equivalently,
+  %
+  %   min[phi(x) = (1/2)<Ax,x> - <b,x>]
+  %
+  % where ``A`` is positive-definite. In the process, we need to
+  % compute a number of search directions ``p`` and optimal step
+  % lengths ``alpha``; i.e.,
+  %
+  %   x_{k+1} = x_{k} + alpha_{k}*p_{k}
+  %
+  % This function computes alpha_{k} in the formula above.
+  %
+  % INPUT:
+  %
+  %   - ``r`` -- The residual, Ax - b, at the current step.
+  %
+  %   - ``A`` -- The matrix ``A`` in the formulation above.
+  %
+  %   - ``p`` -- The current search direction.
+  %
+  % OUTPUT:
+  %
+  %   - ``alpha`` -- The minimizer of ``f(x) = x + alpha*p`` along ``p`.
+  %
+  % NOTES:
+  %
+  % All vectors are assumed to be *column* vectors.
+  %

-  ## A simple calculation should convince you that the gradient of
-  ## phi(x) above is Ax - b == r.
+  % A simple calculation should convince you that the gradient of
+  % phi(x) above is Ax - b == r.
alpha = step_length_positive_definite(r, A, p);
end
@@ -1,36 +1,36 @@
function alpha = step_length_positive_definite(g, Q, p)
-  ##
-  ## Find the minimizer alpha of,
-  ##
-  ##   phi(alpha) = f(x + alpha*p)
-  ##
-  ## where ``p`` is a descent direction,
-  ##
-  ##   f(x) = (1/2)<Qx,x> - <b,x>
-  ##
-  ## and ``Q`` is positive-definite.
-  ##
-  ## The closed-form solution to this problem is given in Nocedal and
-  ## Wright, (3.55).
-  ##
-  ## INPUT:
-  ##
-  ##   - ``g`` -- The gradient of f at x.
-  ##
-  ##   - ``Q`` -- The positive-definite matrix in the definition of
-  ##     ``f`` above.
-  ##
-  ##   - ``p`` -- The direction in which ``f`` decreases. The line
-  ##     along which we minimize f(x + alpha*p).
-  ##
-  ## OUTPUT:
-  ##
-  ##   - ``alpha`` -- The value which causes ``f`` to decrease the
-  ##     most.
-  ##
-  ## NOTES:
-  ##
-  ## All vectors are assumed to be *column* vectors.
-  ##
+  %
+  % Find the minimizer alpha of,
+  %
+  %   phi(alpha) = f(x + alpha*p)
+  %
+  % where ``p`` is a descent direction,
+  %
+  %   f(x) = (1/2)<Qx,x> - <b,x>
+  %
+  % and ``Q`` is positive-definite.
+  %
+  % The closed-form solution to this problem is given in Nocedal and
+  % Wright, (3.55).
+  %
+  % INPUT:
+  %
+  %   - ``g`` -- The gradient of f at x.
+  %
+  %   - ``Q`` -- The positive-definite matrix in the definition of
+  %     ``f`` above.
+  %
+  %   - ``p`` -- The direction in which ``f`` decreases. The line
+  %     along which we minimize f(x + alpha*p).
+  %
+  % OUTPUT:
+  %
+  %   - ``alpha`` -- The value which causes ``f`` to decrease the
+  %     most.
+  %
+  % NOTES:
+  %
+  % All vectors are assumed to be *column* vectors.
+  %
alpha = -(g' * p)/(p' * Q * p);
end
@@ -1,17 +1,17 @@
function f = extended_powell1(x)
-  ##
-  ## The extended Powell function. See Dennis & Schnabel, Appendix B,
-  ## problem #2.
-  ##
-  ## This function has a minimum at x=(0,0,...,0) with f(x) == 0. The
-  ## suggested starting point is x0=(3, -1, 0, 1,..., 3, -1, 0, 1).
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The extended Powell function. See Dennis & Schnabel, Appendix B,
+  % problem #2.
+  %
+  % This function has a minimum at x=(0,0,...,0) with f(x) == 0. The
+  % suggested starting point is x0=(3, -1, 0, 1,..., 3, -1, 0, 1).
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
f = NA;
return;
end
index 7cc7aabf2980c2b3e65f791d631644f6db240b10..f60d5a8b246dc3abb07551d2cb5273083e2b2329 100644 (file)
@@ -1,15 +1,15 @@
-  ##
-  ## The gradient of the extended Powell function. See
-  ## extended_powell1.m for more information.
-  ##
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The gradient of the extended Powell function. See
+  % extended_powell1.m for more information.
+  %
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
g = NA;
return;
end
@@ -1,15 +1,15 @@
function H = extended_powell_hessian1(x)
-  ##
-  ## The Hessian of the extended Powell function. See
-  ## extended_powell1.m for more information.
-  ##
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The Hessian of the extended Powell function. See
+  % extended_powell1.m for more information.
+  %
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
H = NA;
return;
end
index 1f545207a0ff0eee5c71a0b820979fff41601382..8d0b4f841761ddbc63c0a730ed3735ba8c93f69b 100644 (file)
@@ -1,17 +1,17 @@
function f = extended_rosenbrock1(x)
-  ##
-  ## The extended Rosenbrock function. See Dennis & Schnabel, Appendix
-  ## B, problem #1.
-  ##
-  ## This function has a minimum at x=(1,1,...,1) with f(x) == 0. The
-  ## suggested starting point is x0=(-1.2, 1,-1.2, 1,...,-1.2, 1).
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The extended Rosenbrock function. See Dennis & Schnabel, Appendix
+  % B, problem #1.
+  %
+  % This function has a minimum at x=(1,1,...,1) with f(x) == 0. The
+  % suggested starting point is x0=(-1.2, 1,-1.2, 1,...,-1.2, 1).
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
f = NA;
return;
end
index 2fd3c8cc63b2c8033fb9ececcb2bb84f1b36bf5c..cbda9a84ce328d81b3845cfa6df3723732915f44 100644 (file)
@@ -1,15 +1,15 @@
-  ##
-  ## The gradient of the extended Rosenbrock function. See
-  ## extended_rosenbrock1.m for more information.
-  ##
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The gradient of the extended Rosenbrock function. See
+  % extended_rosenbrock1.m for more information.
+  %
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
g = NA;
return;
end
index 9c3415d6d254d1d6ab821791754ce5ab61421656..541613d560ffa56d3e613e932b50bb12f75d7b4c 100644 (file)
@@ -1,15 +1,15 @@
function H = extended_rosenbrock_hessian1(x)
-  ##
-  ## The Hessian of the extended Rosenbrock function. See
-  ## extended_rosenbrock1.m for more information.
-  ##
-  ## Since the number of arguments is variable, we take a vector
-  ## instead of its individual components.
-  ##
+  %
+  % The Hessian of the extended Rosenbrock function. See
+  % extended_rosenbrock1.m for more information.
+  %
+  % Since the number of arguments is variable, we take a vector
+  % instead of its individual components.
+  %
n = length(x);

if (odd(n))
-    ## 'm' below must be an integer.
+    % 'm' below must be an integer.
H = NA;
return;
end
index c973a866ceb71428965abed7f9ba79764c6a6942..c151ed616026c0de428f21e36ff757665cd13263 100644 (file)
@@ -1,17 +1,17 @@
function f = himmelblau(x1, x2)
-  ##
-  ## The eponymous function defined by Himmelblau in his Applied
-  ## Nonlinear Programming, 1972.
-  ##
-  ## It has four identical local minima,
-  ##
-  ##   * f(3,2)                                    == 0
-  ##
-  ##   * f(-2.805118086952745,  3.131312518250573) == 0
-  ##
-  ##   * f(-3.779310253377747, -3.283185991286170) == 0
-  ##
-  ##   * f(3.584428340330492, -1.848126526964404)  == 0
-  ##
+  %
+  % The eponymous function defined by Himmelblau in his Applied
+  % Nonlinear Programming, 1972.
+  %
+  % It has four identical local minima,
+  %
+  %   * f(3,2)                                    == 0
+  %
+  %   * f(-2.805118086952745,  3.131312518250573) == 0
+  %
+  %   * f(-3.779310253377747, -3.283185991286170) == 0
+  %
+  %   * f(3.584428340330492, -1.848126526964404)  == 0
+  %
f = (x1^2 + x2 - 11)^2 + (x1 + x2^2 - 7)^2;
end
@@ -1,9 +1,9 @@
function f = himmelblau1(x)
-  ##
-  ## A version of the Himmelblau function which takes a column
-  ## 2-vector instead of two distinct arguments. See himmelblau.m for
-  ## more information.
-  ##
+  %
+  % A version of the Himmelblau function which takes a column
+  % 2-vector instead of two distinct arguments. See himmelblau.m for
+  % more information.
+  %
if (length(x) == 2)
f = himmelblau(x(1), x(2));
else
@@ -1,8 +1,8 @@
-  ##
-  ## The gradient of the Himmelblau function. See himmelblau.m for
-  ## more information.
-  ##
+  %
+  % The gradient of the Himmelblau function. See himmelblau.m for
+  % more information.
+  %
f_x1 = 4*(x1^2 + x2 - 11)*x1 + 2*x2^2 + 2*x1 - 14;
f_x2 = 4*(x2^2 + x1 - 7)*x2 + 2*x1^2 + 2*x2 - 22;
g = [f_x1; f_x2];
index 91db28d23e06484dd82756b0207ce05a99267457..5637fa8dcd2a1612f747a761204b175b66fbe4f6 100644 (file)
@@ -1,9 +1,9 @@
-  ##
-  ## A version of the himmelblau_gradient() function which takes a
-  ## column 2-vector instead of two distinct arguments. See
-  ## himmelblau_gradient.m for more information.
-  ##
+  %
+  % A version of the himmelblau_gradient() function which takes a
+  % column 2-vector instead of two distinct arguments. See
+  % himmelblau_gradient.m for more information.
+  %
if (length(x) == 2)
else
index 24b277ac914958e283283396d1b7b841d9bb1b41..b846a6ac17048323e4151df7115804960bef71af 100644 (file)
@@ -1,8 +1,8 @@
function H = himmelblau_hessian(x1, x2)
-  ##
-  ## The Hessian of the Himmelblau function. See himmelblau.m for more
-  ## information.
-  ##
+  %
+  % The Hessian of the Himmelblau function. See himmelblau.m for more
+  % information.
+  %
H = zeros(2,2);
H(1,1) = 12*x1^2 + 4*x2 - 42;
H(1,2) = 4*x1 + 4*x2;
@@ -1,9 +1,9 @@
function H = himmelblau_hessian1(x)
-  ##
-  ## A version of the himmelblau_hessian() function which takes a column
-  ## 2-vector instead of two distinct arguments. See himmelblau_hessian.m
-  ## for more information.
-  ##
+  %
+  % A version of the himmelblau_hessian() function which takes a column
+  % 2-vector instead of two distinct arguments. See himmelblau_hessian.m
+  % for more information.
+  %
if (length(x) == 2)
H = himmelblau_hessian(x(1), x(2));
else
@@ -1,11 +1,11 @@
function f = powell(x1,x2,x3,x4)
-  ## The Powell function. See Dennis & Schnabel, Appendix B, problem
-  ## #1. (The "regular" Powell function is simply the Extended Powell
-  ## with m=1).
-  ##
-  ## This function has a minimum at x=(0,0,0,0) with f(x) == 0. The
-  ## suggested starting point is x0=(3,-1,0,1).
-  ##
+  % The Powell function. See Dennis & Schnabel, Appendix B, problem
+  % #1. (The "regular" Powell function is simply the Extended Powell
+  % with m=1).
+  %
+  % This function has a minimum at x=(0,0,0,0) with f(x) == 0. The
+  % suggested starting point is x0=(3,-1,0,1).
+  %
f = (x1 + 10*x2)^2 + 5*(x3 - x4)^2;
f = f + (x2 - 2*x3)^4 + 10*(x1 - x4)^4;
end
index b5e1771471658e5c89959107dacc6a3f8808d68d..52937b579fcc9b0a6303ec16bee5fc89cb4a7278 100644 (file)
@@ -1,8 +1,8 @@
function f = powell1(x)
-  ##
-  ## A version of the Powell function which takes a column 4-vector as
-  ## an argument instead of four distinct arguments.
-  ##
+  %
+  % A version of the Powell function which takes a column 4-vector as
+  % an argument instead of four distinct arguments.
+  %
if (length(x) == 4)
f = powell(x(1), x(2), x(3), x(4));
else
index 7cea917937e5e8f2cb176c63e5161c5ed33e85f2..835f3faec60cc7d3f139035e0decfc8be43c6bdb 100644 (file)
@@ -1,8 +1,8 @@
-  ##
-  ## The gradient of the Powell function. See powell.m for more
-  ## information.
-  ##
+  %
+  % The gradient of the Powell function. See powell.m for more
+  % information.
+  %
f_x1 = 40*(x1 - x4)^3 + 2*x1 + 20*x2;
f_x2 = 4*(x2 - 2*x3)^3 + 20*x1 + 200*x2;
f_x3 = -8*(x2 - 2*x3)^3 + 10*x3 - 10*x4;
index e21b315cfac37bccca2133a7eee8655b23998a84..63868cd7dd81f177a774dd5fddf25820cdc39eb5 100644 (file)
@@ -1,9 +1,9 @@
-  ##
-  ## A version of the powell_gradient() function which takes a column
-  ## 4-vector instead of four distinct arguments. See
-  ## powell_gradient.m for more information.
-  ##
+  %
+  % A version of the powell_gradient() function which takes a column
+  % 4-vector instead of four distinct arguments. See
+  % powell_gradient.m for more information.
+  %
if (length(x) == 4)
g = powell_gradient(x(1), x(2), x(3), x(4));
else
@@ -1,8 +1,8 @@
function H = powell_hessian(x1, x2, x3, x4)
-  ##
-  ## The Hessian of the Powell function. See powell.m for more
-  ## information.
-  ##
+  %
+  % The Hessian of the Powell function. See powell.m for more
+  % information.
+  %
H = zeros(4,4);

H(1,1) = 120*(x1 - x4)^2 + 2;
index b11d7f8f338436bfa80198ed20c025811bb1005a..f734b07681d22c5f3d8f3e8513b8cc3eeb70d9a2 100644 (file)
@@ -1,9 +1,9 @@
function H = powell_hessian1(x)
-  ##
-  ## A version of the powell_hessian() function which takes a column
-  ## 4-vector instead of four distinct arguments. See powell_hessian.m
-  ## for more information.
-  ##
+  %
+  % A version of the powell_hessian() function which takes a column
+  % 4-vector instead of four distinct arguments. See powell_hessian.m
+  % for more information.
+  %
if (length(x) == 4)
H = powell_hessian(x(1), x(2), x(3), x(4));
else
@@ -1,11 +1,11 @@
function f = rosenbrock(x1, x2)
-  ##
-  ## The Rosenbrock function. See Dennis & Schnabel, Appendix B, problem #1.
-  ## (The "regular" Rosenbrock function is simply the Extended Rosenbrock
-  ## with m=1).
-  ##
-  ## This function has a minimum at x=(1,1) with f(x) == 0. The
-  ## suggested starting point is x0=(-1.2, 1).
-  ##
+  %
+  % The Rosenbrock function. See Dennis & Schnabel, Appendix B, problem #1.
+  % (The "regular" Rosenbrock function is simply the Extended Rosenbrock
+  % with m=1).
+  %
+  % This function has a minimum at x=(1,1) with f(x) == 0. The
+  % suggested starting point is x0=(-1.2, 1).
+  %
f = 100*(x2 - x1^2)^2 + (1 - x1)^2;
end
index 5e4df05ed8c7d6e776bd858c4ee91538b707bbf5..9aa69c63718ae0d0b72d24a65263673149f49596 100644 (file)
@@ -1,9 +1,9 @@
function f = rosenbrock1(x)
-  ##
-  ## A version of the Rosenbrock function which takes a column
-  ## 2-vector instead of two distinct arguments. See rosenbrock.m for
-  ## more information.
-  ##
+  %
+  % A version of the Rosenbrock function which takes a column
+  % 2-vector instead of two distinct arguments. See rosenbrock.m for
+  % more information.
+  %
if (length(x) == 2)
f = rosenbrock(x(1), x(2));
else
index 50c4546ffffb5560891c1a3dff4abfe08423832d..2b741a407fc0d7c857fb6b1fd59c165bbe0c4d79 100644 (file)
@@ -1,8 +1,8 @@
-  ##
-  ## The gradient of the Rosenbrock function. See rosenbrock.m for more
-  ## information.
-  ##
+  %
+  % The gradient of the Rosenbrock function. See rosenbrock.m for more
+  % information.
+  %
f_x1 = -400*x1*(x2 - x1^2) + 2*x1 - 2;
f_x2 = 200*(x2 - x1^2);
g = [f_x1; f_x2];
@@ -1,9 +1,9 @@
-  ##
-  ## A version of the rosenbrock_gradient() function which takes a
-  ## column 2-vector instead of two distinct arguments. See
-  ## rosenbrock_gradient.m for more information.
-  ##
+  %
+  % A version of the rosenbrock_gradient() function which takes a
+  % column 2-vector instead of two distinct arguments. See
+  % rosenbrock_gradient.m for more information.
+  %
if (length(x) == 2)
else
index 066b5aeaa5dfdf4cd44d6d0202a04d035c0a3ae8..ccbb84d4ba9e6880a4613ace099210ddbf1a509f 100644 (file)
@@ -1,8 +1,8 @@
function H = rosenbrock_hessian(x1, x2)
-  ##
-  ## The Hessian of the Rosenbrock function. See rosenbrock.m for more
-  ## information.
-  ##
+  %
+  % The Hessian of the Rosenbrock function. See rosenbrock.m for more
+  % information.
+  %
H = zeros(2,2);
H(1,1) = 1200*x1^2 - 400*x2 + 2;
H(1,2) = -400*x1;
index 63a105f9757cd204dccc900cce5a4647240e9a36..804098ecfb1535e524b805e1be9dfe3a41768492 100644 (file)
@@ -1,9 +1,9 @@
function H = rosenbrock_hessian1(x)
-  ##
-  ## A version of the rosenbrock_hessian() function which takes a column
-  ## 2-vector instead of two distinct arguments. See rosenbrock_hessian.m
-  ## for more information.
-  ##
+  %
+  % A version of the rosenbrock_hessian() function which takes a column
+  % 2-vector instead of two distinct arguments. See rosenbrock_hessian.m
+  % for more information.
+  %
if (length(x) == 2)
H = rosenbrock_hessian(x(1), x(2));
else
index f7fd625e7bd9297b9fd2e7310a9f0a8b80805889..e6725ff4c9a808be09dd8186f10597f6e6b3b7a3 100644 (file)
@@ -1,11 +1,11 @@
function f = trigonometric1(x)
-  ##
-  ## The trigonometric function. See More, Garbow, and Hillstrom,
-  ## function #26.
-  ##
-  ## This function has a minimum with f(x) == 0. The suggested
-  ## starting point is x0=(1/n, 1/n,...).
-  ##
+  %
+  % The trigonometric function. See More, Garbow, and Hillstrom,
+  % function #26.
+  %
+  % This function has a minimum with f(x) == 0. The suggested
+  % starting point is x0=(1/n, 1/n,...).
+  %
n = length(x);
f = 0;

index 134ce8699bebb2f850f427df344a56b6f086d14f..845dc3290d41fd1fab2efcf87ab163d3e88233c4 100644 (file)
@@ -1,8 +1,8 @@
-  ##
-  ## The gradient of the trigonometric function. See trigonometric1.m
-  ## for more information.
-  ##
+  %
+  % The gradient of the trigonometric function. See trigonometric1.m
+  % for more information.
+  %
n  = length(x);
g  = zeros(n,1);

@@ -12,13 +12,13 @@ function g = trigonometric_gradient1(x)
f_k = n - cos_sum + k*(1 - cos(x(k))) - sin(x(k));

for j = [ 1 : n ]
-      ## Add to the jth component of g the partial of f^2 with
-      ## respect to x(j). The first term that we add here exists
-      ## regardless of j.
+      % Add to the jth component of g the partial of f^2 with
+      % respect to x(j). The first term that we add here exists
+      % regardless of j.
g(j) = g(j) + 2*f_k*sin(x(j));

if (j == k)
-       ## But this term vanishes when j != k.
+       % But this term vanishes when j != k.
g(j) = g(j) + 2*f_k*k*sin(x(k)) - 2*f_k*cos(x(k));
end
end
index f3e40a412b53d9f58305e1b2137e53fa632f686a..8435a3a98cde72f3297be61b7c697275a512e435 100644 (file)
@@ -1,9 +1,9 @@
function H = trigonometric_hessian1(x)
-  ##
-  ## The Hessian of the Trigonometric function. See trigonometric1.m
-  ## for more information. Not my implementation. Too ugly to
-  ## recreate.
-  ##
+  %
+  % The Hessian of the Trigonometric function. See trigonometric1.m
+  % for more information. Not my implementation. Too ugly to
+  % recreate.
+  %
n  = length(x);
H  = zeros(n,n);

index 3e11a40f0f915c88f52da01637c638e6dbdbfd7e..eb8b2e7983b7864565249ee743ab0a72526275c9 100644 (file)
@@ -1,9 +1,9 @@
function f = wood(x1, x2, x3, x4)
-  ##
-  ## The Wood function. See Dennis & Schnabel, Appendix B, problem #5.
-  ## This function has a global minimum at x=(1,1,1,1) with f(x) == 0.
-  ## The suggested starting point is x0=(-3, -1, -3, -1).
-  ##
+  %
+  % The Wood function. See Dennis & Schnabel, Appendix B, problem #5.
+  % This function has a global minimum at x=(1,1,1,1) with f(x) == 0.
+  % The suggested starting point is x0=(-3, -1, -3, -1).
+  %
f = 100*(x1^2 - x2)^2 + (x1 - 1)^2 + (x3 - 1)^2;
f = f + 90*(x3^2 - x4)^2;
f = f + 10.1*((x2 - 1)^2 + (x4 - 1)^2);
index c34de26a9dfe738f59f595b5a5f3731cef472f32..c0bfa201df4edded6be7e3e21c0598d523dd8056 100644 (file)
@@ -1,9 +1,9 @@
function f = wood1(x)
-  ##
-  ## A version of the Wood function which takes a column 4-vector
-  ## instead of four distinct arguments. See wood.m for more
-  ## information.
-  ##
+  %
+  % A version of the Wood function which takes a column 4-vector
+  % instead of four distinct arguments. See wood.m for more
+  % information.
+  %
if (length(x) == 4)
f = wood(x(1), x(2), x(3), x(4));
else
@@ -1,8 +1,8 @@
function g = wood_gradient(x1, x2, x3, x4)
-  ##
-  ## The gradient of the Wood function. See wood.m for more
-  ## information.
-  ##
+  %
+  % The gradient of the Wood function. See wood.m for more
+  % information.
+  %
f_x1 = 400*(x1^2 - x2)*x1 + 2*x1 - 2;
f_x2 = -200*x1^2 + 220.2*x2 + 19.8*x4 - 40;
f_x3 = 360*(x3^2 - x4)*x3 + 2*x3 - 2;
index eabccc94175b38fc9cb6c50757f78e4b66749122..bfe86f32ffda7576100a85b0d24dcc3b0d20113d 100644 (file)
@@ -1,9 +1,9 @@
-  ##
-  ## A version of the wood_gradient() function which takes a column
-  ## for more information.
-  ##
+  %
+  % A version of the wood_gradient() function which takes a column
+  % for more information.
+  %
if (length(x) == 4)
g = wood_gradient(x(1), x(2), x(3), x(4));
else
@@ -1,8 +1,8 @@
function H = wood_hessian(x1, x2, x3, x4)
-  ##
-  ## The Hessian of the Wood function. See wood.m for more
-  ## information.
-  ##
+  %
+  % The Hessian of the Wood function. See wood.m for more
+  % information.
+  %
H = zeros(4,4);
H(1,1) = 1200*x(1)^2 - 400*x(2) + 2;
H(1,2) = -400*x(1);
@@ -1,9 +1,9 @@
function H = wood_hessian1(x)
-  ##
-  ## A version of the wood_hessian() function which takes a column
-  ## 4-vector instead of four distinct arguments. See wood_hessian.m
-  ## for more information.
-  ##
+  %
+  % A version of the wood_hessian() function which takes a column
+  % 4-vector instead of four distinct arguments. See wood_hessian.m
+  % for more information.
+  %
if (length(x) == 4)
H = wood_hessian(x(1), x(2), x(3), x(4));
else
@@ -1,27 +1,27 @@
function [p,delta] = partition(integerN, a, b)
-  ## Partition the interval [a,b] into integerN subintervals. We do not
-  ## require 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.
-  ##
-  ##
+  % Partition the interval [a,b] into integerN subintervals. We do not
+  % require 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`.
+  % 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];
@@ -1,29 +1,29 @@
function PMs = permutation_matrices(integerN)
-  ## Generate all permutation matrices of size ``integerN``.
-  ##
-  ## INPUT:
-  ##
-  ##   - ``integerN`` -- The dimension of the resulting matrices.
-  ##
-  ## OUTPUT:
-  ##
-  ##   - ``PMs`` -- A cell array of permutation matrices.
-  ##
+  % Generate all permutation matrices of size ``integerN``.
+  %
+  % INPUT:
+  %
+  %   - ``integerN`` -- The dimension of the resulting matrices.
+  %
+  % OUTPUT:
+  %
+  %   - ``PMs`` -- A cell array of permutation matrices.
+  %

if (integerN < 1)
PMs = NA;
return;
end

-  ## Append to this as we generate them.
+  % Append to this as we generate them.
PMs = {};

-  ## Generate all permutations of [1,2,...,integerN].
+  % Generate all permutations of [1,2,...,integerN].
permutations = perms([1:integerN]);

for idx = [ 1 : factorial(integerN) ]
sigma = permutations(idx,:);
-    ## Create a permutation matrix from the permutation, sigma.
+    % Create a permutation matrix from the permutation, sigma.
P = eye(integerN) (sigma,:);
PMs{end+1} = P;
end
index 4171374f0cb03799d3312bc29758218101717173..3ed1c0c673948c5ae44217cb78d0d45559335956 100644 (file)
@@ -1,37 +1,37 @@
function A = poisson_matrix(integerN, x0, xN)
-  ##
-  ## In the numerical solution of the poisson equation,
-  ##
-  ##   -u''(x) = f(x)
-  ##
-  ## in one dimension, subject to the boundary conditions,
-  ##
-  ##   u(x0) = 0
-  ##   u(xN) = 1
-  ##
-  ## over the interval [x0,xN], we need to compute a matrix A which
-  ## is then multiplied by the vector of u(x0), ..., u(xN). The entries
-  ## of A are determined from the second order finite difference formula,
-  ## i.e. the second order forward Euler method.
-  ##
-  ## INPUTS:
-  ##
-  ##   * ``integerN`` - An integer representing the number of
-  ##     subintervals we should use to approximate `u`. Must be
-  ##     greater than or equal to 2, since we have at least two
-  ##     values for u(x0) and u(xN).
-  ##
-  ##   * ``f`` - The function on the right hand side of the poisson
-  ##     equation.
-  ##
-  ##   * ``x0`` - The initial point.
-  ##
-  ##   * ``xN`` - The terminal point.
-  ##
-  ## OUTPUTS:
-  ##
-  ##   * ``A`` - The (N+1)x(N+1) matrix of coefficients for u(x0),
-  ##             ..., u(xN).
+  %
+  % In the numerical solution of the poisson equation,
+  %
+  %   -u''(x) = f(x)
+  %
+  % in one dimension, subject to the boundary conditions,
+  %
+  %   u(x0) = 0
+  %   u(xN) = 1
+  %
+  % over the interval [x0,xN], we need to compute a matrix A which
+  % is then multiplied by the vector of u(x0), ..., u(xN). The entries
+  % of A are determined from the second order finite difference formula,
+  % i.e. the second order forward Euler method.
+  %
+  % INPUTS:
+  %
+  %   * ``integerN`` - An integer representing the number of
+  %     subintervals we should use to approximate `u`. Must be
+  %     greater than or equal to 2, since we have at least two
+  %     values for u(x0) and u(xN).
+  %
+  %   * ``f`` - The function on the right hand side of the poisson
+  %     equation.
+  %
+  %   * ``x0`` - The initial point.
+  %
+  %   * ``xN`` - The terminal point.
+  %
+  % OUTPUTS:
+  %
+  %   * ``A`` - The (N+1)x(N+1) matrix of coefficients for u(x0),
+  %             ..., u(xN).

if (integerN < 2)
A = NA
@@ -40,27 +40,27 @@ function A = poisson_matrix(integerN, x0, xN)

[xs,h] = partition(integerN, x0, xN);

-  ## We cannot evaluate u_xx at the endpoints because our
-  ## differentiation algorithm relies on the points directly to the left
-  ## and right of `x`.
+  % 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);

-  ## These are the coefficient vectors for the u(x0) and u(xn)
-  ## constraints. There should be N zeros and a single 1.
+  % 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_xN_coeffs = cat(2, the_rest_zeros, 1);

-  ## Start with the u(x0) row.
+  % 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.
+    % 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 xN and negate the whole thing (see
-  ## the definition of the poisson problem).
+  % Finally, append the last row for xN and negate the whole thing (see
+  % the definition of the poisson problem).
A = - cat(1, A, u_xN_coeffs);
end