85 views

Skip to first unread message

Jan 4, 2024, 4:21:31 AMJan 4

to sage-support

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.

Jan 4, 2024, 5:29:56 AMJan 4

to sage-support

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,

Jan 4, 2024, 6:04:14 AMJan 4

to sage-support

You can get the same errors from pure Maxima if you set domain to "complex",

no?

Jan 4, 2024, 7:15:07 AMJan 4

to sage-support

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,

Jan 4, 2024, 7:57:49 AMJan 4

to sage-s...@googlegroups.com

On Thu, 4 Jan 2024 at 12:15, Emmanuel Charpentier

<emanuel.c...@gmail.com> wrote:

>

> 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…).

The problem is with Sage's subs method:
<emanuel.c...@gmail.com> wrote:

>

> 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…).

sage: e1 = (-40/3)**(2/3)

sage: e2 = (y**(2/3)).subs({y:-40/3})

sage: e1

4*(-5/3)^(2/3)

sage: e2

4*(5/3)^(2/3)

sage: e1.n()

-2.81144221767250 + 4.86956076355288*I

sage: e2.n()

5.62288443534499

When the number to be substituted has an exact cube factor subs loses

the sign under the radical:

sage: (y^(2/3)).subs({y:-12})

(-12)^(2/3)

sage: (y^(2/3)).subs({y:-24})

4*3^(2/3)

This is most likely due to confusing y^(n/m) as (y^n)^(1/m) when it

should be (y^(1/m))^n. These are not equivalent if y^(1/m) means the

principal mth root.

--

Oscar

Jan 4, 2024, 9:37:53 AMJan 4

to sage-s...@googlegroups.com

On 4 January 2024 12:15:07 WET, Emmanuel Charpentier <emanuel.c...@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…

>>>>

>>>> Sage’s solution of the first system is wrong, and 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

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu