Nearest neighbour spacing distribution in the GUE

Power series expansion of the d.e. about the        origin

de=-4*((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)) +
  (s*Derivative[2][y][s])^2;
  
  approx=Series[de,{s,0,12}];
  
  equations=LogicalExpand[approx==0];
  
  y'''[0]=-1/(2*Pi); y[0]=0;y'[0]=0;y''[0]=0;
  
  values = Solve[equations];
  
  ya[s_] = Normal[Series[y[s],{s,0,11}] /. values[[1]]]

Power series expansion of nn(s) about s=0

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}]]

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

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]

Here the numerical solution is integrated to
         determine the mean of the distrubution

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]

Next the numerical solution is compared to the power series solution

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]}]

The function nn(s) can be tabulated

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}];

Comparison with data from the GUE

Code to generate GUE matrices  and plot empirical value of nn(s) (output suppressed)

<<Statistics`ContinuousDistributions`

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]}])

barplot[n_, g_, int_] :=                                                                                    
  (lis = Flatten[frlis2[n, g, int], 1]; Graphics[{Line[lis],AbsoluteThickness[.5]}])

Calculation of theoretical curve for nn(s) and comparison with numerical data

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)"}])

double[15,10000,20]

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.

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]];

dplot = (a1=ggue;a2=g1;a3=g106;a4=g1020;Show[a1,a2,a3,a4,
AxesLabel -> {"t","nn(t)"},PlotRange -> {0,1.7},AspectRatio -> 1])