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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-support/9ca217da-7198-4e6f-9c0d-d4783baaee0bn%40googlegroups.com.
