% File: dss008.m % function [ux]=dss008(xl,xu,n,u) % % Routine dss008 computes the first derivative, u , of a % x % variable u over the spatial domain xl le x le xu from classical % nine-point, eighth-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 following Taylor series (or % polynomials) are given in routines dss002, dss004 and % dss006. % % Nine-point formulas % % (1) Left end, point i = 1 % % 2 3 4 % a(u2 = u1 + u1 ( dx) + u1 ( dx) + u1 ( dx) + u1 ( dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 ( dx) + u1 ( dx) + u1 ( dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % b(u3 = u1 + u1 (2dx) + u1 (2dx) + u1 (2dx) + u1 (2dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (2dx) + u1 (2dx) + u1 (2dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % c(u4 = u1 + u1 (3dx) + u1 (3dx) + u1 (3dx) + u1 (3dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (3dx) + u1 (3dx) + u1 (3dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % d(u5 = u1 + u1 (4dx) + u1 (4dx) + u1 (4dx) + u1 (4dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (4dx) + u1 (4dx) + u1 (4dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % e(u6 = u1 + u1 (5dx) + u1 (5dx) + u1 (5dx) + u1 (5dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (5dx) + u1 (5dx) + u1 (5dx) + ...) % 5x 5f 6x 6f 7x 7f % 2 3 4 % f(u7 = u1 + u1 (6dx) + u1 (6dx) + u1 (6dx) + u1 (6dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (6dx) + u1 (6dx) + u1 (6dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % g(u8 = u1 + u1 (7dx) + u1 (7dx) + u1 (7dx) + u1 (7dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (7dx) + u1 (7dx) + u1 (7dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % h(u9 = u1 + u1 (8dx) + u1 (8dx) + u1 (8dx) + u1 (8dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u1 (8dx) + u1 (8dx) + u1 (8dx) + ...) % 5x 5f 6x 6f 7x 7f % % Constants a, b, c, d, e, f, g and h are selected so that the % coefficients of the u1 terms sum to one and the coefficients of % x % the u1 , u1 , u1 , u1 , u1 , u1 and u1 sum to zero. % 2x 3x 4x 5x 6x, 7x 8x % % 1 1 1 1 1 1 1 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 1 % % 2 2 2 2 2 2 2 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 3 3 3 3 3 3 3 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 4 4 4 4 4 4 4 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 5 5 5 5 5 5 5 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 6 6 6 6 6 6 6 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 7 7 7 7 7 7 7 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % 8 8 8 8 8 8 8 % a + 2 b + 3 c + 4 d + 5 e + 6 f + 7 g + 8 h = 0 % % Simultaneous solution for a, b, c, d, e, f, g and h followed by % the solution of the preceding Taylor series, truncated after the % u terms, for u1 gives the following nine-point approximation % 8x x % % u1 = (1/8f)(-109584u1 + 322560u2 - 564480u3 + 752640u4 % x % -705600u5 + 451584u6 - 188160u7 + 46080u8 % 8 % - 5040u9) + O(dx ) (1) % % The preceding analysis can be repeated to produce nine-point % approximations for the first derivatives u2 , u3 , u4 , ui , % x x x x % un-3 , un-2 , un-1 and un . The results can be summarized by % x x x x % the following Bickley matrix for n = 8, m = 1, p = 0 to 8, % (Bickley, W. G., Formulae for Numerical Differentiation, Math. % Gaz., vol. 25, 1941) % % -109584 322560 -564480 752640 -705600 % % -5040 -64224 141120 -141120 117600 % % 720 -11520 -38304 80640 -50400 % % -240 2880 -20160 -18144 50400 % % 1/8f 144 -1536 8064 -32256 0 % % -144 1440 -6720 20160 -50400 % % 240 -2304 10080 -26880 50400 % % -720 6720 -28224 70560 -117600 % % 5040 -46080 188160 -451584 705600 % % 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,..., i,..., n-3, n-2, n-1, n (taking into account % the symmetry properties of the Bickley matrix). % % Compute the spatial increment dx=(xu-xl)/(n-1); r8fdx=1./(40320.*dx); nm4=n-4; % % Grid point 1 ux( 1)=r8fdx*... (-109584. *u( 1)... +322560. *u( 2)... -564480. *u( 3)... +752640. *u( 4)... -705600. *u( 5)... +451584. *u( 6)... -188160. *u( 7)... +46080. *u( 8)... -5040. *u( 9)); % % Grid point 2 ux( 2)=r8fdx*... (-5040. *u( 1)... -64224. *u( 2)... +141120. *u( 3)... -141120. *u( 4)... +117600. *u( 5)... -70560. *u( 6)... +28224. *u( 7)... -6720. *u( 8)... +720. *u( 9)); % % Grid point 3 ux( 3)=r8fdx*... (+720. *u( 1)... -11520. *u( 2)... -38304. *u( 3)... +80640. *u( 4)... -50400. *u( 5)... +26880. *u( 6)... -10080. *u( 7)... +2304. *u( 8)... -240. *u( 9)); % % Grid point 4 ux( 4)=r8fdx*... (-240. *u( 1)... +2880. *u( 2)... -20160. *u( 3)... -18144. *u( 4)... +50400. *u( 5)... -20160. *u( 6)... +6720. *u( 7)... -1440. *u( 8)... +144. *u( 9)); % % Grid point i for i=5:nm4 ux( i)=r8fdx*... (+144. *u(i-4)... -1536. *u(i-3)... +8064. *u(i-2)... -32256. *u(i-1)... +0. *u(i )... +32256. *u(i+1)... -8064. *u(i+2)... +1536. *u(i+3)... -144. *u(i+4)); end % % Grid point n-3 ux(n-3)=r8fdx*... (-144. *u(n-8)... +1440. *u(n-7)... -6720. *u(n-6)... +20160. *u(n-5)... -50400. *u(n-4)... +18144. *u(n-3)... +20160. *u(n-2)... -2880. *u(n-1)... +240. *u(n )); % % Grid point n-2 ux(n-2)=r8fdx*... (+240. *u(n-8)... -2304. *u(n-7)... +10080. *u(n-6)... -26880. *u(n-5)... +50400. *u(n-4)... -80640. *u(n-3)... +38304. *u(n-2)... +11520. *u(n-1)... -720. *u(n )); % % Grid point n-1 ux(n-1)=r8fdx*... (-720. *u(n-8)... +6720. *u(n-7)... -28224. *u(n-6)... +70560. *u(n-5)... -117600. *u(n-4)... +141120. *u(n-3)... -141120. *u(n-2)... +64224. *u(n-1)... +5040. *u(n )); % % Grid point n ux(n )=r8fdx*... (+5040. *u(n-8)... -46080. *u(n-7)... +188160. *u(n-6)... -451584. *u(n-5)... +705600. *u(n-4)... -752640. *u(n-3)... +564480. *u(n-2)... -322560. *u(n-1)... +109584. *u(n ));