Rational parametrization for B, z and t, and the associated ideal¶
In [1]:
Pol.<z, t, B, V, W> = QQ[]
In [2]:
sys = [V*(1+W)*(1-2*V-V^2*W)-B*(V-2*V^2)*(1+V*W),V-2*V^2-t,W*(1-2*V-V^2*W)-z*(1+W)]
In [3]:
Id = ideal(sys)
Groebner basis to eliminate V and W¶
In [4]:
El = Id.elimination_ideal([V,W])
El
Out[4]:
Ideal (16*z^2*t^5*B^4 - 32*z^3*t^4*B^3 - 32*z^2*t^5*B^3 + 16*z^2*t^4*B^4 - 32*z*t^5*B^4 + 16*z^4*t^3*B^2 + 32*z^3*t^4*B^2 + 24*z^2*t^5*B^2 - 32*z^3*t^3*B^3 + 32*z^2*t^4*B^3 + 64*z*t^5*B^3 - 32*z*t^4*B^4 + 16*t^5*B^4 - 8*z^3*t^4*B - 8*z^2*t^5*B + 16*z^4*t^2*B^2 - 60*z^2*t^4*B^2 - 40*z*t^5*B^2 + 64*z^2*t^3*B^3 + 32*z*t^4*B^3 - 32*t^5*B^3 + 16*t^4*B^4 + z^2*t^5 + 12*z^3*t^3*B + 20*z^2*t^4*B + 8*z*t^5*B - 32*z^3*t^2*B^2 - 69*z^2*t^3*B^2 + 12*z*t^4*B^2 + 16*t^5*B^2 - 32*z*t^3*B^3 - 32*t^4*B^3 + 4*z^2*t^4 + 21*z^3*t^2*B + 9*z^2*t^3*B - 12*z*t^4*B + 15*z^2*t^2*B^2 + 53*z*t^3*B^2 + 16*t^4*B^2 + 6*z^2*t^3 + z^3*t*B - 20*z^2*t^2*B - 21*z*t^3*B + z*t^2*B^2 + 4*z^2*t^2 - z^2*t*B - z*t^2*B + z^2*t) of Multivariate Polynomial Ring in z, t, B, V, W over Rational Field
In [5]:
factor(El.gen(0))
Out[5]:
t * (t + 1) * (16*z^2*t^3*B^4 - 32*z^3*t^2*B^3 - 32*z^2*t^3*B^3 - 32*z*t^3*B^4 + 16*z^4*t*B^2 + 32*z^3*t^2*B^2 + 24*z^2*t^3*B^2 + 64*z^2*t^2*B^3 + 64*z*t^3*B^3 + 16*t^3*B^4 - 8*z^3*t^2*B - 8*z^2*t^3*B - 32*z^3*t*B^2 - 84*z^2*t^2*B^2 - 40*z*t^3*B^2 - 32*z*t^2*B^3 - 32*t^3*B^3 + z^2*t^3 + 20*z^3*t*B + 28*z^2*t^2*B + 8*z*t^3*B + 15*z^2*t*B^2 + 52*z*t^2*B^2 + 16*t^3*B^2 + 3*z^2*t^2 + z^3*B - 19*z^2*t*B - 20*z*t^2*B + z*t*B^2 + 3*z^2*t - z^2*B - z*t*B + z^2)
We see that there is a global factor t*(t+1) that we can remove
In [6]:
pol = El.gen(0)//(t*(t+1))
pol
Out[6]:
16*z^2*t^3*B^4 - 32*z^3*t^2*B^3 - 32*z^2*t^3*B^3 - 32*z*t^3*B^4 + 16*z^4*t*B^2 + 32*z^3*t^2*B^2 + 24*z^2*t^3*B^2 + 64*z^2*t^2*B^3 + 64*z*t^3*B^3 + 16*t^3*B^4 - 8*z^3*t^2*B - 8*z^2*t^3*B - 32*z^3*t*B^2 - 84*z^2*t^2*B^2 - 40*z*t^3*B^2 - 32*z*t^2*B^3 - 32*t^3*B^3 + z^2*t^3 + 20*z^3*t*B + 28*z^2*t^2*B + 8*z*t^3*B + 15*z^2*t*B^2 + 52*z*t^2*B^2 + 16*t^3*B^2 + 3*z^2*t^2 + z^3*B - 19*z^2*t*B - 20*z*t^2*B + z*t*B^2 + 3*z^2*t - z^2*B - z*t*B + z^2
Change of variables to simplify expressions
In [7]:
trans = pol(t=(1-t)/8).numerator()
trans
Out[7]:
-16*z^2*t^3*B^4 - 256*z^3*t^2*B^3 + 32*z^2*t^3*B^3 + 48*z^2*t^2*B^4 + 32*z*t^3*B^4 - 1024*z^4*t*B^2 + 256*z^3*t^2*B^2 - 24*z^2*t^3*B^2 + 512*z^3*t*B^3 + 416*z^2*t^2*B^3 - 64*z*t^3*B^3 - 48*z^2*t*B^4 - 96*z*t^2*B^4 - 16*t^3*B^4 - 64*z^3*t^2*B + 8*z^2*t^3*B + 1024*z^4*B^2 + 1536*z^3*t*B^2 - 600*z^2*t^2*B^2 + 40*z*t^3*B^2 - 256*z^3*B^3 - 928*z^2*t*B^3 - 64*z*t^2*B^3 + 32*t^3*B^3 + 16*z^2*B^4 + 96*z*t*B^4 + 48*t^2*B^4 - z^2*t^3 - 1152*z^3*t*B + 200*z^2*t^2*B - 8*z*t^3*B - 1792*z^3*B^2 + 312*z^2*t*B^2 + 296*z*t^2*B^2 - 16*t^3*B^2 + 480*z^2*B^3 + 320*z*t*B^3 - 96*t^2*B^3 - 32*z*B^4 - 48*t*B^4 + 27*z^2*t^2 + 1728*z^3*B + 792*z^2*t*B - 136*z*t^2*B + 312*z^2*B^2 - 776*z*t*B^2 + 48*t^2*B^2 - 192*z*B^3 + 96*t*B^3 + 16*B^4 - 243*z^2*t - 1512*z^2*B + 360*z*t*B + 440*z*B^2 - 48*t*B^2 - 32*B^3 + 729*z^2 - 216*z*B + 16*B^2
Expressing B as a solution of the simplified polynomial equation on B, z, t¶
In [8]:
z = polygen(QQ, 'z')
trans = Frac(QQ['z'])['t']['B'](trans)
Using DynPuiseux to expand B as a Puiseux series in t, with z as a parameter
In [9]:
from puiseux import puiseux
In [10]:
puiseux(trans, 4)
Out[10]:
[alg0 + (((z^2 - 5/8*z)/(z^2 - 1/2*z + 1/16))*alg0 + 9/16*z/(z^2 - 1/2*z + 1/16))*t + alg1*t^(3/2) + (((z^4 - 9/8*z^3 + 47/128*z^2 - 7/128*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*alg0 + (19/32*z^3 - 55/256*z^2 + 11/256*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*t^2 + O(t^(5/2))]
In [11]:
res = _
In [12]:
parent(res[0])
Out[12]:
Puiseux Series Ring in t over Univariate Quotient Polynomial Ring in alg1 over Univariate Quotient Polynomial Ring in alg0 over Fraction Field of Univariate Polynomial Ring in z over Rational Field with modulus (z - 1)*alg0^2 + (-8*z^2 + 7*z + 1)*alg0 - 27/4*z with modulus (1024*z^4 - 1792*z^3 + 960*z^2 - 208*z + 16)*alg1^2 - 64*z^2
In [13]:
res[0]
Out[13]:
alg0 + (((z^2 - 5/8*z)/(z^2 - 1/2*z + 1/16))*alg0 + 9/16*z/(z^2 - 1/2*z + 1/16))*t + alg1*t^(3/2) + (((z^4 - 9/8*z^3 + 47/128*z^2 - 7/128*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*alg0 + (19/32*z^3 - 55/256*z^2 + 11/256*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*t^2 + O(t^(5/2))
Switching to symbolic variables to solve the equations in alg0 and alg1, so we can select the solution with nonnegative coefficients
In [14]:
var('z t alg0 alg1')
polas=alg0 + (((z^2 - 5/8*z)/(z^2 - 1/2*z + 1/16))*alg0 + 9/16*z/(z^2 - 1/2*z + 1/16))*t + alg1*t^(3/2) + (((z^4 - 9/8*z^3 + 47/128*z^2 - 7/128*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*alg0 + (19/32*z^3 - 55/256*z^2 + 11/256*z)/(z^4 - z^3 + 3/8*z^2 - 1/16*z + 1/256))*t^2
In [15]:
sols=solve( [(z - 1)*alg0^2 + (-8*z^2 + 7*z + 1)*alg0 - 27/4*z, (1024*z^4 - 1792*z^3 + 960*z^2 - 208*z + 16)*alg1^2 - 64*z^2],alg0,alg1)
In [16]:
sols
Out[16]:
[[alg0 == 1/2*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), alg1 == -2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))], [alg0 == 1/2*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), alg1 == -2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))], [alg0 == 1/2*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), alg1 == 2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))], [alg0 == 1/2*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), alg1 == 2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))]]
In [17]:
polsols=[polas.substitute(sols[i]) for i in range(4)]
polsols
Out[17]:
[t^2*((152*z^3 - 55*z^2 + 11*z)/(256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1) + (128*z^4 - 144*z^3 + 47*z^2 - 7*z)*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/((256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1)*(z - 1))) + t*((8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)*(8*z^2 - 5*z)/((16*z^2 - 8*z + 1)*(z - 1)) + 9*z/(16*z^2 - 8*z + 1)) - 2*t^(3/2)*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1)) + 1/2*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), t^2*((152*z^3 - 55*z^2 + 11*z)/(256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1) + (128*z^4 - 144*z^3 + 47*z^2 - 7*z)*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/((256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1)*(z - 1))) + t*((8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)*(8*z^2 - 5*z)/((16*z^2 - 8*z + 1)*(z - 1)) + 9*z/(16*z^2 - 8*z + 1)) - 2*t^(3/2)*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1)) + 1/2*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), t^2*((152*z^3 - 55*z^2 + 11*z)/(256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1) + (128*z^4 - 144*z^3 + 47*z^2 - 7*z)*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/((256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1)*(z - 1))) + t*((8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)*(8*z^2 - 5*z)/((16*z^2 - 8*z + 1)*(z - 1)) + 9*z/(16*z^2 - 8*z + 1)) + 2*t^(3/2)*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1)) + 1/2*(8*z^2 - sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1), t^2*((152*z^3 - 55*z^2 + 11*z)/(256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1) + (128*z^4 - 144*z^3 + 47*z^2 - 7*z)*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/((256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1)*(z - 1))) + t*((8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)*(8*z^2 - 5*z)/((16*z^2 - 8*z + 1)*(z - 1)) + 9*z/(16*z^2 - 8*z + 1)) + 2*t^(3/2)*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1)) + 1/2*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1)]
We look at the expansion in z of 2 coefficients in t, to select the only nonnegative series
In [18]:
ps=[polsols[i].coefficient(t,0).series(z,5) for i in range(4)]
ps
Out[18]:
[27/4*z + (-27/16)*z^2 + (-81/32)*z^3 + (-1107/256)*z^4 + Order(z^5), 1 + 5/4*z + 27/16*z^2 + 81/32*z^3 + 1107/256*z^4 + Order(z^5), 27/4*z + (-27/16)*z^2 + (-81/32)*z^3 + (-1107/256)*z^4 + Order(z^5), 1 + 5/4*z + 27/16*z^2 + 81/32*z^3 + 1107/256*z^4 + Order(z^5)]
In [19]:
pps=[polsols[i].coefficient(t,3/2).series(z,5) for i in range(4)]
pps
Out[19]:
[2*z + 13*z^2 + 267/4*z^3 + 2521/8*z^4 + Order(z^5), 2*z + 13*z^2 + 267/4*z^3 + 2521/8*z^4 + Order(z^5), (-2)*z + (-13)*z^2 + (-267/4)*z^3 + (-2521/8)*z^4 + Order(z^5), (-2)*z + (-13)*z^2 + (-267/4)*z^3 + (-2521/8)*z^4 + Order(z^5)]
In [20]:
sol=polsols[1]
sol
Out[20]:
t^2*((152*z^3 - 55*z^2 + 11*z)/(256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1) + (128*z^4 - 144*z^3 + 47*z^2 - 7*z)*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/((256*z^4 - 256*z^3 + 96*z^2 - 16*z + 1)*(z - 1))) + t*((8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)*(8*z^2 - 5*z)/((16*z^2 - 8*z + 1)*(z - 1)) + 9*z/(16*z^2 - 8*z + 1)) - 2*t^(3/2)*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1)) + 1/2*(8*z^2 + sqrt(4*z^2 - 5*z + 1)*(4*z - 1) - 7*z - 1)/(z - 1)
Computing the asymptotics of B in t, then in z¶
We first look at the asymptotics in t, near the dominant singularity t=1/8
In [21]:
sol.coefficient(t,3/2)
Out[21]:
-2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))
Careful about the convention for 'SingularityAnalysis'!!! (t^(3/2) : argument=-3/2)
In [22]:
asymptotic_expansions.SingularityAnalysis('n', 1/8, -3/2, precision=1)
Out[22]:
3/4/sqrt(pi)*8^n*n^(-5/2) + O(8^n*n^(-7/2))
In [23]:
asymptT=_.truncate(1).exact_part()
asymptTz=asymptT*sol.coefficient(t,3/2)
asymptTz
Out[23]:
(-3/2*z/(sqrt(pi)*sqrt(4*z^2 - 5*z + 1)*(4*z - 1)))*8^n*n^(-5/2)
We now make use of AsymptoticRing to compute the asymptotics in z (not very efficient but simple to write)
In [24]:
A.<v> = AsymptoticRing(growth_group='v^QQ', coefficient_ring=SR); A
Out[24]:
Asymptotic Ring <v^QQ> over Symbolic Ring
In [25]:
Rat = Frac(QQ['z'])
In [26]:
f=sol.coefficient(t,3/2)
f
Out[26]:
-2*z/(sqrt(4*z^2 - 5*z + 1)*(4*z - 1))
In [27]:
asy = Rat(f^2)(1/4-1/v)^(1/2)
asy
Out[27]:
1/8*sqrt(1/3)*v^(3/2) - 7/12*sqrt(1/3)*v^(1/2) + 5/12*sqrt(1/3)*v^(-1/2) - 23/54*sqrt(1/3)*v^(-3/2) + 155/324*sqrt(1/3)*v^(-5/2) - 91/162*sqrt(1/3)*v^(-7/2) + 329/486*sqrt(1/3)*v^(-9/2) - 605/729*sqrt(1/3)*v^(-11/2) + 1001/972*sqrt(1/3)*v^(-13/2) - 50765/39366*sqrt(1/3)*v^(-15/2) + 192049/118098*sqrt(1/3)*v^(-17/2) - 121771/59049*sqrt(1/3)*v^(-19/2) + 2792335/1062882*sqrt(1/3)*v^(-21/2) - 5356309/1594323*sqrt(1/3)*v^(-23/2) + 6871825/1594323*sqrt(1/3)*v^(-25/2) - 8840510/1594323*sqrt(1/3)*v^(-27/2) + 136805035/19131876*sqrt(1/3)*v^(-29/2) - 9821575/1062882*sqrt(1/3)*v^(-31/2) + 3089867495/258280326*sqrt(1/3)*v^(-33/2) - 6010286975/387420489*sqrt(1/3)*v^(-35/2) + O(v^(-37/2))
In [28]:
SR(asy.truncate(1).exact_part())
Out[28]:
1/8*sqrt(1/3)*v^(3/2)
In [29]:
asymptotic_expansions.SingularityAnalysis('m', 1/4, 3/2, precision=1)
Out[29]:
2/sqrt(pi)*4^m*m^(1/2) + O(4^m*m^(-1/2))
In [30]:
simplify(2/sqrt(pi)*sqrt(1/3)*3/4/sqrt(pi))
Out[30]:
1/2*sqrt(3)/pi
In [ ]:
In [ ]: