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