#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
Satoshi Adachi
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)
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) batch(solve_system_of_algebraic_equations_OK.maxima)

batching #p/Volumes/HFS+2/home/adachi/work/285/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;

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

The result is
-------------------------------------------------------------------------------
Maxima 5.14.0cvs http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
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) batch(solve_system_of_algebraic_equations_NG.maxima)

batching #p/Volumes/HFS+2/home/adachi/work/285/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)
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) batch(solve_system_of_algebraic_equations_with_algebraic.maxima)

batching #p/Volumes/HFS+2/home/adachi/work/285/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,
Satoshi Adachi

Discussion

  • Barton Willis
    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:

    (%i92) load(topoly)$

    (%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
    true_sol : adjoin(si, true_sol);
    (%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

     
  • Satoshi Adachi
    Satoshi Adachi
    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.
    Your information is very helpful!

    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:
    *
    * load("topoly"); // for elim()
    * algebraic:true;
    *
    * S.Adachi 2008/04/10
    */

    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().
    *
    * S.Adachi 2008/04/10
    */

    load(topoly);
    load("algsys_willisbl.maxima");

    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:
    -------------------------------------------------------------------------------
    [Macintosh:~/work/288:1] adachi% maxima -b test_algsys_willisbl.maxima
    Maxima 5.14.0cvs http://maxima.sourceforge.net
    Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
    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) batch(test_algsys_willisbl.maxima)

    batching #p/Volumes/HFS+2/home/adachi/work/288/test_algsys_willisbl.maxima
    (%i2) load(topoly)
    (%o2) /usr/local/share/maxima/5.14.0cvs/share/contrib/topoly.lisp
    (%i3) load(algsys_willisbl.maxima)
    (%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"
    [Macintosh:~/work/288:2] adachi%
    -------------------------------------------------------------------------------

    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,
    Satoshi Adachi

     
  • Robert Dodier
    Robert Dodier
    2008-06-21

    • labels: --> Lisp Core - Solving equations