A problem with the Sage solvers

[Emmanuel Charpentier]

ask.sagemat.org question 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).

​

[Emmanuel Charpentier]

On the Sympy list, Cris Smith points out 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…

​

[Dima Pasechnik]

sympy 1.9 is coming: https://trac.sagemath.org/ticket/32542

> --
> 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...@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.