dss050=function(xl,xu,n,u,ux,nl,nu) { # # An extensive set of documentation comments detailing the derivation # of the following tenth order finite differences (FDs) is not given # here to conserve space. The comments are available in the Matlab # version of dss050 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* ( 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]); if(nu==1) uxx[n]=rdxs* ( 7.561626984126982*u[ n] -44.437301587301590*u[n-1]+ 138.593253968253980*u[n-2]-295.519841269841320*u[n-3]+ 457.029761904761930*u[n-4]-521.113333333333460*u[n-5]+ 439.394444444444330*u[n-6]-271.261904761904700*u[n-7]+ 119.413690476190470*u[n-8] -35.551587301587304*u[n-9]+ 6.423730158730158*u[n-10] -0.532539682539683*u[n-11]); if(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); if(nu==2) uxx[n]=rdxs* (-10.128622763920379*u[n ]+20.000000000000000*u[n-1]- 22.499999999999972*u[n-2]+26.666666666666572*u[n-3]- 26.249999999999943*u[n-4]+20.159999999999968*u[n-5]- 11.666666666666686*u[n-6] +4.897959183673493*u[n-7]- 1.406250000000000*u[n-8] +0.246913580246911*u[n-9]- 0.020000000000001*u[n-10]+5.857936507936508*ux[n ]*dx); # # dx in from boundaries (x=xl+dx,x=xu-dx) 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]); uxx[n-1]=rdxs* ( 0.532539682539683*u[n ] +1.171150793650794*u[n-1]- 9.289682539682540*u[n-2]+21.434523809523807*u[n-3]- 31.912698412698408*u[n-4]+35.258333333333333*u[n-5]- 29.046666666666663*u[n-6]+17.623015873015873*u[n-7]- 7.654761904761902*u[n-8] +2.254960317460317*u[n-9]- 0.403968253968254*u[n-10]+0.033253968253968*u[n-11]); # # 2*dx in from boundaries (x=xl+2*dx,x=xu-2*dx) 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]); uxx[n-2]=rdxs* ( -0.033253968253968*u[n ]+0.931587301587302*u[n-1]- 1.023611111111112*u[n-2]-1.973809523809523*u[n-3]+ 4.973809523809524*u[n-4]-5.575555555555556*u[n-5]+ 4.531666666666666*u[n-6]-2.709523809523809*u[n-7]+ 1.162301587301588*u[n-8]-0.338888888888889*u[n-9]+ 0.060198412698413*u[n-10]-0.004920634920635*u[n-11]); # # 3*dx in from boundaries (x=xl+3*dx,x=xu-3*dx) 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]); uxx[n-3]=rdxs* ( 0.004920634920635*u[n ]-0.092301587301587*u[n-1]+ 1.256349206349206*u[n-2]-2.106150793650793*u[n-3]+ 0.461904761904762*u[n-4]+1.076666666666667*u[n-5]- 1.028888888888889*u[n-6]+0.634523809523810*u[n-7]- 0.273809523809524*u[n-8]+0.079761904761905*u[n-9]- 0.014126984126984*u[n-10]+0.001150793650794*u[n-11]); # # 4*dx in from boundaries (x=xl+4*dx,x=xu-4*dx) 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]); uxx[n-4]=rdxs* ( -0.001150793650794*u[n ]+0.018730158730159*u[n-1]- 0.168253968253968*u[n-2]+1.509523809523809*u[n-3]- 2.675793650793650*u[n-4]+1.373333333333333*u[n-5]+ 0.013333333333334*u[n-6]-0.117460317460318*u[n-7]+ 0.064880952380952*u[n-8]-0.020634920634921*u[n-9]+ 0.003809523809524*u[n-10]-0.000317460317460*u[n-11]); # # Remaining interior points (x=xl+5*dx,...,x=xu-5*dx) for(i in 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]+ 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 ]); # # All points concluded (x=xl,...,x=xu) return(c(uxx)); }