
# Define the rings of differential and shift operators
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)


# The generating function of these walks satisfies the differential equation
diff = z^2*(4*z-1)*(4*z+1)*Dz^3+2*z*(4*z+1)*(16*z-3)*Dz^2+(2*(112*z^2+14*z-3))*Dz+64*z+12
diff

(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


# The counting sequence satisfies the recurrence defined by the following operator
rec = diff.to_S(Shift)
rec

(-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


# We verify the first terms match what we expect and check with OEIS
LST = rec.to_list([1,2],10)
LST

[1, 2, 6, 18, 60, 200, 700, 2450, 8820, 31752]


# Search the inital sequence terms in the OEIS *this requires an internet connection -- it is not required in the following*
oeis(LST)

0: A005566: Number of walks of length n on square lattice, starting at origin, staying in first quadrant.


# Asymptotically, the counting sequence is a C-linear combination of the following
show(rec.generalized_series_solutions(n=2))


# There is a basis of solutions for the C-vector space defined by the differential operator, with the following expansions at the origin
diff.local_basis_expansions(0, order = 4)

[z^(-2) - 4*z^(-1)*log(z) - 8*log(z) - 16*z*log(z) - 8*z, z^(-1), 1 + 2*z]


# Note the final basis element is the generating function, so in this basis the generating function has coordinates
ini = [0,0,1]


# There is a basis of solutions for the C-vector space defined by the differential operator, with the following expansions at z=1/4
diff.local_basis_expansions(1/4, order = 4)

[log(z - 1/4) - 6*(z - 1/4)*log(z - 1/4) + 31*(z - 1/4)^2*log(z - 1/4) - 150*(z - 1/4)^3*log(z - 1/4) - 2/3*(z - 1/4)^3,
 1 - 14*(z - 1/4)^2 + 108*(z - 1/4)^3,
 (z - 1/4) - 15/2*(z - 1/4)^2 + 43*(z - 1/4)^3]


# We can compute the change of basis matrix between the basis with series at the origin and the basis with series at 1/4
# We take a low numeric accuracy here to fit on screen
show(diff.numerical_transition_matrix([0,1/4],1e-1))


# Thus the generating function has a singular expansion at z=1/4 which begins
bas = diff.numerical_transition_matrix([0,1/4]) * vector(ini)
loc = diff.local_basis_expansions(1/4,2)
add([a*b for [a,b] in zip(bas,loc)])

([-1.2732395447351627+/-2.54e-17]+[+/-2.65e-22]*I)*log(z - 1/4) + ([7.6394372684109761+/-5.21e-17]+[+/-1.59e-21]*I)*(z - 1/4)*log(z - 1/4) + ([-2.3906971441245563+/-1.99e-17]+[4.0000000000000000+/-2.78e-17]*I)*1 + ([12.890661954217663+/-2.55e-16]+[-24.000000000000000+/-1.12e-16]*I)*(z - 1/4)


# The dominant singular behaviour is given by the term with log(z-1/4)
# Using the coefficient extraction [z^n]log(1-4*z) = -4^n/n dominant asymptotic behaviour is
var('n')
-bas[0].real() * 4^n/n

([1.2732395447351627 +/- 2.54e-17])*4^n/n


# We will determine the leading constant to be 4/pi
(4/pi).n()

1.27323954473516

Example 2.29 (Simple Lattice Walks in $\mathbb{N}^2$)¶