First shot at a solve_poisson function.

[octave.git] / solve_poisson.m

diff --git a/solve_poisson.m b/solve_poisson.m
index be1b1ec0faa5449ab27cc536f34f7e4c1e312c75..946d1d8aa582a39da300d7bf226046d82025714e 100644 (file)
--- 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