## Example 7.5 (A Non-Combinatorial Series Expansion)¶

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)

In :
# 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

In :
# Set basic info for this example
var('x,y,z,λ,t,u')
H = 1 - (x+y+z) + 81/8*x*y*z

In :
# 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()))

The ideal generated by [a1 - 1/2*lambdaR - 1/2, a2 - 1/2*lambdaR - 1/2, a3 - 1/2*lambdaR - 1/2, b1 - 1/2*lambdaI, b2 - 1/2*lambdaI, b3 - 1/2*lambdaI, lambdaR^3 + 3*lambdaR^2 + 49/27*lambdaR - 27/8*lambdaI^6 - 5*lambdaI^4 - 32/27*lambdaI^2 + 49/81, lambdaR*lambdaI - 9/8*lambdaI^5 - 5/3*lambdaI^3 + 49/81*lambdaI, lambdaI^7 + 16/9*lambdaI^5 + 64/81*lambdaI^3 - 3136/19683*lambdaI] has dimension 0.

In :
# 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])

In :
# 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()

In :
# As expected, there are a finite number of solutions
# At these solutions, the t variable is a root of the following
GBxy[-1].factor()

Out:
(1/10460353203) * (t - 3) * (t - 1) * (3*t - 1) * (t^2 + 3) * (t^2 - 3*t + 9) * (t^2 - t + 1) * (t^2 + t + 1) * (t^2 + 3*t + 3) * (3*t^2 + 1) * (3*t^2 + 3*t + 1) * (9*t^2 - 3*t + 1) * (9*t^3 + 12*t^2 + 4*t - 27) * (81*t^3 + 36*t^2 + 4*t - 9) * (81*t^6 + 108*t^5 + 108*t^4 + 534*t^3 + 340*t^2 - 108*t + 729) * (729*t^6 - 972*t^5 + 540*t^4 - 144*t^3 - 308*t^2 + 108*t + 243)
In :
# 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)))

Out:
[0.3333333333333334?, 0.3451548307798085?]
In :
# 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')
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])])

[-2/3, -2/3, -2/3]
[-2/3, -2/3, -2/3]
[-2/3, -2/3, -2/3]
[-2/3, -2/3, -2/3]
[-2/3, -2/3, -2/3]

In :
# 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

Out:
u^3 - 8/27*u + 8/81
In :
# 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)))

Out:
1
In :
# 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, " is represented by\t",Q[vars]/Pd)

In :
# 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")

The root u = -0.6666666666666667? gives the critical point
[x,y] = [-0.6666666666666667?, -0.6666666666666667?]

The root u = 0.3333333333333334? - 0.1924500897298753?*I gives the critical point
[x,y] = [0.3333333333333334? - 0.1924500897298753?*I, 0.3333333333333334? - 0.1924500897298753?*I]

The root u = 0.3333333333333334? + 0.1924500897298753?*I gives the critical point
[x,y] = [0.3333333333333334? + 0.1924500897298753?*I, 0.3333333333333334? + 0.1924500897298753?*I]


In :
# Select the minimal critical points
[u1,u2] = filter(lambda U: not sys.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)

In :
# 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)

Asymptotics are given by adding the evaluations of

At u = -1/3*I*sqrt(1/3) + 1/3 and u = 1/3*I*sqrt(1/3) + 1/3

In :
# We break out the steps to help algebraic simplification
maxima_calculus('algebraic: true;')

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)

Explicitly, the dominant asymptotic behaviour is

In :
# 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)]

In :
# 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/ASMn.subs(n=10)).n())
print("At n = 100, the ratio of asympt. to actual sequence is",(LSTlong/ASMn.subs(n=100)).n())
print("At n = 1000, the ratio of asympt. to actual sequence is",(LSTlong/ASMn.subs(n=1000)).n())

At n = 10, the ratio of asympt. to actual sequence is 0.961192243551694
At n = 100, the ratio of asympt. to actual sequence is 0.996111746232913
At n = 1000, the ratio of asympt. to actual sequence is 0.999611117301680

In [ ]: