We run our algorithm and compute the Kronecker representation "by hand" using Sage's capability to compute Gröbner bases.
Requirements: ore_algebra (for optional heuristic check of asymptotics at end)
# We will need our procedure to get the Hessian matrix appearing in asymptotics
# This code copied from Chapter 5
def getHes(H,vars):
dd = len(vars)
V = zero_vector(SR,dd)
U = matrix(SR,dd)
M = matrix(SR,dd-1)
for j in range(dd):
V[j] = 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
# Set basic info for this example
var('x,y,z,λ,t,u')
H = 1 - (x+y+z) + 81/8*x*y*z
# Because we no longer have a combinatorial example, we must split the
# variables into real/imaginary components and define the extended critical point system
var('a1,a2,a3,b1,b2,b3,lambdaR,lambdaI')
assume(a1,'real'); assume(a2,'real'); assume(a3,'real')
assume(b1,'real'); assume(b2,'real'); assume(b3,'real')
Hc = H.subs(x=a1+I*b1, y=a2+I*b2, z=a3+I*b3)
Hr = Hc.real()
Hi = Hc.imag()
CPab1 = [diff(Hr,v)*v + diff(Hr,w)*w - lambdaR for (v,w) in zip([a1,a2,a3],[b1,b2,b3])]
CPab2 = [diff(Hi,v)*v + diff(Hi,w)*w - lambdaI for (v,w) in zip([a1,a2,a3],[b1,b2,b3])]
# Compute a Gröbner basis for the critical point equations in the a and b variables
R = PolynomialRing(QQ,8,'a1,a2,a3,b1,b2,b3,lambdaR,lambdaI',order='lex')
J = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2])
GB = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2]).groebner_basis()
# Note there are a finite number of (complex) solutions!
print("The ideal generated by {} has dimension {}.".format(GB,J.dimension()))
# To determine minimality, we add the extra equations using x and y variables
# First, define xj + I*yj as a root of H
# Then encode that the coordinate-wise moduli of xj + I*yj is a multiple of crit pt aj + I*bj
# Finally, add the extra equations derived from critical point method
var('x1,x2,x3,y1,y2,y3,ν')
Hrxy = Hr.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)
Hixy = Hi.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)
mod = [X^2+Y^2 - t*(A^2+B^2) for (X,Y,A,B) in zip([x1,x2,x3],[y1,y2,y3],[a1,a2,a3],[b1,b2,b3])]
nuxy = [(Y - ν*X)*diff(Hrxy,X) - (X + ν*Y)*diff(Hrxy,Y) for (X,Y) in zip([x1,x2,x3],[y1,y2,y3])]
# Define the ideal with all these equations
S = PolynomialRing(QQ,16,'ν,lambdaR,lambdaI,x1,x2,x3,y1,y2,y3,b1,b2,b3,a1,a2,a3,t',order='lex')
Jxy = S.ideal([S(k) for k in list(GB) + [Hrxy,Hixy] + mod + nuxy])
# Compute the corresponding Gröbner basis. This can take several minutes.
# (Note: the same Gröbner computation in Maple is almost instantaneous)
GBxy = Jxy.groebner_basis()
# As expected, there are a finite number of solutions
# At these solutions, the t variable is a root of the following
GBxy[-1].factor()
# There are two values between 0 and 1 -- corresponding to non-minimal critical points
tPol = QQbar[t](SR(GBxy[-1]))
badt = list(filter(lambda r: SR(r).is_real() and r>0 and r<1, tPol.roots(multiplicities=False)))
badt
# Here we list all the non-minimal critical points
# In fact, there is only one non-minimal critical point
# This implies all other critical points are minimal
S2 = PolynomialRing(QQbar,15,'ν,lambdaR,lambdaI,x1,x2,x3,y1,y2,y3,b1,b2,b3,a1,a2,a3',order='lex')
for T in badt:
J = S2.ideal([S2(SR(k.subs(t=T))) for k in GBxy])
V = J.variety()
for v in V:
print([(v[j] + QQbar(I)*v[k]).radical_expression() for (j,k) in zip([a1,a2,a3],[b1,b2,b3])])
# Now we find asymptotics in terms of the minimal critical points
# Form the extended critical point system
vars = [x,y,z,λ]
C = [H] + [v * diff(H,v) - λ for v in vars[0:-1]]
# Define the multivariate polynomial ring and compute a Gröbner basis
R = PolynomialRing(QQbar,6,'x,y,z,λ,t,u',order='lex')
C = C + [u - x]
GB = R.ideal([R(k) for k in C]).groebner_basis()
# The Kronecker is representation parametrized by the roots of the following polynomial
P = QQ[u](GB[-1])
P
# The polynomial P is square-free: its gcd with its derivative is zero
# (Note: we need to cast between polynomial and symbolic rings to compute desired quantities)
QQ[u](P).gcd(QQ[u](SR(P).diff(u)))
# Compute the Q polynomials for the Kronecker representation from the Gröbner basis
Pd = diff(P,u)
Q = {}
for v in vars:
[Qv] = list(filter(lambda p: SR(p).has(v), GB))
_,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)
sys = [SR(Q[v]/Pd) for v in vars]
# For instance, here is the first variable in the Kronecker representation
show("In the Kronecker rep. the variable \t",vars[0], " is represented by\t",Q[vars[0]]/Pd)
# There are three minimal critical points: the non-minimal point (-2/3,-2/3,-2/3)
# and two minimal critical points
rts = P.roots(QQbar, multiplicities=False)
for uval in rts:
print("The root u = {} gives the critical point".format(uval))
print("[x,y] =",[p.subs(u == uval) for p in sys[:-2]],"\n")
# Select the minimal critical points
[u1,u2] = filter(lambda U: not sys[0].subs(u==U)== -2/3, P.roots(QQbar, multiplicities=False))
# Compute the Hessian determinant Htilde for final asymptotic formula
hes = vars[-2]*diff(H,vars[-2])*getHes(H,vars[0:-1])
hesdet = hes.determinant().factor()
# Add the Hessian, the product of the variables, and the negation of the numerator
# to the Kronecker rep. We use a straightforward but inefficient approach by
# computing another Gröbner basis!
var('T,h,g')
PolySys = [P, T - prod(vars[:-1]), h - hesdet, g+1] + C
R2 = PolynomialRing(QQ,8,'g,T,h,x,y,z,λ,u',order='lex')
GB2 = R2.ideal([R2(str(k)) for k in PolySys]).groebner_basis()
for v in [T,h,g]:
[Qv] = list(filter(lambda p: SR(p).has(v), GB2))
_,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)
# Get final asymptotic constants
d = 3
A = Q[g]/Q[λ]
B = Q[λ]^(d-1)/Q[h]*Pd^(2-d)
C = Pd/Q[T]
# Get result of algorithm
var('n')
ASM = sqrt(n*2*pi)^((1-d)) * A * sqrt(B) * SR(C)^n
print("Asymptotics are given by adding the evaluations of")
show(ASM)
print("At u = {} and u = {}".format(u1.radical_expression(),u2.radical_expression()))
# We break out the steps to help algebraic simplification
maxima_calculus('algebraic: true;')
A1 = QQbar(SR(A).subs(u==u1)).radical_expression()
B1 = QQbar(SR(B).subs(u==u1)).radical_expression()
C1 = QQbar(SR(C).subs(u==u1)).radical_expression()
A2 = QQbar(SR(A).subs(u==u2)).radical_expression()
B2 = QQbar(SR(B).subs(u==u2)).radical_expression()
C2 = QQbar(SR(C).subs(u==u2)).radical_expression()
ASMn = sqrt(n*2*pi)^((1-d))* QQbar(A1*sqrt(B1)).radical_expression() * C1^n
ASMn = ASMn + sqrt(n*2*pi)^((1-d)) * QQbar(A2*sqrt(B2)).radical_expression() * C2^n
# Print result
print("Explicitly, the dominant asymptotic behaviour is")
show(ASMn)
# Using the ore_algebra package we can heuristically check our derived asymptotics
# First, generate the first 20 terms of the diagonal (takes ~20 seconds)
cfs = QQ[x,y,z]((1/H).taylor((x,0),(y,0),(z,0),58))
LST = [cfs[k,k,k] for k in range(20)]
# Next, use ore_algebra to guess a recurrence
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
rec = guess(LST,Shift)
# Use this recurrence to quickly generate 1000 terms of diagonal
LSTlong = rec.to_list(LST[:5],1001)
# Print results
print("At n = 10, the ratio of asympt. to actual sequence is",(LSTlong[10]/ASMn.subs(n=10)).n())
print("At n = 100, the ratio of asympt. to actual sequence is",(LSTlong[100]/ASMn.subs(n=100)).n())
print("At n = 1000, the ratio of asympt. to actual sequence is",(LSTlong[1000]/ASMn.subs(n=1000)).n())