X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=poisson_matrix.m;h=3ed1c0c673948c5ae44217cb78d0d45559335956;hp=4171374f0cb03799d3312bc29758218101717173;hb=b12c6c2a4bf4cef29b2e08b743c92889505c7ed9;hpb=e1b71b4ca7cfa08ac76744a17a3778d4ccfaa7e2 diff --git a/poisson_matrix.m b/poisson_matrix.m index 4171374..3ed1c0c 100644 --- a/poisson_matrix.m +++ b/poisson_matrix.m @@ -1,37 +1,37 @@ 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 @@ -40,27 +40,27 @@ function A = poisson_matrix(integerN, x0, xN) [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