# Method of lines/example implementation/dss004

% File: dss004.m % function [ux]=dss004(xl,xu,n,u) % % Function dss004 computes the first derivative, u , of a % x % variable u over the spatial domain xl le x le xu from classical % five-point, fourth-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 routine dss002. % % Five-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 % % Constants a, b, c and d are selected so that the coefficients % of the u1 terms sum to one and the coefficients of the u1 , % x 2x % u1 and u1 terms sum to zero % 3x 4x % % a + 2b + 3c + 4d = 1 % % a + 4b + 9c + 16d = 0 % % a + 8b + 27c + 64d = 0 % % a + 16b + 81c + 256d = 0 % % Simultaneous solution for a, b, c and d followed by the solu- % tion of the preceding Taylor series, truncated after the u % 4x % terms, for u1 gives the following five-point approximation % x % 4 % u1 = (1/12dx)(-25u1 + 48u2 - 36u3 + 16u4 - 3u5) + O(dx ) (1) % x % % (2) Interior point, i = 2 % % 2 3 4 % a(u1 = u2 + u2 (-dx) + u2 (-dx) + u2 (-dx) + u2 (-dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u2 (-dx) + u2 (-dx) + u2 (-dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % b(u3 = u2 + u2 ( dx) + u2 ( dx) + u2 ( dx) + u2 ( dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u2 ( dx) + u2 ( dx) + u2 ( dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % c(u4 = u2 + u2 (2dx) + u2 (2dx) + u2 (2dx) + u2 (2dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u2 (2dx) + u2 (2dx) + u2 (2dx) + ...) % 5x 5f 6x 6f 7x 7f % % 2 3 4 % d(u5 = u2 + u2 (3dx) + u2 (3dx) + u2 (3dx) + u2 (3dx) % x 1f 2x 2f 3x 3f 4x 4f % % 5 6 7 % + u2 (3dx) + u2 (3dx) + u2 (3dx) + ...) % 5x 5f 6x 6f 7x 7f % % -a + b + 2c + 3d = 1 % % a + b + 4c + 9d = 0 % % -a + b + 8c + 27d = 0 % % a + b + 16c + 81d = 0 % % Simultaneous solution for a, b, c and d followed by the solu- % tion of the preceding Taylor series, truncated after the u % 4x % terms, for u1 gives the following five-point approximation % x % 4 % u2 = (1/12dx)(-3u1 - 10u2 + 18u3 - 6u4 + u5) + O(dx ) (2) % x % % (3) Interior point i, i ne 2, n-1 % % 2 3 % a(ui-2 = ui + ui (-2dx) + ui (-2dx) + ui (-2dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + ui (-2dx) + ui (-2dx) + ui (-2dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % b(ui-1 = ui + ui ( -dx) + ui ( -dx) + ui ( -dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + ui ( -dx) + ui ( -dx) + ui ( -dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % c(ui+1 = ui + ui ( dx) + ui ( dx) + ui ( dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + ui ( dx) + ui ( dx) + ui ( dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % d(ui+2 = ui + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) + ...) % 4x 4f 5x 5f 6x 6f % % -2a - b + c + 2d = 1 % % 4a + b + c + 4d = 0 % % -8a - b + c + 8d = 0 % % 16a + b + c + 16d = 0 % % Simultaneous solution for a, b, c and d followed by the solu- % tion of the preceding Taylor series, truncated after the u % 4x % terms, for u1 gives the following five-point approximation % x % 4 % ui = (1/12dx)(ui-2 - 8ui-1 + 0ui + 8ui+1 - ui+2) + O(dx ) (3) % x % % (4) Interior point, i = n-1 % % 2 3 % a(un-4 = un-1 + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) + ... % 4x 4f 5x 5f 6x 6f % % 2 3 % b(un-3 = un-1 + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) + ... % 4x 4f 5x 5f 6x 6f % % 2 3 % c(un-2 = un-1 + un-1 ( -dx) + un-1 (- -x) + un-1 ( -dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un-1 ( -dx) + un-1 ( -dx) + un-1 ( -dx) + ... % 4x 4f 5x 5f 6x 6f % % 2 3 % d(un = un-1 + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) + ... % 4x 4f 5x 5f 6x 6f % % -3a - 2b - c + d = 1 % % 9a + 4b + c + d = 0 % % -27a - 8b - c + d = 0 % % 81a + 16b + c + d = 0 % % Simultaneous solution for a, b, c and d followed by the solu- % tion of the preceding Taylor series, truncated after the u % 4x % terms, for u1 gives the following five-point approximation % x % 4 % un-1 = (1/12dx)(-un-4 + 6un-3 - 18un-2 + 10un-1 + 3un) + O(dx ) % x % (4) % % (5) Right end, point i = n % % 2 3 % a(un-4 = un + un (-4dx) + un (-4dx) + un (-4dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un (-4dx) + un (-4dx) + un (-4dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % b(un-3 = un + un (-3dx) + un (-3dx) + un (-3dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un (-3dx) + un (-3dx) + un (-3dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % c(un-2 = un + un (-2dx) + un (-2dx) + un (-2dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un (-2dx) + un (-2dx) + un (-2dx) + ...) % 4x 4f 5x 5f 6x 6f % % 2 3 % d(un-1 = un + un ( -dx) + un ( -dx) + un ( -dx) % x 1f 2x 2f 3x 3f % % 4 5 6 % + un ( -dx) + un ( -dx) + un ( -dx) + ...) % 4x 4f 5x 5f 6x 6f % % -4a - 3b - 2c - d = 1 % % 16a + 9b + 4c + d = 0 % % -64a - 27b - 8c - d = 0 % % 256a + 81b + 16c + d = 0 % % Simultaneous solution for a, b, c and d followed by the solu- % tion of the preceding Taylor series, truncated after the u % 4x % terms, for u1 gives the following five-point approximation % x % 4 % un = (1/12dx)(3un-4 - 16un-3 + 36un-2 - 48un-1 + 25un) + O(dx ) % x % (5) % % The weighting coefficients for equations (1) to (5) can be % summarized as % % -25 48 -36 16 -3 % % -3 -10 18 -6 1 % % 1/12 1 -8 0 8 -1 % % -1 6 -18 10 3 % % 3 -16 36 -48 25 % % which are the coefficients reported by Bickley for n = 4, m = % 1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical % Differentiation, Math. Gaz., vol. 25, 1941. Note - the Bickley % coefficients have been divided by a common factor of two). % % Equations (1) to (5) can now be programmed to generate the % derivative u (x) of function u(x). % x % % Compute the spatial increment dx=(xu-xl)/(n-1); r4fdx=1./(12.*dx); nm2=n-2; % % Equation (1) (note - the rhs of equations (1), (2), (3), (4) % and (5) have been formatted so that the numerical weighting % coefficients can be more easily associated with the Bickley % matrix above) ux( 1)=r4fdx*... ( -25.*u( 1) +48.*u( 2) -36.*u( 3) +16.*u( 4) -3.*u( 5)); % % Equation (2) ux( 2)=r4fdx*... ( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5)); % % Equation (3) for i=3:nm2 ux( i)=r4fdx*... ( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2)); end % % Equation (4) ux(n-1)=r4fdx*... ( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n)); % % Equation (5) ux( n)=r4fdx*... ( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n)); %