Nearest neighbour spacing distribution in the GUE
This is an unevaluated Mathematica notebook. If you do not
have a copy of Mathematica handy, you may want to also view the
evaluated version
Power series expansion of the d.e. about the origin
In[]:=
Power series expansion of nn(s) about s=0
In[]:=
py[s_] := Simplify[ya[s]/s]
pin[s_] := 2*Integrate[py[2*t],{t,0,s}]
pnn1[s_] := Series[Exp[pin[Pi*s]],{s,0,10}]
pnn2[s_] := -2*Pi*py[2*Pi*s]
pnn[s_] = Normal[Series[pnn1[s]*pnn2[s],{s,0,9}]]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Numerical solution of the d.e., calculation
of the mean and comparison with power series
expansion
Here the power series solution at s=1.1 is used as the
initial condition
In[]:=
de=2*((1 + y[s] - s*Derivative[1][y][s])*
(Derivative[1][y][s]^2 -
(1 - (1 + y[s] - s*Derivative[1][y][s])^(1/2))^2))^(1/2) +
s*Derivative[2][y][s];
sol = NDSolve[{de == 0, y[1.1]==N[ya[1.1]],
y'[1.1]==N[ya'[1.1]]},y,{s,1.1,16},AccuracyGoal -> 10,
PrecisionGoal -> 10, WorkingPrecision -> 20, MaxSteps -> 2000]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Here the numerical solution is integrated to
determine the mean of the distrubution
In[]:=
ty[s_] = y[s] /. sol; tty[s_]= ty[s][[1]];
vy[s_] := tty[s]/s /; 1 <= s <=13
vy[s_] := Simplify[ya[s]/s] /; 0 <= s < 1
in[s_] := 2*NIntegrate[vy[2*t],{t,0,s}]
mean=NIntegrate[Exp[in[s]],{s,0,2*N[Pi]}]/N[Pi]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Next the numerical solution is compared to the power series solution
In[]:=
ty[s_] = y[s] /. sol; tty[s_]= ty[s][[1]];
vy[s_] := tty[s]/s /; 1 <= s <=13
vy[s_] := Simplify[ya[s]/s] /; 0 <= s < 1
in[s_] := 2*NIntegrate[vy[2*t],{t,0,s}]
nn[s_] := -2*N[Pi]*vy[2*N[Pi]*s]*Exp[in[Pi*s]]
Plot[{nn[s],pnn[s]},{s,.001,.5},PlotStyle->{
GrayLevel[0],GrayLevel[.6]}]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
The function nn(s) can be tabulated
In[]:=
ty[s_] = y[s] /. sol; tty[s_]= ty[s][[1]];
vy[s_] := tty[s]/s /; 1.1 <= s <=16
vy[s_] := N[Simplify[ya[s]/s]] /; 0 <= s < 1.1
in[s_] := 2*NIntegrate[vy[2*t],{t,0,s}]
nn[s_] := -2*N[Pi]*vy[2*N[Pi]*s]*Exp[in[Pi*s]]
t = Table[{j*.01,nn[j*.01]},{j,1,170}];
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Comparison with data from the GUE
Code to generate GUE matrices and plot empirical value of nn(s) (output
suppressed)
In[]:=
<<Statistics`ContinuousDistributions`
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
upperdiag[n_] :=
Table[If[i < j, Random[ndis2] + I*Random[ndis2], 0], {i, n}, {j, n}]
ndis2 = NormalDistribution[0, 1/2];
diagmat[n_] := DiagonalMatrix[Table[Random[ndis1], {n}]]
ndis1 = NormalDistribution[0, 2^(-1/2)];
herm[n_] := (s = upperdiag[n]; diagmat[n] + s + Conjugate[Transpose[s]])
evlist[n_, g_] :=
(ev = {}; Do[evall = Sort[Re[Eigenvalues[herm[n]]]];
a = evall[[(n + 1)/2]]; b = evall[[(n + 1)/2 + 1]];
c = evall[[(n + 1)/2 - 1]];
If[a - c < b - a, AppendTo[ev, a - c], AppendTo[ev, b - a]], {g}];
Flatten[ev])
nev[n_, g_, int_] := Round[N[(int*evlist[n, g])/(Pi/Sqrt[2*n])]]/int
frlis2[n_, g_, int_] :=
(qq = nev[n, g, int]; Table[{{N[(i-1/2 )/int],
N[(int*Count[qq, i/int])/g]}, {N[(i+1/2)/int], N[(int*Count[qq, i/int])/g]}}
, {i, 1, Round[1.7*int]}])
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
barplot[n_, g_, int_] :=
(lis = Flatten[frlis2[n, g, int], 1]; Graphics[{Line[lis],AbsoluteThickness[.5]}])
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Calculation of theoretical curve for nn(s) and comparison with numerical data
In[]:=
nnplot = ListPlot[t,PlotJoined -> True,PlotStyle -> AbsoluteThickness[.5]]
double[n_,g_,int_] := (a1=nnplot;a2=barplot[n,g,int];Show[a1,a2,
AspectRatio -> 1,AxesLabel -> {"t","nn(t)"}])
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
double[15,10000,20]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Comparison with data from zeros of Riemann
zeta function for 10^6 zeros about the 10^20
zero, 10^6 zero and first zero
Here we read in the data as a list. The data gives the number
of zeros in the intervals [j*.05, (j+1)*.05) (j=0,...40). The third,
fourth and fifth columns contain the data for zeros about
zero 1, zero 10^6 and zero 10^20 respectively. The theoretical
value of nn(s) is plotted from the Table t.
In[]:=
data1 = ReadList["data.txt",Number];
num106 = Table[ Text["*",{(j-.5)*.05,20*data1[[5*(j-1)+4]]/1000000}],{j,1,34}];
num1 = Table[ Text["o",{(j-.5)*.05,20*data1[[5*(j-1)+3]]/1000000}],{j,1,34}];
num1020= Table[{(j-.5)*.05,20*data1[[5*(j-1)+5]]/1000000},{j,1,34}];
g1020 = ListPlot[num1020, PlotStyle -> AbsolutePointSize[4]];
g106 = Show[ Graphics [num106]];
g1 = Show[ Graphics [num1]];
ggue = ListPlot[t,PlotJoined -> True,PlotStyle -> AbsoluteThickness[.5]];
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
In[]:=
dplot = (a1=ggue;a2=g1;a3=g106;a4=g1020;Show[a1,a2,a3,a4,
AxesLabel -> {"t","nn(t)"},PlotRange -> {0,1.7},AspectRatio -> 1])
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above