{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Higher Asymptotics of Lattice Path Walks on {N,S,E,W} Steps in a Quarter Plane\n", "By Keith Ritchie and Stephen Melczer\n", "\n", "Here we compute high-order asymptotic expansions for the number of lattice walks beginning at the origin of $\\mathbb{Z}^2$, taking the steps $\\{(\\pm1,0),(0,\\pm1)\\}$, and staying in the non-negative quadrant $\\mathbb{N}^2$. Details on why this computes the desired expansion can be found in Chapter 6 of [**An Invitation to Analytic Combinatorics**](https://melczer.ca/textbook/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we introduce necessary preliminaries, and define our function to compute asymptotic expansions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Define Weyl algebra (consisting of differentials and variables)\n", "W. = DifferentialWeylAlgebra(QQ);\n", "dt1,dt2 = W.differentials()\n", "var('T1 T2')\n", "\n", "#Definition of application of differential operator.\n", "def applyOp(Dop,f):\n", " return sum([i[1]*T1^(i[0][0][0])*T2^(i[0][0][1])*diff(f,T1,i[0][1][0],T2,i[0][1][1]) for i in Dop.list()])\n", "\n", "#Define Psi function\n", "Ps = -log(cos(T1)+cos(T2)) + log(2) - (T1^2+T2^2)/4\n", "\n", "#Define P functions for rho = 1 (P1) and rho = -1 (P2)\n", "P1 = (1 + e^(i*T1))*(1 + e^(i*T2))\n", "P2 = (1 - e^(i*T1))*(1 - e^(i*T2))\n", "\n", "#Define differential operator Epsilon.\n", "Epsilon = (-2)*(dt1^2+dt2^2)\n", "\n", "# Define variables and assumptions\n", "var('l j n')\n", "assume(l>0, j>0, n>0)\n", "assume(l, 'integer', j, 'integer', n, 'integer')\n", "\n", "#Sequence of powers of Epsilon applied to f*P1 for arbitrary l at (0,0)\n", "def a1(p,l):\n", " return applyOp(Epsilon^p, Ps^l*P1).subs(T1=0,T2=0)\n", "\n", "#Sequence of powers of Epsilon applied to f*P2 for arbitrary l at (0,0)\n", "def a2(p,l):\n", " return applyOp(Epsilon^p, Ps^l*P2).subs(T1=0,T2=0)\n", "\n", "#Define Kj terms for rho = 1 and rho = -1 as a sequence in j following the general construction on page 234.\n", "def Kj1(j):\n", " return (-1)^j *sum(a1(l+j,l)/(2^(j+l)*factorial(l)*factorial(j+l)) for l in (0..2*j))\n", "\n", "def Kj2(j):\n", " return (-1)^j *sum(a2(l+j,l)/(2^(j+l)*factorial(l)*factorial(j+l)) for l in (0..2*j))\n", "\n", "#Define expression for asymptotics of Cn as a function of M, a natural number.\n", "def Cn(M):\n", " return 4^n/(pi*n)*sum((Kj1(j)+Kj2(j)*(-1)^n)*n^(-j) for j in (0..M))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we give some examples" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "4*4^n/(pi*n)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "-2*4^n*(3/n - 2)/(pi*n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Generate first few expansions\n", "ASM0 = Cn(0)\n", "ASM1 = Cn(1)\n", "show(ASM0)\n", "show(ASM1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "1/2*4^n*((2*(-1)^n + 19)/n^2 - 12/n + 8)/(pi*n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Take order two expansion\n", "ASM2 = Cn(2)\n", "show(ASM2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We check our formulas by comparing them to numerically generated data." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's check our formula\n", "# Code to generate number of walks of length n\n", "@CachedFunction \n", "def WalksIJ(i,j,n):\n", " if i<0 or j<0:\n", " return 0\n", " elif n==0 and i==0 and j==0:\n", " return 1\n", " elif n==0:\n", " return 0\n", " else:\n", " 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)\n", " \n", "def Walks(n): \n", " return sum([WalksIJ(i,j,n) for i in range(n+1) for j in range(n+1)])\n", "\n", "# Test code\n", "[Walks(k) for k in range(10)]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "116424\n", "116987.123578877\n" ] } ], "source": [ "# Pretty close already at n=10\n", "print(Walks(10))\n", "print(ASM2.subs(n=10).n())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each successive expansion adds about one decimal place of accuracy when $n=20$, as expected." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.07412849449210\n", "0.993568857405193\n", "1.00061782565030\n" ] } ], "source": [ "print(ASM0.subs(n=20).n()/Walks(20))\n", "print(ASM1.subs(n=20).n()/Walks(20))\n", "print(ASM2.subs(n=20).n()/Walks(20))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.1", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }