Ask Your Question
1

Finding the minimal polynomial for a solution of a trigonometric system of quadratic equations

asked 2023-08-16 18:17:59 +0100

aobara gravatar image

updated 2023-08-18 21:45:11 +0100

Following the answers provided in a previous Post https://ask.sagemath.org/question/706..., I tried to apply the proposed methods to a similar trigonometric system of quadratic equations.

$ 8 r^2 =((1+\sqrt{3})\cos(\alpha) r+2\cos(\beta)r-\frac{1}{2})^2 +((1+\sqrt{3})\sin(\alpha)r+2 \sin(\beta) r-\frac{1}{2})^2 $

$ 4 r^2 = (2(1+\sqrt{3})\cos(\alpha)r+2\cos(\beta)r-1)^2 +(2(1+\sqrt{3})\sin(\alpha)r+2\sin(\beta)r)^2 $

$ 4 r^2 = (((1+\sqrt{3})\cos(\alpha)r+2\cos(\beta)r +(1+\sqrt{3})\sin(\alpha)r+2\sin(\beta)r)\frac{1}{2}-\frac{1}{2})^2\\ \qquad+((1-(1+\sqrt{3})\cos(\alpha)r-2\cos(\beta)r +(1+\sqrt{3})\sin(\alpha)r+2 \sin(\beta)r)\frac{1}{2}-\frac{1}{2})^2 $

I am interested in finding the minimal polynomial of the root $r$ between 0 and 1 for which a numerical approximation is $r \approx 0.1225582838364$, $\alpha \approx 0.174507496606164$ and $\beta \approx 0.417629154722126$.

Below is my (unsuccessfull) attempt to use the variety() function to solve this problem

