Highly Symmetric Walk Asymptotics

We give a selection of asymptotics for highly symmetric walks. Note: In the textbook we derive explicit asymptotic formulas which can be computed by hand. Here we verify those formulas on the examples from the textbook by computing asymptotics directly using the smooth point asymptotic results from Chapter 5.

Requirements: ore_algebra (partially)

# This cell imports the code for smooth asymptotics from Chapter 5

# Set a parameter to help simplify some algebraic numbers
maxima_calculus('algebraic: true;')

# Procedure to get Hessian appearing in asymptotics
# Input: H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
#        CP, a dictionary mapping elements of vars
# Output: The Hessian H defined in Lemma 5.5 of the textbook at the point w defined by CP
def getHes(H,r,vars,CP):
    dd = len(vars)
    V = zero_vector(SR,dd)
    U = matrix(SR,dd)
    M = matrix(SR,dd-1)

    for j in range(dd):
        V[j] = r[j]/r[-1]
        for i in range(dd):
            U[i,j] = vars[i]*vars[j]*diff(H,vars[i],vars[j])/vars[-1]/diff(H,vars[-1])
    for i in range(dd-1):
        for j in range(dd-1):
            M[i,j] = V[i]*V[j] + U[i,j] - V[j]*U[i,-1] - V[i]*U[j,-1] + V[i]*V[j]*U[-1,-1]
            if i == j: M[i,j] = M[i,j] + V[i]
    return M(CP)

# Procedure to apply differential operator to f and set all variables to zero
# Input: dop, element of a DifferentialWeylAlgebra over a polynomial ring
#        f, an element of the base polynomial ring of dop
# Output: dop(f) evaluated when all variables are zero
def eval_op(dop, f):
    if len(f.parent().gens()) == 1:
        return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[E[0][1][0]] for E in dop])
    else:
        return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[(v for v in E[0][1])] for E in dop])

# Procedure to get critical points of rational function with denominator H, in direction r
# Input: H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
# Output: Solutions (if found by solve) of the smooth critical point equations of H in the direction r
def critpt(H,r,vars):
    d = len(vars)
    criteqs = [H] + [r[j]*vars[0]*diff(H,vars[0]) - r[0]*vars[j]*diff(H,vars[j]) for j in range(1,d)]
    return solve(criteqs,vars,solution_dict=true)

# Procedure to compute asymptotic contribution of a strictly minimal contributing point
# Input: G, member of the symbolic ring
#        H, member of the symbolic ring
#        r, direction vector (which can contain symbolic entries)
#        vars, vector of variables
#        CP, a dictionary mapping elements of vars to coordinates of a strictly minimal contributing point
#        M, positive integer describing the number of terms in the asymptotic expansion to compute
#        g, parametrization of variable vars[-1] near CP, in terms of the remaining variables
# Output: ex, pw, se such that ex*pw*(se+O(n^(M-1)) gives an asymptotic expansion of the r-diagonal of 
#         G/H in the variables vars, to order M.
# NOTE: Unlike the textbook, M here refers to the number of terms in the expansion
#       (not the order of the expansion, so M should be at least 1)
def smoothContrib(G,H,r,vars,CP,M,g):
    # Preliminary definitions
    dd = len(vars)
    field = SR
    tvars = list(var('t%d'%i) for i in range(dd-1))
    dvars = list(var('dt%d'%i) for i in range(dd-1))

    # Define differential Weyl algebra and set variable names
    W = DifferentialWeylAlgebra(PolynomialRing(field,tvars))
    WR = W.base_ring()
    T = PolynomialRing(field,tvars).gens()
    D = list(W.differentials())

    # Compute Hessian matrix and differential operator Epsilon
    HES = getHes(H,r,vars,CP)
    HESinv = HES.inverse()
    v = matrix(W,[D[k] for k in range(dd-1)])
    Epsilon = -(v * HESinv.change_ring(W) * v.transpose())[0,0]

    # Define quantities for calculating asymptotics
    tsubs = [v == v.subs(CP)*exp(I*t) for [v,t] in zip(vars,tvars)]
    tsubs += [vars[-1]==g.subs(tsubs)]
    P = (-G/g/diff(H,vars[-1])).subs(tsubs)
    psi = log(g.subs(tsubs)/g.subs(CP)) + I * add([r[k]*tvars[k] for k in range(dd-1)])/r[-1]
    v = matrix(SR,[tvars[k] for k in range(dd-1)])
    psiTilde = psi - (v * HES * v.transpose())[0,0]/2

    # Recursive function to convert symbolic expression to polynomial in t variables
    def to_poly(p,k):
        if k == 0:
            return add([a*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])
        return add([to_poly(a,k-1)*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])

    # Compute Taylor expansions to sufficient orders
    N = 2*M
    PsiSeries = to_poly(taylor(psiTilde,*((v,0) for v in tvars), N),dd-2)
    PSeries = to_poly(taylor(P,*((v,0) for v in tvars), N),dd-2)

    # Precompute products used for asymptotics
    EE = [Epsilon^k for k in range(3*M-2)]
    PP = [PSeries] + [0 for k in range(2*M-2)]
    for k in range(1,2*M-1):
        PP[k] = PP[k-1]*PsiSeries
    
    # Function to compute constants appearing in asymptotic expansion
    def Clj(l,j):
        return (-1)^j*SR(eval_op(EE[l+j],PP[l]))/(2^(l+j)*factorial(l)*factorial(l+j))
    
    # Compute different parts of asymptotic expansion
    var('n')
    ex = (prod([1/v^k for (v,k) in zip(vars,r)]).subs(CP).canonicalize_radical())^n
    pw = (r[-1]*n)^((1-dd)/2)
    se = sqrt((2*pi)^(1-dd)/HES.det()) * add([add([Clj(l,j) for l in range(2*j+1)])/(r[-1]*n)^j for j in range(M)])
    
    return ex, pw, se.canonicalize_radical()

