Path: Blog/PUBLIC/1 - NSEW-Walks-PUBLIC.ipynb
Description:

Sage worksheet for a blog post on melczer.ca/blog/

Views: 177
Visibility: Unlisted (only visible to those who know the link)
Project: Main Files
In [1]:
```# This worksheet corresponds to a blogpost available on melczer.ca/blog/
# It explores the enumeration of walks on the steps N,S,E,W restricted to the first quadrant
# Much of the worksheet needs the ore_algebra package, available from https://github.com/mkauers/ore_algebra/
```
In [2]:
```# To begin, this function computes the number of walks on {N,S,E,W}, restricted to the first quadrant,
# which start at the origin and take n steps ending at (i,j)
@CachedFunction
def WalksIJ(i,j,n):
if i<0 or j<0:
return 0
elif n==0 and i==0 and j==0:
return 1
elif n==0:
return 0
else:
return WalksIJ(i-1,j,n-1) + WalksIJ(i+1,j,n-1) + WalksIJ(i,j-1,n-1) + WalksIJ(i,j+1,n-1)

# Compute the number of walks of length n, now ending anywhere in the quadrant
def Walks(n):
return sum([WalksIJ(i,j,n) for i in range(n+1) for j in range(n+1)])
```
In [3]:
```# First terms of the counting sequence
LST = [Walks(k) for k in range(11)]
print(LST)
```
Out[3]:
[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752, 116424]
In [4]:
```# Search OEIS for sequence with same terms up to index 10
OE = oeis(LST)[0]
print(OE)
```
Out[4]:
A005566: Number of walks of length n on square lattice, starting at origin, staying in first quadrant. 0: a(n) is the number of involutions of length 2n which are invariant under the reverse-complement map and have no decreasing subsequences of length 5. - _Eric S. Egge_, Oct 21 2008
In [5]:
```# The ore_algebra package is needed from here down
from ore_algebra import *

# Create the shift ore_algebra with rational coefficients
Ind.<n> = PolynomialRing(QQ)
Shift.<Sn> = OreAlgebra(Ind)

# Guess a recurrence satisfied by the sequence (need about 30 terms to find the recurrence)
LST = [Walks(k) for k in range(30)]
rec = guess(LST,Shift)
print(rec)
```
Out[5]:
(-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32
In [6]:
```# Use recurrence to generate many terms
LongLST = rec.to_list(LST,301)
print(LongLST[300])
```
Out[6]:
17523545967408829999363642806657506800280885492603311381994686521023183995818494526343136704823884594622982400835490350268164727785035412382829989246361352225213220805387866514176
In [7]:
```# Probability that a walk of length 300 stays in the first quadrant
(LongLST[300]/4^(300)).n()
```
Out[7]:
0.00422303415339021
In [8]:
```# Compute the first two asymptotic terms of a basis of the recurrence
rec.generalized_series_solutions(2)
```
Out[8]:
[4^n*n^(-1)*(1 - 3/2*n^(-1) + O(n^(-2))), (-4)^n*n^(-3)*(1 - 9/2*n^(-1) + O(n^(-2)))]
In [9]:
```# Approximate the constant lambda_1 such that w_n = lambda_1 * 4^n/n
LongLST = rec.to_list(LST,10001)
print((LongLST[100]/(4^100/100)).n())
print((LongLST[1000]/(4^1000/1000)).n())
print((LongLST[10000]/(4^10000/10000)).n())
```
Out[9]:
1.25446885783361 1.27133302123899 1.27304859221955
In [10]:
```# Create the differential algebra to encode linear differential equation
# (and shift algebra which expresses linear recurrence for coefficients)
Pols.<z> = PolynomialRing(QQ)
Diff.<Dz> = OreAlgebra(Pols)

# Guess an annihilating linear differential equation for generating function
diffWalk = guess(LST,Diff)
print(diffWalk)
```
Out[10]:
(16*z^4 - z^2)*Dz^3 + (128*z^3 + 8*z^2 - 6*z)*Dz^2 + (224*z^2 + 28*z - 6)*Dz + 64*z + 12
In [11]:
```# Converting from the differential equation to a recurrence gives rec (up to a left multiple)
print(Diff(diffWalk).to_S(Shift))
print((n + 2)*rec)
```
Out[11]:
(-n^3 - 9*n^2 - 26*n - 24)*Sn^2 + (8*n^2 + 36*n + 40)*Sn + 16*n^3 + 80*n^2 + 128*n + 64 (-n^3 - 9*n^2 - 26*n - 24)*Sn^2 + (8*n^2 + 36*n + 40)*Sn + 16*n^3 + 80*n^2 + 128*n + 64
In [12]:
```# Singularities of the generating function are roots of the leading differential term
```
Out[12]:
(4*z - 1) * (4*z + 1) * z^2
In [13]:
```# Get expansions at 0 of a basis of solutions for the differential equation
print(diffWalk.local_basis_expansions(0, order = 4))
```
Out[13]:
[z^(-2) - 4*z^(-1)*log(z) - 8*log(z) - 16*z*log(z) - 8*z, z^(-1), 1 + 2*z]
In [14]:
```# Get expansions at 1/4 of a basis of solutions for the differential equation
print(diffWalk.local_basis_expansions(1/4, order = 3))
```
Out[14]:
[log(z - 1/4) - 6*(z - 1/4)*log(z - 1/4) + 31*(z - 1/4)^2*log(z - 1/4), 1 - 14*(z - 1/4)^2, (z - 1/4) - 15/2*(z - 1/4)^2]
In [15]:
```# Find connection coefficients
ini = [0,0,1]
bas = diffWalk.numerical_transition_matrix([0,1/4]) * vector(ini)
print(bas)
```
Out[15]:
([-1.2732395447351627 +/- 2.54e-17] + [+/- 2.65e-22]*I, [-2.3906971441245563 +/- 1.99e-17] + [4.0000000000000000 +/- 2.78e-17]*I, [12.890661954217663 +/- 2.55e-16] + [-24.000000000000000 +/- 1.12e-16]*I)
In [16]:
```# Find the first coefficient to 1000 decimal places (takes around 2 seconds)
(diffWalk.numerical_transition_matrix([0,1/4],eps=1e-1000)* vector(ini))[0]
```
Out[16]:
[-1.2732395447351626861510701069801148962756771659236515899813387524711743810738122807209104221300246876485827418140636642951432947689166292230237392858535987138339197354992726205846219711754024615097391431697391812935468991219337890320946590409137815980457275236951206952213916489639152874955027526842297847995715472019893787819156768871865847742645992493358917024375955892175030525982925357139283119669931147888798709447935996995430046813694308674489401501372843723802959040779156829180746701444744199559730824426172542040257625982225311773281835739849567853267248481344242879850712958999906622932354822380405601299254205151107965704870409775950091814510221179031250645443716638278540904994185125596862001960023820788567124552223742810522016801305419681673984928499244916496251718727139876735314816926032604496069722128544177372731260597966617807819702831990012426351265118541792748660383765866297523255998072612616627947763148718624697387405122935160933003656475466210501492002089817437692256900796035 +/- 1.23e-1001] + [+/- 2.93e-1016]*I
In [17]:
```# The first basis element has leading singular term log(t-1/4)
# The extraction [z^n]log(1-4z) = -4^n/n implies our sequence grows as 4^n/n * C where C is
-bas[0]
```
Out[17]:
[1.2732395447351627 +/- 2.54e-17] + [+/- 2.65e-22]*I
In [18]:
```# The constant turns out to be 4/pi
(4/pi).n()
```
Out[18]:
1.27323954473516
In [0]:
```
```
In [0]:
```################################################################################################
# The following code is used to generate animations of random walks in the first quadrant
################################################################################################

# Generate endpoints of random walk using steps in Steps
def walkSum(Steps):
wk = [[0, 0]]
for k in Steps:
wk.append([wk[-1][0]+k[0],wk[-1][1]+k[1]])
return wk

# Find index walk wk leaves non-negative quadrant
def leaveIndex(wk):
for k in range(len(wk)):
if (wk[k][0]<0) or (wk[k][1]<0):
return k;
return len(wk)+1;

# Plot the state of the random walk wk at time k
# bd denotes the plot boundaries
def singlePlot(wk, k, bd):
if k < leaveIndex(wk):
oldSteps = line(wk[0:k], color='black', thickness=2)
newSteps = line(wk[k-1:k+1], color='red', thickness=2)
else:
oldSteps = line(wk[0:k], alpha=0.03, color='gray')
newSteps = line(wk[k-1:k+1], alpha=0.03, color='gray', thickness=2)
return oldSteps + newSteps

# Plot the state of all random walks in WK at time k
# bd denotes the plot boundaries
def singleListPlot(WK, k, bd):
plt = line([])
for wk in WK:
plt += singlePlot(wk,k,bd)

plt.set_axes_range(xmin=-bd,xmax=bd, ymin=-bd,ymax=bd)
numLeft = len(list(filter(lambda x:leaveIndex(x)>k,WK)))
txt = text("Number of steps: %s" %(k),(0,1.05),axis_coords=True, horizontal_alignment="left")
txt += text("Walks in first quadrant: %s" %(numLeft),(0,1),axis_coords=True,horizontal_alignment="left")

# Add plots and text, and set boundaries and aspect ratio
plt = plt + txt
plt.set_axes_range(xmin=-bd,xmax=bd, ymin=-bd,ymax=bd)
return plt

# Sequence of plots of first N steps of the random walks in WK
# bd denotes the plot boundaries
def plotWalk(WK, N, bd):
return [singleListPlot(WK,k,bd) for k in range(1,N+1)]
```
In [0]:
```# Define set of steps
ST = [[1,0],[-1,0],[0,1],[0,-1]]

# Define parameters K and N
K = 100
N = 300

# Generate steps for K random walks of length N
StepSeq = [[choice(ST) for k in range(N)] for j in range(K)]

# Generate the random walks and find maximum coordinate
WK = [walkSum(Steps) for Steps in StepSeq]
bd = max(flatten(WK))

# Get sequence of plots and convert to Sage Animation class
ani = animate(plotWalk(WK, N, bd+1))

# To display the animation in Sage, need imagemagick installed
# On a Mac with Homebrew, run "brew install imagemagick" in the terminal
# Then run the following to see the animation
os.environ["PATH"]+=":/usr/local/bin"
ani
```
In [0]:
```
```
In [12]:
```
```
In [0]:
```
```
In [0]:
```
```
In [0]:
```
```
In [0]:
```
```
In [0]:
```
```
In [0]:
```
```
In [0]:
```
```