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). if (integerN < 2) A = NA return end [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`. 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. 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. A = u_x0_coeffs; for x = differentiable_points % Append each row obtained from the forward Euler method to A. u_row = forward_euler(2, xs, x); A = cat(1, A, u_row); end % Finally, append the last row for xN and negate the whole thing (see % the definition of the poisson problem). A = - cat(1, A, u_xN_coeffs); end