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