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

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

Out[2]:
[(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]
In [3]:
# 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

Out[3]:
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
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 one less than the degree of P
[Q[v].degree() for v in vars]

Out[6]:
[7, 7, 7, 7]
In [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]

Out[7]:
[5, 4, 5, 4]
In [8]:
# Note that P factors
P.factor()

Out[8]:
(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)
In [9]:
# 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[9]:
u^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234
In [10]:
# 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]


In [11]:
# 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?]

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

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

Out[13]:
[0.548232473567013?,
9.997110519821022?,
0.5648805044366720?,
0.5648805044366720?]
In [14]:
# 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^4 - 1211/78*u^3 + 172525/3042*u^2 - 242111/3042*u + 9161/234=0.

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

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

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

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

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