DSS.rar_dss_run

% 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) % %