Loading [MathJax]/jax/output/HTML-CSS/jax.js
Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 1 year ago

aobara gravatar image

I finally managed to solve the various issues I encountered with this problem. After a while I realized that the RR.ideal(...) function does not work properly if the system of equations contains references to either sqrt(3) of AA(cos(pi/6)). Instead, everything works perfectly with the introduction of additional 2 variables 'c2, s2', and 2 additionnal purely polynomial equations c2-1/2, c2^2+s2^2-1.

Previously, with the use of sqrt(3) of AA(cos(pi/6)), J1=RR.ideal(...) took a very long time but did not crash. Then the call Sols = J1.variety() returned Solution dimension = 0. Well it turns out this is not correct... After substituting with additional variables c2 and s2 it returned Solution dimension = 1, which is an indication that the initial system of equation is not solvable. Indeed, it turns out that equation 1) and 3) are equivalent. So I had to replace equation 3) with another independant equation as follows:

8r2=((1+3)cos(α)r+2cos(β)r12)2+((1+3)sin(α)r+2sin(β)r12)2

4r2=(2(1+3)cos(α)r+2cos(β)r1)2+(2(1+3)sin(α)r+2sin(β)r)2

4r2=((1+3)cos(α)r+2(cos(β)/23/2sin(β))r1/2)2+((1+3)sin(α)r+2(3/2cos(β))+sin(β))/2)r1/2)2

Updating the initial code yields

reset()
from time import time as stime
RR.<r, s0, c0, s1, c1, c2, s2>=AA[]
Sys2 = [ # eq1
    -8*r^2+((1+2*s2)*c0*r+2*c1*r-1/2)^2 +((1+2*s2)*s0*r+2*s1*r-1/2)^2,
    # eq2
    -4*r^2+(2*(1+2*s2)*c0*r+2*c1*r-1)^2+(2*(1+2*s2)*s0*r+2*s1*r)^2,
    # eq3
    -4*r^2+((1+2*s2)*c0*r+2*(c1*c2-s1*s2)*r-1/2)^2
          +((1+2*s2)*s0*r+2*(c1*s2+s1*c2)*r-1/2)^2,
    # Basic
    s0^2 + c0^2 - 1,
    s1^2 + c1^2 -1, 
    c2 - 1/2,
    s2^2+c2^2-1]
J1 = RR.ideal(Sys2)
print("Solution dimension = ", J1.dimension())
t0 = stime()
Sols = J1.variety()
t1 = stime()
print("Number of solutions = ", len(Sols), ", found in %f6.3 seconds."%(t1-t0))
Errs = [abs(s[r]-0.1225582838364) for s in Sols]
Sol = Sols[Errs.index(min(Errs))]
show(Sols)
t2 = stime()
print(Sol[r])
P = Sol[r].minpoly()
t3 = stime()
print("Minimal polynomial (found in %f6.3 seconds) :"%(t3-t2))
print(P)

The runtime is incredibly fast and returns the expected minimal polynomial

Solution dimension =  0
Number of solutions =  8 , found in 1.0070516.3 seconds.
Minimal polynomial (found in 0.0939966.3 seconds) :
x^8 - 396/853*x^6 + 597/6824*x^4 - 27/6824*x^2 + 9/218368

Which turns out to be biquadratic of degree 8, which means: I also get a closed form solution for r!!!

r=121853(3961713126516+7440630.1225582838364

And that gets far beyond what I was expecting!!!! ;)

Thanks to everybody for your advises and efforts!!!

@Emmannuel charpentier: you were right all along, there is nothing wrong with variety(), if you convert trigonometric expressions or radical expressions into purely polynomial expressions.

click to hide/show revision 2
No.2 Revision

I finally managed to solve the various issues I encountered with this problem. After a while I realized that the RR.ideal(...) function does not work properly if the system of equations contains references to either sqrt(3) of AA(cos(pi/6)). Instead, everything works perfectly fine with the introduction of additional 2 variables 'c2, s2', and 2 additionnal purely polynomial equations c2-1/2, c2^2+s2^2-1.

Previously, with the use of sqrt(3) of AA(cos(pi/6)), J1=RR.ideal(...) took a very long time but did not crash. Then the call Sols = J1.variety() returned Solution dimension = 0. Well it turns out this is not correct... After substituting with additional variables c2 and s2 it returned Solution dimension = 1, which is an indication that the initial system of equation is not solvable. Indeed, it turns out that equation 1) and 3) are equivalent. So I had to replace equation 3) with another independant equation as follows:

8r2=((1+3)cos(α)r+2cos(β)r12)2+((1+3)sin(α)r+2sin(β)r12)2

4r2=(2(1+3)cos(α)r+2cos(β)r1)2+(2(1+3)sin(α)r+2sin(β)r)2

4r2=((1+3)cos(α)r+2(cos(β)/23/2sin(β))r1/2)2+((1+3)sin(α)r+2(3/2cos(β))+sin(β))/2)r1/2)2

Updating the initial code yields

reset()
from time import time as stime
RR.<r, s0, c0, s1, c1, c2, s2>=AA[]
Sys2 = [ # eq1
    -8*r^2+((1+2*s2)*c0*r+2*c1*r-1/2)^2 +((1+2*s2)*s0*r+2*s1*r-1/2)^2,
    # eq2
    -4*r^2+(2*(1+2*s2)*c0*r+2*c1*r-1)^2+(2*(1+2*s2)*s0*r+2*s1*r)^2,
    # eq3
    -4*r^2+((1+2*s2)*c0*r+2*(c1*c2-s1*s2)*r-1/2)^2
          +((1+2*s2)*s0*r+2*(c1*s2+s1*c2)*r-1/2)^2,
    # Basic
    s0^2 + c0^2 - 1,
    s1^2 + c1^2 -1, 
    c2 - 1/2,
    s2^2+c2^2-1]
J1 = RR.ideal(Sys2)
print("Solution dimension = ", J1.dimension())
t0 = stime()
Sols = J1.variety()
t1 = stime()
print("Number of solutions = ", len(Sols), ", found in %f6.3 seconds."%(t1-t0))
Errs = [abs(s[r]-0.1225582838364) for s in Sols]
Sol = Sols[Errs.index(min(Errs))]
show(Sols)
t2 = stime()
print(Sol[r])
P = Sol[r].minpoly()
t3 = stime()
print("Minimal polynomial (found in %f6.3 seconds) :"%(t3-t2))
print(P)

The runtime is incredibly fast and returns the expected minimal polynomial

Solution dimension =  0
Number of solutions =  8 , found in 1.0070516.3 seconds.
Minimal polynomial (found in 0.0939966.3 seconds) :
x^8 - 396/853*x^6 + 597/6824*x^4 - 27/6824*x^2 + 9/218368

Which turns out to be biquadratic of degree 8, which means: I also get a closed form solution for r!!!

r=121853(3961713126516+7440630.1225582838364

And that gets far beyond what I was expecting!!!! expecting ;)

Thanks to everybody for your advises and efforts!!!efforts!

@Emmannuel charpentier: you were right all along, there is nothing wrong with variety(), if you convert trigonometric expressions or radical expressions into purely polynomial expressions.

expressions.