From f32a60e5aceefce5b4b7497b9d295c8175297110 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 14 Sep 2012 23:52:13 -0400 Subject: [PATCH] First shot at a solve_poisson function. --- solve_poisson.m | 46 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/solve_poisson.m b/solve_poisson.m index be1b1ec..946d1d8 100644 --- a/solve_poisson.m +++ b/solve_poisson.m @@ -12,6 +12,50 @@ function u = solve_poisson(integerN, f) ## over the interval [0,1]. It is assumed that the function `f` is at ## least continuous on [0,1]. ## + ## INPUTS: + ## + ## * ``integerN`` - An integer representing the number of subintervals + ## we should use to approximate `u`. + ## + ## * ``f`` - The function on the right hand side of the poisson + ## equation. + ## + ## OUTPUTS: + ## + ## * u - The function `u` that approximately solves -u_xx(x) = f(x). + ## + + [xs,h] = partition(integerN, 0, 1); + + ## 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) + + ## Prepend and append a zero row to f(xs). + ## Oh, and negate the whole thing (see the formula). + F = - cat(2, 0, f(differentiable_points), 0) + + ## 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_x4_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 x4. + A = cat(1, A, u_x4_coeffs); + + ## Solve the system. This gives us the value of u(x) at xs. + U = A \ F'; - [xs,h] = partition(integerN, a, b); + plot(xs,U) end -- 2.43.2