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/9bb67b28-ec60-484f-98f0-95181cedb5b0n%40googlegroups.com.
