Method of lines/example implementation/dss042
From Scholarpedia
% File: dss042.m
%
function [uxx]=dss042(xl,xu,n,u,ux,nl,nu)
%
% Function dss042 computes a second-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 is for a set of second-order, four-point
% approximations for a second derivative that can be used at the
% boundaries of a spatial domain. These approximations have the
% features
%
% (1) Only interior and boundary points are used (i.e., no
% fictitious points are used)
%
% (2) The normal derivative at the boundary is included as part
% of the approximation for the second derivative
%
% (3) Approximations for the boundary conditions are not used.
%
% The derivation is by Professor Gilbert A. Stengle, Department of
% Mathematics, Lehigh University, Bethlehem, pa 18015, and was done
% on December 7, 1985.
%
% For an approximation at the left boundary, involving the points
% i, i+1, i+2 and i+3, consider the following Taylor series expan-
% sions
%
% ux(i)( dx) uxx(i)( dx)**2 uxxx(i)( dx)**3
% u(i+1) = u(i) + ---------- + -------------- + --------------- +...
% 1 2 6
%
%
% ux(i)(2dx) uxx(i)(2dx)**2 uxxx(i)(2dx)**3
% u(i+2) = u(i) + ---------- + -------------- + --------------- +...
% 1 2 6
%
% If we now form the following linear combination, involving con-
% stants a, b, c and d to be determined, and use the preceding two
% Taylor series,
%
% a*u(i) + b*ux(i) + c*u(i+1) + d*u(i+2)
%
% we have
%
% a*u(i) + b*ux(i) + c*u(i+1) + d*u(i+2) =
%
% (a + b + c + d)*u(i) +
%
% (b + dx*c + 2*dx*d)*ux(i) +
%
% (c*(dx**2)/2 + d*((2*dx)**2)/2)*uxx(i) +
%
% (c*(dx**3)/6 + d*((2*dx)**3)/6)*uxxx(i) + O(dx**4)
%
% The third derivative, uxxx(i), can be dropped by taking
%
% c = -8*d
%
% The second derivative, uxx(i), can be retained by taking
%
% (dx**2)(c/2 + 2*d) = 1
%
% which, when combined with the preceding result gives
%
% d = -1/(2*(dx**2))
%
% c = 4/(dx**2)
%
% The first derivative, ux(i), can be dropped by taking
%
% b + dx*c + 2*dx*d = 0
%
% or
%
% b = -dx*c - 2*dx*d = -4/dx - 2*dx*(-1/(2*(dx**2))) = -3/dx
%
% Finally, u(i), can be dropped by taking
%
% a = - c - d = 8*d - d = -7*d = -7/(2*(dx**2))
%
% If we now solve for the derivative of interest, uxx(i),
%
% uxx(i) = -7/(2(dx**2))*u(i) - 3/dx*ux(i)
%
% + 8/(dx**2)*u(i+1) - 1/(2*(dx**2))u(i+2) + O(dx**2)
%
% = (1/(2*(dx**2)))*(-u(i+2) + 8*u(i+1) - 7*u(i) - 6*dx*ux(i))
%
% + O(dx**2)
%
% which is the four-point, second-order approximation for the second
% derivative, uxx(i), including the first derivative, ux(i).
%
% Four checks of this approximation can easily be made for u(i) =
% 1, u(i) = x, u(i) = x**2 and u(i) = x**3
%
% uxx(i) = (1/(2*(dx**2)))*(-1 + 8*1 - 7*1 - 6*dx*0) = 0
%
% uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx) + 8*(x + dx)
%
% -7*x - 6*dx*1) = 0
%
% uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx)**2 + 8*(x + dx)**2
%
% - 7*(x**2) - 6*dx*(2*x))
%
% = (- x**2 - 4*x*dx - 4*dx**2
%
% + 8*x**2 + 16*x*dx + 8*dx**2
%
% - 7*x**2 - 12*x*dx)/(2*(dx**2)) = 2
%
% uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx)**3 + 8*(x + dx)**3
%
% - 7*(x**3) - 6*dx*(3*x**2))
%
% = (1/(2*(dx**2)))*(- x**3 - 6*dx*x**2 - 12*x*dx**2
%
% - 8*dx**3 + 8*x**3 + 24*dx*x**2 + 24*x*dx**2 + 8*dx**3
%
% - 7*x**3 - 18*dx*x**2)
%
% = (1/(2*(dx**2)))*(12*x*dx**2) = 6*x
%
% The preceding approximation for uxx(i) can be applied at the
% left boundary value of x by taking i = 1. An approximation at
% the right boundary is obtained by taking dx = -dx and reversing
% the subscripts in the preceding approximation, with i = n
%
% uxx(i)
%
% = (1/(2*(dx**2)))*(-u(i-2) + 8*u(i-1) - 7*u(i) + 6*dx*ux(i))
%
% + O(dx**2)
%
% To obtain approximations of the second derivative which do not
% involve the first derivative, we take as the linear combination
%
% a*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3)
%
% we have
%
% a*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) =
%
% (a + b + c + d)*u(i)+
%
% (dx*b + 2*dx*c + 4*dx*d)*ux(i)+
%
% (b*(dx**2)/2 + c*((2*dx)**2)/2 + d*((3*dx)**2)/2)*uxx(i) +
%
% (b*(dx**3)/6 + c*((2*dx)**3)/6 + d*((3*dx)**3)/6)*uxx(i) +
%
% O(dx**4)
%
% The third derivative, uxxx(i), can be dropped by taking
%
% b + 8*c + 27*d = 0
%
% The second derivative, uxx(i), can be retained by taking
%
% (dx**2)*(b/2 + 2*c + (9/2)*d) = 1
%
% The first derivative can be dropped by taking
%
% b + 2*c + 3*d = 0
%
% Solution of the preceding equations for c and d by elimination of
% b gives
%
% 6*c + 24*d = 0
%
% 4*c + 18*d = -2/(dx**2)
%
% Then, eliminating c, gives
%
% (18 - 16)*d = -2/(dx**2)
%
% or
%
% d = -1/(dx**2)
%
% c = (24/6)/(dx**2) = 4/(dx**2)
%
% b = -8/(dx**2) + 3/(dx**2) = -5/(dx**2)
%
% u(i) can be dropped by taking
%
% a + b + c + d = 0
%
% or
%
% a = (5 - 4 + 1)/(dx**2) = 2/(dx**2)
%
% If we now solve for the derivative of interest, uxx(i),
%
% uxx(i) = (1/dx**2)*(2*u(i) - 5*u(i+1) + 4*u(i+2) - 1*u(i+3))
%
% + O(dx**2)
%
% Which is the four-point, second-order approximation for the second
% derivative, uxx(i), without the first derivative, ux(i).
%
% Four checks of this approximation can easily be made for u(i) =
% 1, u(i) = x, u(i) = x**2 and u(i) = x**3
%
% uxx(i) = (1/dx**2)*(2 - 5 + 4 - 1) = 0
%
% uxx(i) = (1/dx**2)*(2*x - 5*(x + dx) + 4*(x + 2*dx)
%
% - 1*(x + 3*dx)) = 0
%
% uxx(i) = (1/dx**2)*(2*x**2 - 5*(x + dx)**2 + 4*(x + 2*dx)**2
%
% - 1*(x + 3*dx)**2) = 2
%
% uxx(i) = (1/dx**2)*(2*x**3 - 5*(x + dx)**3 + 4*(x + 2*dx)**3
%
% - 1*(x + 3*dx)**3)
%
% = (1/dx**2)*(2*x**3 - 5*x**3 - 15*x*dx**2
%
% - 15*dx*x**2 - 5*dx**3 + 4*x**3 + 24*dx*x**2
%
% + 48*x*dx**2 + 32*dx**3 - x**3 - 9*dx*x**2
%
% - 27*x*dx**2 - 27dx**3)
%
% = (1/dx**2)*(6*x*dx**2) = 6*x
%
% The preceding approximation for uxx(i) can be applied at the
% left boundary value of x by taking i = 1. An approximation at
% the right boundary is obtained by taking dx = -dx and reversing
% the subscripts in the preceding approximation, with i = n
%
% uxx(i) = (1/dx**2)*(2*u(i) - 5*u(i-1) + 4*u(i-2) - 1*u(i-3))
%
% + O(dx**2)
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% uxx at the left boundary, without ux
if nl==1
uxx(1)=(( 2.)*u( 1)...
+(-5.)*u( 2)...
+( 4.)*u( 3)...
+(-1.)*u( 4))/(dx^2);
%
% uxx at the left boundary, including ux
elseif nl==2
uxx(1)=((-7.)*u( 1)...
+( 8.)*u( 2)...
+(-1.)*u( 3))/(2.*dx^2)...
+(-6.)*ux( 1) /(2.*dx);
end
%
% uxx at the right boundary, without ux
if nu==1
uxx(n)=(( 2.)*u(n )...
+(-5.)*u(n-1)...
+( 4.)*u(n-2)...
+(-1.)*u(n-3))/(dx^2);
%
% uxx at the right boundary, including ux
elseif nu==2
uxx(n)=((-7.)*u(n )...
+( 8.)*u(n-1)...
+(-1.)*u(n-2))/(2.*dx^2)...
+( 6.)*ux(n ) /(2.*dx);
end
%
% uxx at the interior grid points
for i=2:n-1
uxx(i)=(u(i+1)-2.*u(i)+u(i-1))/dx^2;
end
%
