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;
[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
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
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
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
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
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
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
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
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 = [];
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 ]
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;
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
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
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;
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;
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
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
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) );
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
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
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
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;
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);
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
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
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
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
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
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
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
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
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
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
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];
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
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;
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
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
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
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;
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
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;
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
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
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
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];
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
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;
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
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;
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);
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
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);
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);
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
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;
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
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);
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
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];
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
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
[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
projects / octave.git / commitdiff