## #1370 (F=0 and G=0) and (F=0 and F-G=0) have different solutions.

open
nobody
5
2008-06-21
2008-03-17
No

Dear Developers of Maxima,

I met a set of two albebraic equations F = 0 and G = 0
of two variables x and y such that (i) solve([F=0,G=0],[x,y]) returns
no solution [] (this is wrong), whereas (ii) solve([F=0,F-G=0],[x,y])
returns true solutions. This is strange, since the first set of
equations F=0 and G=0 should have the same solutions of the second set of
equations F=0 and F-G=0.

(1) The Maxima program to solve F=0 and F-G=0 is
-------------------------------------------------------------------------------
/*
* solve_system_of_algebraic_equations_OK.maxima:
*
* F is a polynomial of x and y, which is degree 3 in x and 2 in y.
* G is a polynomial of x and y, which is degree 3 in x and 2 in y.
*
* solve() and algsys() can solve the system equations F=0 and F-G=0.
*/

display2d:false;

F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;

G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;

solutions:solve([F=0,F-G=0],[x,y]); /* OK: solutions are found. */
-------------------------------------------------------------------------------

The result is
-------------------------------------------------------------------------------
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) batch(solve_system_of_algebraic_equations_OK.maxima)

(%i2) display2d : false
(%o2) false
(%i3) F:1+x+2*a*x*y+3*x^2*y+2*a*x^2*y^2+2*x^3*y^2
(%o3) 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1
(%i4) G:1+a^2*y+2*a*x*y+x^2*y+a^2*x*y^2+2*a*x^2*y^2+x^3*y^2
(%o4) x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1
(%i5) solutions:solve([F = 0,F-G = 0],[x,y])
(%o5) [[x = (sqrt(a)*(a^2-4*a)+sqrt(a-4)*(2*a-a^2))/(2*sqrt(a-4)),
y = -1/(a*sqrt(a^2-4*a))],
[x = -(sqrt(a-4)*(a^2-2*a)+sqrt(a)*(a^2-4*a))/(2*sqrt(a-4)),
y = 1/(a*sqrt(a^2-4*a))]]
(%o6) "solve_system_of_algebraic_equations_OK.maxima"
------------------------------------------------------------------------------

===============================================================================

(2) The Maxima program to solve F=0 and G=0 is
-------------------------------------------------------------------------------
/*
* solve_system_of_algebraic_equations_NG.maxima:
*
* F is a polynomial of x and y, which is degree 3 in x and 2 in y.
* G is a polynomial of x and y, which is degree 3 in x and 2 in y.
*
* solve() and algsys() can NOT solve the system equations F=0 and G=0.
*/

display2d:false;

F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;

G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;

-------------------------------------------------------------------------------

The result is
-------------------------------------------------------------------------------
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) batch(solve_system_of_algebraic_equations_NG.maxima)

(%i2) display2d : false
(%o2) false
(%i3) F:1+x+2*a*x*y+3*x^2*y+2*a*x^2*y^2+2*x^3*y^2
(%o3) 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1
(%i4) G:1+a^2*y+2*a*x*y+x^2*y+a^2*x*y^2+2*a*x^2*y^2+x^3*y^2
(%o4) x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1
(%i5) solutions:solve([F = 0,G = 0],[x,y])
(%o5) []
(%o6) "solve_system_of_algebraic_equations_NG.maxima"
-------------------------------------------------------------------------------

===============================================================================

(3) I also tried algebraic:true to solve F=0 and G=0. The program is
-------------------------------------------------------------------------------
/*
* solve_system_of_algebraic_equations_with_algebraic_NG.maxima:
*
* F is a polynomial of x and y, which is degree 3 in x and 2 in y.
* G is a polynomial of x and y, which is degree 3 in x and 2 in y.
*
* solve() and algsys() can NOT solve the system equations F=0 and G=0.
*/

display2d:false;

algebraic:true;

F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;

G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;

solutions:solve([F=0,G=0],[x,y]); /* NG: internal error! */
-------------------------------------------------------------------------------

The result is further strange for me as
-------------------------------------------------------------------------------
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) batch(solve_system_of_algebraic_equations_with_algebraic.maxima)

(%i2) display2d : false
(%o2) false
(%i3) algebraic:true
(%o3) true
(%i4) F:1+x+2*a*x*y+3*x^2*y+2*a*x^2*y^2+2*x^3*y^2
(%o4) 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1
(%i5) G:1+a^2*y+2*a*x*y+x^2*y+a^2*x*y^2+2*a*x^2*y^2+x^3*y^2
(%o5) x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1
(%i6) solutions:solve([F = 0,G = 0],[x,y])
Polynomial quotient is not exact
-- an error. To debug this try debugmode(true);
-------------------------------------------------------------------------------

Sincerely yours,

## Discussion

• Barton Willis
2008-03-17

Logged In: YES
user_id=895922
Originator: NO

Thanks for reporting this bug. The function algsys isn't nearly as good as we would like.
The best I can offer is my favorite workaround (that might work); it's something like:

(%i93) algebraic : true\$

Use 'elim' to triangularize the equations---this can create spurious solutions:

(%i94) elim([F,G],[x,y]);
(%o94) [[],[x^4+a^2*x^3-a*x^3+a^3*x^2-a^2*x^2+a^3*x,x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1]]

