Method of lines/example implementation/dss048

From Scholarpedia
Jump to: navigation, search
    % File: dss048.m
    %
      function [uxx]=dss048(xl,xu,n,u,ux,nl,nu)
    %
    %  Function dss048 computes a eighth-order approximation of a
    %  second-order derivative, with or without the normal derivative
    %  at the boundary. 
    %
    %  Argument list
    %
    %     xl      Left value of the spatial independent variable (input)
    %
    %     xu      Right value of the spatial independent variable (input)
    %
    %     n       Number of spatial grid points, including the end
    %             points (input)
    %
    %     u       One-dimensional array of the dependent variable to be
    %             differentiated (input)
    %
    %     ux      One-dimensional array of the first derivative of u.
    %             The end values of ux, ux(1) and ux(n), are used in
    %             Neumann boundary conditions at x = xl and x = xu,
    %             depending on the arguments nl and nu (see the de-
    %             scription of nl and nu below)
    %
    %     uxx     One-dimensional array of the second derivative of u
    %             (output)
    %
    %     nl      Integer index for the type of boundary condition at
    %             x = xl (input).  the allowable values are
    %
    %                1 - Dirichlet boundary condition at x = xl
    %                    (ux(1) is not used)
    %
    %                2 - Neumann boundary condition at x = xl
    %                    (ux(1) is used)
    %
    %     nu      Integer index for the type of boundary condition at
    %             x = xu (input).  The allowable values are
    %
    %                1 - Dirichlet boundary condition at x = xu
    %                    (ux(n) is not used)
    %
    %                2 - Neumann boundary condition at x = xu
    %                    (ux(n) is used)
    %
    %  The following eighth-order, ten-point approximations for a
    %  second derivative 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.
    %
    %  Weighting coefficients for finite difference approximations, computed
    %  by an algorithm of B. Fornberg (1,2), are used in the following 
    %  approximations.  The combination of these coefficients to give the 
    %  second derivative approximations at the boundaries in x, including 
    %  the first derivative, was suggested by Professor Fornberg.
    %
    %  (1)  Fornberg, B., Fast Generation of Weights in Finite Difference
    %       Formulas, In Recent Developments In Numerical Methods And
    %       Software for ODEs/DAEs/PDEs, G. Byrne and W. E. Schiesser (eds), 
    %       World Scientific, River Edge, NJ, 1992
    %
    %  (2)  Fornberg, B., Calculation of Weights In Finite Difference
    %       Formulas, SIAM Review, Vol. 40, No. 3, Pp 685-691, September, 
    %       1999
    %
    %  Grid spacing
       dx=(xu-xl)/(n-1);
       rdxs=1.0/dx^2;
    %
    %  uxx at the left boundary
    %
    %     Without ux 
             if nl==1
             uxx(1)=rdxs*...
                    ( 6.463293650793652*u(1)...
                    -32.921428571428578*u(2)...
                    +83.842857142857156*u(3)...
                   -139.755555555555555*u(4)...
                   +162.474999999999999*u(5)...
                   -132.500000000000000*u(6)...
                    +74.544444444444444*u(7)...
                    -27.628571428571430*u(8)...
                     +6.083928571428572*u(9)...
                     -0.603968253968254*u(10));
    %
    %     With ux
             elseif nl==2
             uxx(1)=rdxs*...
                   ( -8.914169501133784*u(1)...
                    +16.000000000000000*u(2)...
                    -14.000000000000000*u(3)...
                    +12.444444444444444*u(4)...
                     -8.750000000000000*u(5)...
                     +4.480000000000000*u(6)...
                     -1.555555555555555*u(7)...
                     +0.326530612244898*u(8)...
                     -0.031250000000000*u(9)... 
                     -5.435714285714286*ux(1)*dx);
             end
    %
    %  uxx at the right boundary
    %
    %     Without ux 
             if nu==1
             uxx(n)=rdxs*...
                   ( -0.603968253968254*u(n-9)...
                     +6.083928571428572*u(n-8)...
                    -27.628571428571430*u(n-7)...
                    +74.544444444444444*u(n-6)...
                   -132.500000000000000*u(n-5)...
                   +162.474999999999999*u(n-4)...
                   -139.755555555555555*u(n-3)...
                    +83.842857142857156*u(n-2)...
                    -32.921428571428578*u(n-1)...
                     +6.463293650793652*u(n  ));
    %
    %     With ux 
             elseif nu==2
             uxx(n)=rdxs*...
                    (-0.031250000000000*u(n-8)...
                     +0.326530612244905*u(n-7)...
                     -1.555555555555555*u(n-6)...
                     +4.480000000000000*u(n-5)...
                     -8.750000000000000*u(n-4)...
                    +12.444444444444444*u(n-3)...
                    -14.000000000000000*u(n-2)...
                    +16.000000000000000*u(n-1)...
                     -8.914169501133790*u(n  )...
                     +5.435714285714287*ux(n  )*dx);
             end
    %
    %  uxx at the interior grid points
    %
    %     i = 2 
             uxx(2)=rdxs*...
                     (0.603968253968254*u(1)...
                     +0.423611111111111*u(2)...
                     -5.742857142857143*u(3)...
                    +11.366666666666665*u(4)...
                    -12.922222222222221*u(5)...
                    +10.274999999999999*u(6)...
                     -5.666666666666665*u(7)...
                     +2.068253968253968*u(8)...
                     -0.450000000000000*u(9)...
                     +0.044246031746032*u(10));
    %
    %     i = 3
             uxx(3)=rdxs*...
                    (-0.044246031746032*u(1)...
                     +1.046428571428571*u(2)...
                     -1.567460317460318*u(3)...
                     -0.433333333333333*u(4)...
                     +2.075000000000000*u(5)...
                     -1.772222222222222*u(6)...
                     +0.983333333333333*u(7)...
                     -0.357142857142857*u(8)...
                     +0.077182539682540*u(9)...
                     -0.007539682539683*u(10));
    %
    %     i = 4
             uxx(4)=rdxs*...
                     (0.007539682539683*u(1)...
                     -0.119642857142857*u(2)...
                     +1.385714285714286*u(3)...
                     -2.472222222222222*u(4)...
                     +1.150000000000000*u(5)...
                     +0.175000000000000*u(6)...
                     -0.188888888888889*u(7)...
                     +0.078571428571429*u(8)...
                     -0.017857142857143*u(9)...
                     +0.001785714285714*u(10));
    %
    %     i = n-3
             uxx(n-3)=rdxs*...
                     (0.001785714285714*u(n-9)...
                     -0.017857142857143*u(n-8)...
                     +0.078571428571429*u(n-7)...
                     -0.188888888888889*u(n-6)...
                     +0.175000000000000*u(n-5)...
                     +1.150000000000000*u(n-4)...
                     -2.472222222222224*u(n-3)...
                     +1.385714285714286*u(n-2)...
                     -0.119642857142857*u(n-1)...
                     +0.007539682539683*u(n  ));
    %
    %     i = n-2 
             uxx(n-2)=rdxs*...
                    (-0.007539682539683*u(n-9)...
                     +0.077182539682540*u(n-8)...
                     -0.357142857142857*u(n-7)...
                     +0.983333333333333*u(n-6)...
                     -1.772222222222222*u(n-5)...
                     +2.075000000000001*u(n-4)...
                     -0.433333333333335*u(n-3)...
                     -1.567460317460318*u(n-2)...
                     +1.046428571428572*u(n-1)...
                     -0.044246031746032*u(n  ));
    %
    %     i = n-1 
             uxx(n-1)=rdxs*...
                     (0.044246031746032*u(n-9)...
                     -0.450000000000000*u(n-8)...
                     +2.068253968253968*u(n-7)...
                     -5.666666666666667*u(n-6)...
                    +10.275000000000000*u(n-5)...
                    -12.922222222222222*u(n-4)...
                    +11.366666666666667*u(n-3)...
                     -5.742857142857142*u(n-2)...
                     +0.423611111111112*u(n-1)...
                     +0.603968253968254*u(n  ));
    %
    %     i = 5, 6,..., n-4
             for i=5:n-4
             uxx(i)=rdxs*...
                    (-0.001785714285714*u(i-4)...
                     +0.025396825396825*u(i-3)...
                     -0.200000000000000*u(i-2)...
                     +1.600000000000000*u(i-1)...
                     -2.847222222222221*u(i  )...
                     +1.600000000000000*u(i+1)...
                     -0.200000000000000*u(i+2)...
                     +0.025396825396825*u(i+3)...
                     -0.001785714285714*u(i+4));
             end
    %
    
    Personal tools
    Namespaces

    Variants
    Actions
    Navigation
    Focal areas
    Activity
    Tools