solve_poisson.m

   1 function u = solve_poisson(integerN, f)
   2   ##
   3   ## Numerically solve the poisson equation,
   4   ##
   5   ##   -u_xx(x) = f(x)
   6   ##
   7   ## in one dimension, subject to the boundary conditions,
   8   ##
   9   ##   u(0) = 0
  10   ##   u(1) = 1
  11   ##
  12   ## over the interval [0,1]. It is assumed that the function `f` is at
  13   ## least continuous on [0,1].
  14   ##
  15   ## INPUTS:
  16   ##
  17   ##   * ``integerN`` - An integer representing the number of subintervals
  18   ##     we should use to approximate `u`.
  19   ##
  20   ##   * ``f`` - The function on the right hand side of the poisson
  21   ##     equation.
  22   ##
  23   ## OUTPUTS:
  24   ##
  25   ##   * u - The function `u` that approximately solves -u_xx(x) = f(x).
  26   ##
  27
  28   [xs,h] = partition(integerN, 0, 1);
  29
  30   ## We cannot evaluate u_xx at the endpoints because our
  31   ## differentiation algorithm relies on the points directly to the left
  32   ## and right of `x`.
  33   differentiable_points = xs(2:end-1)
  34
  35   ## Prepend and append a zero row to f(xs).
  36   ## Oh, and negate the whole thing (see the formula).
  37   F = - cat(2, 0, f(differentiable_points), 0)
  38
  39   ## These are the coefficient vectors for the u(x0) and u(xn)
  40   ## constraints. There should be N zeros and a single 1.
  41   the_rest_zeros = zeros(1, integerN);
  42   u_x0_coeffs = cat(2, 1, the_rest_zeros);
  43   u_x4_coeffs = cat(2, the_rest_zeros, 1);
  44
  45   ## Start with the u(x0) row.
  46   A = u_x0_coeffs;
  47
  48   for x = differentiable_points
  49     ## Append each row obtained from the forward Euler method to A.
  50     u_row = forward_euler(2, xs, x);
  51     A = cat(1, A, u_row);
  52   end
  53
  54   ## Finally, append the last row for x4.
  55   A = cat(1, A, u_x4_coeffs);
  56
  57   ## Solve the system. This gives us the value of u(x) at xs.
  58   U = A \ F';
  59
  60   plot(xs,U)
  61 end