Add diffusion_matrix_sparse() and its tests.
[octave.git] / 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