Method of lines/example implementation/dss010

From Scholarpedia
Jump to: navigation, search
    %  File: dss010.m
    %
       function [ux]=dss010(xl,xu,n,u)
    %
    %  Function dss010 computes the first derivative, u , of a
    %                                                 x
    %  variable u over the spatial domain xl le x le xu from classical
    %  eleven-point, tenth-order finite difference approximations
    %
    %  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)
    %
    %  The mathematical details of the finite difference approxima-
    %  tions are given in routines dss002, dss004, dss006 and
    %  dss008, and 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/10 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
       dx=(xu-xl)/(n-1);
       r10fdx=1./(3628800.*dx);
       nm5=n-5;
    %
    %  Grid point 1
       ux(  1)=r10fdx*...
         (-10628640.      *u(   1)...
          +36288000.      *u(   2)...
          -81648000.      *u(   3)...
          +145152000.     *u(   4)...
          -190512000.     *u(   5)...
          +182891520.     *u(   6)...
          -127008000.     *u(   7)...
          +62208000.      *u(   8)...
          -20412000.      *u(   9)...
          +4032000.       *u(  10)...
          -362880.        *u(  11));
    %
    %  Grid point 2
       ux(  2)=r10fdx*...
         (-362880.        *u(   1)...
          -6636960.       *u(   2)...
          +16329600.      *u(   3)...
          -21772800.      *u(   4)...
          +25401600.      *u(   5)...
          -22861440.      *u(   6)...
          +15240960.      *u(   7)...
          -7257600.       *u(   8)...
          +2332800.       *u(   9)...
          -453600.        *u(  10)...
          +40320.         *u(  11));
    %
    %  Grid point 3
       ux(  3)=r10fdx*...
         (+40320.         *u(   1)...
          -806400.        *u(   2)...
          -4419360.       *u(   3)...
          +9676800.       *u(   4)...
          -8467200.       *u(   5)...
          +6773760.       *u(   6)...
          -4233600.       *u(   7)...
          +1935360.       *u(   8)...
          -604800.        *u(   9)...
          +115200.        *u(  10)...
          -10080.         *u(  11));
    %
    %  Grid point 4
       ux(  4)=r10fdx*...
         (-10080.         *u(   1)...
          +151200.        *u(   2)...
          -1360800.       *u(   3)...
          -2756160.       *u(   4)...
          +6350400.       *u(   5)...
          -3810240.       *u(   6)...
          +2116800.       *u(   7)...
          -907200.        *u(   8)...
          +272160.        *u(   9)...
          -50400.         *u(  10)...
          +4320.          *u(  11));
    %
    %  Grid point 5
       ux(  5)=r10fdx*...
         (+4320.          *u(   1)...
          -57600.         *u(   2)...
          +388800.        *u(   3)...
          -2073600.       *u(   4)...
          -1330560.       *u(   5)...
          +4354560.       *u(   6)...
          -1814400.       *u(   7)...
          +691200.        *u(   8)...
          -194400.        *u(   9)...
          +34560.         *u(  10)...
          -2880.          *u(  11));
    %
    %  Grid point i, i ne 1 to 5, n-4 to n
       for i=6:nm5
         ux(  i)=r10fdx*...
         (-2880.          *u( i-5)...
          +36000.         *u( i-4)...
          -216000.        *u( i-3)...
          +864000.        *u( i-2)...
          -3024000.       *u( i-1)...
          +0.             *u( i  )...
          +3024000.       *u( i+1)...
          -864000.        *u( i+2)...
          +216000.        *u( i+3)...
          -36000.         *u( i+4)...
          +2880.          *u( i+5));
       end
    %
    %  Grid point n-4
       ux(n-4)=r10fdx*...
         (+2880.          *u(n-10)...
          -34560.         *u( n-9)...
          +194400.        *u( n-8)...
          -691200.        *u( n-7)...
          +1814400.       *u( n-6)...
          -4354560.       *u( n-5)...
          +1330560.       *u( n-4)...
          +2073600.       *u( n-3)...
          -388800.        *u( n-2)...
          +57600.         *u( n-1)...
          -4320.          *u( n  ));
    %
    %  Grid point n-3
       ux(n-3)=r10fdx*...
         (-4320.          *u(n-10)...
          +50400.         *u( n-9)...
          -272160.        *u( n-8)...
          +907200.        *u( n-7)...
          -2116800.       *u( n-6)...
          +3810240.       *u( n-5)...
          -6350400.       *u( n-4)...
          +2756160.       *u( n-3)...
          +1360800.       *u( n-2)...
          -151200.        *u( n-1)...
          +10080.         *u( n  ));
    %
    %  Grid point n-2
       ux(n-2)=r10fdx*...
         (+10080.         *u(n-10)...
          -115200.        *u( n-9)...
          +604800.        *u( n-8)...
          -1935360.       *u( n-7)...
          +4233600.       *u( n-6)...
          -6773760.       *u( n-5)...
          +8467200.       *u( n-4)...
          -9676800.       *u( n-3)...
          +4419360.       *u( n-2)...
          +806400.        *u( n-1)...
          -40320.         *u( n  ));
    %
    %  Grid point n-1
       ux(n-1)=r10fdx*...
         (-40320.         *u(n-10)...
          +453600.        *u( n-9)...
          -2332800.       *u( n-8)...
          +7257600.       *u( n-7)...
          -15240960.      *u( n-6)...
          +22861440.      *u( n-5)...
          -25401600.      *u( n-4)...
          +21772800.      *u( n-3)...
          -16329600.      *u( n-2)...
          +6636960.       *u( n-1)...
          +362880.        *u( n  ));
    %
    %  Grid point n
       ux(n  )=r10fdx*...
         (+362880.        *u(n-10)...
          -4032000.       *u( n-9)...
          +20412000.      *u( n-8)...
          -62208000.      *u( n-7)...
          +127008000.     *u( n-6)...
          -182891520.     *u( n-5)...
          +190512000.     *u( n-4)...
          -145152000.     *u( n-3)...
          +81648000.      *u( n-2)...
          -36288000.      *u( n-1)...
          +10628640.      *u( n  ));
    %
    
    Personal tools
    Namespaces

    Variants
    Actions
    Navigation
    Focal areas
    Activity
    Tools