Determining asymptotics of the central binomial coefficients squared.
Requirements: ore_algebra package
# 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)
# Start with the previously calculated shift and differential operators
rec = (n+1)^2*Sn - 4*(2*n+1)^2
diff = z*(1-16*z)*Dz^2+(1-32*z)*Dz-4
# Asymptotic expansions of a basis to the recurrence
rec.generalized_series_solutions()
[16^n*n^(-1)*(1 - 1/4*n^(-1) + 1/32*n^(-2) + 1/128*n^(-3) - 5/2048*n^(-4) + O(n^(-5)))]
# The generating function has singularities at the roots of the leading coefficient polynomial
diff.leading_coefficient().factor()
(-16) * (z - 1/16) * z
# There is a basis with the following series expansions near the origin
origin_basis = diff.local_basis_expansions(0, order = 3)
show(origin_basis)
# At the origin, the power series element is the generating function
# So the generating function is expressed in this basis with coordinates
ini = [0,1]
print("The generating function has series expansion {} + ... at the origin".format(add([origin_basis[k] * ini[k] for k in range(2)])))
The generating function has series expansion 1 + 4*z + 36*z^2 + ... at the origin
# Expansions of a basis near the generating function singularity z=1/16
show(diff.local_basis_expansions(1/16, order = 2))
# We can compute the change of basis matrix between the basis with series at the origin and the basis with series at 1/6
# We take a low numeric accuracy here to fit on screen
show(diff.numerical_transition_matrix([0,1/16],1e-10))
# Thus the generating function has a singular expansion at z=1/16 which begins
bas = diff.numerical_transition_matrix([0,1/16]) * vector(ini)
loc = diff.local_basis_expansions(1/16,2)
add([a*b for [a,b] in zip(bas,loc)])
([-0.31830988618379067+/-2.14e-18]+[+/-2.98e-38]*I)*log(z - 1/16) + ([1.2732395447351627+/-2.54e-17]+[+/-1.19e-37]*I)*(z - 1/16)*log(z - 1/16) + ([2.5464790894703253+/-7.71e-17]+[-4.0000000000000000+/-4.17e-17]*I)*(z - 1/16) + ([+/-2.37e-30]+[1.0000000000000000+/-6.94e-18]*I)*1
# The extraction [z^n]log(1-16z) = -16^n/n then implies the central binomial coefficients have asymptotic expansion
var('n')
bas[0].real() * (-16^n/n)
([0.31830988618379067 +/- 2.14e-18])*16^n/n
# The connection constant is 1/pi
(1/pi).n()
0.318309886183791