X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=advection_matrix.m;h=591d012a3388ae30a40a9692bad97c3722bc6c31;hp=5e76ee9ea19891250f262fcdd1675c334acfe971;hb=b12c6c2a4bf4cef29b2e08b743c92889505c7ed9;hpb=e1b71b4ca7cfa08ac76744a17a3778d4ccfaa7e2 diff --git a/advection_matrix.m b/advection_matrix.m index 5e76ee9..591d012 100644 --- a/advection_matrix.m +++ b/advection_matrix.m @@ -1,82 +1,82 @@ 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] + % + % 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; @@ -85,28 +85,28 @@ function S = advection_matrix(integerN, x0, xN) [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. + % 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. + % 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. + % 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. + % 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. + % Finally, append the last row for xN. S = cat(1, S, u_xN_coeffs); end