diff2_neumann_bcs_4thorder.m

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

      function [xzz_z0,xzz_zL]=diff2_neumann_bcs_4thorder(z0,zL,x,xz_z0,xz_zL,n0,nL)
%...
%...  Authors: A. Vande Wouwer, P. Saucez, W.E. Schiesser (2002)
%...
%...  function diff2_neumann_bcs_4thorder computes a six-point approximation of
%...  the second derivative at the boundaries, taking the value of the first 
%...  derivative at the boundary (Neumann boundary condition) into account.
%...
%...  argument list
%...
%...     x       one-dimensional array of the dependent variable
%...             to be differentiated (input)
%...
%...     xz_z0   value of the first derivative at the left boundary 
%...             z = z0 (input). (ignored if nl different from 1)
%...
%...     xz_zL   value of the first derivative at the right boundary 
%...             z = zL (input). (ignored if nu different from 1)
%...
%...     n0      integer flag for Neumann boundary condition at
%...             z = z0 (input).  (uxl=ux(1) is used if nl = 1)
%...
%...     nL      integer flag for Neumann of boundary condition at
%...             z = zL (input).  (uxu=ux(n) is used if nu = 1)
%...
%...     xzz_z0  value of the second derivative at the left boundary 
%...             z = z0 (output)
%...
%...     xzz_zL  value of the second derivative at the right boundary 
%...             z = zL (output)
%...
%...  origin of the approximations
%...
%...  the following derivation is for a set of fourth-order, six-point
%...  approximations for a second derivative that can be 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)  the normal derivative at the boundary is included as part
%...          of the approximation for the second derivative
%...
%...     (3)  approximations for the boundary conditions are not used.
%...
%...  the following derivation was completed by W.E. Schiesser, depts
%...  of che and csee, Lehigh University, Bethlehem, PA 18015, USA, on
%...  December 15, 1986.
%...
%...  ******************************************************************
%...
%...  for grid point 1, an approximation with a neumann boundary
%...  condition of the form
%...
%...     a*x(i+1) + b*x(i+2) + c*x(i+3) + d*x(i+4) + e*xz(i) + f*x(i)
%...
%...  will be used.  the corresponding algebraic equations are
%...
%...     to drop xz
%...
%...        a + 2*b + 3*c + 4*d + e = 0                        (1)
%...
%...     to retain xzz
%...
%...        a + 4*b + 9*c + 16*d = 2                           (2)
%...
%...     to drop xzzz
%...
%...        a + 8*b + 27*c + 64*d = 0                          (3)
%...
%...     to drop xzzzz
%...
%...        a + 16*b + 81*c + 256*d = 0                        (4)
%...
%...     to drop xzzzzz
%...
%...        a + 32*b + 243*c + 1024*d = 0                      (5)
%...
%...  equations (1) to (5) can be solved for a, b, c, d and e.  if
%...
%...  equation (2) is subtracted from equations (3), (4) and (5),
%...
%...     4*b + 18*c + 48*d = -2                                (6)
%...
%...     12*b + 72*c + 240*d = -2                              (7)
%...
%...     28*b + 234*c + 1008*d = -2                            (8)
%...
%...  equations (6), (7) and (8) can be solved for b, c and d
%...
%...     18*c + 96*d = 4                                       (9)
%...
%...     108*c + 672*d = 12                                    (10)
%...
%...  equations (9) and (10) can be solved for c and d, c = 8/9, d =
%...  -1/8.
%...
%...  from equation (6), b = -3.  from equation (2), a = 8.  from
%...  equation (1), e = -25/6.
%...
%...  the final differentiation formula is then obtained as
%...
%...  8*x(i+1) - 3*x(i+2) + (8/9)*x(i+3) - (1/8)*x(i+4)
%...
%...  - (25/6)*xz(i)*dz
%...
%...  = (8 - 3 + (8/9) - (1/8))*x(i) + xzz*(dz**2) + o(dz**6)
%...
%...  or
%...
%...  xzz(i) = (1/12*dz**2)*((-415/6)*x(i) + 96*x(i+1) - 36*x(i+2)
%...                                                                (35)
%...  + (32/3)*x(i+3) - (3/2)*x(i+4) - 50*xz(i)*dz) + o(dz**4)
%...
%...  equation (35) will be applied at i = 1 and i = n
%...
%...  xzz(1) = (1/12*dz**2)*((-415/6)*x(1) + 96*x(2) - 36*x(3)
%...                                                                (36)
%...  + (32/3)*x(4) - (3/2)*x(5) - 50*xz(1)*dz) + o(dz**4)
%...
%...  xzz(n) = (1/12*dz**2)*((-415/6)*x(n) + 96*x(n-1) - 36*x(n-2)
%...                                                                (37)
%...  + (32/3)*x(n-3) - (3/2)*x(n-4) + 50*xz(n)*dz) + o(dz**4)
%...
%...  ******************************************************************
%...
%...  grid spacing
      n=length(x);
      dz=(z0-zL)/(n-1);
%...
%...  calculate xzz at the left boundary, including xz
      if n0==1
      xzz_z0 = (-415/6*x(1)+96*x(2)-36*x(3)+32/3*x(4)-3/2*x(5))/(12*dz^2)-50*xz_z0/(12*dz);
      end
%...
%...  calculate xzz at the right boundary, including xz
      if nL==1
      xzz_zL =(-415/6*x(n)+96*x(n-1)-36*x(n-2)+32/3*x(n-3)-3/2*x(n-4))/(12*dz^2)+50*xz_zL/(12*dz);
      end