From: Michael Orlitzky
Date: Wed, 6 Feb 2013 19:47:35 +0000 (-0500)
Subject: Add the advection_matrix() function and a test for it.
X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=commitdiff_plain;h=b39006334c16dcd7256cb7c36e88996229863e6b
Add the advection_matrix() function and a test for it.
---
diff --git a/advection_matrix.m b/advection_matrix.m
new file mode 100644
index 0000000..5e76ee9
--- /dev/null
+++ b/advection_matrix.m
@@ -0,0 +1,112 @@
+function S = advection_matrix(integerN, x0, xN)
+ ##
+ ## The numerical solution of the advection-diffusion equation,
+ ##
+ ## -d*u''(x) + v*u'(x) + r*u = f(x)
+ ##
+ ## in one dimension, subject to the boundary conditions,
+ ##
+ ## u(x0) = u(xN)
+ ##
+ ## u'(x0) = u'(xN)
+ ##
+ ## over the interval [x0,xN] gives rise to a linear system:
+ ##
+ ## AU = h^2*F
+ ##
+ ## where h = 1/n, and A is given by,
+ ##
+ ## A = d*K + v*h*S + r*h^2*I.
+ ##
+ ## We will call the matrix S the "advection matrix," and it will be
+ ## understood that the first row (corresponding to j=0) is to be
+ ## omitted; since we have assumed that when j=0, u(xj) = u(x0) =
+ ## u(xN) and likewise for u'. ignored (i.e, added later).
+ ##
+ ## INPUTS:
+ ##
+ ## * ``integerN`` - An integer representing the number of
+ ## subintervals we should use to approximate `u`. Must be greater
+ ## than or equal to 2, since we have at least two values for u(x0)
+ ## and u(xN).
+ ##
+ ## * ``f`` - The function on the right hand side of the poisson
+ ## equation.
+ ##
+ ## * ``x0`` - The initial point.
+ ##
+ ## * ``xN`` - The terminal point.
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``S`` - The NxN matrix of coefficients for the vector [u(x1),
+ ## ..., u(xN)].
+ ##
+ ## EXAMPLES:
+ ##
+ ## For integerN=4, x0=0, and x1=1, we will have four subintervals:
+ ##
+ ## [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]
+ ##
+ ## The first row of the matrix 'S' should compute the "derivative"
+ ## at x1=0.25. By the finite difference formula, this is,
+ ##
+ ## u'(x1) = (u(x2) - u(x0))/2
+ ##
+ ## = (u(x2) - u(x4))/2
+ ##
+ ## Therefore, the first row of 'S' should look like,
+ ##
+ ## 2*S1 = [0, 1, 0, -1]
+ ##
+ ## and of course we would have F1 = [0] on the right-hand side.
+ ## Likewise, the last row of S should correspond to,
+ ##
+ ## u'(x4) = (u(x5) - u(x3))/2
+ ##
+ ## = (u(x1) - u(x3))/2
+ ##
+ ## So the last row of S will be,
+ ##
+ ## 2*S4 = [1, 0, -1, 0]
+ ##
+ ## Each row 'i' in between will have [-1, 0, 1] beginning at column
+ ## (i-1). So finally,
+ ##
+ ## 2*S = [0, 1, 0, -1]
+ ## [-1, 0, 1, 0]
+ ## [0, -1, 0, 1]
+ ## [1, 0, -1, 0]
+
+ if (integerN < 2)
+ S = NA;
+ return
+ end
+
+ [xs,h] = partition(integerN, x0, xN);
+
+ ## We cannot evaluate u_xx at the endpoints because our
+ ## differentiation algorithm relies on the points directly to the
+ ## left and right of `x`. Since we're starting at j=1 anyway, we cut
+ ## off two from the beginning.
+ differentiable_points = xs(3:end-1);
+
+ ## 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 - 3);
+ u_x0_coeffs = cat(2, the_rest_zeros, [0.5, 0, -0.5]);
+ u_xN_coeffs = cat(2, [0.5, 0, -0.5], the_rest_zeros);
+
+ ## Start with the u(x0) row.
+ S = u_x0_coeffs;
+
+ for x = differentiable_points
+ ## Append each row obtained from the forward Euler method to S.
+ ## Chop off x0 first.
+ u_row = central_difference(xs(2:end), x);
+ S = cat(1, S, u_row);
+ end
+
+ ## Finally, append the last row for xN.
+ S = cat(1, S, u_xN_coeffs);
+end
diff --git a/tests/advection_matrix_tests.m b/tests/advection_matrix_tests.m
new file mode 100644
index 0000000..7509b07
--- /dev/null
+++ b/tests/advection_matrix_tests.m
@@ -0,0 +1,13 @@
+## Try a 4x4 advection matrix. See the docstring in advection_matrix.m
+## for an explanation of the expected result.
+
+expected = (1/2) * [ 0, 1, 0, -1; ...
+ -1, 0, 1, 0; ...
+ 0, -1, 0, 1; ...
+ 1, 0, -1, 0 ];
+
+actual = advection_matrix(4, 0, 1);
+
+unit_test_equals("4x4 advection matrix is correct", ...
+ true, ...
+ expected == actual);