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