## Example 7.3 (Apéry Revisited)¶

We run our algorithm and compute the Kronecker representation "by hand" using Sage's capability to compute Gröbner bases. Note that Example 7.1 uses the Maple package of Melczer and Salvy to automatically compute this.

Requirements: None

In [1]:
# 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 [2]:
# Set basic info for the Apéry example
var('w,x,y,z,λ,t,u')
vars = [w,x,y,z,λ,t]
H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)

# Form the extended critical point system
C = [H] + [v * diff(H,v) - λ for v in vars[0:-2]]
S = C + [H.subs([v==t*v for v in vars])]
S

Out[2]:
[-(x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1)*(z + 1) + 1,
-(x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1)*(z + 1) - λ,
-(w*(x + 1)*(y + 1)*y*(z + 1)*z + (x*y*z + y*z + y + z + 1)*w*(y + 1)*(z + 1))*x - λ,
-((x*z + z + 1)*w*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*w*(x + 1)*(z + 1))*y - λ,
-((x*y + y + 1)*w*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1))*z - λ,
-(t^3*x*y*z + t^2*y*z + t*y + t*z + 1)*(t*x + 1)*(t*y + 1)*(t*z + 1)*t*w + 1]
In [3]:
# Define the multivariate polynomial ring and compute a Gröbner basis
R = PolynomialRing(QQbar,7,'w,x,y,z,λ,t,u',order='lex')
S = S + [u - (w+t)]
GB = R.ideal([R(k) for k in S]).groebner_basis()

# The Kronecker is representation parametrized by the roots of the following polynomial
P = QQ[u](GB[-1])
P

Out[3]:
u^14 + 1144*u^13 + 561458*u^12 + 153334416*u^11 + 25196311153*u^10 + 2498103299744*u^9 + 139498034620672*u^8 + 3517543585569916*u^7 + 11318250382599528*u^6 + 16555760259324420*u^5 + 19429225239935936*u^4 + 19058949700404944*u^3 + 9535423903452608*u^2 - 2094865315255056*u - 85832780580609892
In [4]:
# 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[4]:
1
In [5]:
# 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]]," divided by\t",Pd)

In [6]:
# All polynomials in the Kronecker representation have degree close to degree(P) = 14
[Q[v].degree() for v in vars]

Out[6]:
[13, 13, 12, 12, 13, 13]
In [7]:
# They have coefficients that are at most 18 digits long
# As expected, this is an *order of magnitude* smaller than the constants
# appearing in the Gröbner basis, which have constants with ~100 digits
[Q[v].norm(infinity).log(10).round() for v in vars]

Out[7]:
[18, 17, 17, 17, 17, 18]
In [8]:
max(map(lambda p: abs(QQ(p).numerator()),GB[0].coefficients())).log(10).round()

Out[8]:
99
In [9]:
# Note that P factors
P.factor()

Out[9]:
(u^2 + 162*u - 167) * (u^12 + 982*u^11 + 402541*u^10 + 88286768*u^9 + 10961079084*u^8 + 737152378392*u^7 + 21909849528196*u^6 + 91252409193628*u^5 + 194304964440524*u^4 + 317508355295408*u^3 + 441800743647348*u^2 + 511124563867704*u + 513968745991676)
In [10]:
# The critical points, where t=1, are given by the values of u which are roots of
Pt = gcd(P,Pd-Q[t])
Pt

Out[10]:
u^2 + 162*u - 167
In [11]:
# There is a critical point with positive coordinates (and one with negative coordinates)
for uval in Pt.roots(QQbar, multiplicities=False):
print("There is a critical point where the coordinates {} equal".format(vars[:-2]))
print([p.subs(u == uval) for p in sys[:-2]])

There is a critical point where the coordinates [w, x, y, z] equal
[-164.0243866176395?, -0.4142135623730951?, -0.7071067811865475?, -0.7071067811865475?]
There is a critical point where the coordinates [w, x, y, z] equal
[0.02438661763951283?, 2.414213562373095?, 0.7071067811865475?, 0.7071067811865475?]

