Spectral method code

LGY, 2017

Intro

Each spectral setup module takes as input the desired interval and total number of grid points, and returns a list containing the grid {(Subscript[x, i])}, a list {(Subscript[w, i])} of quadrature integration weights, first and second order differentiation matrices (d^( ))and dd, and a vector-valued cardinal function χ:

(Subscript[f, i]):= f((Subscript[x, i]))
		f'((Subscript[x, i]))≈(∑_(j=1)^n)(Subscript[d, ij])(Subscript[f, j])
		f''((Subscript[x, i]))≈(∑_(j=1)^n)(Subscript[dd, ij])(Subscript[f, j])
		f(x) ≈ (∑_(i=1)^n)(Subscript[f, i])(Subscript[χ(x), i])
		(∫_(Subscript[x, min])^(Subscript[x, max]))dx w(x) f(x) ≈ (∑_(i=1)^n)(Subscript[w, i])(Subscript[f, i])

Modules for finite, non-periodic intervals take an additional endpoint parameter which is list of two logical values: 
	{True,True}:	 include both endpoints in grid (Gauss-Lobatto), 
	{True,False}:	 include (Subscript[x, min]) but not (Subscript[x, max])  in grid (Gauss-Radau),
	{False,True}:	 include (Subscript[x, max]) but not (Subscript[x, min]) in grid (Gauss-Radau),
	{False,False}:	include neither endpoint  in grid (Gauss).

Modules for periodic or semi-infinite intervals take an additional endpoint parameter which is a single logical value:
	True:	  include endpoint in grid,
	False:	 do not include endpoint in grid.

Modules which must solve for grid points numerically take a final optional parameter specifying the working precision (default, machine precision).  Interval endpoints may also be passed 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.  Cardinal functions are carefully coded in forms which remain accurate (avoid nasty precision loss) when evaluated very close to grid points.

Refs: 
	Jean-Paul Berrut & Lloyd N. Trefethen, Barycentric Lagrange Interpolation, SIAM review 46 (2004) 501-517
	Walter Gautschi, Gauss-Radau formulae for Jacobi and Laguerre weight functions, MathCompSim 54 (2000) 403-412
	Z.S. Zheng Guanghui Huang, Computation of Jacobi-Gauss Lobatto Quadrature Based on Derivative Relation

SetOptions[Notebooks["SpectralCode2.nb"][[1]],AutoGeneratedPackage->Automatic]

Barycentric Lagrange interpolation

