From: Michael Orlitzky Date: Mon, 15 Oct 2012 19:00:55 +0000 (-0400) Subject: Factor out the poisson_matrix function. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=6de25e4c4273789fd2f98e309db4ef3f0902b185;p=octave.git Factor out the poisson_matrix function. --- diff --git a/poisson_matrix.m b/poisson_matrix.m new file mode 100644 index 0000000..696185f --- /dev/null +++ b/poisson_matrix.m @@ -0,0 +1,65 @@ +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. + A = cat(1, A, u_xN_coeffs); +end