advection_matrix.m

   1 function S = advection_matrix(integerN, x0, xN)
   2   %
   3   % The numerical solution of the advection-diffusion equation,
   4   %
   5   %   -d*u''(x) + v*u'(x) + r*u = f(x)
   6   %
   7   % in one dimension, subject to the boundary conditions,
   8   %
   9   %   u(x0)  = u(xN)
  10   %
  11   %   u'(x0) = u'(xN)
  12   %
  13   % over the interval [x0,xN] gives rise to a linear system:
  14   %
  15   %   AU = h^2*F
  16   %
  17   % where h = 1/n, and A is given by,
  18   %
  19   %   A = d*K + v*h*S + r*h^2*I.
  20   %
  21   % We will call the matrix S the "advection matrix," and it will be
  22   % understood that the first row (corresponding to j=0) is to be
  23   % omitted; since we have assumed that when j=0, u(xj) = u(x0) =
  24   % u(xN) and likewise for u'. ignored (i.e, added later).
  25   %
  26   % INPUTS:
  27   %
  28   %   * ``integerN`` - An integer representing the number of
  29   %   subintervals we should use to approximate `u`. Must be greater
  30   %   than or equal to 2, since we have at least two values for u(x0)
  31   %   and u(xN).
  32   %
  33   %   * ``f`` - The function on the right hand side of the poisson
  34   %   equation.
  35   %
  36   %   * ``x0`` - The initial point.
  37   %
  38   %   * ``xN`` - The terminal point.
  39   %
  40   % OUTPUTS:
  41   %
  42   %   * ``S`` - The NxN matrix of coefficients for the vector [u(x1),
  43   %   ..., u(xN)].
  44   %
  45   % EXAMPLES:
  46   %
  47   % For integerN=4, x0=0, and x1=1, we will have four subintervals:
  48   %
  49   %   [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]
  50   %
  51   % The first row of the matrix 'S' should compute the "derivative"
  52   % at x1=0.25. By the finite difference formula, this is,
  53   %
  54   %   u'(x1) = (u(x2) - u(x0))/2
  55   %
  56   %          = (u(x2) - u(x4))/2
  57   %
  58   % Therefore, the first row of 'S' should look like,
  59   %
  60   %   2*S1 = [0, 1, 0, -1]
  61   %
  62   % and of course we would have F1 = [0] on the right-hand side.
  63   % Likewise, the last row of S should correspond to,
  64   %
  65   %   u'(x4) = (u(x5) - u(x3))/2
  66   %
  67   %          = (u(x1) - u(x3))/2
  68   %
  69   % So the last row of S will be,
  70   %
  71   %   2*S4 = [1, 0, -1, 0]
  72   %
  73   % Each row 'i' in between will have [-1, 0, 1] beginning at column
  74   % (i-1). So finally,
  75   %
  76   %   2*S = [0, 1, 0, -1]
  77   %          [-1, 0, 1, 0]
  78   %         [0, -1, 0, 1]
  79   %         [1, 0, -1, 0]
  80
  81   if (integerN < 2)
  82     S = NA;
  83     return
  84   end
  85
  86   [xs,h] = partition(integerN, x0, xN);
  87
  88   % We cannot evaluate u_xx at the endpoints because our
  89   % differentiation algorithm relies on the points directly to the
  90   % left and right of `x`. Since we're starting at j=1 anyway, we cut
  91   % off two from the beginning.
  92   differentiable_points = xs(3:end-1);
  93
  94   % These are the coefficient vectors for the u(x0) and u(xn)
  95   % constraints. There should be N zeros and a single 1.
  96   the_rest_zeros = zeros(1, integerN - 3);
  97   u_x0_coeffs = cat(2, the_rest_zeros, [0.5, 0, -0.5]);
  98   u_xN_coeffs = cat(2, [0.5, 0, -0.5], the_rest_zeros);
  99
 100   % Start with the u(x0) row.
 101   S = u_x0_coeffs;
 102
 103   for x = differentiable_points
 104     % Append each row obtained from the forward Euler method to S.
 105     % Chop off x0 first.
 106     u_row = central_difference(xs(2:end), x);
 107     S = cat(1, S, u_row);
 108   end
 109
 110   % Finally, append the last row for xN.
 111   S = cat(1, S, u_xN_coeffs);
 112 end