SpectralInterp[grid_]:=
Module[{np=Length[grid],p,d,dd,dx,card},
p=1/Array[Times@@(grid[[#]]-Drop[grid,{#}])&,np];
card[x_]:=Module[{tmp},
Quiet[tmp=p/(x-grid);
tmp/=Plus@@tmp;
tmp/.Indeterminate->1]];
dx=IdentityMatrix[np];
dx-=1/(dx+Outer[Subtract,grid,grid]);
d=dx Outer[Divide,p,p];
dd=-dx+Plus@@d;
d=Transpose[d];
dd*=-2d;
d-=DiagonalMatrix[Plus@@@d];
dd-=DiagonalMatrix[Plus@@@dd];
{card,d,dd}]

Periodic interval:

Fourier:			w(x)=1

SpectralFourier[xmin_,xmax_,n_?IntegerQ,endpt_?BooleanQ]:=
Module[{m,grid,w,d,dd,card},
m=(xmax-xmin)/(2π);
w=ConstantArray[(xmax-xmin)/n,n];
If[endpt,
	grid=2π#/n&/@Range[0,n-1];
	card=func$[x,Sinc[n/2(x-grid[[#]])]/Sinc[1/2(x-grid[[#]])]
	If[EvenQ[n],Cos[1/2(x-grid[[#]])],1]&/@Range[1,n]]
		/.func$->Function;
	d=Table[If[i==j,
			0,
			1/2(((-1))^(i+j))If[EvenQ[n],Cot,Csc][(i-j)π/n]],
		{i,1,n},{j,1,n}];
	dd=Table[If[i== j,
			-1/12If[EvenQ[n],(n^2)+2,(n^2)-1],
			-1/2 (((-1))^(i+j))Csc[(i-j)π/n]If[EvenQ[n],Csc,Cot][(i-j)π/n]],
		{i,1,n},{j,1,n}];,
	grid=(2#-1)π/n&/@Range[1,n];
	card=func$[x,Sinc[n/2(x-grid[[#]])]/Sinc[1/2(x-grid[[#]])]
	If[EvenQ[n],Cos[1/2(x-grid[[#]])],1]&/@Range[1,n]]
		/.func$->Function;
	d=Table[If[i==j,
			0,
			1/2(((-1))^(i+j))If[EvenQ[n],Cot,Csc][(i-j)π/n]],
		{i,1,n},{j,1,n}];
	dd=Table[If[i== j,
			-1/12If[EvenQ[n],(n^2)+2,(n^2)-1],
			-1/2 (((-1))^(i+j))Csc[(i-j)π/n]If[EvenQ[n],Csc,Cot][(i-j)π/n]],
			{i,1,n},{j,1,n}];];
{m grid+xmin,w,d/m,dd/(m^2),card[(#-xmin)/m]&}];

Finite interval:

Chebyshev:		w(x)= ((((Subscript[x, max])-x))^(-1/2))(((x-(Subscript[x, min])))^(-1/2))

SpectralChebyshev[xmin_,xmax_,n_?IntegerQ,endpts_]:=
Module[{nm1=n-1,b,m,w,z,grid,d,dd,card},
b=(xmax+xmin)/2;
m=(xmin-xmax)/2;
Which[endpts[[1]]&&endpts[[2]],
	grid=Array[Cos[# π/nm1]&,n,0];
	w=π/nm1 Prepend[1/2]@Append[1/2]@ConstantArray[1,n-2],
endpts[[1]]&&!endpts[[2]],
	grid=Array[Cos[2# π/(2n-1)]&,n,0];
	w=π/(2n-1)Prepend[1]@ConstantArray[2,nm1],
!endpts[[1]]&&endpts[[2]],
	grid=Array[Cos[(2# +1)π/(2n-1)]&,n,0];
	w=π/(2n-1)Append[1]@ConstantArray[2,nm1],
!endpts[[1]]&&!endpts[[2]],
	grid=Array[Cos[(2# +1) π/(2n)]&,n,0];
	w=ConstantArray[π/n,n]];
{card,d,dd}=SpectralInterp[grid];
{m grid+b,w,d/m,dd/(m^2),card[(#-b)/m]&}];

Legendre:		w(x)= 1

SpectralGegenbauer[α_][xmin_,xmax_,n_?IntegerQ,endpts_,prec_:MachinePrecision]:=
SpectralJacobi[α-1/2,α-1/2][xmin,xmax,n,endpts,prec];

Gegenbauer:		w(x)= ((((Subscript[x, max])-x))^α)(((x-(Subscript[x, min])))^α)

SpectralGegenbauer[α_][xmin_,xmax_,n_?IntegerQ,endpts_,prec_:MachinePrecision]:=
SpectralJacobi[α-1/2,α-1/2][xmin,xmax,n,endpts,prec];

Jacobi:			w(x)=((((Subscript[x, max])-x))^α)(((x-(Subscript[x, min])))^β)

SpectralJacobi[α_,β_][xmin_,xmax_,n_?IntegerQ,endpts_,prec_:MachinePrecision]:=
Module[{np1=n+1,nm1=n-1,b,m,w,z,grid,d,dd,card},
b=(xmax+xmin)/2;
m=(xmax-xmin)/2;
Which[endpts[[1]]&&endpts[[2]],
	grid=z/.NSolve[JacobiP[n-2,α+1,β+1,z]==0,z,Reals,WorkingPrecision->prec];
	w=((JacobiP[nm1,α,β,#])^(-2))&/@grid;
	w*=(2^(α+β+1))Gamma[n+α]Gamma[n+β]/Gamma[n]/Gamma[np1+α+β]/nm1;
	w=Prepend[(2^(α+β+1))Gamma[α+2]Gamma[β+1]/Gamma[α+β+3]
	Binomial[nm1+α,n-2]/Binomial[nm1+β,n-2]/Binomial[n+α+β,n-2]][w];
	w=Append[(2^(α+β+1))Gamma[β+2]Gamma[α+1]/Gamma[α+β+3]
	Binomial[nm1+β,n-2]/Binomial[nm1+α,n-2]/Binomial[n+α+β,n-2]][w];
	grid=Prepend[-1]@Append[1]@grid,
endpts[[1]]&&!endpts[[2]],
	grid=z/.NSolve[JacobiP[nm1,α,β+1,z]==0,z,Reals,WorkingPrecision->prec];
	w=(1-#)((JacobiP[nm1,α,β,#])^(-2))&/@grid;
	w*=(2^(α+β))Gamma[n+α]Gamma[n+β]/Gamma[n]/Gamma[np1+α+β]/(n+β);
	w=Prepend[(2^(α+β+1))Gamma[β+1]Gamma[n+α]/Gamma[np1+α+β]/Binomial[n+β,nm1]][w];
	grid=Prepend[-1]@grid,
!endpts[[1]]&&endpts[[2]],
	grid=z/.NSolve[JacobiP[n-1,α+1,β,z]==0,z,Reals,WorkingPrecision->prec];
	w=(1+#)((JacobiP[nm1,α,β,#])^(-2))&/@grid;
	w*=(2^(α+β))Gamma[n+β]Gamma[n+α]/Gamma[n]/Gamma[np1+α+β]/(n+α);
	w=Append[(2^(α+β+1))Gamma[α+1]Gamma[n+β]/Gamma[np1+α+β]/Binomial[n+α,nm1]][w];
	grid=Append[1]@grid,
!endpts[[1]]&&!endpts[[2]],
	grid=z/.NSolve[JacobiP[n,α,β,z]==0,z,Reals,WorkingPrecision->prec];
	w=(((JacobiP[nm1,α+1,β+1,#] JacobiP[np1,α,β,#]))^(-1))&/@grid;
	w*=-(2^(α+β+1)) (2 np1+α+β)(((np1+α+β))^(-2)) 
	Gamma[np1+α] Gamma[np1+β]/Gamma[n+2]/Gamma[np1+α+β]];
{card,d,dd}=SpectralInterp[grid];
{m grid+b,w (m^(1+α+β)),d/m,dd/(m^2),card[(#-b)/m]&}];

Semi-infinite interval:

(Associated) Laguerre:	w(x) (e^(-μ x))(x^α)

SpectralLaguerre[α_][μ_,np_?IntegerQ,endpt_?BooleanQ,prec_:MachinePrecision]:=
Module[{npp1=np+1,npm1=np-1,z,w,grid,d,dd,card},
If[endpt,
grid=z/.NSolve[LaguerreL[npm1,α+1,z]==0,z,Reals,WorkingPrecision->prec];
w=((LaguerreL[npm1,α,#])^(-2))&/@grid;
w*=Gamma[α+1]Binomial[npm1+α,npm1]/(np+α);
w=Prepend[Gamma[α+1]/Binomial[np+α,npm1]][w];
grid=Prepend[0][grid];,
grid=z/.NSolve[LaguerreL[np,α,z]==0,z,Reals,WorkingPrecision->prec];
w=# (((npp1 LaguerreL[npp1,α,#]))^(-2))&/@grid;
w*=Gamma[npp1+α]/(np!);];
{card,d,dd}=SpectralInterp[grid];
{grid/μ,w (μ^(-α-1)),d μ,dd (μ^2),card[μ #]&}]

Infinite interval:

Hermite:		 w(x)=(e^(-(μ^2) (x^2)))

SpectralHermite[μ_,np_?IntegerQ,prec_:MachinePrecision]:=
Module[{npm1=np-1,z,w,grid,d,dd,card},
grid=z/.NSolve[HermiteH[np,z]==0,z,Reals,WorkingPrecision->prec];
w=((HermiteH[npm1,#])^(-2))&/@grid;
w*=(2^npm1)npm1!Sqrt[π]/np;
{card,d,dd}=SpectralInterp[grid];
{grid/μ,w (μ^(-1)),d μ,dd (μ^2),card[μ #]&}]

Multiple identical subdomains

MultiDomain[xmin_,xmax_,nd_,np_,approx_]:=
Module[{i,x,xp,grid,w,d,dd,card,subcard,dx=(xmax-xmin)/nd},
xp=ConstantArray[0,{nd,np}];
w=ConstantArray[0,{nd,np}];
subcard=ConstantArray[0,nd];
For[i=1;x=xmin,i≤nd,i++;x+=dx,
{xp[[i]],w[[i]]d,dd,subcard[[i]]}=approx[x,x+dx,np,{True,True}]];
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],Flatten[w],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}]];

Consistency tests

fouriertest[f_,xmin_,xmax_,np_,endpt_]:=Module[{},
{grid,w,d,dd,card}=SpectralFourier[xmin,xmax,np,endpt];
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,Integrate[f[x],{x,xmin,xmax}]-fv.w}/.x->RandomReal[{xmin,xmax}]];

chebyshevtest[f_,xmin_,xmax_,np_,endpts_]:=Module[{},
{grid,w,d,dd,card}=SpectralChebyshev[xmin,xmax,np,endpts];
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,Integrate[(((xmax-x))^(-1/2))(((x-xmin))^(-1/2))f[x]//Simplify,{x,xmin,xmax}]-fv.w}/.x->RandomReal[{xmin,xmax}]];

jacobitest[f_,α_,β_,xmin_,xmax_,np_,endpts_]:=Module[{},
{grid,w,d,dd,card}=SpectralJacobi[α,β][xmin,xmax,np,endpts];
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,Integrate[(((xmax-x))^α)(((x-xmin))^β)f[x],{x,xmin,xmax}]-fv.w}/.x->RandomReal[{xmin,xmax}]];

laguerretest[f_,μ_,α_,np_,endpt_]:=Module[{},
{grid,w,d,dd,card}=SpectralLaguerre[α][μ,np,endpt];
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,Integrate[(E^(-μ x))(x^α)f[x],{x,0,Infinity}]-fv.w}/.x->RandomReal[2/μ]];

hermitetest[f_,μ_,np_]:=Module[{},
{grid,w,d,dd,card}=SpectralHermite[μ,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,Integrate[(E^(-(μ^2) (x^2)))f[x],{x,-Infinity,Infinity}]-fv.w}/.x->RandomReal[2/μ]];

fouriertest[Cos[2π#]&,0,2,5,#]&/@{True,False}//FullSimplify//TableForm

chebyshevtest[((#)^8)&,0,1.0,11,#]&/@{{True,True},{True,False},{False,True},{False,False}}//TableForm

jacobitest[((#)^3)&,3/2,2,0,1,21,#]&/@{{True,True},{True,False},{False,True},{False,False}}//TableForm

laguerretest[((#)^5)&,2,3,11,#]&/@{True,False}//TableForm

hermitetest[((#)^4)&,1/2,11]