three_point_centered_uni_d2.m

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

      function [D]=three_point_centered_uni_D2(z0,zL,n)
%...
%...  A. Vande Wouwer, P. Saucez and W.E. Schiesser (2002)
%...
%...  function three_point_centered_uni_D2 returns the differentiation matrix  
%...  for computing the second derivative, xzz, of a variable x over the spatial domain
%...  z0 < x < zL from classical three-point, second-order finite difference 
%...  approximations (this function can be used in combination with four_point_uni_neumann_bc, 
%...  and replaces dss042)
%...
%...  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)
%...
%...  four-point approximations of the second derivative are used at the
%...  boundaries of a spatial domain.  these approximations have the
%...  features
%...
%...     (1)  only interior and boundary points are used (i.e., no
%...          fictitious points are used)
%...
%...     (2)  approximations for the boundary conditions are not used.
%...
%...  the derivation is by professor Gilbert A. Stengle, Department of
%...  Mathematics, Lehigh University, Bethlehem, PA 18015, and was done
%...  on december 7, 1985.
%...
%...  for an approximation at the left boundary, involving the points
%...  i, i+1, i+2 and i+3, consider the following taylor series expan-
%...  sions
%...
%...                  xz(i)( dz)   xzz(i)( dz)**2   xzzz(i)( dz)**3
%...  x(i+1) = x(i) + ---------- + -------------- + --------------- +...
%...                       1             2                 6
%...
%...
%...                  xz(i)(2dz)   xzz(i)(2dz)**2   xzzz(i)(2dz)**3
%...  x(i+2) = x(i) + ---------- + -------------- + --------------- +...
%...                       1             2                 6
%...
%...
%...  to obtain approximations of the second derviavtive which do not
%...  involve the first derivative, we take as the linear combination
%...
%...     a*x(i) + b*x(i+1) + c*x(i+2) + d*x(i+3) =
%...
%...  we have
%...
%...     a*x(i) + b*x(i+1) + c*x(i+2) + d*x(i+3) =
%...
%...     (a + b + c + d)*x(i)+
%...
%...     (dz*b + 2*dz*c + 4*dz*d)*xz(i)+
%...
%...     (b*(dz**2)/2 + c*((2*dz)**2)/2 + d*((3*dz)**2)/2)*xzz(i) +
%...
%...     (b*(dz**3)/6 + c*((2*dz)**3)/6 + d*((3*dz)**3)/6)*xzz(i) +
%...
%...     o(dz**4)
%...
%...  the third derivative, xzzz(i), can be dropped by taking
%...
%...     b + 8*c + 27*d = 0
%...
%...  the second derivative, xzz(i), can be retained by taking
%...
%...     (dz**2)*(b/2 + 2*c + (9/2)*d) = 1
%...
%...  the first derivative can be dropped by taking
%...
%...     b + 2*c + 3*d = 0
%...
%...  solution of the preceding equations for c and d by elimination of
%...  b gives
%...
%...     6*c + 24*d = 0
%...
%...     4*c + 18*d = -2/(dz**2)
%...
%...  then, eliminating c, gives
%...
%...     (18 - 16)*d = -2/(dz**2)
%...
%...  or
%...
%...     d = -1/(dz**2)
%...
%...     c = (24/6)/(dz**2) = 4/(dz**2)
%...
%...     b = -8/(dz**2) + 3/(dz**2) = -5/(dz**2)
%...
%...  x(i) can be dropped by taking
%...
%...     a + b + c + d = 0
%...
%...  or
%...
%...     a = (5 - 4 + 1)/(dz**2) = 2/(dz**2)
%...
%...  if we now solve for the derivative of interest, xzz(i),
%...
%...     xzz(i) = (1/dz**2)*(2*x(i) - 5*x(i+1) + 4*x(i+2) - 1*x(i+3))
%...
%...            + o(dz**2)
%...
%...  which is the four-point, second-order approximation for the second
%...  derivative, xzz(i), without the first derivative, xz(i).
%...
%...  four checks of this approximation can easily be made for x(i) =
%...  1, x(i) = z, x(i) = z**2 and x(i) = z**3
%...
%...     xzz(i) = (1/dz**2)*(2 - 5 + 4 - 1) = 0
%...
%...     xzz(i) = (1/dz**2)*(2*z - 5*(z + dz) + 4*(z + 2*dz)
%...
%...              - 1*(z + 3*dz)) = 0
%...
%...     xzz(i) = (1/dz**2)*(2*z**2 - 5*(z + dz)**2 + 4*(z + 2*dz)**2
%...
%...              - 1*(z + 3*dz)**2) = 2
%...
%...     xzz(i) = (1/dz**2)*(2*z**3 - 5*(z + dz)**3 + 4*(z + 2*dz)**3
%...
%...            - 1*(z + 3*dz)**3)
%...
%...             = (1/dz**2)*(2*z**3 - 5*z**3 - 15*z*dz**2
%...
%...             - 15*dz*z**2 - 5*dz**3 + 4*z**3 + 24*dz*z**2
%...
%...             + 48*z*dz**2 + 32*dz**3 - z**3 - 9*dz*z**2
%...
%...             - 27*z*dz**2 - 27dz**3)
%...
%...             = (1/dz**2)*(6*z*dz**2) = 6*z
%...
%...  the preceding approximation for xzz(i) can be applied at the
%...  left boundary value of z by taking i = 1.  an approximation at
%...  the right boundary is obtained by taking dz = -dz and reversing
%...  the subscripts in the preceding approximation, with i = n
%...
%...     xzz(i) = (1/dz**2)*(2*x(i) - 5*x(i-1) + 4*x(i-2) - 1*x(i-3))
%...
%...            + o(dz**2)
%...                    
%...
%...  compute the spatial increment
      dz=(zL-z0)/(n-1);
      rdzs=1/(dz^2);
%...
%...  sparse discretization matrix      
%...
%...  interior points      
      D=diag(+1*ones(n-1,1),-1)+diag(-2*ones(n,1),0)+diag(+1*ones(n-1,1),1);      
%...
%...  boundary points      
      D(1,1:4) = [+2 -5 +4 -1];
      D(n,(n-3):n) = [-1 +4 -5 +2];
%...      
      D=rdzs*D;
      D=sparse(D);