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