dss006=function(xl,xu,n,u) { # # An extensive set of documentation comments detailing the derivation # of the following second order finite differences (FDs) is not given # here to conserve space. The comments are available in the Matlab # version of dss006 in http://www.pdecomp.net. The derivation is also # detailed in Schiesser, W. E., The Numerical Method of Lines Integration # of Partial Differential Equations, Academic Press, San Diego, 1991. # # Preallocate arrays ux=rep(0,n); # # Grid spacing dx=(xu-xl)/(n-1); # # 1/(6!*dx) for subsequent use r6fdx=1./(720*dx); # # ux vector # # Boundaries (x=xl,x=xu) ux[1]=r6fdx* (-1764*u[ 1]+4320*u[ 2]-5400*u[ 3]+4800*u[ 4]- 2700*u[ 5] +864*u[ 6] -120*u[ 7]); ux[n]=r6fdx* ( 1764*u[ n]-4320*u[n-1]+5400*u[n-2]-4800*u[n-3]+ 2700*u[n-4] -864*u[n-5] +120*u[n-6]); # # dx in from boundaries (x=xl+dx,x=xu-dx) ux[ 2]=r6fdx* (-120*u[ 1]-924*u[ 2]+1800*u[ 3]-1200*u[ 4]+ 600*u[ 5]-180*u[ 6] +24*u[ 7]); ux[n-1]=r6fdx* ( 120*u[ n]+924*u[n-1]-1800*u[n-2]+1200*u[n-3]- 600*u[n-4]+180*u[n-5] -24*u[n-6]); # # 2*dx in from boundaries (x=xl+2*dx,x=xu-2*dx) ux[ 3]=r6fdx* ( 24*u[ 1]-288*u[ 2]-420*u[ 3]+960*u[ 4]- 360*u[ 5] +96*u[ 6] -12*u[ 7]); ux[n-2]=r6fdx* (-24*u[ n]+288*u[n-1]+420*u[n-2]-960*u[n-3]+ 360*u[n-4] -96*u[n-5] +12*u[n-6]); # # Interior points (x=xl+3*dx,...,x=xu-3*dx) for(i in 4:(n-3)) ux[i]=r6fdx* (-12*u[i-3]+108*u[i-2]-540*u[i-1]+ 12*u[i+3]-108*u[i+2]+540*u[i+1]); # # All points concluded (x=xl,...,x=xu) return(c(ux)); }