poisson_matrix.m

   1 function A = poisson_matrix(integerN, x0, xN)
   2   ##
   3   ## In the numerical solution of the poisson equation,
   4   ##
   5   ##   -u''(x) = f(x)
   6   ##
   7   ## in one dimension, subject to the boundary conditions,
   8   ##
   9   ##   u(x0) = 0
  10   ##   u(xN) = 1
  11   ##
  12   ## over the interval [x0,xN], we need to compute a matrix A which
  13   ## is then multiplied by the vector of u(x0), ..., u(xN). The entries
  14   ## of A are determined from the second order finite difference formula,
  15   ## i.e. the second order forward Euler method.
  16   ##
  17   ## INPUTS:
  18   ##
  19   ##   * ``integerN`` - An integer representing the number of
  20   ##     subintervals we should use to approximate `u`. Must be
  21   ##     greater than or equal to 2, since we have at least two
  22   ##     values for u(x0) and u(xN).
  23   ##
  24   ##   * ``f`` - The function on the right hand side of the poisson
  25   ##     equation.
  26   ##
  27   ##   * ``x0`` - The initial point.
  28   ##
  29   ##   * ``xN`` - The terminal point.
  30   ##
  31   ## OUTPUTS:
  32   ##
  33   ##   * ``A`` - The (N+1)x(N+1) matrix of coefficients for u(x0),
  34   ##             ..., u(xN).
  35
  36   if (integerN < 2)
  37     A = NA
  38     return
  39   end
  40
  41   [xs,h] = partition(integerN, x0, xN);
  42
  43   ## We cannot evaluate u_xx at the endpoints because our
  44   ## differentiation algorithm relies on the points directly to the left
  45   ## and right of `x`.
  46   differentiable_points = xs(2:end-1);
  47
  48   ## These are the coefficient vectors for the u(x0) and u(xn)
  49   ## constraints. There should be N zeros and a single 1.
  50   the_rest_zeros = zeros(1, integerN);
  51   u_x0_coeffs = cat(2, 1, the_rest_zeros);
  52   u_xN_coeffs = cat(2, the_rest_zeros, 1);
  53
  54   ## Start with the u(x0) row.
  55   A = u_x0_coeffs;
  56
  57   for x = differentiable_points
  58     ## Append each row obtained from the forward Euler method to A.
  59     u_row = forward_euler(2, xs, x);
  60     A = cat(1, A, u_row);
  61   end
  62
  63   ## Finally, append the last row for xN and negate the whole thing (see
  64   ## the definition of the poisson problem).
  65   A = - cat(1, A, u_xN_coeffs);
  66 end