CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Project: ALEA ACSV
Views: 17
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004
Tweet
Share
Share
Kernel: SageMath 9.5

Exercises: An Invitation to Analytic Combinatorics in Several Variables

Created by Stephen Melczer

This notebook complements the exercises found on the ALEA ACSV course webpage

A quick Sage tutorial can be found here (try an interactive version in your browser here).

See https://melczer.ca/textbook for further Sage notebooks solving problems in analytic combinatorics in several varibles. In particular, this notebook (also available as a static HTML page) gives an algorithm to compute asymptotic terms. Don't use it to solve these exercises, but use it to further check your answers! (Or to compute asymptotics for other problems!)

# Helper function to compute the Hessian matrix M from the function H(vars)
# in the direction R at the point CP, specified by a list of substitutions.
# Copied from melczer.ca/textbook/
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.subs(CP)


# Helper function to compute leading asymptotics of the R-diagonal of G(vars)/H(vars)
# determined by the Main Asymptotic Theorem of Smooth ACSV at the point CP, specified
# by a list of substitutions. We take det(M) as an input that can be computed by the
# above function.
var('n')
def leadingASM(G,H,detM,R,vars,CP):
    dd = len(R)
    lcoeff = -G/vars[-1]/H.diff(vars[-1])
    exp = 1/mul([vars[k]^R[k] for k in range(dd)])^n

    ASM = exp * (2*pi*n*R[-1])^((1-dd)/2) / sqrt(detM) * lcoeff
    return ASM.subs(CP)

These functions can be used to compute the matrix M appearing in asymptotics, as well as the leading asymptotic term in an asymptotic expansion.

Here is an example of their use to find asymptotics for the main diagonal of 1/(1−x−y−z)1/(1-x-y-z).

# Introduce variables x, y, and z
var('x y z')
H = 1 - x - y - z

# In the main diagonal direction this has a critical point at (1/3,1/3,1/3)
CPs = solve([H,x*diff(H,x) - y*diff(H,y), x*diff(H,x) - z*diff(H,z)],[x,y,z])
show(CPs)
[[x=(13),y=(13),z=(13)]]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\left[x = \left(\frac{1}{3}\right), y = \left(\frac{1}{3}\right), z = \left(\frac{1}{3}\right)\right]\right]
# Let CP be the critical point, defined as a list [x == 1/3, y == 1/3, z == 1/3]
CP = CPs[0]

# Get the matrix M 
M = getHes(H,[1,1,1],[x,y,z],CP)
show("M = ", M)
M =(2112)\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|M|\verb| |\verb|=| \left(\begin{array}{rr} 2 & 1 \\ 1 & 2 \end{array}\right)
# Get and print leading asymptotics
ASM = leadingASM(1,H,M.determinant(),[1,1,1],[x,y,z],CP)
print("The dominant asymptotic behaviour of the main diagonal is")
show(ASM)
The dominant asymptotic behaviour of the main diagonal is
32 π(127)nn\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\sqrt{3}}{2 \, \pi \left(\frac{1}{27}\right)^{n} n}

We can check our asymptotic approximation by computing series terms.

# First, define the ring of formal power series (more efficient for computations)
# We use capital letters to denote the formal power series variables
P.<X,Y,Z> =  QQ[['X,Y,Z']] 

# Computes the series expansion up to precision 3*N
N = 10
ser = 1/(1-X-Y-Z + O(X^(3*N+1)))

# Check ratio of asymptotic formula to actual coefficients -- this should go to 1!
(ser.coefficients()[X^N*Y^N*Z^N]/ASM.subs(n=N)).n()
0.978030899006393

Question 1: Delannoy Numbers

The Delannoy number da,bd_{a,b} is the number of paths from the origin (0,0)(0,0) to the point (a,b)(a,b) using only the steps N=(0,1)\textsf{N}=(0,1), E=(1,0)\textsf{E} = (1,0), and NE=(1,1)\textsf{NE}=(1,1).

(a) Prove the recurrence $$ d_{a,b} = \begin{cases} 1 &: \text{ if $a=0$ or b=0b=0} \ d_{a-1,b} + d_{a,b-1} + d_{a-1,b-1} &:\text{ otherwise} \end{cases} Concludethat Conclude that D(x,y) = \sum_{a,b\geq0}d_{a,b}x^ay^b = \frac{1}{1-x-y-xy}. $$

(b) Use the Main Theorem of Smooth ACSV to find asymptotics of dn,nd_{n,n} as the (1,1)(1,1)-diagonal of D(x,y)D(x,y). What are the critical points in the (1,1)(1,1) direction? Which are minimal?

(c) Use the Main Theorem of Smooth ACSV to find asymptotics of the (r,s)(r,s)-diagonal of D(x,y)D(x,y) for any r,s>0r,s>0.

