dss010
From Scholarpedia
< Method of lines | example implementation
| Samir Hamdi et al. (2007), Scholarpedia, 2(10):2859. | doi:10.4249/scholarpedia.5398 | revision #25160 [link to/cite this article] | |||||||||||||||||||
Curator: Samir Hamdi, GENIVAR (Hydropower & Hydraulics) and Solid Mechanics Laboratory, Ecole Polytechnique, Paris, France
Curator: William E. Schiesser, Lehigh University and University of Pennsylvania, USA
Curator: Graham W Griffiths, School of Engineering and Mathematical Sciences, City University, London, UK
% 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 ));
%
