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/336c7823-9d04-4d62-b76f-8aea7779f2d2n%40googlegroups.com.
