Method of lines/example implementation/dss008
From Scholarpedia
% 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 ));
%