# Procedure to aid in printing an asymptotic expansion
# Procedure to get critical points of rational function with denominator H, in direction r
# Input: ex,pw,se as returned by smoothContrib(G,H,r,vars,CP,M,g)
# Output: None (function pretty prints the asymptotic expression defined by ex,pw,se, and M)
def disp_asm(ex,pw,se,M):
    show(ex*pw,LatexExpr("\\Bigg("), se, LatexExpr("+ O\\Bigg("), n^(-M), LatexExpr("\\Bigg)\\Bigg)"))

# Test if all coordinates in a list of numbers are real and positive
def pos_coords(L):
    return all(map(lambda c: (c>0) and c.is_real(), L))

# Find critical points with real positive coordinates
def walk_pos_critpt(L,vars):
    return filter(lambda l: pos_coords([v.subs(l) for v in vars]), L)

Example 6.2 (Highly Symmetric Models in Two Dimensions)

# Enter the short step sets defining highly symmetric quadrant models
N = (0,1); SS = (0,-1); E = (1,0); W = (-1,0);
NE = (1,1); NW = (-1,1); SE = (1,-1); SW = (-1,-1);
HS1  = [N,SS,E,W]
HS2  = [NE,SE,NW,SW]
HS3  = [N,SS,NE,SE,NW,SW]
HS4  = [N,SS,E,W,NW,SW,SE,NE]

# Define basic quantities
var('x,y,t')
vars = [x,y,t]
r = [1,1,1]
M = 1

# Loop through the four models under consideration
for ST in [HS1,HS2,HS3,HS4]:
    # Define the rational function for the model
    S = add([x^i*y^j for (i,j) in ST])
    G = (1+x)*(1+y)
    H = 1-t*x*y*S

    # Get the minimal critical point with positive coordinates (and prove it is unique)
    [cp] = walk_pos_critpt(critpt(H,r,vars),vars)
    
    # Compute and print asymptotics
    g = solve(H,vars[-1])[0].rhs()
    ex,pw,se = smoothContrib(G,H,r,vars,cp,M,g)
    print("The number of quadrant models on the highly-symmetric step set {} has asymptotic expansion".format(ST))
    disp_asm(ex,pw,se,M)

The number of quadrant models on the highly-symmetric step set [(0, 1), (0, -1), (1, 0), (-1, 0)] has asymptotic expansion

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{4^{n}}{n} \Bigg( \frac{4}{\pi} + O\Bigg( \frac{1}{n} \Bigg)\Bigg)$$

