dss050
From Scholarpedia
< Method of lines | example implementation
| Samir Hamdi et al. (2007), Scholarpedia, 2(10):2859. | revision #25165 [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: dss050.m
%
function [uxx]=dss050(xl,xu,n,u,ux,nl,nu)
%
% Function dss050 computes a tenth-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 tenth-order, twelve-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);
%
% 1/dx**2 for subsequent use
rdxs=1.0/dx^2;
%
% uxx at the left boundary
%
% Without ux
if nl==1
uxx(1)=rdxs*...
( 7.561626984126982*u(1)...
-44.437301587301590*u(2)...
+138.593253968253980*u(3)...
-295.519841269841320*u(4)...
+457.029761904761930*u(5)...
-521.113333333333460*u(6)...
+439.394444444444330*u(7)...
-271.261904761904700*u(8)...
+119.413690476190470*u(9)...
-35.551587301587304*u(10)...
+6.423730158730158*u(11)...
-0.532539682539683*u(12));
%
% With ux
elseif nl==2
uxx(1)=rdxs*...
(-10.128622763920379*u(1)...
+20.000000000000000*u(2)...
-22.499999999999972*u(3)...
+26.666666666666572*u(4)...
-26.249999999999943*u(5)...
+20.159999999999968*u(6)...
-11.666666666666686*u(7)...
+4.897959183673493*u(8)...
-1.406250000000000*u(9)...
+0.246913580246911*u(10)...
-0.020000000000001*u(11)...
-5.857936507936508*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux
if nu==1
uxx(n)=rdxs*...
( -0.532539682539683*u(n-11)...
+6.423730158730161*u(n-10)...
-35.551587301587304*u(n-9)...
+119.413690476190480*u(n-8)...
-271.261904761904760*u(n-7)...
+439.394444444444390*u(n-6)...
-521.113333333333340*u(n-5)...
+457.029761904761930*u(n-4)...
-295.519841269841270*u(n-3)...
+138.593253968253980*u(n-2)...
-44.437301587301590*u(n-1)...
+7.561626984126984*u(n ));
%
% With ux
elseif nu==2
uxx(n)=rdxs*...
( -0.019999999999999*u(n-10)...
+0.246913580246904*u(n-9)...
-1.406249999999986*u(n-8)...
+4.897959183673493*u(n-7)...
-11.666666666666686*u(n-6)...
+20.159999999999968*u(n-5)...
-26.249999999999943*u(n-4)...
+26.666666666666629*u(n-3)...
-22.499999999999972*u(n-2)...
+20.000000000000000*u(n-1)...
-10.128622763920383*u(n )...
+5.857936507936508*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2
uxx(2)=rdxs*...
( 0.532539682539683*u(1)...
+1.171150793650794*u(2)...
-9.289682539682540*u(3)...
+21.434523809523807*u(4)...
-31.912698412698408*u(5)...
+35.258333333333333*u(6)...
-29.046666666666663*u(7)...
+17.623015873015873*u(8)...
-7.654761904761902*u(9)...
+2.254960317460317*u(10)...
-0.403968253968254*u(11)...
+0.033253968253968*u(12));
%
% i = 3
uxx(3)=rdxs*...
( -0.033253968253968*u(1)...
+0.931587301587302*u(2)...
-1.023611111111112*u(3)...
-1.973809523809523*u(4)...
+4.973809523809524*u(5)...
-5.575555555555556*u(6)...
+4.531666666666666*u(7)...
-2.709523809523809*u(8)...
+1.162301587301588*u(9)...
-0.338888888888889*u(10)...
+0.060198412698413*u(11)...
-0.004920634920635*u(12));
%
% i = 4
uxx(4)=rdxs*...
( 0.004920634920635*u(1)...
-0.092301587301587*u(2)...
+1.256349206349206*u(3)...
-2.106150793650793*u(4)...
+0.461904761904762*u(5)...
+1.076666666666667*u(6)...
-1.028888888888889*u(7)...
+0.634523809523810*u(8)...
-0.273809523809524*u(9)...
+0.079761904761905*u(10)...
-0.014126984126984*u(11)...
+0.001150793650794*u(12));
%
% i = 5
uxx(5)=rdxs*...
( -0.001150793650794*u(1)...
+0.018730158730159*u(2)...
-0.168253968253968*u(3)...
+1.509523809523809*u(4)...
-2.675793650793650*u(5)...
+1.373333333333333*u(6)...
+0.013333333333334*u(7)...
-0.117460317460318*u(8)...
+0.064880952380952*u(9)...
-0.020634920634921*u(10)...
+0.003809523809524*u(11)...
-0.000317460317460*u(12));
%
% i = n-4
uxx(n-4)=rdxs*...
(-0.000317460317460*u(n-11)...
+0.003809523809524*u(n-10)...
-0.020634920634921*u(n-9)...
+0.064880952380952*u(n-8)...
-0.117460317460317*u(n-7)...
+0.013333333333333*u(n-6)...
+1.373333333333333*u(n-5)...
-2.675793650793652*u(n-4)...
+1.509523809523810*u(n-3)...
-0.168253968253968*u(n-2)...
+0.018730158730159*u(n-1)...
-0.001150793650794*u(n ));
%
% i = n-3
uxx(n-3)=rdxs*...
( 0.001150793650794*u(n-11)...
-0.014126984126984*u(n-10)...
+0.079761904761905*u(n-9)...
-0.273809523809524*u(n-8)...
+0.634523809523809*u(n-7)...
-1.028888888888889*u(n-6)...
+1.076666666666667*u(n-5)...
+0.461904761904762*u(n-4)...
-2.106150793650793*u(n-3)...
+1.256349206349206*u(n-2)...
-0.092301587301587*u(n-1)...
+0.004920634920635*u(n ));
%
% i = n-2
uxx(n-2)=rdxs*...
(-0.004920634920635*u(n-11)...
+0.060198412698413*u(n-10)...
-0.338888888888889*u(n-9)...
+1.162301587301588*u(n-8)...
-2.709523809523809*u(n-7)...
+4.531666666666666*u(n-6)...
-5.575555555555556*u(n-5)...
+4.973809523809526*u(n-4)...
-1.973809523809526*u(n-3)...
-1.023611111111110*u(n-2)...
+0.931587301587302*u(n-1)...
-0.033253968253968*u(n ));
%
% i = n-1
uxx(n-1)=rdxs*...
( 0.033253968253968*u(n-11)...
-0.403968253968254*u(n-10)...
+2.254960317460318*u(n-9)...
-7.654761904761904*u(n-8)...
+17.623015873015873*u(n-7)...
-29.046666666666670*u(n-6)...
+35.258333333333340*u(n-5)...
-31.912698412698411*u(n-4)...
+21.434523809523810*u(n-3)...
-9.289682539682541*u(n-2)...
+1.171150793650794*u(n-1)...
+0.532539682539683*u(n ));
%
% i = 6, 7,..., n-5
for i=6:n-5
uxx(i)=rdxs*...
( 0.000317460317460*u(i-5)...
-0.004960317460317*u(i-4)...
+0.039682539682540*u(i-3)...
-0.238095238095238*u(i-2)...
+1.666666666666667*u(i-1)...
-2.927222222222222*u(i )...
+1.666666666666667*u(i+1)...
-0.238095238095238*u(i+2)...
+0.039682539682540*u(i+3)...
-0.004960317460317*u(i+4)...
+0.000317460317460*u(i+5));
end
%
