On the Sympy list, Cris Smith points out <https://groups.google.com/g/sympy/c/EB_Z6h3ZRDg/m/DuAIQkquBAAJ> that factoring the orignal equations is enough to allow Sympy’s solve go get a correct solution. It turns out that he’s right for Sympy 1.9 (current), but not for Sympy 1.8 (current in Sage 9.5.beta7).
Sorry for the noise… Le dimanche 5 décembre 2021 à 20:55:32 UTC+1, Emmanuel Charpentier a écrit : > ask.sagemat.org question > <https://ask.sagemath.org/question/59063/weird-c-values-from-solving-system-of-equations/> > > demonstrating a problem common to all free equation solvers : solve > > $$ > \begin{align > > > > > *}-a{1}^{3} a{2} + a{1} a{2}^{2} \ -3 \, a{1}^{2} a{2} b{1} + 2 \, a{1} > a{2} b{2} + a{2}^{2} b{2} - a{1} b{2} \ -a{1}^{2} a{2}^{2} + a{2}^{3} \ -2 > \, a{1}^{2} a{2} b{2} - 2 \, a{2}^{2} b{1} + 3 \, a{2}^{2} b{2}\end{align* > } > $$ > > Unk = var('a1 a2 b1 b2') > eq1 = a1 * a2^2 - a2 * a1^3 > eq2 = 2*a1*a2*b2 + b2*a2^2 - 3*a2*a1^2*b1 - a1*b2 > eq3 = a2^3 - a2^2*a1^2 > eq4 = 3*a2^2*b2 - 2*a2*a1^2*b2 - 2*a2^2*b1 > Sys = [eq1, eq2, eq3, eq4] > > Solving eq1 for a2 : > > S1 = eq1.factor().solve(a2, solution_dict=True) ;print ("S1 = ", S1) > > S1 = [{a2: a1^2}, {a2: 0}] > > gives us a set of solutions which are *also* solutions of eq3 : > > [eq3.subs(s) for s in S1] > > [0, 0] > > It easy to check that there are no ther roots to eq3: > > flatten([eq3.solve(v, algorithm="sympy") for v in eq3.variables()]) > > [a1 == sqrt(a2), a1 == -sqrt(a2), a2 == 0, a2 == a1^2] > > S1[1] turns out to be *also* a solution of eq4: > > eq4.subs(S1[1]) > > 0 > > As above, we can substitute S1[0] in eq4 and merge the resulting > solutions of eq4 to S1[0] suitably updated: > > E4=eq4.subs(S1[0]) > S4=flatten([E4.solve(v, solution_dict=True) for v in E4.variables()]) > S134t=S1[1:]for s in S4: > S0={u:S1[0][u].subs(s) if "subs" in dir(S1[0][u]) else S1[0][u] for u in > S1[0].keys()} > S134t+=[S0.copy()|s] > S134t > > [{a2: 0}, {a2: 0, a1: 0}, {a2: a1^2, b1: 1/2*b2}, {a2: a1^2, b2: 2*b1}] > > This set being redundant, we trim it manually : > > S134=[S134t[u] for u in (0, 3)] > S134 > > [{a2: 0}, {a2: a1^2, b2: 2*b1}] > > Substituting these solutions in eq2 gives us a set of equations, whose > solutions for their variables complete the partial solutions of S134, thus > giving the set of solutions of the complete system : > > S1234=[]for s in S134: > E=eq2.subs(s) > S=flatten([E.solve(v, solution_dict=True) for v in E.variables()]) > for s1 in S: > s0={u:s[u].subs(s1) if "subs" in dir(s[u]) else s[u] for u in > s.keys()} > S1234+=[s0.copy()|s1] > S1234 > > [{a2: 0, a1: 0}, > {a2: 0, b2: 0}, > {a2: 1/324*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + > 16*(-I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^2, > b2: 2*b1, > a1: -1/2*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) - > 8/9*(-I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 4/3}, > {a2: 1/324*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) + > 16*(I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^2, > b2: 2*b1, > a1: -1/2*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) - > 8/9*(I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 4/3}, > {a2: 1/81*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + > 16/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 12)^2, > b2: 2*b1, > a1: (1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 16/9/(1/9*I*sqrt(101)*sqrt(3) > + 37/27)^(1/3) + 4/3}, > {a2: 0, b2: 2*b1, a1: 0}, > {a2: a1^2, b2: 0, b1: 0}] > > Let’s check these solution by substitution in the original system: > > Chk=[[e.subs(s) for e in Sys] for s in S1234] ; Chk > > [[0, 0, 0, 0], > [0, 0, 0, 0], > [0, > -1/104976*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + > 16*(-I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^4*b1 - > 1/1458*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + > 16*(-I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^3*b1 + > 1/9*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + > 16*(-I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)*b1, > 0, > 0], > [0, > -1/104976*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) + > 16*(I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^4*b1 - > 1/1458*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) + > 16*(I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)^3*b1 + > 1/9*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) + > 16*(I*sqrt(3) + 1)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 24)*b1, > 0, > 0], > [0, > -1/6561*b1*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + > 16/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 12)^4 + > 4/729*b1*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + > 16/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 12)^3 - > 2/9*b1*(9*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + > 16/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 12), > 0, > 0], > [0, 0, 0, 0], > [0, 0, 0, 0]] > > The fly in the ointment is that once substituted, for some solutions, eq2 > is a first-degree monomial in b1, whose coefficient is in each case > numerically quite close to 0 : > > K=[Chk[u][1].coefficient(b1) for u in range(2,5)] > [u.n() for u in K] > > [6.66133814775094e-16 + 4.44089209850063e-16*I, > -8.88178419700125e-16 - 4.44089209850063e-16*I, > -2.84217094304040e-14] > > but cannot be proven to be 0 ([u.is_zero() for u in K] never returns). > > Furthermore [u.factor().n() for u in K] aborts (on Sage 9.5.beta7). > > But Sympy *seems* to be able to prove that these coefficients are 0 : > > [u._sympy_().factor()._sage_() for u in K] > > [0, 0, 0] > > Using algorithm="sympy" to get the solutions to eq2 leads to *different* > solutions for the quadrinomial in a1, expressed as a trigonometric > expression : > > S1234s=[]for s in S134: > E=eq2.subs(s) > S=flatten([E.solve(v, solution_dict=True, algorithm="sympy") for v in > E.variables()]) > for s1 in S: > s0={u:s[u].subs(s1) if "subs" in dir(s[u]) else s[u] for u in > s.keys()} > S1234s+=[s0.copy()|s1] > S1234s > > [{a2: 0, a1: 0}, > {a2: 0, b2: 0}, > {a2: 0, b2: 2*b1, a1: 0}, > {a2: 16/9*(2*cos(1/3*arctan(3/37*sqrt(303))) + 1)^2, > b2: 2*b1, > a1: 8/3*cos(1/3*arctan(3/37*sqrt(303))) + 4/3}, > {a2: (-2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) + > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^2, > b2: 2*b1, > a1: -2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) + > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3}, > {a2: (2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) - > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^2, > b2: 2*b1, > a1: 2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) - > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3}, > {a2: a1^2, b2: 0, b1: 0}] > > However, these expressions lead againt to first-order monomials in b1, > again unprovably null : > > Chks=[[e.subs(s) for e in Sys] for s in S1234s] ; Chks > > [[0, 0, 0, 0], > [0, 0, 0, 0], > [0, 0, 0, 0], > [0, > -256/81*b1*(2*cos(1/3*arctan(3/37*sqrt(303))) + 1)^4 + > 256/27*b1*(2*cos(1/3*arctan(3/37*sqrt(303))) + 1)^3 - > 8/3*b1*(2*cos(1/3*arctan(3/37*sqrt(303))) + 1), > 0, > 0], > [0, > -(-2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) + > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^4*b1 + > 4*(-2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) + > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^3*b1 + > 4/9*(3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) - > 3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) - > 4*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 4*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) > + 37/27)^(1/3) + 4*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 4*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 3*cos(1/3*arctan(3/37*sqrt(303))) + > 3*I*sin(1/3*arctan(3/37*sqrt(303))) - 6)*b1, > 0, > 0], > [0, > -(2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) - > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^4*b1 + > 4*(2/3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) - > 8/9*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - > 8/9*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 8/9*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 8/9*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 2/3*cos(1/3*arctan(3/37*sqrt(303))) - > 2/3*I*sin(1/3*arctan(3/37*sqrt(303))) + 4/3)^3*b1 + > 4/9*(-3*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303))) + > 3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 4*I*sqrt(3)*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 4*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) > + 37/27)^(1/3) + 4*cos(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) - 4*I*sin(1/3*arctan(3/37*sqrt(303)))/abs(1/9*I*sqrt(303) + > 37/27)^(1/3) + 3*cos(1/3*arctan(3/37*sqrt(303))) + > 3*I*sin(1/3*arctan(3/37*sqrt(303))) - 6)*b1, > 0, > 0], > [0, 0, 0, 0]] > > Again, the numerical values are close to 0 : > > Ks=[Chks[u][1].coefficient(b1) for u in range(2,5)] > [u.n() for u in Ks] > > [0.000000000000000, 0.000000000000000, -8.88178419700125e-16] > > But various attempts to prove the nullity at least partially fail : > > print([u._sympy_().simplify()._sage_() for u in Ks]) > print([u._sympy_().factor().simplify().is_zero for u in Ks]) > > [0, 0, -256/81*(2*sin(1/6*pi - 1/3*arctan(3/37*sqrt(303))) - 1)^4 - > 256/27*(2*sin(1/6*pi - 1/3*arctan(3/37*sqrt(303))) - 1)^3 - > 8/3*sqrt(3)*sin(1/3*arctan(3/37*sqrt(303))) + > 8/3*cos(1/3*arctan(3/37*sqrt(303))) - 8/3] > [True, True, None] > > For what it’s worth, Mathematica expresses its results without expressing > the roots of the quadrinomial, but denoting them as such roots : > > mathematica.Reduce([u==0 for u in Sys],[a1, a2, b1, b2]) > > (a1 == 0 && a2 == 0) || (a2 == 0 && a1 != 0 && b2 == 0) || > ((a1 == Root[2 - 4*#1^2 + #1^3 & , 1, 0] || > a1 == Root[2 - 4*#1^2 + #1^3 & , 2, 0] || > a1 == Root[2 - 4*#1^2 + #1^3 & , 3, 0]) && a2 == a1^2 && b2 == 2*b1) || > (a2 == a1^2 && a1*(2 - 4*a1^2 + a1^3) != 0 && b1 == 0 && b2 == 0) > > And, curiously, Sympy left to its own devices leads to an *erroneous* > solution, further explored on the |Sympyn list]( > https://groups.google.com/g/sympy/c/EB_Z6h3ZRDg). > > -- You received this message because you are subscribed to the Google Groups "sage-support" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/50fba344-ea3c-4ed0-adf2-a7b4403f4fcan%40googlegroups.com.
