five_point_centered_uni_d2.m

有关the method of line即MOL算法的matlab代码

      function [D]=five_point_centered_uni_D2(z0,zL,n)
%...
%...  A. Vande Wouwer, P. Saucez and W.E. Schiesser (2002)
%...
%...  function five_point_centered_uni_D2 returns the differentiation matrix  
%...  for computing the second derivative, xzz, of a variable x over the spatial 
%...  z0 < z < zL from classical five-point, fourth-order finite difference 
%...  domain approximations (this function can be used in combination with  
%...  six_point_uni_neumann_bc and replaces dss044)
%...
%...  argument list
%...
%...     z0      left value of the spatial independent variable (input)
%...
%...     zL      right value of the spatial independent variable (input)
%...
%...     n       number of spatial grid points, including the end
%...             points (input)
%...
%...  the following derivation was completed by W.E. Schiesser, Depts
%...  of CHE and CSEE, Lehigh University, Bethlehem, PA 18015, USA, on
%...  December 15, 1986.
%...
%...  ******************************************************************
%...
%...  (1)  xzz at the interior points 3, 4,..., n-2
%...
%...  to develop a set of fourth-order correct differentiation formulas
%...  for the second derivative xzz, we consider first the interior
%...  grid points at which a symmetric formula can be used.
%...
%...  if we consider a formula of the form
%...
%...     a*x(i-2) + b*x(i-1) + e*x(i) + c*x(i+1) + d*x(i+2)
%...
%...  taylor series expansions of x(i-2), x(i-1), x(i+1) and x(i+2)
%...  can be substituted into this formula.  we then consider the
%...  linear albegraic equations relating a, b, c and d which will
%...  retain certain terms, i.e., xzz, and drop others, e.g., xzzz,
%...  xzzzz and xzzzzz.
%...
%...  thus, for grid points 3, 4,..., n-2
%...
%...     to retain xzz
%...
%...        4*a + b + c +  4*d = 2                              (1)
%...
%...     to drop xzzz
%...
%...       -8*a - b + c +  8*d = 0                              (2)
%...
%...     to drop xzzzz
%...
%...       16*a + b + c + 16*d = 0                              (3)
%...
%...     to drop xzzzzz
%...
%...      -32*a - b + c + 32*d = 0                              (4)
%...
%...  equations (1) to (4) can be solved for a, b, c and d.  if equa-
%...  tion (1) is added to equation (2)
%...
%...        -4*a + 2*c + 12*d = 2                               (5)
%...
%...  if equation (1) is subtracted from equation (3)
%...
%...        12*a + 12*d = -2                                    (6)
%...
%...  if equation (1) is added to equation (4)
%...
%...        -28*a + 2*c + 36*d = 2                              (7)
%...
%...  equations (5) to (7) can be solved for a, c and d.  if equation
%...  (5) is subtracted from equation (7), and the result combined
%...  with equation (6)
%...
%...         12*a + 12*d = -2                                   (6)
%...
%...        -24*a + 24*d = 0                                    (8)
%...
%...  equations (6) and (8) can be solved for a and d.  from (8), a
%...  = d.  from equation (6), a = -1/12 and d = -1/12.  then, from
%...  equation (5), c = 4/3, and from equation (1), b = 4/3.
%...
%...  the final differentiation formula is then obtained as
%...
%...     (-1/12)*x(i-2) +   (4/3)*x(i-1) +
%...
%...       (4/3)*x(i+1) + (-1/12)*x(i+2)
%...
%...     (-1/12 + 4/3 - 1/12 + 4/3)*x(i) + xzz(i)*(dz**2) + o(dz**6)
%...
%...  or
%...
%...     xzz(i) = (1/(12*dz**2))*(-1*x(i-2) + 16*x(i-1)
%...
%...                  - 30*x(i) + 16*x(i+1) -  1*x(i+2)         (9)
%...
%...                  + o(dz**4)
%...
%...  note that the xz term drops out, i.e., the basic equation is
%...
%...        -2*a - b + c + 2*d =
%...
%...        -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0
%...
%...  equation (9) was obtained by dropping all terms in the underlying
%...  taylor series up to and including the fifth derivative, xzzzzz.
%...  thus, equation (9) is exact for polynomials up to and including
%...  fifth order.  this can be checked by substituting the functions
%...  1, z, z**2, z**3, z**4 and z**5 in equation (9) and computing the
%...  corresponding derivatives for comparison with the known second
%...  derivatives.  this is done for 1 merely by summing the weighting
%...  coefficients in equation (9), which should sum to zero, i.e.,
%...  -1 + 16 - 30 + 16 -1 = 0.
%...
%...  for the remaining functions, the algebra is rather involved, but
%...  these functions can be checked numerically, i.e., numerical values
%...  of z**2, z**3, z**4 and z**5 can be substituted in equation (9)
%...  and the computed derivatives can be compared with the know numeri-
%...  cal second derivatives.  this is not a proof of correctness of
%...  equation (9), but would likely detect any errors in equation (9).
%...
%...
%...  ******************************************************************
%...
%...  (2)  xzz at the interior points i = 2 and n-1
%...
%...  for grid point 2, we consider a formula of the form
%...
%...     a*x(i-1) + f*x(i) + b*x(i+1) + c*x(i+2) + d*x(i+3) + e*x(i+4)
%...
%...  taylor series expansions of x(i-1), x(i+1), x(i+2), x(i+3) and
%...  x(i+4) when substituted into this formula give linear algebraic
%...  equations relating a, b, c, d and e.
%...
%...     to drop xz
%...
%...        -a + b + 2*c + 3*d + 4*e = 0                       (10)
%...
%...     to retain xzz
%...
%...        a + b + 4*c + 9*d + 16*e = 2                       (11)
%...
%...     to drop xzzz
%...
%...       -a + b + 8*c + 27*d + 64*e = 0                      (12)
%...
%...     to drop xzzzz
%...
%...        a + b + 16*c + 81*d + 256*e = 0                    (13)
%...
%...     to drop xzzzzz
%...
%...       -a + b + 32*c + 243*d + 1024*e = 0                  (14)
%...
%...  equations (11), (12), (13) and (14) can be solved for a, b, c,
%...  d and e.  if equation (10) is added to equation (11)
%...
%...     2*b + 6*c + 12*d +20*e = 2                            (15)
%...
%...  if equation (10) is subtracted from equation (12)
%...
%...     6*c + 24*d + 60*e = 0                                 (16)
%...
%...  if equation (10) is added to equation (13)
%...
%...     2*b + 18*c + 84*d + 260*e = 0                         (17)
%...
%...  if equation (10) is subtracted from equation (14)
%...
%...     30*c + 240*d + 1020*e = 0                             (18)
%...
%...  equations (15), (16), (17) and (18) can be solved for b, c, d
%...  and e.
%...
%...     6*c + 24*d + 60*e = 0                                 (16)
%...
%...  if equation (15) is subtracted from equation (17)
%...
%...     12*c + 72*d + 240*e = -2                              (19)
%...
%...     30*c + 240*d + 1020*e = 0                             (18)
%...
%...  equations (16), (18) and (19) can be solved for c, d and e.  if
%...  two times equation (16) is subtracted from equation (19),
%...
%...     24*d + 120*e = -2                                     (20)
%...
%...  if five times equation (16) is subtracted from equation (18),
%...
%...     120*d + 720*e = 0                                     (21)
%...
%...  equations (20) and (21) can be solved for d and e.  from (21),
%...  e = (-1/6)*d.  substitution in equation (20) gives d = -1/2.
%...  thus, e = 1/12.  from equation (16), c = 7/6.  from equation
%...  (15), b = -1/3.  from equation (10), a = 5/6.
%...
%...  the final differentiation formula is then obtained as
%...
%...  (5/6)*x(i-1) + (-1/3)*x(i+1) + (7/6)*x(i+2) + (-1/2)*x(i+3)
%...
%...  + (1/12)*x(i+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*x(i)
%...
%...  + xzz*(dz**2) + o(dz**6)
%...
%...  or
%...
%...  xzz(i) = (1/12*dz**2)*(10*x(i-1) - 15*x(i) - 4*x(i+1)
%...                                                           (22)
%...         + 14*x(i+2) - 6*x(i+3) + 1*x(i+4)) + o(dz**4)
%...
%...  equation (22) will be applied at i = 2 and n-1.  thus
%...
%...  xzz(2) = (1/12*dz**2)*(10*x(1) - 15*x(2) - 4*x(3)
%...                                                           (23)
%...         + 14*x(4) - 6*x(5) + 1*x(6)) + o(dz**4)
%...
%...  xzz(n-1) = (1/12*dz**2)*(10*x(n) - 15*x(n-1) - 4*x(n-2)
%...                                                           (24)
%...           + 14*x(n-3) - 6*x(n-4) + 1*x(n-5)) + o(dz**4)
%...
%...
%...  ******************************************************************
%...
%...  (3)  xzz at the boundary points 1 and n
%...
%...  finally, for grid point 1, an approximation with a dirichlet
%...  boundary condition of the form
%...
%...  a*x(i+1) + b*x(i+2) + c*x(i+3) + d*x(i+4) + e*x(i+5) + f*x(i)
%...
%...  can be used.  the corresponding algebraic equations are
%...
%...     to drop xz
%...
%...        a + 2*b + 3*c + 4*d + 5*e = 0                      (25)
%...
%...     to retain xzz
%...
%...        a + 4*b + 9*c + 16*d + 25*e = 2                    (26)
%...
%...     to drop xzzz
%...
%...        a + 8*b + 27*c + 64*d + 125*e = 0                  (27)
%...
%...     to drop xzzzz
%...
%...        a + 16*b + 81*c + 256*d + 625*e = 0                (28)
%...
%...     to drop xzzzzz
%...
%...        a + 32*b + 243*c + 1024*d + 3125*e = 0             (29)
%...
%...  equations (25), (26), (27), (28) amd (29) can be solved for a,
%...  b, c, d and e.
%...
%...        2*b + 6*c + 12*d + 20*e = 2                        (30)
%...
%...        6*b + 24*c + 60*d + 120*e = 0                      (31)
%...
%...        14*b + 78*c + 252*d + 620*e = 0                    (32)
%...
%...        30*b + 240*c + 1020*d + 3120*e = 0                 (33)
%...
%...  equations (30), (31), (32) and (33) can be solved for b, c, d
%...  and e
%...
%...        6*c + 24*d + 60*e = -6                             (34)
%...
%...        36*c + 168*d + 480*e = -14                         (35)
%...
%...        150*c + 840*d + 2820*e = -30                       (36)
%...
%...  equations (34), (35) and (36) can be solved for c, d and e
%...
%...        24*d + 120*e = 22                                  (37)
%...
%...        240*d + 1320*e = 120                               (38)
%...
%...  from equations (37) and (38), d = 61/12, e = -5/6.  from equation
%...  (34), c = -13.  from equation (30), b = 107/6.  from equation
%...  (25), a = -77/6.
%...
%...  the final differentiation formula is then obtained as
%...
%...  (-77/6)*x(i+1) + (107/6)*x(i+2) - 13*x(i+3) + (61/12)*x(i+4)
%...
%...  - (5/6)*x(i+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*x(i) +
%...
%...  xzz(i)*(dz**2) + o(dz**6)
%...
%...  or
%...
%...  xzz(i) = (1/12*dz**2)*(45*x(i) - 154*x(i+1) + 214*x(i+2)
%...                                                                (39)
%...         - 156*x(i+3) + 61*x(i+4) - 10*x(i+5)) + o(dz**4)
%...
%...  equation (52) will be applied at i = 1 and i = n
%...
%...  xzz(1) = (1/12*dz**2)*(45*x(1) - 154*x(2) + 214*x(3)
%...                                                                (40)
%...         - 156*x(4) + 61*x(5) - 10*x(6)) + o(dz**4)
%...
%...  xzz(n) = (1/12*dz**2)*(45*x(n) - 154*x(n-1) + 214*x(n-2)
%...                                                                (41)
%...         -156*x(n-3) + 61*x(n-4) - 10*x(n-5)) + o(dz**4)
%...
%...
%...  ******************************************************************
%...
%...  compute the spatial increment
      dz=(zL-z0)/(n-1);
      r12dzs=1/(12*dz^2);
%...
%...  sparse discretization matrix      
%...
%...  interior points (i = 3, 4,..., n-2 (equation 9))     
      D=diag(-1*ones(n-2,1),-2)+diag(+16*ones(n-1,1),-1)+diag(-30*ones(n,1),0)+diag(+16*ones(n-1,1),1)+diag(-1*ones(n-2,1),2);      
%...
%...  i = 2 (equation 23) and i = n-1 (equation 24)
      D(2,1:6) = [+10 -15 -4 +14 -6 +1];
      D(n-1,(n-5):n) = [+1 -6 +14 -4 -15 +10];
%...
%...  boundary points      
      D(1,1:6) = [+45 -154 +214 -156 +61 -10];
      D(n,(n-5):n) = [-10 +61 -156 +214 -154 +45];
%...      
      D=r12dzs*D;
      D=sparse(D);