Replace ##-style comments with %-style comments in all non-test code.
authorMichael Orlitzky <michael@orlitzky.com>
Fri, 22 Mar 2013 17:39:31 +0000 (13:39 -0400)
committerMichael Orlitzky <michael@orlitzky.com>
Fri, 22 Mar 2013 17:39:31 +0000 (13:39 -0400)
56 files changed:
advection_matrix.m
advection_matrix_sparse.m
c_inner_product.m
c_norm.m
central_difference.m
diffusion_matrix_sparse.m
divided_difference.m
divided_difference_coefficients.m
envelope.m
fixed_point_method.m
forward_euler.m
forward_euler1.m
is_upper_triangular.m
legendre_p.m
legendre_p_tilde.m
misc/even.m
misc/odd.m
newtons_method.m
optimization/step_length_cgm.m
optimization/step_length_positive_definite.m
optimization/test_functions/extended_powell1.m
optimization/test_functions/extended_powell_gradient1.m
optimization/test_functions/extended_powell_hessian1.m
optimization/test_functions/extended_rosenbrock1.m
optimization/test_functions/extended_rosenbrock_gradient1.m
optimization/test_functions/extended_rosenbrock_hessian1.m
optimization/test_functions/himmelblau.m
optimization/test_functions/himmelblau1.m
optimization/test_functions/himmelblau_gradient.m
optimization/test_functions/himmelblau_gradient1.m
optimization/test_functions/himmelblau_hessian.m
optimization/test_functions/himmelblau_hessian1.m
optimization/test_functions/powell.m
optimization/test_functions/powell1.m
optimization/test_functions/powell_gradient.m
optimization/test_functions/powell_gradient1.m
optimization/test_functions/powell_hessian.m
optimization/test_functions/powell_hessian1.m
optimization/test_functions/rosenbrock.m
optimization/test_functions/rosenbrock1.m
optimization/test_functions/rosenbrock_gradient.m
optimization/test_functions/rosenbrock_gradient1.m
optimization/test_functions/rosenbrock_hessian.m
optimization/test_functions/rosenbrock_hessian1.m
optimization/test_functions/trigonometric1.m
optimization/test_functions/trigonometric_gradient1.m
optimization/test_functions/trigonometric_hessian1.m
optimization/test_functions/wood.m
optimization/test_functions/wood1.m
optimization/test_functions/wood_gradient.m
optimization/test_functions/wood_gradient1.m
optimization/test_functions/wood_hessian.m
optimization/test_functions/wood_hessian1.m
partition.m
permutation_matrices.m
poisson_matrix.m

index 5e76ee9ea19891250f262fcdd1675c334acfe971..591d012a3388ae30a40a9692bad97c3722bc6c31 100644 (file)
@@ -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 = [0] 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 = [0] 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 @@
 function S = advection_matrix_sparse(integerN)
-  ##
-  ## Sparse version of the advection_matrix function. See
-  ## advection_matrix.m for details.
-  ##
+  %
+  % Sparse version of the advection_matrix function. See
+  % advection_matrix.m for details.
+  %
 
   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
index 1d7bfd1e67122823724ebafb8ee919fe58c9980e..b53deef48adeba65439e81961ddf865dce4d7e60 100644 (file)
@@ -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);
   ip = quad(integrand, a, b);
 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
index ad58fd44546fcd3031ef2d5e7514ecd5d0c6c08b..ec9c7ed45fd41c4ee3b5f056e09f9bd7118baf98 100644 (file)
@@ -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
index a40f5597d6b47712ad7e0b397585d71641b80ec8..92e991575b54460cbf3475ebb92e68d7faa8253f 100644 (file)
@@ -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
index f6248f67b1dd40bf22ebdfed0947ddd6edbc47fa..ad63956cfa603c097084d15291962c552d6c9d7a 100644 (file)
@@ -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
index 2b1da25003990e21bb91657f1dae8ce6ad085590..2a53f8945cc72b44ca6280e04a4e58b0c662e44d 100644 (file)
@@ -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 @@
 function g = extended_powell_gradient1(x)
-  ##
-  ## 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
index 972e690da5d83c0d667fff69a2169b3ade8b12bf..4d5f50f96a4544c356368a9b73e5559757a3567d 100644 (file)
@@ -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 @@
 function g = extended_rosenbrock_gradient1(x)
-  ##
-  ## 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
index c975af25713f405b3c83b076b1a4759825bc9c2e..2598db7c9dbb4b3cefa02e29ec99fcad40e35945 100644 (file)
@@ -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
index b50e7b6265879d06677b64259ad2912bb88a8eef..44b59f6e7cc3ad1bcfbf7aa7453125f4a1767ab3 100644 (file)
@@ -1,8 +1,8 @@
 function g = himmelblau_gradient(x1,x2)
-  ##
-  ## 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 @@
 function g = himmelblau_gradient1(x)
-  ##
-  ## 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)
     g = himmelblau_gradient(x(1), 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;
index 84efb3100c4582e58b96c9755f231de2fcd835b3..bf4b56162a46d3673dead735ada273656d379cc9 100644 (file)
@@ -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
index 5031316539e587600f61bcc220adec9771086c89..5ea74e388f10094ba981b188f6c5541ca8c701a0 100644 (file)
@@ -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 @@
 function g = powell_gradient(x1,x2,x3,x4)
-  ##
-  ## 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 @@
 function g = powell_gradient1(x)
-  ##
-  ## 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
index 48cc390df91b3ac1436bb3ead7c2190c9a3a500a..248149ceef34c14e75913c601b2ba051d1a4ab49 100644 (file)
@@ -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
index 8073a0ae10adfe92ff703abc25abbed24567de97..22b2730726846aaa8f3b5b99f8f433c0a1596a4b 100644 (file)
@@ -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 @@
 function g = rosenbrock_gradient(x1, x2)
-  ##
-  ## 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];
index e8a83b6c34a878badc6cea642f26bc09219eb480..29fa3938c0c4411470dc7a484e88fe1a73384b99 100644 (file)
@@ -1,9 +1,9 @@
 function g = rosenbrock_gradient1(x)
-  ##
-  ## 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)
     g = rosenbrock_gradient(x(1), 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 @@
 function g = trigonometric_gradient1(x)
-  ##
-  ## 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
index eb9f3e643051d2d751faff487fda218934e5a1bf..adf37c2a60f04d4ce48e78187866fae11b757d99 100644 (file)
@@ -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 @@
 function g = wood_gradient1(x)
-  ##
-  ## A version of the wood_gradient() function which takes a column
-  ## 4-vector instead of four distinct arguments. See wood_gradient.m
-  ## for more information.
-  ##
+  %
+  % A version of the wood_gradient() function which takes a column
+  % 4-vector instead of four distinct arguments. See wood_gradient.m
+  % for more information.
+  %
   if (length(x) == 4)
     g = wood_gradient(x(1), x(2), x(3), x(4));
   else
index f8e38c4456404f4fbcfe7acf14fd7f88f03ac9fb..e71f180bb0600fad5360a5e8e1b8da608fb7fe21 100644 (file)
@@ -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);
index 106cb403aef60337b946767d7b82df8ab5adf19e..4af006a785994a697242ab3083576de7158126ed 100644 (file)
@@ -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
index a08423671716d741302f67944344e3fbf76236db..99e38921d8d5c103b026ad9f9b0e9d98608bf3ac 100644 (file)
@@ -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];
index 9032d2e8eef5528355985d040133201a5421dea2..1ad1e57dc93c07995ee09d870b9100cf7384e749 100644 (file)
@@ -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