sympy 1.9 is coming: https://trac.sagemath.org/ticket/32542
On Mon, Dec 6, 2021 at 11:19 AM Emmanuel Charpentier
<emanuel.charpent...@gmail.com> wrote:
>
> 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…
>
> Le dimanche 5 décembre 2021 à 20:55:32 UTC+1, Emmanuel Charpentier a écrit :
>>
>> 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).
>
> --
> 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.
--
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/CAAWYfq35Rut28%2B_8eg%3DaNxEz2QimvZGLVVNi_YTUUexVmF3ZCA%40mail.gmail.com.
