Pólya's Random Walk Constants

Author

Eric W. Weisstein
February 3, 2015

This notebook downloaded from http://mathworld.wolfram.com/notebooks/Statistics/PolyasRandomWalkConstants.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/PolyasRandomWalkConstants.html.

©2015 Wolfram Research, Inc. except for portions noted otherwise

Watson's Triple Integral

Numerical

(3/(((2π))^3))NIntegrate[1/(3-Cos[x]-Cos[y]-Cos[z]),{x,-Pi,0,Pi},{y,-Pi,0,Pi},{z,-Pi,0,Pi}]//Timing

Analytic Forms

(((Sqrt[(3/2)]) Γ((1/24)) Γ((5/24)) Γ((7/24)) Γ((11/24)))/(16 (π^3)))

u3=Sqrt[6]/(32Pi^3)Times@@(Gamma[#/24]&/@{1,5,7,11});

N[u3,20]

N[1-(1/u3),20]

((12 (12 (Sqrt[2])-10 (Sqrt[3])-7 (Sqrt[6])+18) ((K(((((2-(Sqrt[3])) ((Sqrt[3])-(Sqrt[2]))))^2)))^2))/(π^2))

N[(12/(π^2))(18+12Sqrt[2]-10Sqrt[3]-7Sqrt[6])EllipticK[((2-Sqrt[3])(Sqrt[3]-Sqrt[2]))^2]^2,20]

3 (18+12 (Sqrt[2])-10 (Sqrt[3])-7 (Sqrt[6])) (((Subscript[ϑ, 3])(0,(e^(-(Sqrt[6]) π))))^4)

N[3(18+12Sqrt[2]-10Sqrt[3]-7Sqrt[6])(1+2Sum[Exp[-k^2 π Sqrt[6]],{k,∞}])^4,20]

Plots

Plot[Table[BesselI[0,t/d]^d Exp[-t],{d,3}],{t,0,20},PlotRange->{Automatic,{0,1}},Evaluated->True]

Analytic Integrals

u(3)

Integrate[BesselI[0,t/3]^3 Exp[-t],{t,0,Infinity}]

u(4)

Integrate[BesselI[0,t/4]^4 Exp[-t],{t,0,Infinity}]

Numerical Integrals

FindRoot[((BesselI[0,(t/5)])^5) (E^(-t))==10^-12,{t,10^5}]

FindRoot[((BesselI[0,(t/5)])^5) (E^(-t))==10^-13,{t,10^6}]

Polya[n_,opts___]:=1-1/NIntegrate[BesselI[0,t/n]^n Exp[-t],{t,0,UpperLimit/.{opts}/.UpperLimit->10^9},
opts,SingularityDepth->30,WorkingPrecision->50,PrecisionGoal->5,MaxRecursion->50]

Table[{n,Polya[n,10^8,PrecisionGoal->11]},{n,4,8}]//TableForm

Table[{n,Polya[n,PrecisionGoal->11]},{n,4,8}]//TableForm

Try to get lots of digits

Michael Trott
Feb. 3, 2015

With NSum and NIntegrate you will never get 1000 digits.

I think take a series expansion for BesselI at 0. Apply Binomial theorem for the ^4, then integrate. The sum from n0 to infinity symbolically and from 1 to n0 explicitly.

(Subscript[p, 4])

u4=With[{n=4},NIntegrate[((BesselI[0,t/n])^n)Exp[-t],{t,0,∞}]]

1-1/u4

u4=With[{n=4},NIntegrate[((BesselI[0,t/n])^n)Exp[-t],{t,0,∞},WorkingPrecision->20]]

1-1/u4

u4=With[{n=4},NIntegrate[((BesselI[0,t/n])^n)Exp[-t],{t,0,1000},WorkingPrecision->20]]

1-1/u4

u4=With[{n=4},NIntegrate[((BesselI[0,t/n])^n)Exp[-t],{t,0,10,∞},WorkingPrecision->20,MaxRecursion->15,Method->"DoubleExponential"]]

1-1/u4

Series

With[{n=4},Series[BesselI[0,t/n],{t,0,12}]]

SeriesCoefficient[((BesselI[0,t/4])^4),{t,0,n}]

With[{n=4},Series[((BesselI[0,t/n])^n),{t,0,12}]]

ser=Normal[%];

Plot[((BesselI[0,t/4])^4)-ser,{t,0,10}]

SeriesCoefficient[BesselI[0,t/4],{t,0,n}]//FullSimplify

Sloane

.1932016732249839373418709733293691605758733864501395

8

sum=Sum[((8^(-n))/((((n/2)!))^2))(t^n),{n,2,8,2}]

bin=Sum[Binomial[4,k] sum^k,{k,0,4}]

Plot[((BesselI[0,t/4])^4)-bin,{t,0,10}]

(u4=Block[{$MaxExtraPrecision=10^2},NIntegrate[((BesselI[0,t/4])^4)Exp[-t],{t,0,1000},WorkingPrecision->200]+Integrate[bin Exp[-t],{t,1000,∞}]])//Timing

p4=1-1/u4

Precision[p4]