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]

81 if (integerN < 2)

82 S = NA;

83 return

84 end

86 [xs,h] = partition(integerN, x0, xN);

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);

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);

100 ## Start with the u(x0) row.

101 S = u_x0_coeffs;

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

110 ## Finally, append the last row for xN.

111 S = cat(1, S, u_xN_coeffs);

112 end