dss046=function(xl,xu,n,u,ux,nl,nu) { # # 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 dss046 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* ( 5.211111111111110*u[ 1]-22.300000000000000*u[ 2]+ 43.950000000000000*u[ 3]-52.722222222222200*u[ 4]+ 41.000000000000000*u[ 5]-20.100000000000000*u[ 6]+ 5.661111111111110*u[ 7] -0.700000000000000*u[ 8]); if(nu==1) uxx[n]=rdxs* ( 5.211111111111110*u[n ]-22.300000000000000*u[n-1]+ 43.950000000000000*u[n-2]-52.722222222222200*u[n-3]+ 41.000000000000000*u[n-4]-20.100000000000000*u[n-5]+ 5.661111111111110*u[n-6] -0.700000000000000*u[n-7]); if(nl==2) uxx[1]=rdxs* ( -7.493888888888860*u[ 1]+12.000000000000000*u[ 2]- 7.499999999999940*u[ 3] +4.444444444444570*u[ 4]- 1.874999999999960*u[ 5] +0.479999999999979*u[ 6]- 0.055555555555568*u[ 7] -4.900000000000000*ux[1]*dx); if(nu==2) uxx[n]=rdxs* ( -7.493888888888860*u[n ]+12.000000000000000*u[n-1]- 7.499999999999940*u[n-2] +4.444444444444570*u[n-3]- 1.874999999999960*u[n-4] +0.479999999999979*u[n-5]- 0.055555555555568*u[n-6] +4.900000000000000*ux[n]*dx); # # dx in from boundaries (x=xl+dx,x=xu-dx) uxx[2]=rdxs* ( 0.700000000000000*u[ 1]-0.388888888888889*u[ 2]- 2.700000000000000*u[ 3]+4.750000000000000*u[ 4]- 3.722222222222220*u[ 5]+1.800000000000000*u[ 6]- 0.500000000000000*u[ 7]+0.061111111111111*u[ 8]); uxx[n-1]=rdxs* ( 0.700000000000000*u[n ]-0.388888888888889*u[n-1]- 2.700000000000000*u[n-2]+4.750000000000000*u[n-3]- 3.722222222222220*u[n-4]+1.800000000000000*u[n-5]- 0.500000000000000*u[n-6]+0.061111111111111*u[n-7]); # # 2*dx in from boundaries (x=xl+2*dx,x=xu-2*dx) uxx[3]=rdxs* ( -0.061111111111111*u[ 1]+1.188888888888890*u[ 2]- 2.100000000000000*u[ 3]+0.722222222222223*u[ 4]+ 0.472222222222222*u[ 5]-0.300000000000000*u[ 6]+ 0.088888888888889*u[ 7]-0.011111111111111*u[ 8]); uxx[n-2]=rdxs* ( -0.061111111111111*u[n ]+1.188888888888890*u[n-1]- 2.100000000000000*u[n-2]+0.722222222222223*u[n-3]+ 0.472222222222222*u[n-4]-0.300000000000000*u[n-5]+ 0.088888888888889*u[n-6]-0.011111111111111*u[n-7]); # # Remaining interior points (x=xl+3*dx,...,x=xu-3*dx) for(i in 4:(n-3)) uxx[i]=rdxs* ( 0.011111111111111*u[i-3]-0.150000000000000*u[i-2]+ 1.500000000000000*u[i-1]+0.011111111111111*u[i+3]- 0.150000000000000*u[i+2]+1.500000000000000*u[i+1]- 2.722222222222220*u[i ]); # # All points concluded (x=xl,...,x=xu) return(c(uxx)); }