Example 7.4 (Two Critical Points with Positive Coordinates)

We run our algorithm and compute the Kronecker representation "by hand" using Sage's capability to compute Gröbner bases.
Requirements: None

# 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,λ,t,u')
vars = [x,y,λ,t]
H = (1-x-y)*(20-x-40*y)-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

[(x + 40*y - 20)*(x + y - 1) - 1,
 (2*x + 41*y - 21)*x - λ,
 (41*x + 80*y - 60)*y - λ,
 (t*x + 40*t*y - 20)*(t*x + t*y - 1) - 1]

# Define the multivariate polynomial ring and compute a Gröbner basis
R = PolynomialRing(QQbar,5,'x,y,λ,t,u',order='lex')
S = S + [u - (x+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

u^8 - 198726679/6563778*u^7 + 173956613387/511974684*u^6 - 73067829295337/39934025352*u^5 + 8519230060835315/1557426988728*u^4 - 2508167123647177/259571164788*u^3 + 655159497358069/64892791197*u^2 - 158493265458143/27323280504*u + 170155274912099/119802076056

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

1

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|In|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|Kronecker|\phantom{\verb!x!}\verb|rep.|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|variable|\phantom{\verb!x!}\verb| | x \phantom{\verb!x!}\verb|is|\phantom{\verb!x!}\verb|represented|\phantom{\verb!x!}\verb|by | \frac{899}{39} u^{7} - \frac{123578457923}{255987342} u^{6} + \frac{1313142566253}{369759494} u^{5} - \frac{10209620955321937}{778713494364} u^{4} + \frac{590457682083167}{21630930399} u^{3} - \frac{335415361956118}{10246230189} u^{2} + \frac{2753820801687041}{129785582394} u - \frac{1502824746811243}{259571164788} \phantom{\verb!x!}\verb|divided|\phantom{\verb!x!}\verb|by | 8 u^{7} - \frac{1391086753}{6563778} u^{6} + \frac{173956613387}{85329114} u^{5} - \frac{365339146476685}{39934025352} u^{4} + \frac{8519230060835315}{389356747182} u^{3} - \frac{2508167123647177}{86523721596} u^{2} + \frac{1310318994716138}{64892791197} u - \frac{158493265458143}{27323280504}$$

# All polynomials in the Kronecker representation have degree one less than the degree of P
[Q[v].degree() for v in vars]

[7, 7, 7, 7]

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

[5, 4, 5, 4]

# Note that P factors
P.factor()

(u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234) * (u^4 - 48409909/3281889*u^3 + 27671511239/511974684*u^2 - 19048879957/255987342*u + 18573875659/511974684)

# The critical points, where t=1, are given by the values of u which are roots of 
Pt = gcd(P,Pd-Q[t])
Pt

u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234

# Now there are two critical point with positive coordinates
for uval in Pt.roots(QQbar, multiplicities=False):
    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 = 1.548232473567013? gives the critical point
[x,y] = [0.548232473567013?, 0.3099773361396642?] 

The root u = 10.997110519821022? gives the critical point
[x,y] = [9.997110519821022?, 0.2527749731647175?] 

The root u = 1.490149016126495? - 0.2807916065035742?*I gives the critical point
[x,y] = [0.4901490161264949? - 0.2807916065035742?*I, 0.5808033325272964? + 0.1623547816338436?*I] 

The root u = 1.490149016126495? + 0.2807916065035742?*I gives the critical point
[x,y] = [0.4901490161264949? + 0.2807916065035742?*I, 0.5808033325272964? - 0.1623547816338436?*I] 

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

# Because one of the positive points is not minimal, there *is* a solution with t in (0,1)
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 [1.000000000000000?, 1.709936749104775?, 0.09218565523266332?, 1.000000000000000?, 0.7114300338237230? - 0.10463158979488784?*I, 0.7114300338237230? + 0.10463158979488784?*I, 1.000000000000000? + 0.?e-27*I, 1.000000000000000? + 0.?e-27*I]
Here are the values lying in (0,1): [0.09218565523266332?]

# Define the critical points with positive real coordinates 
[u1,u2] = filter(lambda r:SR(r).is_real(), Pt.roots(QQbar, multiplicities=False))
cp1 = [r.subs(u==u1) for r in sys[:2]]
cp2 = [r.subs(u==u2) for r in sys[:2]]

# Find the solutions of the extended system which represent the same critical points
alg_sys = [[p.subs(u=uval) for p in sys] for uval in P.roots(QQbar, multiplicities=False)]
cp1t = list(filter(lambda pt: cp1 == pt[:2], alg_sys))
cp2t = list(filter(lambda pt: cp2 == pt[:2], alg_sys))

# Print results
print("The ray from the origin to the critical point cp1 = {} contains {}*cp1 and {}*cp1 \n".format(cp1,cp1t[0][-1],cp1t[1][-1]))
print("The ray from the origin to the critical point cp2 = {} contains {}*cp2 and {}*cp2".format(cp2,cp2t[0][-1],cp2t[1][-1]))

The ray from the origin to the critical point cp1 = [0.548232473567013?, 0.3099773361396642?] contains 1.000000000000000?*cp1 and 1.709936749104775?*cp1 

The ray from the origin to the critical point cp2 = [9.997110519821022?, 0.2527749731647175?] contains 0.09218565523266332?*cp2 and 1.000000000000000?*cp2

# No other critical point has the same coordinate-wise modulus as the minimal critical point
[abs(sys[0].subs(u==uval)) for uval in  Pt.roots(QQbar, multiplicities=False)]

[0.548232473567013?,
 9.997110519821022?,
 0.5648805044366720?,
 0.5648805044366720?]

# 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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{35061 \, u^{3} - 203363 \, u^{2} + 381283 \, u - 234261}{12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111}, y = \frac{20982 \, u^{3} - 263923 \, u^{2} + 644528 \, u - 446047}{4 \, {\left(12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111\right)}}\right]$$

where u satisfies u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234=0.

# 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 - prod(vars[:-2]), h - hesdet, g+1] + S
R2 = PolynomialRing(QQ,8,'g,T,h,x,y,λ,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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[T = \frac{796653 \, u^{3} - 5489818 \, u^{2} + 11297717 \, u - 7320092}{78 \, {\left(12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111\right)}}, h = \frac{21576789 \, u^{3} - 92722412 \, u^{2} + 118817977 \, u - 40527214}{39 \, {\left(12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111\right)}}, g = \left(-1\right)\right]$$

where u satisfies u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234=0.

# Get final asymptotic constants
d = 2
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)==cp1[0], newP.roots(QQbar, multiplicities=False))

# Get result of algorithm
var('n')
ASM = ((n*2*pi)^((1-d)/2) * A * sqrt(B) * SR(C)^n).simplify()

# Parse Sage output to get numerical approximation of asymptotics
def numeval(e):
    E = SR(e)
    if E.is_constant():
        return E.n()
    else:
        return E
numASM = prod([numeval(k) for k in ASM.subs(u=posu).expand().operands()])

# Print result
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(numASM)

Asymptotics are given by evaluating

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\sqrt{39} {\left(12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111\right)} \left(\frac{78 \, {\left(12168 \, u^{3} - 141687 \, u^{2} + 345050 \, u - 242111\right)}}{796653 \, u^{3} - 5489818 \, u^{2} + 11297717 \, u - 7320092}\right)^{n} \sqrt{\frac{588627 \, u^{3} - 2845362 \, u^{2} + 4562963 \, u - 2409968}{21576789 \, u^{3} - 92722412 \, u^{2} + 118817977 \, u - 40527214}}}{{\left(588627 \, u^{3} - 2845362 \, u^{2} + 4562963 \, u - 2409968\right)} \sqrt{\pi n}}$$

At the root 1.548232473567013? of the polynomial u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234.
This gives dominant asymptotic behaviour

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{0.0545997615256984 \cdot 5.884442204019508?^{n}}{\sqrt{n}}$$

# Because u is the root of a quadric polynomial, we could express asymptotics in radicals
# Those brave enough to view the output can uncomment and run the next line of code
# print(ASM.subs(u=posu.radical_expression()).simplify_full())