dss048=function(xl,xu,n,u,ux,nl,nu) { # # An extensive set of documentation comments detailing the derivation # of the following eighth order finite differences (FDs) is not given # here to conserve space. The comments are available in the Matlab # version of dss048 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); uxx=rep(0,n); # # Grid spacing dx=(xu-xl)/(n-1); # # 1/(dx**2) for subsequent use rdxs=1/dx^2; # # uxx vector # # Boundaries (x=xl,x=xu) 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]); if(nu==1) uxx[n]=rdxs* ( 6.463293650793652*u[n ] -32.921428571428578*u[n-1]+ 83.842857142857156*u[n-2]-139.755555555555555*u[n-3]+ 162.474999999999999*u[n-4]-132.500000000000000*u[n-5]+ 74.544444444444444*u[n-6] -27.628571428571430*u[n-7]+ 6.083928571428572*u[n-8]- 0.603968253968254*u[n-9]); if(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); if(nu==2) uxx[n]=rdxs* ( -8.914169501133790*u[n ]+16.000000000000000*u[n-1]- 14.000000000000000*u[n-2]+12.444444444444444*u[n-3]- 8.750000000000000*u[n-4] +4.480000000000000*u[n-5]- 1.555555555555555*u[n-6] +0.326530612244905*u[n-7]- 0.031250000000000*u[n-8]+ 5.435714285714287*ux[n ]*dx); # # dx in from boundaries (x=xl+dx,x=xu-dx) 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]); uxx[n-1]=rdxs* ( 0.603968253968254*u[n ] +0.423611111111112*u[n-1]- 5.742857142857142*u[n-2]+11.366666666666667*u[n-3]- 12.922222222222222*u[n-4]+10.275000000000000*u[n-5]- 5.666666666666667*u[n-6] +2.068253968253968*u[n-7]- 0.450000000000000*u[n-8] +0.044246031746032*u[n-9]); # # 2*dx in from boundaries (x=xl+2*dx,x=xu-2*dx) 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]); uxx[n-2]=rdxs* (-0.044246031746032*u[n ]+1.046428571428572*u[n-1]- 1.567460317460318*u[n-2]-0.433333333333335*u[n-3]+ 2.075000000000001*u[n-4]-1.772222222222222*u[n-5]+ 0.983333333333333*u[n-6]-0.357142857142857*u[n-7]+ 0.077182539682540*u[n-8]-0.007539682539683*u[n-9]); # # 3*dx in from boundaries (x=xl+3*dx,x=xu-3*dx) 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]); uxx[n-3]=rdxs* (0.007539682539683*u[n ]-0.119642857142857*u[n-1]+ 1.385714285714286*u[n-2]-2.472222222222222*u[n-3]+ 1.150000000000000*u[n-4]+0.175000000000000*u[n-5]- 0.188888888888889*u[n-6]+0.078571428571429*u[n-7]- 0.017857142857143*u[n-8]+0.001785714285714*u[n-9]); # # Remaining interior points (x=xl+4*dx,...,x=xu-4*dx) for(i in 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]- 0.001785714285714*u[i+4]+0.025396825396825*u[i+3]- 0.200000000000000*u[i+2]+1.600000000000000*u[i+1]- 2.847222222222221*u[i ]); # # All points concluded (x=xl,...,x=xu) return(c(uxx)); }