In [12]:
# We can compute the algebraic values of t at solutions of the system
tvals = [sys[-1].subs(u == uval) for uval in P.roots(QQbar, multiplicities=False)]
print("At the solutions of the system, t takes the values {}".format(tvals))

# And we can verify none lie in (0,1)
# This means the point with positive coordinates is *minimal*
num_in_zo = list(filter(lambda c: c.is_real() and c>0 and c<1, tvals))
print("Here are the values lying in (0,1): {}".format(num_in_zo))

At the solutions of the system, t takes the values [-0.005979808060975540?, 1.000000000000000?, 1.139687851637440?, 1.580884722469457?, 2.409071027725915?, 1.000000000000000?, 0.7665952278602726? - 2.211004338122884?*I, 0.7665952278602726? + 2.211004338122884?*I, -2.040073239787960? - 0.6346682615488171?*I, -2.040073239787960? + 0.6346682615488171?*I, -0.7916404858532640? - 1.499007434480210?*I, -0.7916404858532640? + 1.499007434480210?*I, 0.5032866008950336? - 1.375523098438217?*I, 0.5032866008950336? + 1.375523098438217?*I]
Here are the values lying in (0,1): []

In [13]:
# We can simplify the system to the points where t = 1
# First, compute reciprocal of the new derivative modulo the new P
newP = Pt
newPd = diff(newP,u)
_,invPd,_ = Pd.xgcd(newP)

# Then update the Q polynomials using polynomial division
newQ = {}
for v in vars[0:-1]:
_,newQ[v] = QQ[u](Q[v] * newPd * invPd).quo_rem(newP)

# Print the new representation
print("The Kronecker representation for the critical points (where t=1 in the extended system) is given by")
show([v == newQ[v]/newPd for v in vars[0:-2]])
print("where u satisfies {}=0.".format(newP))

The Kronecker representation for the critical points (where t=1 in the extended system) is given by

where u satisfies u^2 + 162*u - 167=0.

In [14]:
# Compute the Hessian determinant Htilde for final asymptotic formula
hes = vars[-3]*diff(H,vars[-3])*getHes(H,vars[0:-2])
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 = [newP, T - w*x*y*z, h - hesdet, g+1] + S
R2 = PolynomialRing(QQ,10,'g,T,h,w,x,y,z,λ,t,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))
_,newQ[v] = QQ[u](newPd * v.subs(solve(SR(Qv),v))).quo_rem(newP)

# Print the new elements of the representation
print("The Kronecker representation for the critical points (where t=1 in the extended system) is given by")
show([v == newQ[v]/newPd for v in [T,h,g]])
print("where u satisfies {}=0.".format(newP))

The Kronecker representation for the critical points (where t=1 in the extended system) is given by

where u satisfies u^2 + 162*u - 167=0.

In [15]:
# Get final asymptotic constants
d = 4
A = newQ[g]/newQ[λ]
B = newQ[λ]^(d-1)/newQ[h]*newPd^(2-d)
C = newPd/newQ[T]

# Find root of P that gives critical point with positive coordinates
[posu] = filter(lambda uval: sys[0].subs(u == uval)>0, newP.roots(QQbar, multiplicities=False))

# Print result of algorithm
var('n')
ASM = ((n*2*pi)^(-3/2) * A * sqrt(B) * SR(C)^n).simplify()
print("Asymptotics are given by evaluating")
show(ASM)
print("At the root {} of the polynomial {}.".format(posu,newP))
print("This gives dominant asymptotic behaviour")
show(ASM.subs(u=posu))

Asymptotics are given by evaluating

At the root 1.024386617639513? of the polynomial u^2 + 162*u - 167.
This gives dominant asymptotic behaviour

In [16]:
# Because u is the root of a quadratic polynomial, we can express asymptotics in radicals
maxima_calculus('algebraic: true;') # This helps simplify roots