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') 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/c6a7f267-dd7b-4f7b-807b-42f59374c6e0n%40googlegroups.com.
