restart:
with(LinearAlgebra):with(Groebner):
Digits := 30:
hesMatrix := proc(H,r,vars)
local d,i,j,V,U,M:
d := nops(vars):
for j from 1 to d do
V[j] := 1:
for i from 1 to d do
U[i,j] := vars[i]*vars[j]*diff(H,vars[i],vars[j])/vars[d]/diff(H,vars[d]):
od:
od:
M := Matrix(d-1,d-1):
for i from 1 to d-1 do for j from 1 to d-1 do
M[i,j] := V[i]*V[j] + U[i,j] - V[j]*U[i,d] - V[i]*U[j,d] + V[i]*V[j]*U[d,d];
if i = j then
M[i,j] := M[i,j] + V[i]:
fi:
od:od:
return M:
end:AperyH := 1-w*(1+x)*(1+y)*(1+z)*(1+y+z+y*z+x*y*z);# Create extended
S := H, seq(j*diff(H,j)-lambda,j=[w,x,y,z]), subs(w=w*t,x=x*t,y=y*t,z=z*t,H):
vars := [w,x,y,z,lambda,t]:
d := 4:
# Find Grobner Basis -- each element aside from the first is a polynomial in u and one additional variable, and is linear in the additional variable. Thus, u is a separating linear form
GB := Groebner[Basis]([u-(w+t),S],plex(op(vars),u)):
P := GB[1];
Pd := diff(GB[1],u):# Get Kronecker representation
for v in vars do
EQ := op(select(has,GB,v));
Q[v] := rem(Pd*solve(EQ,v), P, u):
od:
sys := [seq(Q[v]/Pd,v=vars)]:seq(degree(Q[v]),v=vars);# The polynomials in the Kronecker representation have coefficients of bitsize ~18 digits
round(log[10](maxnorm(P)));
seq(round(log[10](maxnorm(Q[v]))),v=vars);# P factors
factor(P);# The critical points, where t=1, are given by the values of u which are roots of
Pt := gcd(P,Pd-Q[t]);# There is a critical point with positive coordinates
for uval in [solve(Pt)] do
[w,x,y,z] = evala([seq(subs(u=uval,Q[v]/Pd),v=[w,x,y,z])]); evalf(%,5);
od;# Thus, there are two values of u with t=1
# Numerically, we can see that there are no solutions with 0<t<1
# This means a point with crit point with positive real coordinates is minimal
seq(subs(u=k,Q[t]/Pd),k=fsolve(P,u));# Refine to the system without t
newP := Pt:
newPd:=diff(newP,u):
gcdex(Pd,newP,u,'invPd'):
for j in remove(has,vars,t) do
Q[j] := rem(Q[j]*newPd*invPd,newP,u):
od:
P := newP:
Pd := newPd: # The new Kronecker rep
P=0;
[w,x,y,z,lambda]=[Q[w]/Pd,Q[x]/Pd,Q[y]/Pd,Q[z]/Pd,Q[lambda]/Pd]; # Hessian Determinant
Htdet := factor(vars[4]^3*diff(H,vars[4])^3*Determinant(hesMatrix(H,[1,1,1,1],vars[1..4])));# Add -G = -1, w*x*y*z, Htdet into the Kronecker representation
GB2 := Groebner[Basis]([P,op(GB),T - w*x*y*z,ht - Htdet, g+1],plex(g,T,ht,lambda,w,x,y,z,t,u)):
for v in [g,T,ht] do
EQ := op(select(has,GB2,v));
Q[v] := rem(Pd*solve(EQ,v), P, u):
od:
Q[ht]; Q[T];# Get final asymptotic constants
A := Q[g]/Q[lambda]:
B := simplify(Q[lambda]^(d-1)/Q[ht]*Pd^(2-d)):
C := simplify(Pd/Q[T]):# Get solution of Pt corresponding to positive coordinate solution
r1 := select(a->evalf(subs(u=a,Q[x]/Pd))>0,[solve(Pt,u)]);
evalf(r1,9);# Output of the algorithm is
ASM := (n*2*Pi)^(-3/2) * A*sqrt(B) * C^n:
ASM, P, r1;simplify(evala(subs(u=op(r1),ASM)));
evalf(%,7);Example with Two Positive Critical PointsH := (1-x-y)*(20-x-40*y)-1;# Create system using linear form u=x+y+t
ff := H, seq(j*diff(H,j)-lambda,j=[x,y]), subs(x=x*t,y=y*t,H):
vars := [x,y,t,lambda]:
d := 2:
# Find Grobner Basis -- each polynomial is linear in the variables aside from u,
# so u is a separating linear form
GB := Groebner[Basis]([u-(x+t),ff],plex(lambda,x,y,t,u)):
P := GB[1];
Pd := diff(GB[1],u):# Find Kronecker representation from Grobner basis
for v in vars do
EQ := op(select(has,GB,v));
Q[v] := rem(Pd*solve(EQ,v), P, u):
od:
sys := simplify([seq(Q[v]/Pd,v=vars)]):# The critical points, where t=1, are given by the values of u which are roots of
Pt := gcd(P,Pd-Q[t]);# There are two crit points with non-negative coordinates
u1,u2 := fsolve(Pt);
for uu in [u1,u2] do
evalf(uu,5),evalf(subs(u=uu,sys[1..2]),10);
od;# Looking at the full set of solutions shows that the second point (x2,y2) is not minimal (there is a point approximately 0.0921*(x2,y2). The first point (x1,y1) is minimal (the only points on the line t(x_1,y_1) have t=1 and t=1.709...
for uu in fsolve(P) do
subs(u=uu,sys);
od;# There are no other points with the same coordinate-wise modulus
for uu in fsolve(Pt,complex) do
map(abs,subs(u=uu,sys[1..2]));
od;# Refine to the system without t
newP := Pt:
newPd:=diff(newP,u):
gcdex(Pd,newP,u,'invPd'):
for j in remove(has,vars,t) do
Q[j] := rem(Q[j]*newPd*invPd,newP,u):
od:
P := newP:
Pd := newPd: # The new Kronecker rep
P=0;
[x,y,lambda]=[Q[x]/Pd,Q[y]/Pd,Q[lambda]/Pd]; # Hessian Determinant
Htdet := factor(vars[d]^(d-1)*diff(H,vars[d])^(d-1)*Determinant(hesMatrix(H,[seq(1,k=1..d)],vars[1..d])));# Add -G = -1, w*x*y*z, Htdet into the Kronecker representation
GB2 := Groebner[Basis]([P,op(GB),T - mul(v,v=vars[1..d]),ht - Htdet, g+1],plex(g,T,ht,op(vars),u)):
for v in [g,T,ht] do
EQ := op(select(has,GB2,v));
Q[v] := rem(Pd*solve(EQ,v), P, u):
od:
Q[ht]; Q[T];# Get final asymptotic constants
A := Q[g]/Q[lambda]:
B := simplify(Q[lambda]^(d-1)/Q[ht]*Pd^(2-d)):
C := simplify(Pd/Q[T]):# Output of the algorithm is the following
# Make sure u1 is the root of P which gives (x1,y1) ~ (.5482, .3099)
ASM := (n*2*Pi)^((1-d)/2) * simplify(A)*sqrt(simplify(B)) * C^n:
'ASM' = ASM; 'P' = P; 'u'=u1;# Numerical approximation of asymptotics
evalf(subs(u=u1,ASM));