On 4 January 2024 12:15:07 WET, Emmanuel Charpentier
<emanuel.charpent...@gmail.com> wrote:
>
>
>Indeed :
>(%i29) domain:complex; (%o29) complex (%i30) Sol1C:solve(Sys1, Unks);
>(%o30) [[x = 8,y = -40/3,l = (2*25^(1/3))/(3*9^(1/3))]] (%i31)
>Sol2C:solve(Sys2, Unks); (%o31) [] (%i32) map(lambda([w], map(lambda([v],
>subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys1))), Sol1C);
>(%o32) [[((8*(-1)^(2/3)*5^(5/3))/3^(2/3)-(40*25^(1/3))/9^(1/3))/12,
>((8*(-1)^(1/3)*3^(2/3)*5^(1/3)*25^(1/3))/9^(1/3)+40)
>/(2*(-1)^(1/3)*3^(2/3)*5^(1/3)),0]] (%i33) %, numer; (%o33)
>[[0.08333333333333333*(56.22884435344994
>*(0.8660254037844387*%i-0.4999999999999998) -56.22884435344994),
>0.1405721108836249*(39.99999999999999
>*(0.8660254037844386*%i+0.5000000000000001)
>+40)*(0.5000000000000001-0.8660254037844386*%i), 0]] (%i34) map(lambda([w],
>map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))),
>Sys2))), Sol2C); (%o34) [] (%i35) %, numer; (%o35) []
>
>But this does not explain whiy Sage is uneble o check (numericalmlmy or
>otherwise) the solutions given by Sympy or Mathematica, which check in
>Sympy (I didn’t yet try to check them in Mathematica, the limitations of
>the current Mathematica interface make this bothersome…).
>
>IMHO, this would deserve a critiocal ticket (possibly a blocker one), but I
>do not know how to report this efficiently. Suggestions welcome…
open an issue with this data there
>
>HTH,
>
>Le jeudi 4 janvier 2024 à 12:04:14 UTC+1, Dima Pasechnik a écrit :
>
>> You can get the same errors from pure Maxima if you set domain to
>> "complex",
>> no?
>>
>> On Thursday, January 4, 2024 at 10:29:56 AM UTC Emmanuel Charpentier wrote:
>>
>>> The problem seems Sage-specific : the same systems solve correctly (up to
>>> numerical noise) in “pure” Maxima :
>>> ;;; Loading #P"/usr/lib/x86_64-linux-gnu/ecl-21.2.1/sb-bsd-sockets.fas"
>>> ;;; Loading #P"/usr/lib/x86_64-linux-gnu/ecl-21.2.1/sockets.fas" Maxima
>>> 5.46.0 https://maxima.sourceforge.io using Lisp ECL 21.2.1 Distributed
>>> under the GNU Public License. See the file COPYING. Dedicated to the memory
>>> of William Schelter. The function bug_report() provides bug reporting
>>> information. (%i1) display2d:false; (%o1) false (%i2) Unks:[x,y,l]; (%o2)
>>> [x,y,l] (%i3) f:10*x^(1/3)*y^(2/3); (%o3) 10*x^(1/3)*y^(2/3) (%i4) g:5*x -
>>> 6*y; (%o4) 5*x-6*y (%i5) h:5*x^2 + 6*y; (%o5) 6*y+5*x^2 (%i6) fx:diff(f,x);
>>> (%o6) (10*y^(2/3))/(3*x^(2/3)) (%i7) fy:diff(f, y); (%o7)
>>> (20*x^(1/3))/(3*y^(1/3)) (%i8) gx:diff(g, x); (%o8) 5 (%i9) gy:diff(g, y);
>>> (%o9) -6 (%i10) hx:diff(h, x); (%o10) 10*x (%i11) hy:diff(h, y); (%o11) 6
>>> (%i12) Sys1:[fx=l*gx,fy=l*gy,g=120]; (%o12) [(10*y^(2/3))/(3*x^(2/3)) =
>>> 5*l,(20*x^(1/3))/(3*y^(1/3)) = -6*l, 5*x-6*y = 120] (%i13)
>>> Sys2:[fx=l*hx,fy=l*hy,h=120]; (%o13) [(10*y^(2/3))/(3*x^(2/3)) =
>>> 10*l*x,(20*x^(1/3))/(3*y^(1/3)) = 6*l, 6*y+5*x^2 = 120] (%i14)
>>> Sol1:solve(Sys1, Unks); (%o14) [[x = 8,y = -40/3,l = (2*5^(2/3))/3^(5/3)]]
>>> (%i15) Sol2:solve(Sys2, Unks); (%o15) [[x = (2*sqrt(6))/sqrt(5),y = 16,l =
>>> 18750^(1/6)/9], [x = -(2*sqrt(6))/sqrt(5),y = 16,l = -18750^(1/6)/9]]
>>> (%i16) map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u],
>>> ratsimp(lhs(u)-rhs(u))), Sys1))), Sol1); (%o16) [[0,0,0]] (%i17)
>>> map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u],
>>> ratsimp(lhs(u)-rhs(u))), Sys2))), Sol2); (%o17)
>>> [[(5^(1/3)*(5*2^(11/3)-(2^(8/3)*5^(1/6)*6^(5/6)*18750^(1/6))/3))
>>> /(3*2^(2/3)*6^(1/3)),
>>> -(2^(7/3)*18750^(1/6)-2^(7/3)*5^(5/6)*6^(1/6))/(3*2^(4/3)),0],
>>> [(5^(1/3)*(5*2^(11/3)-(2^(8/3)*5^(1/6)*6^(5/6)*18750^(1/6))/3))
>>> /(3*2^(2/3)*6^(1/3)),
>>> -(2^(7/3)*5^(5/6)*6^(1/6)-2^(7/3)*18750^(1/6))/(3*2^(4/3)),0]] (%i18) %,
>>> numer; (%o18) [[1.404069277432862e-15,0.0,0],[1.404069277432862e-15,0.0,0]]
>>>
>>> HTH,
>>>
>>> Le jeudi 4 janvier 2024 à 10:21:31 UTC+1, Emmanuel Charpentier a écrit :
>>>
>>>> Motivation : see [this post], which shows a case where Sage fails to
>>>> find the roots of a three equations system.
>>>>
>>>> I signalled in this thread that Sympy was able to find these roots. But
>>>> I stumbled on a difficulty checking these solutions.
>>>>
>>>> Set up the systems :
>>>> # Pretext : https://groups.google.com/g/sage-support/c/Nw12vYR0L0U
>>>> Unks=var('x,y,l
>>>> <https://groups.google.com/g/sage-support/c/Nw12vYR0L0UUnks=var('x,y,l>')
>>>> f(x, y) = 10*x^(1/3)*y^(2/3) g(x, y) = 5*x - 6*y h(x, y) = 5*x^2 + 6*y fx
>>>> =
>>>> diff(f,x) fy = diff(f, y) gx = diff(g, x) gy = diff(g, y) hx = diff(h, x)
>>>> hy = diff(h, y) Sys1 = [fx(x, y)==l*gx(x, y),fy(x, y)==l*gy(x, y),g(x,
>>>> y)==120] Sys2 = [fx(x, y)==l*hx(x, y),fy(x, y)==l*hy(x, y),h(x, y)==120]
>>>>
>>>> Sage’s (default) solution of the first system
>>>> DSol1 = solve(Sys1, Unks, solution_dict=True)
>>>>
>>>> As already shown in the original thread, Sage (slowly) fails to solve
>>>> the second system. Sympy can solve it (anfd the first one too…) :
>>>> SSol1 = solve(Sys1, Unks, algorithm="sympy") SSol2 = solve(Sys2, Unks,
>>>> algorithm="sympy")
>>>>
>>>> For what it’s worth, chech the “competition” :
>>>> MSol1 = [{u[1].sage():u[2].sage() for u in s} for s in
>>>> mathematica(Sys1).Solve(Unks)] MSol2 = [{u[1].sage():u[2].sage() for u in
>>>> s} for s in mathematica(Sys2).Solve(Unks)]
>>>>
>>>> Both Sympy and Mathematica find one solution to the first system and two
>>>> for the second. These solutions have different ex^pressins, but that is
>>>> not
>>>> problematic. The fly in the ointment is that Sage has trouble checking
>>>> these solutions. A simple numerical check of these solutions is :
>>>> sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s in
>>>> DSol1] [[-0.000977, 8.44 - 4.87*I, 0.000]] sage:
>>>> [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s in SSol1]
>>>> [[7.03 - 4.06*I, 0.000488*I, 0.000]] sage:
>>>> [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys2] for s in SSol2]
>>>> [[0.000977, 0.000, 0.000], [-0.000977 + 0.00195*I, 0.000122 - 0.000244*I,
>>>> 0.000]] sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s
>>>> in MSol1] [[7.03 - 4.06*I, 0.000, 0.000]] sage:
>>>> [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys2] for s in MSol2]
>>>> [[-0.00146 - 0.000977*I, 0.000122 + 0.000244*I, 0.000], [0.00195,
>>>> 0.000244,
>>>> 0.000]]
>>>>
>>>> None of these solutions seems to check. Uh oh…
>>>>
>>>> However, using sympy to compute the same numerical checks gives
>>>> different results :
>>>> sage:
>>>> [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e
>>>> in Sys1] for s in DSol1] [[-7.03 + 4.06*I, 8.43 - 4.87*I, 0]] sage:
>>>> [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e
>>>> in Sys1] for s in SSol1] [[0, 0, 0]] sage:
>>>> [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e
>>>> in Sys2] for s in SSol2] [[0, 0, 0], [0, 0, 0]] sage:
>>>> [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e
>>>> in Sys1] for s in MSol1] [[0, 0, 0]] sage:
>>>> [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e
>>>> in Sys2] for s in MSol2] [[0, 0, 0], [0, 0, 0]]
>>>>
>>>> Sage’s solution to the first system does not check ; Both Sympy’s and
>>>> Mathematica’s solutions of both systems check.
>>>>
>>>> This tells us that
>>>>
>>>> -
>>>>
>>>> Sage’s solution of the first system is wrong, and that
>>>> -
>>>>
>>>> Sage’s computaion of the numerical checks of knows solutins is *also*
>>>> wrong.
>>>>
>>>> This suggests a bug in Sage’s computations ; this problem *may* be
>>>> bound to the use of fractional powers.
>>>>
>>>> Any hint on the “right” way to report this issue eficiently welcome.
>>>>
>>>>
>>>
>
--
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/5407328D-6FC6-44E8-9AC7-7D31A29F3E03%40gmail.com.