Function to numerically compute terms in the expansion

Check your computed asymptotics against this function!

# First, define the ring of formal power series (more efficient for computations)
P.<X,Y> =  QQ[['X,Y']] 

# Code to compute the coefficient of x^(N*R) * y^(N*S) where R, S, and N are positive integers
R, S, N = 1, 2, 10
N2 = (R+S)*N
ser = 1/(1-X-Y-X*Y + O(X^(N2+1)))
coef = ser.coefficients()[X^(R*N)*Y^(S*N)]
print("The coefficient [x^{}y^{}]D(x,y) = {}".format(N*R,N*S,coef))
The coefficient [x^10y^20]D(x,y) = 4354393801

Question 2: Apéry Asymptotics

Recall from lecture that a key step in Apéry's proof of the irrationality of ζ(3)\zeta(3) is determining the exponential growth of the sequence that can be encoded as the main diagonal of F(x,y,z,t)=11−t(1+x)(1+y)(1+z)(1+y+z+yz+xyz). F(x,y,z,t) = \frac{1}{1 - t(1+x)(1+y)(1+z)(1+y+z+yz+xyz)}. Use the Main Theorem of Smooth ACSV to find dominant asymptotics of this sequence.

Function to numerically compute terms in the expansion

Check your computed asymptotics against this function!

# Numerically compute coefficient of (xyzt)^n (can take a long time for large N)
P.<X,Y,Z> =  QQ['X,Y,Z'] 
N = 30
ser = ((1+X)*(1+Y)*(1+Z)*(1+Y+Z+Y*Z+X*Y*Z))^N
coef = ser[X^N*Y^N*Z^N]
print("The coefficient [(xyzt)^({})]F(x,y,z,t) = {}".format(N,coef))
The coefficient [(xyzt)^(30)]F(x,y,z,t) = 11320115195385966907843180411829810312080825

Question 3: Pathological Directions

(a) Find asymptotics of the (r,s)(r,s)-diagonal of F(x,y)=11−x−xyF(x,y) = \frac{1}{1-x-xy} for any 0<s<r0<s<r.

(b) What are the critical points of F(x,y)=11−x−xyF(x,y) = \frac{1}{1-x-xy} in the (r,s)(r,s) direction when 0<r≤s0<r \leq s? Which are minimal? Characterize the behaviour of the (r,s)(r,s) diagonal when 0<r≤s0<r \leq s.

Function to numerically compute terms in the expansion

Check your computed asymptotics against this function!

# First, define the ring of formal power series (more efficient for computations)
P.<X,Y> =  QQ[['X,Y']] 

# Code to compute the coefficient of x^(N*R) * y^(N*S) where R, S, and N are positive integers
R, S, N = 2, 1, 10
N2 = (R+S)*N
ser = 1/(1-X-X*Y + O(X^(N2+1)))
coef = ser.coefficients()[X^(R*N)*Y^(S*N)]
print("The coefficient [x^{}y^{}]F(x,y) = {}".format(N*R,N*S,coef))
The coefficient [x^20y^10]F(x,y) = 184756

Question 4: A Composition Limit Theorem

An integer composition of size n∈Nn\in\mathbb{N} is an ordered tuple of positive integers (of any length) that sum to nn. A composition can be viewed as an integer partition where the order of the summands matters. Let ck,nc_{k,n} denote the number of compositions of size nn that contain kk ones.

(a) If you know the symbolic method, species theory, or similar enumerative constructions, prove that C(u,x)=∑n,k≥0ck,nukxn=1−x1−2x−(u−1)x(1−x). C(u,x) = \sum_{n,k\geq0}c_{k,n}u^kx^n = \frac{1-x}{1-2x-(u-1)x(1-x)}.

(b) Prove that the distribution for the number of ones in a composition of size nn satisfies a local central limit theorem as n→∞n\rightarrow\infty. More precisely, find a constant t>0t>0 and normal density νn(s)\nu_n(s) such that sup⁡s∈Z∣tncs,n−νn(s)∣→0 \sup_{s \in \mathbb{Z}} |t^nc_{s,n} - \nu_n(s)| \rightarrow 0 as n→∞n\rightarrow\infty.

Function to plot series coefficients

Check your computed distribution against this function!

# Plot series terms versus computed density
K.<U> =  QQ['U']
P.<X> =  K[['X']]

# Set the value of n to test
N = 200
mser = (1 - X)/(1 - 2*X - (U-1)*X*(1-X) + O(X^(N+1)))
uvals = mser[N]

plt = point([])
for k in range(N):
    plt += point([k,uvals[k]])

print("The following plot shows the coefficients of [x^({})]C(u,x)".format(N))
plt
The following plot shows the coefficients of [x^(200)]C(u,x)
Image in a Jupyter notebook