The number of quadrant models on the highly-symmetric step set [(1, 1), (1, -1), (-1, 1), (-1, -1)] has asymptotic expansion

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{4^{n}}{n} \Bigg( \frac{2}{\pi} + O\Bigg( \frac{1}{n} \Bigg)\Bigg)$$

The number of quadrant models on the highly-symmetric step set [(0, 1), (0, -1), (1, 1), (1, -1), (-1, 1), (-1, -1)] has asymptotic expansion

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{6^{n}}{n} \Bigg( \frac{\sqrt{3} \sqrt{2}}{\pi} + O\Bigg( \frac{1}{n} \Bigg)\Bigg)$$

The number of quadrant models on the highly-symmetric step set [(0, 1), (0, -1), (1, 0), (-1, 0), (-1, 1), (-1, -1), (1, -1), (1, 1)] has asymptotic expansion

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{8^{n}}{n} \Bigg( \frac{8}{3 \, \pi} + O\Bigg( \frac{1}{n} \Bigg)\Bigg)$$

Example 6.3 (A Weighted Model)

# Define basic quantities
var('x,y,t,a,b,c')
assume(a>0,b>0,c>0)
vars = [x,y,t]
r = [1,1,1]

# Define the rational function for the model
S = (y+1/y)*(a*x+b+a/x) + (c*x+c/x)
G = (1+x)*(1+y)
H = 1-t*x*y*S

# Get the minimal critical point with positive coordinates (and prove it is unique)
[cp] = list(walk_pos_critpt(critpt(H,r,vars),vars))
cp

# Compute and print asymptotics
M = 1
g = solve(H,vars[-1])[0].rhs()
ex,pw,se = smoothContrib(G,H,r,vars,cp,M,g)
print("The number of weighted quadrant models has asymptotic expansion")
disp_asm(ex,pw,se.factor(),M)

The number of weighted quadrant models has asymptotic expansion

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(4 \, a + 2 \, b + 2 \, c\right)}^{n}}{n} \Bigg( \frac{2 \, {\left(2 \, a + b + c\right)}}{\pi \sqrt{2 \, a + b} \sqrt{2 \, a + c}} + O\Bigg( \frac{1}{n} \Bigg)\Bigg)$$

Example 6.5 (Higher Asymptotics Terms of Simple Walks)

Note: this matches the computations for Example 5.6 in Chapter 5.

# Define basic quantities
var('x,y,t')
r = [1,1,1]
vars = [x,y,t]

# Define the rational function for the model
S = add([x^i*y^j for (i,j) in HS1])
G = (1+x)*(1+y)
H = 1-t*x*y*S

# Find both (minimal) critical points
[sigma, tau] = critpt(H,r,vars)

# Determine the asymptotic contribution of each critical point
M = 4
g = solve(H,vars[-1])[0].rhs()
ex1,pw1,se1 = smoothContrib(G,H,r,vars,sigma,M,g)
ex2,pw2,se2 = smoothContrib(G,H,r,vars,tau,M,g)

# Print the asymptotic contribution of each critical point
print("The asymptotic contribution of {} to diagonal asymptotics is".format(sigma))
disp_asm(ex1,pw1,se1.expand(),M)
print("The asymptotic contribution of {} to diagonal asymptotics is".format(tau))
disp_asm(ex2,pw2,se2.expand(),M)

The asymptotic contribution of {x: 1, y: 1, t: 1/4} to diagonal asymptotics is

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{4^{n}}{n} \Bigg( \frac{4}{\pi} - \frac{6}{\pi n} + \frac{19}{2 \, \pi n^{2}} - \frac{63}{4 \, \pi n^{3}} + O\Bigg( \frac{1}{n^{4}} \Bigg)\Bigg)$$

The asymptotic contribution of {x: -1, y: -1, t: -1/4} to diagonal asymptotics is

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\left(-4\right)^{n}}{n} \Bigg( \frac{1}{\pi n^{2}} - \frac{9}{2 \, \pi n^{3}} + O\Bigg( \frac{1}{n^{4}} \Bigg)\Bigg)$$

