dss002

%  File: dss002.m
%
   function [ux]=dss002(xl,xu,n,u)
%
%  Function dss002 computes the first derivative, u , of a
%                                                  x
%  variable u over the spatial domain xl le x le xu.
%
%  Argument list
%
%     xl      Lower boundary value of x (input)
%
%     xu      Upper boundary value of x (input)
%
%     n       Number of grid points in the x domain including the
%             boundary points (input)
%
%     u       One-dimensional array containing the values of u at
%             the n grid point points for which the derivative is
%             to be computed (input)
%
%     ux      One-dimensional array containing the numerical
%             values of the derivatives of u at the n grid points
%             (output)
%
%  Function dss002 computes the first derivative, u , of a
%                                                  x
%  variable u over the spatial domain xl le x le xu from the
%  classical three-point, second-order finite difference approxi-
%  tions
%
%                                       2
%  u1  = (1/2dx)(-3u1 + 4u2 - u3) + O(dx ) (left boundary,     (1)
%    x                                         x = xl)
%
%                                   2
%  ui  = (1/2dx)(ui+1 - ui-1) + O(dx ) (interior point,        (2)
%    x                                   x ne xl, xu)
%
%                                          2
%  un  = (1/2dx)(3un - 4un-1 + un-2) + O(dx ) (right boundary, (3)
%    x                                            x = xu)
%
%  Equations (1) to (3) apply over a grid in x with corresponding
%  values of the function u(x) represented as
%
%   u1      u2       u3         ui        un-2      un-1    un
%
%  x=xl  x=xl+dx  x=xl+2dx ... X=xi ... X=xu-2dx  x=xu-dx  x=xu
%
%  The origin of equations (1) to (3) is outlined below.
%
%  Consider the following polynomial in x of arbitrary order
%
%                                     2             3
%  u(x) = a0 + a1(x - x0) + a2(x - x0)  + a3(x - x0)  + ....   (4)
%
%  We seek the values of the coefficients a0, a1, a2, ... for a
%  particular function u(x).  If x = x0 is substituted in equation
%  (4), we have immediately a0 = u(x0).  Next, if equation (4) is
%  differentiated with respect to x,
%
%                                                   2
%  du(x)/dx = u (x) = a1 + 2a2(x - x0) + 3a3(x - x0)  + ...    (5)
%              x
%
%  Again, with x = x0, a1 = du(x0)/dx = u (x0).  Differentiation
%                                        x
%  of equation (5) in turn gives
%
%  d2u(x)/dx2 = u  (x) = 2a2 + 6a3(x - x0) + ...
%                2x
%
%  And for x = x0, a2 = u  (x0)/2f (2f = 1*2, i.e., 2 factorial).
%                        2x
%
%  We can continue this process of differentiation followed by the
%  substitution x = x0 to obtain the successive coefficients in
%  equation (4), a3, a4, ...  Finally, substitution of these co-
%  efficients in equation (4) gives
%
%                                                 2
%  u(x) = u(x0) + u (x0)(x - x0) + u  (x0)(x - x0)  +
%                  x       1f       2x       2f
%                                                              (6)
%                                3                  4
%                 U  (x0)(x - x0)  + u  (x0)(x - x0)  + ...
%                  3x       3f        4x       4f
%
%  The correspondence between equation (6) and the well-known
%  Taylor series should be clear.  Thus the expansion of a
%  function, u(x), around a neighboring point x0 in terms of u(x0)
%  and the derivatives of u(x) at x = x0 is equivalent to approxi-
%  mating u(x) near x0 by a polynomial.
%
%  Equation (6) is the starting point for the derivation of the
%  classical finite difference approximations of derivatives such
%  as the three-point formulas of equations (1), (2) and (3).  We
%  will now consider the derivation of these three-point formulas
%  in a standard format that can then be extended to higher
%  multi-point formulas in other subroutines, e.g., five-point
%  formulas in routine dss004.
%
%  Three-point formulas
%
%     (1)  Left end, point i = 1
%
%  If equation (6) is written around the points x = xl for x = xl +
%  dx and x = xl + 2dx, for which the corresponding values of u(x)
%  are u1, u2 and u3 (u1 and u2 are separated with respect to x by
%  distance dx as are u2 and u3, i.e., we assume a uniform grid
%  spacing, dx, for independent variable x)
%
%                                2            3
%  u2 = u1 + u1 ( dx) + u1  ( dx)  + u1  ( dx)  + ...          (7)
%              x  1f      2x  2f       3x  3f
%
%                                2            3
%  u3 = u1 + u1 (2dx) + u1  (2dx)  + u1  (2dx)  + ...          (8)
%              x  1f      2x  2f       3x  3f
%
%  We can now take a linear combination of equations (7) and (8)
%  by first multiplying equation (7) by a constant, a, and equa-
%  tion (8) by constant b
%
%                                  2           3
%  a(u2 = u1 + u1 ( dx) + u1  ( dx) + u1  ( dx) + ...)         (9)
%                x  1f      2x  2f      3x  3f
%
%                                  2           3
%  b(u3 = u1 + u1 (2dx) + u1  (2dx) + u1  (2dx) + ...)        (10)
%                x  1f      2x  2f      3x  3f
%
%  Constants a and b are then selected so that the coefficients of
%  the u1  terms sum to one (since we are interested in obtaining
%        x
%  a finite difference approximation for this first derivative).
%  Also, we select a and b so that the coefficients of the u1
%                                                            2x
%  terms sum to zero in order to drop out the contribution of this
%  second derivative (the basic idea is to drop out as many of the
%  derivatives as possible in the Taylor series beyond the deri-
%  vative of interest, in this case u1 , in order to produce a
%                                     x
%  finite difference approximation for the derivative of maximum
%  accuracy).  In this case we have only two constants, a and b,
%  to select so we can drop out only the second derivative, u1  ,
%                                                             2x
%  in the Taylor series (in addition to retaining the first deri-
%  vative).  This procedure leads to two linear algebraic equa-
%  tions in the two constants
%
%  a + 2b = 1
%
%  a + 4b = 0
%
%  Solution of these equations for a and b gives
%
%  a = 2, b = -1/2
%
%  Solution of equations (9) and (10) for u1  with these values of
%  a and b gives equation (1)               x
%
%                                      2
%  u1 = (1/2dx)(-3u1 + 4u2 - u3) + O(dx )                      (1)
%    x
%               2
%  The term O(dx ) indicates a principal error term due to trunca-
%                                                2
%  tion of the Taylor series which is of order dx .  This term in
%                    2
%  fact equals u1  dx /3f, which is easily obtained in deriving
%                3x
%  equation (1).
%
%  This same basic procedure can now be applied to the derivation
%  of equations (2) and (3).
%
%     (2)  Interior point i
%
%                                    2           3
%  a(ui-1 = ui + ui (-dx) + ui  (-dx) + ui  (-dx) + ...)
%                  x  1f      2x  2f      3x  3f
%
%                                    2           3
%  b(ui+1 = ui + ui ( dx) + ui  ( dx) + ui  ( dx) + ...)
%                  x  1f      2x  2f      3x  3f
%
%  -a + b = 1
%
%   a + b = 0
%
%  a = 1/2, b = -1/2
%                                   2
%  ui  = (1/2dx)(ui+1 - ui-1) + O(dx )                         (2)
%    x
%
%     (3)  Right end, point i = n
%
%                                      2            3
%  a(un-2 = un + un (-2dx) + un  (-2dx) + un  (-2dx) + ...)
%                  X   1f      2x   2f      3x   3f
%
%                                      2            3
%  b(un-1 = un + un ( -dx) + un  ( -dx) + un  ( -dx) + ...)
%                  x   1f      2x   2f      3x   3f
%
%  -2a - b = 1
%
%   4a + b = 0
%
%   a = -2, b = 1/2
%                                          2
%  un  = (1/2dx)(3un - 4un-1 + un-2) + O(dx )                  (3)
%    x
%
%  The weighting coefficients for equations (1), (2) and (3) can
%  be summarized as
%
%          -3   4  -1
%
%     1/2  -1   0   1
%
%           1  -4   3
%
%  Which are the coefficients reported by Bickley for n = 2, m =
%  1, p = 0, 1, 2 (Bickley, W. G., Formulae for Numerical Differ-
%  entiation, Math. Gaz., vol. 25, 1941).
%
%  Equations (1), (2) and (3) can now be programmed to generate
%  the derivative u (x) of u(x).
%                  x
%
%  Compute the spatial increment
   dx=(xu-xl)/(n-1);
   r2fdx=1./(2.*dx);
   nm1=n-1;
%
%  Equation (1) (note - the rhs of the finite difference approxi-
%  tions, equations (1), (2) and (3) have been formatted so that
%  the numerical weighting coefficients can be more easily associ-
%  ated with the Bickley matrix listed above)
   ux(1)=r2fdx*...
     (     -3.   *u(  1)     +4.   *u(  2)     -1.   *u(  3));
%
%  Equation (2)
   for i=2:nm1
     ux(i)=r2fdx*...
     (     -1.   *u(i-1)     +0.   *u(  i)     +1.   *u(i+1));
   end
%
%  Equation (3)
   ux(n)=r2fdx*...
     (      1.   *u(n-2)     -4.   *u(n-1)     +3.   *u(  n));
%