Method of lines/example implementation/dss044
From Scholarpedia
% File: ds044.m
%
function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%
% Function dss044 computes a fourth-order approximation of a
% second-order derivative, with or without the normal derivative
% at the boundary.
%
% Argument list
%
% xl Left value of the spatial independent variable (input)
%
% xu Right value of the spatial independent variable (input)
%
% n Number of spatial grid points, including the end
% points (input)
%
% u One-dimensional array of the dependent variable to be
% differentiated (input)
%
% ux One-dimensional array of the first derivative of u.
% The end values of ux, ux(1) and ux(n), are used in
% Neumann boundary conditions at x = xl and x = xu,
% depending on the arguments nl and nu (see the de-
% scription of nl and nu below)
%
% uxx One-dimensional array of the second derivative of u
% (output)
%
% nl Integer index for the type of boundary condition at
% x = xl (input). The allowable values are
%
% 1 - Dirichlet boundary condition at x = xl
% (ux(1) is not used)
%
% 2 - Neumann boundary condition at x = xl
% (ux(1) is used)
%
% nu Integer index for the type of boundary condition at
% x = xu (input). The allowable values are
%
% 1 - Dirichlet boundary condition at x = xu
% (ux(n) is not used)
%
% 2 - Neumann boundary condition at x = xu
% (ux(n) is used)
%
% The following derivation was completed by W. E. Schiesser, Depts
% of CHE and Math, Lehigh University, Bethlehem, PA 18015, USA, on
% December 15, 1986. Additional details are given in function
% dss042.
%
% ******************************************************************
%
% (1) uxx at the interior points 3, 4,..., n-2
%
% To develop a set of fourth-order correct differentiation formulas
% for the second derivative uxx, we consider first the interior
% grid points at which a symmetric formula can be used.
%
% If we consider a formula of the form
%
% a*u(i-2) + b*u(i-1) + e*u(i) + c*u(i+1) + d*u(i+2)
%
% Taylor series expansions of u(i-2), u(i-1), u(i+1) and u(i+2)
% can be substituted into this formula. We then consider the
% linear albegraic equations relating a, b, c and d which will
% retain certain terms, i.e., uxx, and drop others, e.g., uxxx,
% uxxxx and uxxxxx.
%
% Thus, for grid points 3, 4,..., n-2
%
% To retain uxx
%
% 4*a + b + c + 4*d = 2 (1)
%
% To drop uxxx
%
% -8*a - b + c + 8*d = 0 (2)
%
% To drop uxxxx
%
% 16*a + b + c + 16*d = 0 (3)
%
% To drop uxxxxx
%
% -32*a - b + c + 32*d = 0 (4)
%
% Equations (1) to (4) can be solved for a, b, c and d. If equa-
% tion (1) is added to equation (2)
%
% -4*a + 2*c + 12*d = 2 (5)
%
% If equation (1) is subtracted from equation (3)
%
% 12*a + 12*d = -2 (6)
%
% If equation (1) is added to equation (4)
%
% -28*a + 2*c + 36*d = 2 (7)
%
% Equations (5) to (7) can be solved for a, c and d. If equation
% (5) is subtracted from equation (7), and the result combined
% with equation (6)
%
% 12*a + 12*d = -2 (6)
%
% -24*a + 24*d = 0 (8)
%
% Equations (6) and (8) can be solved for a and d. From (8), a
% = d. From equation (6), a = -1/12 and d = -1/12. Then, from
% equation (5), c = 4/3, and from equation (1), b = 4/3.
%
% The final differentiation formula is then obtained as
%
% (-1/12)*u(i-2) + (4/3)*u(i-1) +
%
% (4/3)*u(i+1) + (-1/12)*u(i+2)
%
% (-1/12 + 4/3 - 1/12 + 4/3)*u(i) + uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/(12*dx**2))*(-1*u(i-2) + 16*u(i-1)
%
% - 30*u(i) + 16*u(i+1) - 1*u(i+2) (9)
%
% + O(dx**4)
%
% Note that the ux term drops out, i.e., the basic equation is
%
% -2*a - b + c + 2*d =
%
% -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0
%
% Equation (9) was obtained by dropping all terms in the underlying
% Taylor series up to and including the fifth derivative, uxxxxx.
% Thus, equation (9) is exact for polynomials up to and including
% fifth order. This can be checked by substituting the functions
% 1, x, x**2, x**3, x**4 and x**5 in equation (9) and computing the
% corresponding derivatives for comparison with the known second
% derivatives. This is done for 1 merely by summing the weighting
% coefficients in equation (9), which should sum to zero, i.e.,
% -1 + 16 - 30 + 16 -1 = 0.
%
% For the remaining functions, the algebra is rather involved, but
% these functions can be checked numerically, i.e., numerical values
% of x**2, x**3, x**4 and x**5 can be substituted in equation (9)
% and the computed derivatives can be compared with the know numeri-
% cal second derivatives. This is not a proof of correctness of
% equation (9), but would likely detect any errors in equation (9).
%
% ******************************************************************
%
% (2) uxx at the interior points i = 2 and n-1
%
% For grid point 2, we consider a formula of the form
%
% a*u(i-1) + f*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) + e*u(i+4)
%
% Taylor series expansions of u(i-1), u(i+1), u(i+2), u(i+3) and
% u(i+4) when substituted into this formula give linear algebraic
% equations relating a, b, c, d and e.
%
% To drop ux
%
% -a + b + 2*c + 3*d + 4*e = 0 (10)
%
% To retain uxx
%
% a + b + 4*c + 9*d + 16*e = 2 (11)
%
% To drop uxxx
%
% -a + b + 8*c + 27*d + 64*e = 0 (12)
%
% To drop uxxxx
%
% a + b + 16*c + 81*d + 256*e = 0 (13)
%
% To drop uxxxxx
%
% -a + b + 32*c + 243*d + 1024*e = 0 (14)
%
% Equations (11), (12), (13) and (14) can be solved for a, b, c,
% d and e. If equation (10) is added to equation (11)
%
% 2*b + 6*c + 12*d +20*e = 2 (15)
%
% If equation (10) is subtracted from equation (12)
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (10) is added to equation (13)
%
% 2*b + 18*c + 84*d + 260*e = 0 (17)
%
% If equation (10) is subtracted from equation (14)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (15), (16), (17) and (18) can be solved for b, c, d
% and e.
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (15) is subtracted from equation (17)
%
% 12*c + 72*d + 240*e = -2 (19)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (16), (18) and (19) can be solved for c, d and e. If
% two times equation (16) is subtracted from equation (19),
%
% 24*d + 120*e = -2 (20)
%
% If five times equation (16) is subtracted from equation (18),
%
% 120*d + 720*e = 0 (21)
%
% Equations (20) and (21) can be solved for d and e. From (21),
% e = (-1/6)*d. Substitution in equation (20) gives d = -1/2.
% thus, e = 1/12. From equation (16), c = 7/6. From equation
% (15), b = -1/3. From equation (10), a = 5/6.
%
% The final differentiation formula is then obtained as
%
% (5/6)*u(i-1) + (-1/3)*u(i+1) + (7/6)*u(i+2) + (-1/2)*u(i+3)
%
% + (1/12)*u(i+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*u(i)
%
% + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(10*u(i-1) - 15*u(i) - 4*u(i+1)
% (22)
% + 14*u(i+2) - 6*u(i+3) + 1*u(i+4)) + O(dx**4)
%
% Equation (22) will be applied at i = 2 and n-1. thus
%
% uxx(2) = (1/12*dx**2)*(10*u(1) - 15*u(2) - 4*u(3)
% (23)
% + 14*u(4) - 6*u(5) + 1*u(6)) + O(dx**4)
%
% uxx(n-1) = (1/12*dx**2)*(10*u(n) - 15*u(n-1) - 4*u(n-2)
% (24)
% + 14*u(n-3) - 6*u(n-4) + 1*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% (3) uxx at the boundary points 1 and n
%
% Finally, for grid point 1, an approximation with a Neumann bound-
% ary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*ux(i) + f*u(i)
%
% Will be used. the corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + e = 0 (25)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d = 2 (26)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d = 0 (27)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d = 0 (28)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d = 0 (29)
%
% Equations (25) to (29) can be solved for a, b, c, d and e. If
%
% Equation (26) is subtracted from equations (27), (28) and (29),
%
% 4*b + 18*c + 48*d = -2 (30)
%
% 12*b + 72*c + 240*d = -2 (31)
%
% 28*b + 234*c + 1008*d = -2 (32)
%
% Equations (30), (31) and (32) can be solved for b, c and d
%
% 18*c + 96*d = 4 (33)
%
% 108*c + 672*d = 12 (34)
%
% Equations (3) and (34) can be solved for c and d, c = 8/9, d =
% -1/8.
%
% From equation (30), b = -3. From equation (26), a = 8. From
% equation (25), e = -25/6.
%
% The final differentiation formula is then obtained as
%
% 8*u(i+1) - 3*u(i+2) + (8/9)*u(i+3) - (1/8)*u(i+4)
%
% - (25/6)*ux(i)*dx
%
% = (8 - 3 + (8/9) - (1/8))*u(i) + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*((-415/6)*u(i) + 96*u(i+1) - 36*u(i+2)
% (35)
% + (32/3)*u(i+3) - (3/2)*u(i+4) - 50*ux(i)*dx) + O(dx**4)
%
% Equation (35) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*((-415/6)*u(1) + 96*u(2) - 36*u(3)
% (36)
% + (32/3)*u(4) - (3/2)*u(5) - 50*ux(1)*dx) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*((-415/6)*u(n) + 96*u(n-1) - 36*u(n-2)
% (37)
% + (32/3)*u(n-3) - (3/2)*u(n-4) + 50*ux(n)*dx) + O(dx**4)
%
% Alternatively, for grid point 1, an approximation with a Dirichlet
% boundary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*u(i+5) + f*u(i)
%
% can be used. The corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + 5*e = 0 (38)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d + 25*e = 2 (39)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d + 125*e = 0 (40)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d + 625*e = 0 (41)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d + 3125*e = 0 (42)
%
% Equations (38), (39), (40), (41) amd (42) can be solved for a,
% b, c, d and e.
%
% 2*b + 6*c + 12*d + 20*e = 2 (43)
%
% 6*b + 24*c + 60*d + 120*e = 0 (44)
%
% 14*b + 78*c + 252*d + 620*e = 0 (45)
%
% 30*b + 240*c + 1020*d + 3120*e = 0 (46)
%
% Equations (43), (44), (45) and (46) can be solved for b, c, d
% and e
%
% 6*c + 24*d + 60*e = -6 (47)
%
% 36*c + 168*d + 480*e = -14 (48)
%
% 150*c + 840*d + 2820*e = -30 (49)
%
% Equations (47), (48) and (49) can be solved for c, d and e
%
% 24*d + 120*e = 22 (50)
%
% 240*d + 1320*e = 120 (51)
%
% From equations (50) and (51), d = 61/12, e = -5/6. From equation
% (47), c = -13. From equation (43), b = 107/6. From equation (38),
% a = -77/6.
%
% The final differentiation formula is then obtained as
%
% (-77/6)*u(i+1) + (107/6)*u(i+2) - 13*u(i+3) + (61/12)*u(i+4)
%
% - (5/6)*u(i+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*u(i) +
%
% uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(45*u(i) - 154*u(i+1) + 214*u(i+2)
% (52)
% - 156*u(i+3) + 61*u(i+4) - 10*u(i+5)) + O(dx**4)
%
% Equation (52) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*(45*u(1) - 154*u(2) + 214*u(3)
% (53)
% - 156*u(4) + 61*u(5) - 10*u(6)) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*(45*u(n) - 154*u(n-1) + 214*u(n-2)
% (54)
% -156*u(n-3) + 61*u(n-4) - 10*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% 1/(12*dx**2) for subsequent use
r12dxs=1./(12.0*dx^2);
%
% uxx at the left boundary
%
% Without ux (equation (53))
if nl==1
uxx(1)=r12dxs*...
( 45.0*u(1)...
-154.0*u(2)...
+214.0*u(3)...
-156.0*u(4)...
+61.0*u(5)...
-10.0*u(6));
%
% With ux (equation (36))
elseif nl==2
uxx(1)=r12dxs*...
(-415.0/6.0*u(1)...
+96.0*u(2)...
-36.0*u(3)...
+32.0/3.0*u(4)...
-3.0/2.0*u(5)...
-50.0*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux (equation (54))
if nu==1
uxx(n)=r12dxs*...
( 45.0*u(n )...
-154.0*u(n-1)...
+214.0*u(n-2)...
-156.0*u(n-3)...
+61.0*u(n-4)...
-10.0*u(n-5));
%
% With ux (equation (37))
elseif nu==2
uxx(n)=r12dxs*...
(-415.0/6.0*u(n )...
+96.0*u(n-1)...
-36.0*u(n-2)...
+32.0/3.0*u(n-3)...
-3.0/2.0*u(n-4)...
+50.0*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2 (equation (23))
uxx(2)=r12dxs*...
( 10.0*u(1)...
-15.0*u(2)...
-4.0*u(3)...
+14.0*u(4)...
-6.0*u(5)...
+1.0*u(6));
%
% i = n-1 (equation (24))
uxx(n-1)=r12dxs*...
( 10.0*u(n )...
-15.0*u(n-1)...
-4.0*u(n-2)...
+14.0*u(n-3)...
-6.0*u(n-4)...
+1.0*u(n-5));
%
% i = 3, 4,..., n-2 (equation (9))
for i=3:n-2
uxx(i)=r12dxs*...
( -1.0*u(i-2)...
+16.0*u(i-1)...
-30.0*u(i )...
+16.0*u(i+1)...
-1.0*u(i+2));
end