reset()
from time import time as stime
RR.<x, y, r, s0, c0, s1, c1>=AA[]
Sys2 = [ # a
    -8*r^2+((1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1/2-x)^2+((1+2*AA(cos(pi/6)))*s0*r+2*s1*r-1/2-y)^2,
    # b
    -4*r^2+(2*(1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1-x)^2+(2*(1+2*AA(cos(pi/6)))*s0*r+2*s1*r-y)^2,
    # eq3
    -4*r^2+(((1+2*AA(cos(pi/6)))*c0*r+2*c1*r+(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)/2-1/2)^2+((1-(1+2*AA(cos(pi/6)))*c0*r- 
2*c1*r+(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)/2-1/2)^2,
    # eq4
    -8*r^2+((1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1/2)^2+((1+2*AA(cos(pi/6)))*s0*r+2*s1*r-1/2)^2,
    # eq5
    -4*r^2+(2*(1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1)^2+(2*(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)^2,
    # eq6
    -4*r^2 +(-1/2*(1 + 2*AA(cos(pi/6)) * r * (c0 - s0) + r * s1 + r * c1)^2 +(1/2*(r*(-(1+2*AA(cos(pi/6)))) * (s0 + c0) + 2 * s1 - 
2 * c1) + 1))^2,
    # Basic
    s0^2 + c0^2 - 1,
    s1^2 + c1^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))]
t2 = stime()
P = Sol[r].minpoly()
t3 = stime()
print("Minimal polynomial (found in %f6.3 seconds) :"%(t3-t2))
print(P)

Unfortunately there seems to be a crash in the Variety() call and no result is returned. It this because this particular system is more complex than in the one in the previous post (one more variable)?

In advance, thanks for your help.

kr,

EDIT

After checking the 3 equations again I am pretty sure we dont need x and y variables so I gave it another try, but still not working...

r = 0.1225582838364
a = 0.1745074966061
b = 0.4176291547221

eq1 = -8*r^2+((1+sqrt(3))*cos(a)*r+2*cos(b)*r-1/2)^2 +((1+sqrt(3))*sin(a)*r+2*sin(b)*r-1/2)^2
eq2 = -4*r^2+(2*(1+sqrt(3))*cos(a)*r+2*cos(b)*r-1)^2+(2*(1+sqrt(3))*sin(a)*r+2*sin(b)*r)^2
eq3 = -4*r^2+(((1+sqrt(3))*cos(a)*r+2*cos(b)*r+(1+sqrt(3))*sin(a)*r+2*sin(b)*r)/2-1/2)^2+((1-(1+sqrt(3))*cos(a)*r- 
2*cos(b)*r+(1+sqrt(3))*sin(a)*r+2*sin(b)*r)/2-1/2)^2

show(eq1.numerical_approx())
show(eq2.numerical_approx())
show(eq3.numerical_approx())


reset()
from time import time as stime
RR.<r, s0, c0, s1, c1>=AA[]
Sys2 = [ # eq1
    -8*r^2+((1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1/2)^2 +((1+AA(cos(-7*pi/12)))*s0*r+2*s1*r-1/2)^2,
    # eq2
    -4*r^2+(2*(1+2*AA(cos(pi/6)))*c0*r+2*c1*r-1)^2+(2*(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)^2,
    # eq3
    -4*r^2+(((1+2*AA(cos(pi/6)))*c0*r+2*c1*r+(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)/2-1/2)^2+((1- 
    (1+2*AA(cos(pi/6)))*c0*r-2*c1*r+(1+2*AA(cos(pi/6)))*s0*r+2*s1*r)/2-1/2)^2,
    # Basic
    s0^2 + c0^2 - 1,
    s1^2 + c1^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))]
t2 = stime()
P = Sol[r].minpoly()
t3 = stime()
print("Minimal polynomial (found in %f6.3 seconds) :"%(t3-t2))
print(P)
edit retag flag offensive close merge delete

Comments

1

Again, your polynomial expressions bear no visible relation to your original problem. How do you derive them ?

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-08-17 12:14:35 +0100 )edit
1

FWIW, on Sage 10.1.rc0, .variety() didn't return in more than 2 hours, but didn't raise an error either...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-08-17 12:16:10 +0100 )edit
1

What did you edit ? I worked on the initial version of yur question, no longer available...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-08-17 12:41:03 +0100 )edit

Regarding my last edit, I only added numerical approximations for $\alpha$ and $\beta$, as I previously only gave the approximation for $r$.

aobara gravatar imageaobara ( 2023-08-17 13:17:07 +0100 )edit

As you suggested, I did compare the original equations to the polynomial expressions in the sage program. Indeed, there were differences and I am pretty sure x and y variables are not needed in this case. I gave it another try and now I am pretty sure the polynomial expressions agree with the initial system of equations. Hoever it still does not produces a result.. After the variety() call nothing is displayed. .

aobara gravatar imageaobara ( 2023-08-18 22:07:43 +0100 )edit

2 Answers

Sort by ยป oldest newest most voted
3

answered 2023-08-17 19:13:18 +0100

Max Alekseyev gravatar image

updated 2023-08-17 21:49:18 +0100

You don't need to use .variety() for finding minimal polynomial here. Instead, you'd need to eliminate all variables except r by calling J1.elimination_ideal( [x, y, s0, c0, s1, c1] ) and then factoring the generator(s) of the resulted elimination ideal. Among the irreducible factors you'd need those that have a root between 0 and 1, which can be tested with the number_of_roots_in_interval() method.

UPD. It turns out that using AA is an overkill (and slow!), while QuadraticField(3) is well enough. There are 2 minimal polynomials each having 4 roots in the interval $(0,1)$:

r^30 + (2545815/1801202*a - 2449853/900601)*r^28 + (338486/900601*a + 603176/900601)*r^27 + (-446648257/152201569*a + 12043728993/2435225104)*r^26 + (1383259681/304403138*a - 1226101327/152201569)*r^25 + (-2015430969/374650016*a + 1914137211/187325008)*r^24 + (-190717225521/31657926352*a + 324665390559/31657926352)*r^23 + (2115188315493437/556419713562752*a - 7222812032062617/1112839427125504)*r^22 + (94112339836993/69552464195344*a - 597799751694513/278209856781376)*r^21 + (-253291049328919143/188069863184210176*a + 212261477857858291/94034931592105088)*r^20 + (-4781584174232475/47017465796052544*a + 9833948490455267/47017465796052544)*r^19 + (37367570590359697/376139726368420352*a - 121023099795610641/752279452736840704)*r^18 + (20999651247970581/188069863184210176*a - 18900003993155813/94034931592105088)*r^17 + (10453083474612541/376139726368420352*a - 65279351862120215/1504558905473681408)*r^16 + (-44064505165782737/752279452736840704*a + 4726239042423265/47017465796052544)*r^15 + (16227635441524291/3009117810947362816*a - 115369507693548305/12036471243789451264)*r^14 + (9663004251123147/752279452736840704*a - 33032622309102601/1504558905473681408)*r^13 + (-49709256497213141/12036471243789451264*a + 84729889946852853/12036471243789451264)*r^12 + (-5521741764471403/6018235621894725632*a + 595756732310263/376139726368420352)*r^11 + (34673497535139875/48145884975157805056*a - 7436862138499025/6018235621894725632)*r^10 + (-345885411536503/3009117810947362816*a + 4709623543156535/24072942487578902528)*r^9 + (-4590502012308313/96291769950315610112*a + 7923331959568715/96291769950315610112)*r^8 + (931535286989029/48145884975157805056*a - 798863982924283/24072942487578902528)*r^7 + (-72056983925195/48145884975157805056*a + 1981603083884441/770334159602524880896)*r^6 + (-22084555538089/48145884975157805056*a + 75321892689733/96291769950315610112)*r^5 + (603181898011/6018235621894725632*a - 16561591112283/96291769950315610112)*r^4 + (-9971674893/6018235621894725632*a + 18538673951/6018235621894725632)*r^3 + (-9264649743/6018235621894725632*a + 31621338863/12036471243789451264)*r^2 + (34988843/188069863184210176*a - 973560803/3009117810947362816)*r - 689105/94034931592105088*a + 9773497/752279452736840704

and

r^30 + (2545815/1801202*a - 2449853/900601)*r^28 + (22486/12337*a - 19960/12337)*r^27 + (-446648257/152201569*a + 12043728993/2435225104)*r^26 + (-928907615/304403138*a + 802016915/152201569)*r^25 + (1796288795/4870450208*a + 118006079/2435225104)*r^24 + (81394169511/411553042576*a - 237170434365/411553042576)*r^23 + (307425189198909/556419713562752*a - 850699208797337/1112839427125504)*r^22 + (37573295226625/69552464195344*a - 246974801001649/278209856781376)*r^21 + (-60239193343637863/188069863184210176*a + 45658652734554995/94034931592105088)*r^20 + (-549608776161751/3616728138157888*a + 1125227431684579/3616728138157888)*r^19 + (27371283258215057/376139726368420352*a - 104755137145182993/752279452736840704)*r^18 + (4560692726127805/188069863184210176*a - 3981105851513519/94034931592105088)*r^17 + (2072252269222557/376139726368420352*a - 6123796422206871/1504558905473681408)*r^16 + (-5679025703445661/752279452736840704*a + 920083549412993/94034931592105088)*r^15 + (1776887757881603/3009117810947362816*a - 5606403805021457/12036471243789451264)*r^14 + (2194126156104539/752279452736840704*a - 7191165002066257/1504558905473681408)*r^13 + (-19953398578164437/12036471243789451264*a + 32145091694720885/12036471243789451264)*r^12 + (-429591257748895/6018235621894725632*a + 56305615331019/376139726368420352)*r^11 + (18329629484328995/48145884975157805056*a - 3895636964469777/6018235621894725632)*r^10 + (-65451127131541/752279452736840704*a + 3523020195864751/24072942487578902528)*r^9 + (-2101099156041305/96291769950315610112*a + 3631888764561931/96291769950315610112)*r^8 + (37361161246969/3703529613473677312*a - 31842123800333/1851764806736838656)*r^7 + (-24177646026627/48145884975157805056*a + 649334834531225/770334159602524880896)*r^6 + (-17869643587017/48145884975157805056*a + 60925199706821/96291769950315610112)*r^5 + (446110416779/6018235621894725632*a - 12203233778267/96291769950315610112)*r^4 + (-180612269/6018235621894725632*a + 1458046367/6018235621894725632)*r^3 + (-9264649743/6018235621894725632*a + 31621338863/12036471243789451264)*r^2 + (34988843/188069863184210176*a - 973560803/3009117810947362816)*r - 689105/94034931592105088*a + 9773497/752279452736840704

where $a = \sqrt{3}$.

edit flag offensive delete link more

Comments

Thanks, and sorry I am really a beginnir with Sage... I tried Sols = J1.elimination_ideal( [x, y, s0, c0, s1, c1] ) show(Sols) but no solution is displayed.

How do you then get the the generators of the ideal and factor them?

aobara gravatar imageaobara ( 2023-08-17 19:45:10 +0100 )edit
1

Strange. I do see that the elimination ideal has 1 generator of degree 60, which factors into 2 polynomials of degree 30 each. See updated answer for details.

Max Alekseyev gravatar imageMax Alekseyev ( 2023-08-17 21:32:41 +0100 )edit
1

Like this:

E = J1.elimination_ideal( [x, y, s0, c0, s1, c1] )
f = E.gen(0).change_ring( QuadraticField(3) )
factor(f)
Max Alekseyev gravatar imageMax Alekseyev ( 2023-08-18 16:09:35 +0100 )edit

I was having some doubts so I first checked the numerical results for the 3 initial equations

var('r a b')

r = 0.1225582838364
a = 0.1745074966061
b = 0.4176291547221

eq1 = -8*r^2+((1+sqrt(3))*cos(a)*r+2*cos(b)*r-1/2)^2 +((1+sqrt(3))*sin(a)*r+2*sin(b)*r-1/2)^2
eq2 = -4*r^2+(2*(1+sqrt(3))*cos(a)*r+2*cos(b)*r-1)^2+(2*(1+sqrt(3))*sin(a)*r+2*sin(b)*r)^2
eq3 = -4*r^2+(((1+sqrt(3))*cos(a)*r+2*cos(b)*r+(1+sqrt(3))*sin(a)*r+2*sin(b)*r)/2-1/2)^2+((1- 
      (1+sqrt(3))*cos(a)*r-2*cos(b)*r+(1+sqrt(3))*sin(a)*r+2*sin(b)*r)/2-1/2)^2

show(eq1.numerical_approx())
show(eq2.numerical_approx())
show(eq3.numerical_approx())

And it did verify pretty ...(more)

aobara gravatar imageaobara ( 2023-08-18 16:21:30 +0100 )edit

I did check the roots in the polynomials you provided, and the closest I found were $0.122113263930323$ for the first one and $0.119342747874449$ for the second one. So it is not quite what is expected since the minimum polynomial should agree with infinite precision to the actural solution $r=0.1225582838364$. Something still is not quite right...

aobara gravatar imageaobara ( 2023-08-18 16:51:13 +0100 )edit
1

answered 2023-08-20 12:39:44 +0100

aobara gravatar image

updated 2023-08-20 12:53:35 +0100

I finally managed to solve the various issues 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 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 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:

$ 8 r^2 =((1+\sqrt{3})\cos(\alpha) r+2\cos(\beta)r-\frac{1}{2})^2 +((1+\sqrt{3})\sin(\alpha)r+2 \sin(\beta) r-\frac{1}{2})^2 $

$ 4 r^2 = (2(1+\sqrt{3})\cos(\alpha)r+2\cos(\beta)r-1)^2 +(2(1+\sqrt{3})\sin(\alpha)r+2\sin(\beta)r)^2 $

$ 4 r^2= ((1+\sqrt{3})\cos(\alpha)r+2(\cos(\beta)/2-\sqrt{3}/2\sin(\beta))r-1/2)^2\\ +((1+\sqrt{3})\sin(\alpha)r+2(\sqrt{3}/2\cos(\beta))+\sin(\beta))/2)r-1/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 = \frac{1}{2} \sqrt{ \frac{1}{853} (396 - 171 \sqrt{3} - \sqrt{-126516 + 74406 \sqrt{3}}} \approx 0.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 or radical expressions into purely polynomial expressions.

edit flag offensive delete link more

Comments

1

Instead of c2-1/2, c2^2+s2^2-1 you can use just s2^2 - 3/4, that is c2 is not really needed here. Also, in this case makes sense to define the ideal over QQ rather than AA. And again, .variety() is not needed here - instead it's more straightforward to use elimination_ideal:

sage: J1.elimination_ideal( [s0, c0, s1, c1, c2, s2] )
Ideal (1901849018368*r^16 - 1110499590144*r^14 + 282845642752*r^12 - 32715190272*r^10 + 2049038848*r^8 - 76393728*r^6 + 1730688*r^4 - 22032*r^2 + 117) of Multivariate Polynomial Ring in r, s0, c0, s1, c1, c2, s2 over Algebraic Real Field
sage: _.gen(0).change_ring(QQ).factor()
(218368*r^8 - 101376*r^6 + 19104*r^4 - 864*r^2 + 9) * (8709376*r^8 - 1042176*r^6 + 49504*r^4 - 1200*r^2 + 13)
Max Alekseyev gravatar imageMax Alekseyev ( 2023-08-20 17:52:20 +0100 )edit
1

The claim the RR.ideal(...) function does not work properly if the system of equations contains references to either sqrt(3) of AA(cos(pi/6)) is not quite correct. It may be faster but I do not any issue with correctness so far. You changed the ideal and thus there is no surprise you get a different solution. The system of equations in your answer is not the same as the one in your question.

Max Alekseyev gravatar imageMax Alekseyev ( 2023-08-20 17:56:24 +0100 )edit

I checked your solution based on J1.elimination_ideal( [s0, c0, s1, c1, c2, s2] ), and it works. I agree that it is more straightforward than .variety(). Regarding my claim about RR.ideal(...), let me try to be a little more precise. Whenever used in conjonction with sqrt(3) or AA(cos(pi/6)), J1.dimension() incorrectly returns 0. This is an issue because it mislead me to believe that my initial system of equations was solvable... However this is not the case.

aobara gravatar imageaobara ( 2023-08-21 10:01:09 +0100 )edit
1

If you think there is a bug in .dimension(), it's worth investigating / reporting. Can you provide a clear minimal example illustrating the issue?

Max Alekseyev gravatar imageMax Alekseyev ( 2023-08-21 17:01:25 +0100 )edit

Hmm I tried to find a simple example to replicate the issue but I could not succeed... Probably a mistake on my part. They might not be any bug in .dimension() after all. Sorry for the false alert...

aobara gravatar imageaobara ( 2023-08-27 17:38:30 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2023-08-16 18:17:59 +0100

Seen: 695 times

Last updated: Aug 20 '23