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;
[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
projects / octave.git / blobdiff