eleven_point_centered_uni_d1.m

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

      function [D]=eleven_point_centered_uni_D1(z0,zL,n)
%...
%...  A. Vande Wouwer, P. Saucez and W.E. Schiesser (2002)
%...
%...  function eleven_point_centered_uni_D1 returns the differentiation matrix 
%...  for computing the first derivative, xz, of a variable x over the spatial domain
%...  z0 < z < zL from classical eleven-point, tenth-order finite difference 
%...  approximations (this function replaces dss0010)
%...
%...
%...  argument list
%...
%...     z0      lower boundary value of z (input)
%...
%...     zL      upper boundary value of z (input)
%...
%...     n       number of grid points in the z domain including the
%...             boundary points (input)
%...
%...  origin of the approximation
%...
%...  the mathematical details of the finite difference approxima-
%...  tions can be summarized by the following bickley matrix
%...  for n = 10, m = 1, p = 0 to 10 (bickley, w. g., formulae for
%...  numerical differentiation, math. gaz., vol. 25, 1941)
%...
%...      -10626840            36288000           -81648000
%...      145152000          -190512000           182891520
%...
%...        -362880            -6636960            16329600
%...      -21772800            25401600           -22861440
%...
%...          40320             -806400            -4419360
%...        9676800            -8467200             6773760
%...
%...         -10080              151200            -1360800
%...       -2756160             6350400            -3810240
%...
%...           4320              -57600              388800
%...       -2073600            -1330560             4354560
%...
%...          -2880               36000             -216000
%...         864000            -3024000                   0
%...
%...           2880              -34560              194400
%...        -691200             1814400            -4354560
%...
%...          -4320               50400             -272160
%...         907200            -2116800             3810240
%...
%...          10080             -115200              604800
%...       -1935360             4233600            -6773760
%...
%...         -40320              453600            -2332800
%...        7257600           -15240960            22861440
%...
%...         362880            -4032000            20412000
%...      -62208000           127008000          -182891520
%...
%...  each entry in this table should be multiplied by 1/10f to
%...  obtain the final weighting coefficients.  from this bickley
%...  matrix, the finite difference approximation of the first
%...  derivative can be programmed for each of the grid points 1, 2,
%...  3, 4, 5,..., i,..., n-4, n-3, n-2, n-1 and n (taking into
%...  account the symmetry properties of the matrix).
%...
%...
%...  compute the spatial increment
      dz=(zL-z0)/(n-1);
      r10fdz=1./(3628800.*dz);
%...
%...  sparse discretization matrix      
%...
%...  interior points      
      D=diag(-2880*ones(n-5,1),-5)+diag(36000*ones(n-4,1),-4)...
       +diag(-216000*ones(n-3,1),-3)+diag(+864000*ones(n-2,1),-2)...
       +diag(-3024000*ones(n-1,1),-1)+diag(+3024000*ones(n-1,1),+1)...
       +diag(-864000*ones(n-2,1),+2)+diag(+216000*ones(n-3,1),+3)...
       +diag(-36000*ones(n-4,1),+4)+diag(+2880*ones(n-5,1),+5);      
%...
%...  boundary points      
      D([1 2 3 4 5],1:11) = [-10628640 [phone] -81648000 [phone] -190512000 [phone] -127008000 [phone] -20412000 [phone];
                             -362880 -6636960 [phone] -21772800 [phone] -22861440 [phone] [phone] +40320;
                             [phone] -4419360 [phone] [phone] [phone] [phone];
                             -10080 [phone] -2756160 [phone] [phone] [phone] +4320;
                             [phone] [phone] -1330560 [phone] [phone] [phone]];
%...                         
      D([n-4 n-3 n-2 n-1 n],(n-10):n) = [[phone] [phone] [phone] +1330560 [phone] [phone];
                                         -4320 [phone] [phone] [phone] +2756160 [phone] +10080;
                                         [phone] [phone] [phone] [phone] +4419360 [phone];
                                         -40320 [phone] [phone] [phone] -25401600 [phone] -16329600 +6636960 +362880; 
                                         [phone] [phone] -62208000 [phone] -182891520 [phone] -145152000 [phone] -36288000 [phone]];
%...      
      D=r10fdz*D;
      D=sparse(D);