Clear["Global`*"];

Comparison of FDMs with spectral methods - derivatives of a function

Periodic functions

Test function and its derivatives

f[x_]=Exp[Sin[x]]
df[x_]=D[f[x],x]
ddf[x_]=D[f[x],x,x]
Plot[{f[x],df[x],ddf[x]},{x,-3π,3π},PlotLegends->{"f(x)","f'(x)","f''(x)"},AxesLabel->{"x","f"}]

Differentiation matrix for the finite-difference method

Clear[h];
getSparseFDMatrix[n_,p_,h_]:=Module[{v,pf2,matrix},
v=ConstantArray[0,n];
pf2=Factorial[p]^2;
Do[
v[[s+1]]=(-1)^(s+1)*pf2/(s*Factorial[p-s]*Factorial[p+s]);
v[[n-s+1]]=-v[[s+1]],
{s,1,p}
];
matrix=SparseArray[Flatten[{Table[Band[{1,i}]->v[[i]],{i,1,n}],
Table[Band[{i,1}]->v[[n+2-i]],{i,2,n}]}],{n,n}];
Return[matrix/h]];
n=10;
getSparseFDMatrix[n,3,h]//MatrixForm

Differentiation matrix of the spectral method

Formula for even number of grid points

Clear[h];
tmp=Flatten[{0,Table[(-1)^i*Cot[i*h/2]/2,{i,1,n-1}]}];
DM=ToeplitzMatrix[tmp,-tmp];
DM//MatrixForm

Formula for odd number of grid points

Clear[h];
tmp=Flatten[{0,Table[(-1)^i*Csc[i*h/2]/2,{i,1,n}]}];
DM=ToeplitzMatrix[tmp,-tmp];
DM//MatrixForm

Numerical derivative using differentiation matrices

ns={15,16,19,20,23,24,28,32,40,48,56,64,72,80,96,112,128};(* number of grid points *)
nm=6; (* number of methods *)
ng=Length[ns]; (* number of grids *)
errorD1=ConstantArray[0.0,{nm,ng}];
errorD2=ConstantArray[0.0,{nm,ng}];
Do[ (* loop over grids *)
n=ns[[i]];
h=N[2π/n];
xs=-π+h*Range[n]; (* notice that point -π is not icluded *)
fxs=N[f[xs]];
eDfxs=N[df[xs]];
eD2fxs=N[ddf[xs]];
Do[ (* loop over methods - the first is spectral method, then FDMs of increasing even order *)
If[m==1, (* spectral method *)
If[Mod[n,2]==0,
tmp=Flatten[{0,Table[0.5*(-1)^i*Cot[i*h/2],{i,1,n-1}]}], (* even number of points *)
tmp=Flatten[{0,Table[0.5*(-1)^i*Csc[i*h/2],{i,1,n-1}]}] (* odd number of points *)
];
DM=ToeplitzMatrix[tmp,-tmp],
(* finite-difference methods *)
DM=getSparseFDMatrix[n,m-1,h]
];
Dfxs=DM.fxs;
D2fxs=DM.Dfxs;
errorD1[[m,i]]=Max[Abs[Dfxs-eDfxs]];
errorD2[[m,i]]=Max[Abs[D2fxs-eD2fxs]],
{m,1,nm}
],
{i,1,Length[ns]}
];

Error of the first derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD1[[i,All]]}],{i,1,nm}],PlotLegends->{"spectral method","FDM order 2","order 4","order 6","order 8","order 10"},AxesLabel->{"n","error"}]

Error of the second derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD2[[i,All]]}],{i,1,nm}],PlotLegends->{"spectral method","FDM order 2","order 4","order 6","order 8","order 10"},AxesLabel->{"n","error"}]

Non-Periodic functions

Problem

Calculate approximate derivatives of a function

f(x)=(e^(-x^2))sin(3x)

using differentiation matrices and study convergence with respect to space grid parameter h.

Test function and its derivatives

f[x_]=Exp[-x^2]Sin[3x]
df[x_]=D[f[x],x]
ddf[x_]=D[f[x],x,x]
Plot[{f[x],df[x],ddf[x]},{x,-2π,2π},PlotRange->All,PlotLegends->{"f(x)","f'(x)","f''(x)"},AxesLabel->{"x","f"}]

Differentiation matrix for the finite-difference method

Clear[h];
getSparseFDBandMatrix[n_,p_,h_]:=Module[{v,pf2,matrix},
v=ConstantArray[0,n];
pf2=Factorial[p]^2;
Do[
v[[s+1]]=(-1)^(s+1)*pf2/(s*Factorial[p-s]*Factorial[p+s]),
{s,1,p}
];
matrix=SparseArray[Flatten[{Table[Band[{1,i}]->v[[i]],{i,1,n}],
Table[Band[{i,1}]->-v[[i]],{i,2,n}]}],{n,n}];
Return[matrix/h]];
n=8;
getSparseFDBandMatrix[n,3,h]//MatrixForm