Use algsys to solve these equations:

(%i95) sol : algsys(first(rest(%)),[x,y])\$

Try to reject the spurious solutions:

(%i96) true_sol : set()\$
(%i97) for si in sol do if every(lambda([a], is(equal(a,0))), ratsimp(subst(si,[F,G]))) then
(%o97) done
(%i98) true_sol;

(%o98) {[x=(a*sqrt(a^2-4*a)-a^2+2*a)/2,y=-sqrt(a^2-4*a)/(a^3-4*a^2)],
[x=-(a*sqrt(a^2-4*a)+a^2-2*a)/2,y=sqrt(a^2-4*a)/(a^3-4*a^2)]}

Two solutions *might* be spurious and two checked out OK. You'll
need to inspect the two rejected solutions to see if they are
really spurious. And a = 0 and a = 4 might be special cases that
need to be checked.

Of course, all this should be automatic; unfortunately, algsys isn't there yet.

Barton

2008-04-11

Logged In: YES
user_id=1953419
Originator: YES

Dear Mr.Barton Willis (willisbl),

Thank you very much for educating me to know how the decifiancy of algsys()
that was described in my previous bug report can be overcome.

At first, I am sorry to be late in writing this letter of thanks.
I had been occupied by duty works until yesterday.

Today, I wrote a wrapper routine of algsys(), which is based fully on your
suggestion and is named algsys_willisbl(), as follows:
-------------------------------------------------------------------------------
/*
* algsys_willisbl.maxima:
*
* better algsys() suggested by willisbl.
*
* [Remark] To use algsys_willisbl(), the following preparation is necessary:
*
* algebraic:true;
*
*/

algsys_willisbl(eqns_lst, vars_lst) :=
block([eqns1_lst,eqns2_lst,sols_cand_lst,sols_lst],
eqns1_lst:map(lambda([e],
if (not atom(e) and operatorp(e,"=")) then lhs(e)-rhs(e) else e),
eqns_lst),
eqns2_lst:first(rest(elim(eqns1_lst, vars_lst))),
sols_cand_lst:algsys(eqns2_lst, vars_lst),
sols_lst:listify(subset(setify(sols_cand_lst),
lambda([sc], every(lambda([a], is(equal(a,0))),
ratsimp(subst(sc, eqns1_lst)))))),
sols_lst);

/* END */
-------------------------------------------------------------------------------

I also wrote a test program for algsys_willisbl() by using the same system of
algebraic equations that was described in my previous bug report of algsys()
as
-------------------------------------------------------------------------------
/*
* test_algsys_willisbl.maxima:
*
* test algsys_willisbl() by applying it to a system of two algebraic
* equations F=0 and G=0, which cannot be solved by algsys().
*
*/

algebraic:true;

display2d:false;

F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;

G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;

sols1:algsys_willisbl([F,G], [x,y]);

sols2:algsys_willisbl([F=0,G=0], [x,y]);

sols3:algsys_willisbl([F+1=1,G+2=2], [x,y]);

/* END */
-------------------------------------------------------------------------------

The result of executation is as follows:
-------------------------------------------------------------------------------
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) batch(test_algsys_willisbl.maxima)

(%o2) /usr/local/share/maxima/5.14.0cvs/share/contrib/topoly.lisp
(%o3) algsys_willisbl.maxima
(%i4) algebraic : true
(%o4) true
(%i5) display2d : false
(%o5) false
(%i6) F:1+x+2*a*x*y+3*x^2*y+2*a*x^2*y^2+2*x^3*y^2
(%o6) 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1
(%i7) G:1+a^2*y+2*a*x*y+x^2*y+a^2*x*y^2+2*a*x^2*y^2+x^3*y^2
(%o7) x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1
(%i8) sols1:algsys_willisbl([F,G],[x,y])
(%o8) [[x = (a*sqrt(a^2-4*a)-a^2+2*a)/2,y = -sqrt(a^2-4*a)/(a^3-4*a^2)],
[x = -(a*sqrt(a^2-4*a)+a^2-2*a)/2,y = sqrt(a^2-4*a)/(a^3-4*a^2)]]
(%i9) sols2:algsys_willisbl([F = 0,G = 0],[x,y])
(%o9) [[x = (a*sqrt(a^2-4*a)-a^2+2*a)/2,y = -sqrt(a^2-4*a)/(a^3-4*a^2)],
[x = -(a*sqrt(a^2-4*a)+a^2-2*a)/2,y = sqrt(a^2-4*a)/(a^3-4*a^2)]]
(%i10) sols3:algsys_willisbl([1+F = 1,2+G = 2],[x,y])
(%o10) [[x = (a*sqrt(a^2-4*a)-a^2+2*a)/2,y = -sqrt(a^2-4*a)/(a^3-4*a^2)],
[x = -(a*sqrt(a^2-4*a)+a^2-2*a)/2,y = sqrt(a^2-4*a)/(a^3-4*a^2)]]
(%o11) "test_algsys_willisbl.maxima"
-------------------------------------------------------------------------------

I did not apply algsys_willisbl() to any other examples, yet.
So, it may have still bugs.

Anyway, I become now able to advance my study again
without the annoyance caused by algsys().
The knowledge that you gave me is very essential to use algsys()
in a better way.
Thank you very much.

Sincerely yours,