# Note that additional terms in the expansion can be computed very efficiently 
# from a P-recurrence satisfied by the diagonal sequence / D-finite equation of the GF
from ore_algebra import *
Pols.<n> = PolynomialRing(QQ)
Shift.<Sn> = OreAlgebra(Pols)

# Import recurrence satisfied by diagonal sequence (obtained from kernel method and creative telescoping)
rec = (-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32

# The asymptotic expansion of the diagonal will be a C-linear combination of these terms
# Looking at the dominant terms allows one to compute the coefficients (here 4/pi and 1/pi)
# This approach can compute about 20 terms in the expansion in one second on a modern laptop
show(rec.generalized_series_solutions(5))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[(4)^{n}n^{-1}\Bigl(1 - \frac{3}{2}n^{-1} + \frac{19}{8}n^{-2} - \frac{63}{16}n^{-3} + \frac{871}{128}n^{-4} + O(n^{-5})\Bigr), (-4)^{n}n^{-3}\Bigl(1 - \frac{9}{2}n^{-1} + \frac{107}{8}n^{-2} - \frac{525}{16}n^{-3} + \frac{9263}{128}n^{-4} + O(n^{-5})\Bigr)\right]$$

Example 6.9 (Higher Order Constants for Simple Walk in Two Dimensions)

# Set basic info
var('x,y,t,a,b')
assume(a,'integer')
assume(b,'integer')
r = [1,1,1]
vars = [x,y,t]

# Define the rational function for the number of quadrant walks on [N,S,E,W]
# which start at the parametrized point (a,b)
S = add([x^i*y^j for (i,j) in HS1])
G = (1-x^(2*a+2))*(1-y^(2*b+2))/x^a/y^b/(1-x)/(1-y)
H = 1-t*x*y*S

# Find both (minimal) critical points
[sigma, tau] = critpt(H,r,vars)

# Determine the asymptotic contribution of each critical point
# (this can take a few minutes to compute the higher order terms)
M = 3
g = solve(H,vars[-1])[0].rhs()
ex1,pw1,se1 = smoothContrib(G,H,r,vars,sigma,M,g)
ex2,pw2,se2 = smoothContrib(G,H,r,vars,tau,M,g)

# Asymptotics for the model with parametrized start is 
# given by adding these two asymptotic contributions

print("The asymptotic contribution of {} to diagonal asymptotics is".format(sigma))
var('N')
se_simp1 = add([k.factor()/n^p for [k,p] in se1.subs(n=1/N).series(N,M).coefficients()])
disp_asm(ex1,pw1,se_simp1,M)

print("The asymptotic contribution of {} to diagonal asymptotics is".format(tau))
disp_asm(ex2,pw2,se2,M)

The asymptotic contribution of {x: 1, y: 1, t: 1/4} to diagonal asymptotics is

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{4^{n}}{n} \Bigg( \frac{4 \, {\left(a + 1\right)} {\left(b + 1\right)}}{\pi} - \frac{2 \, {\left(2 \, a^{2} + 2 \, b^{2} + 4 \, a + 4 \, b + 9\right)} {\left(a + 1\right)} {\left(b + 1\right)}}{3 \, \pi n} + \frac{{\left(36 \, a^{4} + 40 \, a^{2} b^{2} + 36 \, b^{4} + 144 \, a^{3} + 80 \, a^{2} b + 80 \, a b^{2} + 144 \, b^{3} + 516 \, a^{2} + 160 \, a b + 516 \, b^{2} + 744 \, a + 744 \, b + 855\right)} {\left(a + 1\right)} {\left(b + 1\right)}}{90 \, \pi n^{2}} + O\Bigg( \frac{1}{n^{3}} \Bigg)\Bigg)$$

The asymptotic contribution of {x: -1, y: -1, t: -1/4} to diagonal asymptotics is

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\left(-4\right)^{n}}{n} \Bigg( \frac{{\left(\left(-1\right)^{a} {\left(a + 1\right)} b + \left(-1\right)^{a} {\left(a + 1\right)}\right)} \left(-1\right)^{b}}{\pi n^{2}} + O\Bigg( \frac{1}{n^{3}} \Bigg)\Bigg)$$