dss048
From Scholarpedia
< Method of lines | example implementation
| Samir Hamdi et al. (2007), Scholarpedia, 2(10):2859. | revision #25176 [link to/cite this article] | |||||||||||||||||||
Curator: Samir Hamdi, GENIVAR (Hydropower & Hydraulics) and Solid Mechanics Laboratory, Ecole Polytechnique, Paris, France
Curator: William E. Schiesser, Lehigh University and University of Pennsylvania, USA
Curator: Graham W Griffiths, School of Engineering and Mathematical Sciences, City University, London, UK
% File: dss048.m
%
function [uxx]=dss048(xl,xu,n,u,ux,nl,nu)
%
% Function dss048 computes a eighth-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 eighth-order, ten-point approximations for a
% second derivative 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.
%
% Weighting coefficients for finite difference approximations, computed
% by an algorithm of B. Fornberg (1,2), are used in the following
% approximations. The combination of these coefficients to give the
% second derivative approximations at the boundaries in x, including
% the first derivative, was suggested by Professor Fornberg.
%
% (1) Fornberg, B., Fast Generation of Weights in Finite Difference
% Formulas, In Recent Developments In Numerical Methods And
% Software for ODEs/DAEs/PDEs, G. Byrne and W. E. Schiesser (eds),
% World Scientific, River Edge, NJ, 1992
%
% (2) Fornberg, B., Calculation of Weights In Finite Difference
% Formulas, SIAM Review, Vol. 40, No. 3, Pp 685-691, September,
% 1999
%
% Grid spacing
dx=(xu-xl)/(n-1);
rdxs=1.0/dx^2;
%
% uxx at the left boundary
%
% Without ux
if nl==1
uxx(1)=rdxs*...
( 6.463293650793652*u(1)...
-32.921428571428578*u(2)...
+83.842857142857156*u(3)...
-139.755555555555555*u(4)...
+162.474999999999999*u(5)...
-132.500000000000000*u(6)...
+74.544444444444444*u(7)...
-27.628571428571430*u(8)...
+6.083928571428572*u(9)...
-0.603968253968254*u(10));
%
% With ux
elseif nl==2
uxx(1)=rdxs*...
( -8.914169501133784*u(1)...
+16.000000000000000*u(2)...
-14.000000000000000*u(3)...
+12.444444444444444*u(4)...
-8.750000000000000*u(5)...
+4.480000000000000*u(6)...
-1.555555555555555*u(7)...
+0.326530612244898*u(8)...
-0.031250000000000*u(9)...
-5.435714285714286*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux
if nu==1
uxx(n)=rdxs*...
( -0.603968253968254*u(n-9)...
+6.083928571428572*u(n-8)...
-27.628571428571430*u(n-7)...
+74.544444444444444*u(n-6)...
-132.500000000000000*u(n-5)...
+162.474999999999999*u(n-4)...
-139.755555555555555*u(n-3)...
+83.842857142857156*u(n-2)...
-32.921428571428578*u(n-1)...
+6.463293650793652*u(n ));
%
% With ux
elseif nu==2
uxx(n)=rdxs*...
(-0.031250000000000*u(n-8)...
+0.326530612244905*u(n-7)...
-1.555555555555555*u(n-6)...
+4.480000000000000*u(n-5)...
-8.750000000000000*u(n-4)...
+12.444444444444444*u(n-3)...
-14.000000000000000*u(n-2)...
+16.000000000000000*u(n-1)...
-8.914169501133790*u(n )...
+5.435714285714287*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2
uxx(2)=rdxs*...
(0.603968253968254*u(1)...
+0.423611111111111*u(2)...
-5.742857142857143*u(3)...
+11.366666666666665*u(4)...
-12.922222222222221*u(5)...
+10.274999999999999*u(6)...
-5.666666666666665*u(7)...
+2.068253968253968*u(8)...
-0.450000000000000*u(9)...
+0.044246031746032*u(10));
%
% i = 3
uxx(3)=rdxs*...
(-0.044246031746032*u(1)...
+1.046428571428571*u(2)...
-1.567460317460318*u(3)...
-0.433333333333333*u(4)...
+2.075000000000000*u(5)...
-1.772222222222222*u(6)...
+0.983333333333333*u(7)...
-0.357142857142857*u(8)...
+0.077182539682540*u(9)...
-0.007539682539683*u(10));
%
% i = 4
uxx(4)=rdxs*...
(0.007539682539683*u(1)...
-0.119642857142857*u(2)...
+1.385714285714286*u(3)...
-2.472222222222222*u(4)...
+1.150000000000000*u(5)...
+0.175000000000000*u(6)...
-0.188888888888889*u(7)...
+0.078571428571429*u(8)...
-0.017857142857143*u(9)...
+0.001785714285714*u(10));
%
% i = n-3
uxx(n-3)=rdxs*...
(0.001785714285714*u(n-9)...
-0.017857142857143*u(n-8)...
+0.078571428571429*u(n-7)...
-0.188888888888889*u(n-6)...
+0.175000000000000*u(n-5)...
+1.150000000000000*u(n-4)...
-2.472222222222224*u(n-3)...
+1.385714285714286*u(n-2)...
-0.119642857142857*u(n-1)...
+0.007539682539683*u(n ));
%
% i = n-2
uxx(n-2)=rdxs*...
(-0.007539682539683*u(n-9)...
+0.077182539682540*u(n-8)...
-0.357142857142857*u(n-7)...
+0.983333333333333*u(n-6)...
-1.772222222222222*u(n-5)...
+2.075000000000001*u(n-4)...
-0.433333333333335*u(n-3)...
-1.567460317460318*u(n-2)...
+1.046428571428572*u(n-1)...
-0.044246031746032*u(n ));
%
% i = n-1
uxx(n-1)=rdxs*...
(0.044246031746032*u(n-9)...
-0.450000000000000*u(n-8)...
+2.068253968253968*u(n-7)...
-5.666666666666667*u(n-6)...
+10.275000000000000*u(n-5)...
-12.922222222222222*u(n-4)...
+11.366666666666667*u(n-3)...
-5.742857142857142*u(n-2)...
+0.423611111111112*u(n-1)...
+0.603968253968254*u(n ));
%
% i = 5, 6,..., n-4
for i=5:n-4
uxx(i)=rdxs*...
(-0.001785714285714*u(i-4)...
+0.025396825396825*u(i-3)...
-0.200000000000000*u(i-2)...
+1.600000000000000*u(i-1)...
-2.847222222222221*u(i )...
+1.600000000000000*u(i+1)...
-0.200000000000000*u(i+2)...
+0.025396825396825*u(i+3)...
-0.001785714285714*u(i+4));
end
%
