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