dss008=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 dss008 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/(8!*dx) for subsequent use r8fdx=1/(40320*dx); # # ux vector # # Boundaries (x=xl,x=xu) ux[1]=r8fdx* (-109584*u[ 1]+322560*u[ 2]-564480*u[ 3]+752640*u[ 4]- 705600*u[ 5]+451584*u[ 6]-188160*u[ 7] +46080*u[ 8]- 5040*u[ 9]); ux[n]=r8fdx* ( 109584*u[n ]-322560*u[n-1]+564480*u[n-2]-752640*u[n-3]+ 705600*u[n-4]-451584*u[n-5]+188160*u[n-6] -46080*u[n-7]+ 5040*u[n-8]); # # dx in from boundaries (x=xl+dx,x=xu-dx) ux[ 2]=r8fdx* ( -5040*u[ 1]-64224*u[ 2]+141120*u[ 3]-141120*u[ 4]+ 117600*u[ 5]-70560*u[ 6] +28224*u[ 7] -6720*u[ 8]+ 720*u[ 9]); ux[n-1]=r8fdx* ( 5040*u[n ]+64224*u[n-1]-141120*u[n-2]+141120*u[n-3]- 117600*u[n-4]+70560*u[n-5] -28224*u[n-6] +6720*u[n-7]- 720*u[n-8]); # # 2*dx in from boundaries (x=xl+2*dx,x=xu-2*dx) ux[ 3]=r8fdx* ( 720*u[ 1]-11520*u[ 2]-38304*u[ 3]+80640*u[ 4]- 50400*u[ 5]+26880*u[ 6]-10080*u[ 7] +2304*u[ 8]- 240*u[ 9]); ux[n-2]=r8fdx* ( -720*u[n ]+11520*u[n-1]+38304*u[n-2]-80640*u[n-3]+ 50400*u[n-4]-26880*u[n-5]+10080*u[n-6] -2304*u[n-7]+ 240*u[n-8]); # # 3*dx in from boundaries (x=xl+3*dx,x=xu-3*dx) ux[ 4]=r8fdx* ( -240*u[ 1] +2880*u[ 2]-20160*u[ 3]-18144*u[ 4]+ 50400*u[ 5]-20160*u[ 6] +6720*u[ 7] -1440*u[ 8]+ 144*u[ 9]); ux[n-3]=r8fdx* ( 240*u[n ] -2880*u[n-1]+20160*u[n-2]+18144*u[n-3]- 50400*u[n-4]+20160*u[n-5] -6720*u[n-6] +1440*u[n-7]- 144*u[n-8]); # # Interior points (x=xl+4*dx,...,x=xu-4*dx) for(i in 5:(n-4)) ux[ i]=r8fdx* ( 144*u[i-4]-1536*u[i-3]+8064*u[i-2]-32256*u[i-1]- 144*u[i+4]+1536*u[i+3]-8064*u[i+2]+32256*u[i+1]); # # All points concluded (x=xl,...,x=xu) return(c(ux)); }