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… 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/07089e45-9430-447e-990c-003a5ed08311n%40googlegroups.com.
