Spectral method code

LGY, 2017

Each spectral setup module takes as input the desired interval and total number of grid points, and returns a list containing the grid, first and second order differentiation matrices, and a vector-valued cardinal function.  (The Laguerre module for semi-infinite intervals also returns an integration weight vector.)  Give interval endpoints as machine precision (e.g., 1.0 π instead of just π) or extended precision (e.g., 1`30 π instead of π) to force evaluation using the specified precision.  Most cardinal functions are carefully coded in forms which remain accurate (avoid precision loss) when evaluated very close to grid points.

Periodic interval, even or odd # points,  endpoint grid

SpectralFourierExterior[xmin_,xmax_,np_?IntegerQ]:=Module[{m,s,sp,x,xp,d,dd,card},
m=(xmax-xmin)/2/π;
s=(#-xmin)/m&;
sp=2π#/np&/@Range[0,np-1];
xp=m sp + xmin;
card=
func$[x,Sinc[np/2(s[x]-sp[[#]])]/Sinc[1/2(s[x]-sp[[#]])]
If[EvenQ[np],Cos[1/2(s[x]-sp[[#]])],1]&/@Range[1,np]]/.func$->Function;
d=Table[If[i==j,0,1/2(((-1))^(i+j))If[EvenQ[np],Cot,Csc][(i-j)π/np]],{i,1,np},{j,1,np}];
dd=Table[If[i== j,-1/12If[EvenQ[np],(np^2)+2,(np^2)-1],-1/2 (((-1))^(i+j))Csc[(i-j)π/np]If[EvenQ[np],Csc,Cot][(i-j)π/np]],{i,1,np},{j,1,np}];
{xp,d/m,dd/(m^2),card}];

Periodic interval, even or odd # points, interior grid

SpectralFourierInterior[xmin_,xmax_,np_?IntegerQ]:=Module[{m,s,sp,x,xp,d,dd,card},
m=(xmax-xmin)/2/π;
s=(#-xmin)/m&;
sp=(2#-1)π/np&/@Range[1,np];
xp=m sp + xmin;
card=
func$[x,Sinc[np/2(s[x]-sp[[#]])]/Sinc[1/2(s[x]-sp[[#]])]
If[EvenQ[np],Cos[1/2(s[x]-sp[[#]])],1]&/@Range[1,np]]/.func$->Function;
d=Table[If[i==j,0,1/2(((-1))^(i+j))If[EvenQ[np],Cot,Csc][(i-j)π/np]],{i,1,np},{j,1,np}];
dd=Table[If[i== j,-1/12If[EvenQ[np],(np^2)+2,(np^2)-1],-1/2 (((-1))^(i+j))Csc[(i-j)π/np]If[EvenQ[np],Csc,Cot][(i-j)π/np]],{i,1,np},{j,1,np}];
{xp,d/m,dd/(m^2),card}];

Bounded interval, Chebyshev endpoint (Gauss-Lobatto) grid

SpectralChebyshevExterior[xmin_,xmax_,np_?IntegerQ]:=
Module[{npm1,m,b,t,p,tp,sp,xp,d,card},
npm1=np-1;
b=(xmin+xmax)/2;
m=(xmin-xmax)/2;
t=ArcCos[(#-b)/m]&;
tp=π #/npm1&/@Range[0,npm1];
sp=Cos[tp];
xp=m sp + b;
p= If[#==1||#==np,2,1]&/@Range[np];
card=func$[tx,
((p[[#]])^(-1))Sin[tx]Csc[1/2(tx+tp[[#]])]Sinc[npm1(tx-tp[[#]])]/Sinc[1/2(tx-tp[[#]])]&/@Range[np]]/.Sin[tx]Csc[tx/2]->2Cos[tx/2]/.Sin[tx]Csc[(π+tx)/2]->2Sin[tx/2]/.func$->Function;
d=Table[Which[i≠j,(((-1))^(i+j))p[[i]]/p[[j]]/(sp[[i]]-sp[[j]]),
i==1&&j==1,(1+2(npm1^2))/6,
i==np&&j==np,-(1+2(npm1^2))/6,
True,-sp[[j]]/2/(1-((sp[[j]])^2))],{i,1,np},{j,1,np}];
{xp,d/m,d.d/(m^2),card[t[#]]&}]

Bounded interval, Chebyshev interior (roots) grid

SpectralChebyshevInterior[xmin_,xmax_,np_?IntegerQ]:=
Module[{m,b,t,tp,sp,xp,d,dd,card},
b=(xmin+xmax)/2;
m=(xmin-xmax)/2;
t=ArcCos[(#-b)/m]&;
tp=π(#-1/2)/np&/@Range[np];
sp=Cos[tp];
xp=m sp + b;
card=func$[tx,(((-1))^(#-1))Csc[1/2(tx+tp[[#]])]Sinc[np(tx-tp[[#]])]/Sinc[1/2(tx-tp[[#]])]Csc[np tp[[#]]]Sin[tp[[#]]]&/@Range[np]]/.func$->Function;
d=Table[If[i≠j,(((-1))^(i+j))Sqrt[(1-((sp[[j]])^2))/(1-((sp[[i]])^2))]/(sp[[i]]-sp[[j]]),
sp[[j]]/2/(1-((sp[[j]])^2))],{i,1,np},{j,1,np}];
dd=Table[If[i≠j,d[[i,j]](sp[[i]]/(1-((sp[[i]])^2))-2/(sp[[i]]-sp[[j]])),
((sp[[j]])^2)/(((1-((sp[[j]])^2)))^2)-((np^2)-1)/3/(1-((sp[[j]])^2))],{i,1,np},{j,1,np}];
{xp,d/m,dd/(m^2),card[t[#]]&}]

Bounded interval, Gegenbauer endpoint (Gauss-Lobatto) grid

SpectralGegenbauerExterior[λ_][xmin_,xmax_,np_?IntegerQ,prec_:$MachinePrecision]:=Module[{npm1,m,b,s,p,d,sp,x,xp,card},
npm1=np-1;
b=(xmax+xmin)/2;
m=(xmax-xmin)/2;
s=(#-b)/m&;
sp=z/.NSolve[GegenbauerC[np-2,λ+1,z]==0,z,WorkingPrecision->prec];
sp=Join[{-1},sp,{1}];
xp=m sp + b;
p= If[#==1||#==np,2,1]&/@Range[np];
card=func$[sx,2λ(1-(sx^2))gegen[np-2,λ+1,sx]/(4λ(sx-sp[[#]]))/Which[#==1,gegen[np-2,λ+1,-1],#==np,-gegen[np-2,λ+1,1],True,(1-((sp[[#]])^2))(λ+1)gegen[np-3,λ+2,sp[[#]]]]&/@Range[np]]/.func$->Function/.gegen->GegenbauerC;
d=Table[Which[
i==1&&j==1,-1/2-(np-2) (np+2 λ)/(3+2 λ),
i==np&&j==np,1/2+(np-2) (np+2 λ)/(3+2 λ),
i==1&&j==np,(((-1))^np)/2,
i==np&&j==1,-(((-1))^np)/2,
i==1,-gegen[np-2,λ+1,-1]/gegen[np-3,λ+2,sp[[j]]]/((1-sp[[j]]) (((1+sp[[j]]))^2) (λ+1)),
i==np,-gegen[np-2,λ+1,1]/gegen[np-3,λ+2,sp[[j]]]/((1+sp[[j]])(((1-sp[[j]]))^2)(λ+1)),
j==1,(1-sp[[i]]) (λ+1) gegen[np-3,λ+2,sp[[i]]]/gegen[np-2,λ+1,-1],
j==np,(1+sp[[i]]) (λ+1) gegen[np-3,λ+2,sp[[i]]]/gegen[np-2,λ+1,1],
i≠j,gegen[np-3,λ+2,sp[[i]]]/gegen[np-3,λ+2,sp[[j]]](1-((sp[[i]])^2))/((1-((sp[[j]])^2))(sp[[i]]-sp[[j]])),
True,-2 sp[[j]]/(1-((sp[[j]])^2))+(λ+2) gegen[np-4,λ+3,sp[[j]]]/gegen[np-3,λ+2,sp[[j]]]]/.gegen->GegenbauerC,{i,1,np},{j,1,np}];
{xp,d/m,d.d/(m^2),card[s[#]]&}]

Semi-infinite interval, Laguerre interior grid

SpectralLaguerreInterior[a_][μ_,np_?IntegerQ,prec_:MachinePrecision]:=
Module[{npp1,sp,w,s,d,dd,id,dx,p,card},
npp1=np+1;
sp=z/.NSolve[LaguerreL[np,a,z]==0,z,Reals,WorkingPrecision->prec];
p=Table[1/Times@@Rest@RotateLeft[sp[[j]]-sp,j-1],{j,np}];
w=# (((npp1 LaguerreL[npp1,a,#]))^(-2))&/@sp;
w*=Gamma[np+a+1]/(np!);
s=μ#&;
id=IdentityMatrix[np];
dx=id-1/(id+Outer[Subtract,sp,sp]);
d=dx Outer[Divide,p,p];
dd=-dx+Apply[Plus,d];
d=Transpose[d];
dd*=-2d;
d-=DiagonalMatrix[Apply[Plus,d,{1}]];
dd-=DiagonalMatrix[Apply[Plus,dd,{1}]];
card[x_]:=Module[{j,tmp},
tmp=s[x]-sp;
Quiet[tmp=p/tmp;
tmp/=Plus@@tmp;
tmp/.Indeterminate->1]];
{sp/μ,d μ,dd (μ^2),card,w (μ^(-a-1))}]

Multiple subdomains

MultiDomain[xmin_,xmax_,nd_,np_,approx_]:=
Module[{i,x,xp,grid,d,dd,card,subcard,dx=(xmax-xmin)/nd},
xp=ConstantArray[0,{nd,np}];
subcard=ConstantArray[0,nd];
For[i=1;x=xmin,i≤nd,i++;x+=dx,
{xp[[i]],d,dd,subcard[[i]]}=approx[x,x+dx,np]];
card=Function[x,Module[{c,j},
c=ConstantArray[0,{nd,np},SparseArray];
j=domain[xmin,xmax,nd][x];
c[[j]]=subcard[[j]][x];
c//Flatten]];
{Flatten[xp],replicate[d,nd],replicate[dd,nd],card}];

span=Span[#1(#2-1)+1,#1 #2]&;
domain[xmin_,xmax_,nd_]=Max[Min[Ceiling[nd(#-xmin)/(xmax-xmin)],nd],1]&;
replicate[m_,n_]:=ArrayFlatten[Table[If[i==j,SparseArray[m],0],{i,1,n},{j,1,n}]];

Basic consistency tests

spectraltestbounded[f_,xmin_,xmax_,np_,approx_]:=Module[{},
{grid,d,dd,card}=approx[xmin,xmax,np];
fv=f/@grid;
dfv=d.fv//Simplify;
ddfv=dd.fv//Simplify;
{f[x]==card[x].fv,f'[x]==card[x].dfv,f''[x]==card[x].ddfv}/.Sinc[z_]->Sin[z]/z//FullSimplify];

spectraltestsemiinf[f_,μ_,np_,approx_]:=Module[{},
{grid,d,dd,card,w}=approx[μ,np];
fv=f/@grid;
dfv=d.fv//Simplify;
ddfv=dd.fv//Simplify;
{f[x]==card[x].fv,f'[x]==card[x].dfv,f''[x]==card[x].ddfv}/.x->RandomReal[2/μ]//FullSimplify];

spectraltestbounded[Cos[#]&,-2π,2π,8,SpectralFourierExterior]

spectraltestbounded[Sin[2#]&,-π,π,5,SpectralFourierExterior]

spectraltestbounded[Cos[2#]&,0,π,5,SpectralFourierInterior]

spectraltestbounded[Sin[2#]&,0,-π,8,SpectralFourierInterior]

spectraltestbounded[((#)^3)&,-1,1,5,SpectralChebyshevExterior]

spectraltestbounded[((#)^3)&,-1,1,5,SpectralChebyshevInterior]

spectraltestbounded[((#)^3)&,-1,1,7,SpectralGegenbauerExterior[2]]

spectraltestsemiinf[((#)^3)&,1,7,SpectralLaguerreInterior[1]]