Differentiation matrix of the spectral method

getSpectralDiffMatrix[x_]:=Module[{dist,matrix,n},
n=Length[x];
dist=ConstantArray[0,{n,n}];
matrix=ConstantArray[0,{n,n}];
(* precompute distances of points *)
Do[Do[ dist[[i,j]]=x[[i]]-x[[j]],{j,1,n}],{i,1,n}];
(* diagonal elements *)
Do[ matrix[[j,j]]=Sum[If[k==j,0,1/dist[[j,k]]],{k,1,n}],{j,1,n}];
(* off-diagonal elements *)
Do[ 
If[i!=j,
matrix[[i,j]]=Product[If[k!=j&&k!=i ,dist[[i,k]]/dist[[j,k]],1],{k,1,n}]/dist[[j,i]]
],
{j,1,n},{i,1,n}
];
Return[matrix]
];
getSpectralDiffMatrix[{1,1/2,0,-1/2,-1}]//MatrixForm

Numerical derivative using differentiation matrices on the equidistant grid 
=> large error similar to error of an interpolating polynomial of high order on the equidistant grid

ns={8,10,12,14,16,20,24,32,40,48,56};(* number of grid points *)
nm=4; (* number of methods *)
ng=Length[ns]; (* number of grids *)
errorD1=ConstantArray[0.0,{nm,ng}];
errorD2=ConstantArray[0.0,{nm,ng}];
xmin=-2π; xmax=-xmin;
Do[ (* loop over grids *)
n=ns[[i]];
h=N[(xmax-xmin)/n];
xs=xmin+h*Range[n];
fxs=N[f[xs]];
eDfxs=N[df[xs]];
eD2fxs=N[ddf[xs]];
Do[ (* loop over methods - the first is spectral method, then FDMs of increasing even order *)
If[m==1, (* spectral method *)
DM=getSpectralDiffMatrix[xs],
(* finite-difference methods *)
DM=getSparseFDBandMatrix[n,m-1,h]
];
Dfxs=DM.fxs;
D2fxs=DM.Dfxs;
errorD1[[m,i]]=Max[Abs[Dfxs-eDfxs]];
errorD2[[m,i]]=Max[Abs[D2fxs-eD2fxs]],
{m,1,nm}
],
{i,1,Length[ns]}
];

Error of the first derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD1[[i,All]]}],{i,1,nm}],PlotRange->All,PlotLegends->{"spectral method","FDM order 2","order 4","order 6"},AxesLabel->{"n","error"}]

Error of the second derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD2[[i,All]]}],{i,1,nm}],PlotRange->All,PlotLegends->{"spectral method","FDM order 2","order 4","order 6"},AxesLabel->{"n","error"}]

Numerical derivative using differentiation matrices for FDM on the equidistant grid and for the Chebyshev spectral method
=> works because an interpolating polynomial on the Chebyshev grid also works

ns={16,20,24,32,40,48,56,64,72,80,96,112,128};(* number of grid points *)
nm=4; (* number of methods *)
ng=Length[ns]; (* number of grids *)
errorD1=ConstantArray[0.0,{nm,ng}];
errorD2=ConstantArray[0.0,{nm,ng}];
xmin=-π; xmax=-xmin;
Do[ (* loop over methods - the first is spectral method, then FDMs of increasing even order *)
Do[ (* loop over grids *)
n=ns[[i]];
If[m==1,
(* Chebyshev extreme grid *)
xs=N[((xmin+xmax)-(xmax-xmin)*Cos[((Range[n]-1)*π/(n-1))])/2],
(* equidistant grid *)
h=N[(xmax-xmin)/n];
xs=xmin+h*Range[n]
];
fxs=N[f[xs]];
eDfxs=N[df[xs]];
eD2fxs=N[ddf[xs]];
If[m==1, (* spectral method *)
DM=getSpectralDiffMatrix[xs],
(* finite-difference methods *)
DM=getSparseFDBandMatrix[n,m-1,h]
];
Dfxs=DM.fxs;
D2fxs=DM.Dfxs;
errorD1[[m,i]]=Max[Abs[Dfxs[[6;;n-6]]-eDfxs[[6;;n-6]]]];
errorD2[[m,i]]=Max[Abs[D2fxs[[6;;n-6]]-eD2fxs[[6;;n-6]]]],
{i,1,Length[ns]}
],
{m,1,nm}
];

Error of the first derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD1[[i,All]]}],{i,1,nm}],PlotLegends->{"Chebyshev spectral method","FDM order 2","order 4","order 6"},AxesLabel->{"n","error"}]

Error of the second derivative

hs=N[2π/ns];
ListLogLogPlot[Table[Transpose[{ns,errorD2[[i,All]]}],{i,1,nm}],PlotLegends->{"Chebyshev spectral method","FDM order 2","order 4","order 6"},AxesLabel->{"n","error"}]