(* ::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:: *)
SpectralFourierExterior[xmin_,xmax_,np_?IntegerQ]:=Module[{m,s,sp,x,xp,d,dd,card},
m=(xmax-xmin)/2/\[Pi];
s=(#-xmin)/m&;
sp=2\[Pi] #/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)\[Pi]/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)\[Pi]/np]If[EvenQ[np],Csc,Cot][(i-j)\[Pi]/np]],{i,1,np},{j,1,np}];
{xp,d/m,dd/m^2,card}];

(* ::Input::Initialization:: *)
SpectralFourierInterior[xmin_,xmax_,np_?IntegerQ]:=Module[{m,s,sp,x,xp,d,dd,card},
m=(xmax-xmin)/2/\[Pi];
s=(#-xmin)/m&;
sp=(2#-1)\[Pi]/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)\[Pi]/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)\[Pi]/np]If[EvenQ[np],Csc,Cot][(i-j)\[Pi]/np]],{i,1,np},{j,1,np}];
{xp,d/m,dd/m^2,card}];


(* ::Input::Initialization:: *)
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=\[Pi] #/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[(\[Pi]+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+2npm1^2)/6,
i==np&&j==np,-(1+2npm1^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[#]]&}]


(* ::Input::Initialization:: *)
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=\[Pi] (#-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[#]]&}]


(* ::Input::Initialization:: *)
SpectralGegenbauerExterior[\[Lambda]_][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,\[Lambda]+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\[Lambda](1-sx^2)gegen[np-2,\[Lambda]+1,sx]/(4\[Lambda](sx-sp[[#]]))/Which[#==1,gegen[np-2,\[Lambda]+1,-1],#==np,-gegen[np-2,\[Lambda]+1,1],True,(1-sp[[#]]^2)(\[Lambda]+1)gegen[np-3,\[Lambda]+2,sp[[#]]]]&/@Range[np]]/.func$->Function/.gegen->GegenbauerC;
d=Table[Which[
i==1&&j==1,-1/2-(np-2) (np+2 \[Lambda])/(3+2 \[Lambda]),
i==np&&j==np,1/2+(np-2) (np+2 \[Lambda])/(3+2 \[Lambda]),
i==1&&j==np,(-1)^np/2,
i==np&&j==1,-(-1)^np/2,
i==1,-gegen[np-2,\[Lambda]+1,-1]/gegen[np-3,\[Lambda]+2,sp[[j]]]/((1-sp[[j]]) (1+sp[[j]])^2 (\[Lambda]+1)),
i==np,-gegen[np-2,\[Lambda]+1,1]/gegen[np-3,\[Lambda]+2,sp[[j]]]/((1+sp[[j]])(1-sp[[j]])^2 (\[Lambda]+1)),
j==1,(1-sp[[i]]) (\[Lambda]+1) gegen[np-3,\[Lambda]+2,sp[[i]]]/gegen[np-2,\[Lambda]+1,-1],
j==np,(1+sp[[i]]) (\[Lambda]+1) gegen[np-3,\[Lambda]+2,sp[[i]]]/gegen[np-2,\[Lambda]+1,1],
i!=j,gegen[np-3,\[Lambda]+2,sp[[i]]]/gegen[np-3,\[Lambda]+2,sp[[j]]](1-sp[[i]]^2)/((1-sp[[j]]^2)(sp[[i]]-sp[[j]])),
True,-2 sp[[j]]/(1-sp[[j]]^2)+(\[Lambda]+2) gegen[np-4,\[Lambda]+3,sp[[j]]]/gegen[np-3,\[Lambda]+2,sp[[j]]]]/.gegen->GegenbauerC,{i,1,np},{j,1,np}];
{xp,d/m,d.d/m^2,card[s[#]]&}]


(* ::Input::Initialization:: *)
SpectralLaguerreInterior[a_][\[Mu]_,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=\[Mu] #&;
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/\[Mu],d \[Mu],dd \[Mu]^2,card,w \[Mu]^(-a-1)}]


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


(* ::Input::Initialization:: *)
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_,\[Mu]_,np_,approx_]:=Module[{},
{grid,d,dd,card,w}=approx[\[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}/.x->RandomReal[2/\[Mu]]//FullSimplify];

