Square CoCalc Logo
FeaturesSoftwarePricingInfoPoliciesShareSupport Try Sign InSign Up
Path: Blog/PUBLIC/1 - NSEW-Walks-PUBLIC.ipynb
Description:

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

Views: 176
License: none
Visibility: Unlisted (only visible to those who know the link)
Project: Main Files
Raw | Embed | Download |
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) print(OE.comments())
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 diffWalk.leading_coefficient().factor()
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) # Define header text 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]: