(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
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}]


(* ::Input::Initialization:: *)
SpectralFourier[xmin_,xmax_,n_?IntegerQ,endpt_?BooleanQ]:=
Module[{m,grid,w,d,dd,card},
m=(xmax-xmin)/(2\[Pi]);
w=ConstantArray[(xmax-xmin)/n,n];
If[endpt,
	grid=2\[Pi] #/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)\[Pi]/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)\[Pi]/n]If[EvenQ[n],Csc,Cot][(i-j)\[Pi]/n]],
		{i,1,n},{j,1,n}];,
	grid=(2#-1)\[Pi]/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)\[Pi]/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)\[Pi]/n]If[EvenQ[n],Csc,Cot][(i-j)\[Pi]/n]],
			{i,1,n},{j,1,n}];];
{m grid+xmin,w,d/m,dd/m^2,card[(#-xmin)/m]&}];


(* ::Input::Initialization:: *)
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[# \[Pi]/nm1]&,n,0];
	w=\[Pi]/nm1 Prepend[1/2]@Append[1/2]@ConstantArray[1,n-2],
endpts[[1]]&&!endpts[[2]],
	grid=Array[Cos[2# \[Pi]/(2n-1)]&,n,0];
	w=\[Pi]/(2n-1)Prepend[1]@ConstantArray[2,nm1],
!endpts[[1]]&&endpts[[2]],
	grid=Array[Cos[(2# +1)\[Pi]/(2n-1)]&,n,0];
	w=\[Pi]/(2n-1)Append[1]@ConstantArray[2,nm1],
!endpts[[1]]&&!endpts[[2]],
	grid=Array[Cos[(2# +1) \[Pi]/(2n)]&,n,0];
	w=ConstantArray[\[Pi]/n,n]];
{card,d,dd}=SpectralInterp[grid];
{m grid+b,w,d/m,dd/m^2,card[(#-b)/m]&}];


(* ::Input::Initialization:: *)
SpectralGegenbauer[\[Alpha]_][xmin_,xmax_,n_?IntegerQ,endpts_,prec_:MachinePrecision]:=
SpectralJacobi[\[Alpha]-1/2,\[Alpha]-1/2][xmin,xmax,n,endpts,prec];


(* ::Input::Initialization:: *)
SpectralGegenbauer[\[Alpha]_][xmin_,xmax_,n_?IntegerQ,endpts_,prec_:MachinePrecision]:=
SpectralJacobi[\[Alpha]-1/2,\[Alpha]-1/2][xmin,xmax,n,endpts,prec];


(* ::Input::Initialization:: *)
SpectralJacobi[\[Alpha]_,\[Beta]_][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,\[Alpha]+1,\[Beta]+1,z]==0,z,Reals,WorkingPrecision->prec];
	w=JacobiP[nm1,\[Alpha],\[Beta],#]^-2&/@grid;
	w*=2^(\[Alpha]+\[Beta]+1) Gamma[n+\[Alpha]]Gamma[n+\[Beta]]/Gamma[n]/Gamma[np1+\[Alpha]+\[Beta]]/nm1;
	w=Prepend[2^(\[Alpha]+\[Beta]+1) Gamma[\[Alpha]+2]Gamma[\[Beta]+1]/Gamma[\[Alpha]+\[Beta]+3]
	Binomial[nm1+\[Alpha],n-2]/Binomial[nm1+\[Beta],n-2]/Binomial[n+\[Alpha]+\[Beta],n-2]][w];
	w=Append[2^(\[Alpha]+\[Beta]+1) Gamma[\[Beta]+2]Gamma[\[Alpha]+1]/Gamma[\[Alpha]+\[Beta]+3]
	Binomial[nm1+\[Beta],n-2]/Binomial[nm1+\[Alpha],n-2]/Binomial[n+\[Alpha]+\[Beta],n-2]][w];
	grid=Prepend[-1]@Append[1]@grid,
endpts[[1]]&&!endpts[[2]],
	grid=z/.NSolve[JacobiP[nm1,\[Alpha],\[Beta]+1,z]==0,z,Reals,WorkingPrecision->prec];
	w=(1-#)JacobiP[nm1,\[Alpha],\[Beta],#]^-2&/@grid;
	w*=2^(\[Alpha]+\[Beta]) Gamma[n+\[Alpha]]Gamma[n+\[Beta]]/Gamma[n]/Gamma[np1+\[Alpha]+\[Beta]]/(n+\[Beta]);
	w=Prepend[2^(\[Alpha]+\[Beta]+1) Gamma[\[Beta]+1]Gamma[n+\[Alpha]]/Gamma[np1+\[Alpha]+\[Beta]]/Binomial[n+\[Beta],nm1]][w];
	grid=Prepend[-1]@grid,
!endpts[[1]]&&endpts[[2]],
	grid=z/.NSolve[JacobiP[n-1,\[Alpha]+1,\[Beta],z]==0,z,Reals,WorkingPrecision->prec];
	w=(1+#)JacobiP[nm1,\[Alpha],\[Beta],#]^-2&/@grid;
	w*=2^(\[Alpha]+\[Beta]) Gamma[n+\[Beta]]Gamma[n+\[Alpha]]/Gamma[n]/Gamma[np1+\[Alpha]+\[Beta]]/(n+\[Alpha]);
	w=Append[2^(\[Alpha]+\[Beta]+1) Gamma[\[Alpha]+1]Gamma[n+\[Beta]]/Gamma[np1+\[Alpha]+\[Beta]]/Binomial[n+\[Alpha],nm1]][w];
	grid=Append[1]@grid,
!endpts[[1]]&&!endpts[[2]],
	grid=z/.NSolve[JacobiP[n,\[Alpha],\[Beta],z]==0,z,Reals,WorkingPrecision->prec];
	w=(JacobiP[nm1,\[Alpha]+1,\[Beta]+1,#] JacobiP[np1,\[Alpha],\[Beta],#])^-1&/@grid;
	w*=-2^(\[Alpha]+\[Beta]+1) (2 np1+\[Alpha]+\[Beta])(np1+\[Alpha]+\[Beta])^-2 
	Gamma[np1+\[Alpha]] Gamma[np1+\[Beta]]/Gamma[n+2]/Gamma[np1+\[Alpha]+\[Beta]]];
{card,d,dd}=SpectralInterp[grid];
{m grid+b,w m^(1+\[Alpha]+\[Beta]),d/m,dd/m^2,card[(#-b)/m]&}];


(* ::Input::Initialization:: *)
SpectralLaguerre[\[Alpha]_][\[Mu]_,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,\[Alpha]+1,z]==0,z,Reals,WorkingPrecision->prec];
w=LaguerreL[npm1,\[Alpha],#]^-2&/@grid;
w*=Gamma[\[Alpha]+1]Binomial[npm1+\[Alpha],npm1]/(np+\[Alpha]);
w=Prepend[Gamma[\[Alpha]+1]/Binomial[np+\[Alpha],npm1]][w];
grid=Prepend[0][grid];,
grid=z/.NSolve[LaguerreL[np,\[Alpha],z]==0,z,Reals,WorkingPrecision->prec];
w=# (npp1 LaguerreL[npp1,\[Alpha],#])^-2&/@grid;
w*=Gamma[npp1+\[Alpha]]/(np!);];
{card,d,dd}=SpectralInterp[grid];
{grid/\[Mu],w \[Mu]^(-\[Alpha]-1),d \[Mu],dd \[Mu]^2,card[\[Mu] #]&}]


(* ::Input::Initialization:: *)
SpectralHermite[\[Mu]_,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[\[Pi]]/np;
{card,d,dd}=SpectralInterp[grid];
{grid/\[Mu],w \[Mu]^-1,d \[Mu],dd \[Mu]^2,card[\[Mu] #]&}]


(* ::Input::Initialization:: *)
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}]];


(* ::Input::Initialization:: *)
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_,\[Alpha]_,\[Beta]_,xmin_,xmax_,np_,endpts_]:=Module[{},
{grid,w,d,dd,card}=SpectralJacobi[\[Alpha],\[Beta]][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)^\[Alpha] (x-xmin)^\[Beta] f[x],{x,xmin,xmax}]-fv.w}/.x->RandomReal[{xmin,xmax}]];

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

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



