+++ /dev/null
-function u = solve_poisson(integerN, f)
- ##
- ## Numerically solve the poisson equation,
- ##
- ## -u_xx(x) = f(x)
- ##
- ## in one dimension, subject to the boundary conditions,
- ##
- ## u(0) = 0
- ## u(1) = 1
- ##
- ## 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';
-
- plot(xs,U)
-end