Menu ▾ ▴

#1106 algsys/solve cannot find solutions

None
open
David Billinghurst
Lisp Core - Solving equations (86)
5
2016-10-07
2007-02-04
Anonymous
No

(%i169) p1:-x*y^3+y^2+x^4-9*x/8$
(%i170) p2:y^4-x^3*y-9*y/8+x^2$

(%i180) algsys([p1,p2],[x,y]);
(%o180) [[x = 1/2,y = 1],[x = 9/8,y = 9/8],[x = 1,y = 1/2],[x = 0,y = 0]]

the system seems to have 4 (real) solutions.

`? algsys' says:
...
The method is as follows:

(1) First the equations are factored and split into subsystems.

(2) For each subsystem <S_i>, an equation <E> and a variable <x>
are selected. The variable is chosen to have lowest nonzero
degree. Then the resultant of <E> and <E_j> with respect to <x>
is computed for each of the remaining equations <E_j> in the
subsystem <S_i>. This yields a new subsystem <S_i'> in one fewer
variables, as <x> has been eliminated. The process now returns to
(1).
...

so `algsys' uses the resultant of the two polynomial with respect to x or y. i do this by hand:

(%i184) factor(resultant(p1,p2,x));
(%o184) 4096*(y-1)*y*(2*y-1)*(8*y-9)*(y^2+y+1)*(4*y^2+2*y+1)*(64*y^2+72*y+81)

(resultant with respect to x gives a similar result)
this gives 6 additional solution for y that can be found by sove. if i substitute such an y in the original polyomial the resulting polynomials in x have degree 4 and are also solvable. why don't `algsys' or `solve' don't find these solution?

Related

Bugs: #1266
Bugs: #3208

Discussion

  • Barton Willis

    Barton Willis - 2007-02-05

    Logged In: YES
    user_id=895922
    Originator: NO

    Here is a workaround for your equations; the workaround might help in general:

    (%i34) load(grobner)$
    Loading maxima-grobner $Revision: 1.2 $ $Date: 2006/11/08 03:40:02 $
    (%i35) p1:-x*y^3+y^2+x^4-9*x/8$
    (%i36) p2:y^4-x^3*y-9*y/8+x^2$
    (%i37) eqs : map('ratnumer, [p1,p2])$
    (%i38) eqs : poly_reduced_grobner(eqs,[x,y])$
    (%i39) algsys(eqs,[x,y]);
    (%o39) [[x=0,y=0],[x=1,y=1/2],[x=9/8,y=9/8],[x=1/2,y=1],[x=(sqrt(3)*%i-1)/4,y=-(sqrt(3)*%i+1)/2],[x=-
    (sqrt(3)*%i+1)/4,y=(sqrt(3)*%i-1)/2],[x=(9*sqrt(3)*%i-9)/16,y=-(9*sqrt(3)*%i+9)/16],[x=-(9*sqrt(3)*%i+9)/16,y=(9*sqrt(3)*%i-9)/16],[x=(sqrt(3)*%i-1)/2,y=-
    (sqrt(3)*%i+1)/4],[x=-(sqrt(3)*%i+1)/2,y=(sqrt(3)*%i-1)/4]]

     

  • Robert Dodier

    Robert Dodier - 2007-02-08
    • labels: --> Lisp Core - Solving equations
     

  • David Billinghurst

    David Billinghurst - 2016-10-02
    • assigned_to: David Billinghurst
    • Group: --> None
     

  • David Billinghurst

    David Billinghurst - 2016-10-02

    Commit [2986fc], which fixed a number of algsys bugs - see [bugs:#3208] - has made this worse. I think this is due to heap overflow in radcan called from function EBAKSUBST1. My changes didn't cause this, but they did allow the solution to procede to the point of failure.

    (%i1) p1:-x*y^3+y^2+x^4-9*x/8$

(%i2) p2:y^4-x^3*y-9*y/8+x^2$

(%i3) algsys([p1,p2],[x,y]);
Heap exhausted during garbage collection: 0 bytes available, 16 requested.
     

    Related

    Bugs: #3208
    Commit: [2986fc]

  • David Billinghurst

    David Billinghurst - 2016-10-04

    Commit [640ca7] has fixed the crash. Still only the 4 real solutions.

    (%i1) p1:-x*y^3+y^2+x^4-9*x/8$

(%i2) p2:y^4-x^3*y-9*y/8+x^2$

(%i3) algsys([p1,p2],[x,y]);

(%o3) [[x = 1/2,y = 1],[x = 9/8,y = 9/8],[x = 1,y = 1/2],[x = 0,y = 0]]
     

    Related

    Commit: [640ca7]

  • David Billinghurst

    David Billinghurst - 2016-10-07

    I have traced this and looked at one missed solution manually. The others are similar.

    algsys:
    1. decided to eliminate x and solve for y. Reasonable.
    2. calculates resultant(eq,eq2,y) and solves for y to find 4 real roots and 3 complex conjugate pairs.
    3. correctly finds solutions for the 4 real roots for y
    4. fails for the 6 complex roots roots for y

    For each complex root y, algsys
    1. substitutes y into eq1 and eq2
    2. calculates the resultant == 0 to confirm a solution exists
    3. solves a cubic polynomial for x, giving messy roots
    4. fails to confirm the the (x,y) pair is a solution as the symbolic calculation explodes.

    There are a number of opportunities for improvement. It may be desirable to
    1. issue a warning about missed roots when the resultant==0 but no corresponding solution is found.
    2. check constant expressions numerically to see if they are close to zero
    3. convert "messy" roots to floats if appropriate. algsys already returns numerical roots in certain circumstances.
    4. rectform can blow out the expression size a lot. A check on this would be prudent. (Done in commit [e49d61])

    (%i2) float2(e):=scanmap(lambda([u],rectform(float(u))),e,bottomup)

(%i3) eq1:(-x)*y^3+y^2+x^4+((-9)*x)/8
(%i4) eq2:y^4-x^3*y+((-9)*y)/8+x^2
(
%i5) eq1:expand(eq1*8)
                                3       2      4
(%o5)                   (- 8 x y ) + 8 y  + 8 x  - 9 x
(%i6) eq2:expand(eq2*(-8))
                              4       3              2
(%o6)                   (- 8 y ) + 8 x  y + 9 y - 8 x

/* eliminate x */

(%i7) r:resultant(eq1,eq2,x)
                                9          6         3
(%o7)             4096 y (4096 y  - 10440 y  + 7073 y  - 729)

(%i8) rf:factor(r)
                                           2              2
(%o8) 4096 (y - 1) y (2 y - 1) (8 y - 9) (y  + y + 1) (4 y  + 2 y + 1)
                                                                 2
                                                            (64 y  + 72 y + 81)

/* The resultant factors.
   Look at first quadratic factor */

(%i9) f:inpart(rf,6)
                                   2
(%o9)                             y  + y + 1

/* Solve for y and look at first root */

(%i10) Y:solve(f,y)
                         sqrt(3) %i + 1      sqrt(3) %i - 1
(%o10)            [y = - --------------, y = --------------]
                               2                   2
(%i11) Y1:Y[1]
                                   sqrt(3) %i + 1
(%o11)                       y = - --------------
                                         2

/* Substitute Y1 into eq1 and eq2.
   Check for factors.  None.  Tidy up */

(%i12) ev(eq11:eq1,Y1)
               4                   3                             2
(%o12)      8 x  + (sqrt(3) %i + 1)  x - 9 x + 2 (sqrt(3) %i + 1)
(%i13) factor(eq11)
                           4
(%o13)                  8 x  - 17 x + 4 sqrt(3) %i - 4
(%i14) ev(eq21:eq2,Y1)
                                                          4
                              3       2   (sqrt(3) %i + 1)
(%o14) (- 4 (sqrt(3) %i + 1) x ) - 8 x  - -----------------
                                                  2
                                                             9 (sqrt(3) %i + 1)
                                                           - ------------------
                                                                     2
(%i15) factor(eq21)
                               3      3       2
                 8 sqrt(3) %i x  + 8 x  + 16 x  + sqrt(3) %i + 1
(%o15)         - -----------------------------------------------
                                        2
(%i16) eq21:inpart(%,2)
                              3      3       2
(%o16)          8 sqrt(3) %i x  + 8 x  + 16 x  + sqrt(3) %i + 1

/* Solve eq11 and eq22 for x.  resultant=0, so solution exists */ 

(%i17) resultant(eq11,eq21,x)
(%o17)                                 0

(%i18) X:solve(eq21,x)
             - 1   sqrt(3) %i
(%o18) [x = (--- - ----------)
              2        2
           sqrt(10)                    8           1  1/3
 (-------------------------- - ----------------- - --)
     3/2                 3/2           1/2     3   16
  8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
                                sqrt(3) %i   - 1
                             4 (---------- + ---)
                                    2         2
 + -------------------------------------------------------------------------
          1/2     2           sqrt(10)                    8           1  1/3
   9 (%i 3    + 1)  (-------------------------- - ----------------- - --)
                        3/2                 3/2           1/2     3   16
                     8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
          2              sqrt(3) %i   - 1
 - ---------------, x = (---------- + ---)
          1/2                2         2
   3 (%i 3    + 1)
           sqrt(10)                    8           1  1/3
 (-------------------------- - ----------------- - --)
     3/2                 3/2           1/2     3   16
  8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
                                - 1   sqrt(3) %i
                             4 (--- - ----------)
                                 2        2
 + -------------------------------------------------------------------------
          1/2     2           sqrt(10)                    8           1  1/3
   9 (%i 3    + 1)  (-------------------------- - ----------------- - --)
                        3/2                 3/2           1/2     3   16
                     8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
          2
 - ---------------, x = 
          1/2
   3 (%i 3    + 1)
          sqrt(10)                    8           1  1/3
(-------------------------- - ----------------- - --)
    3/2                 3/2           1/2     3   16
 8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
                                       4
 + -------------------------------------------------------------------------
          1/2     2           sqrt(10)                    8           1  1/3
   9 (%i 3    + 1)  (-------------------------- - ----------------- - --)
                        3/2                 3/2           1/2     3   16
                     8 3    (sqrt(3) %i + 1)      27 (%i 3    + 1)
          2
 - ---------------]
          1/2
   3 (%i 3    + 1)

(%i19) Xf:float2(X)
(%o19) [x = 0.4330127018922192 %i - 0.2500000000000001, 
        x = 0.7006292692220368 %i - 0.4045084971874738, 
        x = 0.1545084971874738 - 0.2676165673298174 %i]

(%i20) ev(float2([eq1,eq2]),Y1,X[1])
(%o20) [6.661338147750939e-16 %i + 2.164934898019055e-15, 
                               1.498801083243961e-15 - 3.33066907387547e-15 %i]

(%i21) ev(float2([eq1,eq2]),Y1,X[2])
(%o21) [1.163118960624631 - 2.014581135048569 %i, 
                               1.77635683940025e-15 - 3.996802888650564e-15 %i]

(%i22) ev(float2([eq1,eq2]),Y1,X[3])
(%o22) [11.54086057667739 %i - 6.66311896062463, 
                              1.207367539279858e-15 - 2.692290834716005e-15 %i]
     

    Related

    Commit: [e49d61]


    Last edit: David Billinghurst 2016-10-08
    bug1106_1.mac

  • David Billinghurst

    David Billinghurst - 2016-10-07

    A reduced test case, isolated from the transcript above, is:

    eq11:8*x^4+(sqrt(3)*%i+1)^3*x-9*x+2*(sqrt(3)*%i+1)^2$
eq21:8*sqrt(3)*%i*x^3+8*x^3+16*x^2+sqrt(3)*%i+1$
algsys([eq11,eq21],[x]);

/* returns [], solution is
[x = 0.4330127018922192 %i - 0.2500000000000001,
 y = -(sqrt(3)*%i+1)/2 */
     

    Last edit: David Billinghurst 2016-10-07

  • Aleksas

    Aleksas - 2016-10-07

    (%i1) p1:-xy^3+y^2+x^4-9x/8$
    p2:y^4-x^3y-9y/8+x^2$
    (%i3) p3:eliminate([p1,p2],[y])[1];
    (%o3) 4096x(4096x^9-10440x^6+7073*x^3-729)
    (%i4) sols:solve([p1,p2,p3],[x,y]);
    (%o4) [[x=0,y=0],[x=1/2,y=1],[x=9/8,y=9/8],[x=1,y=1/2],[x=(sqrt(3)%i-1)/2,y=-(sqrt(3)%i+1)/4],[x=-(sqrt(3)%i+1)/2,y=(sqrt(3)%i-1)/4],[x=(3^(5/2)%i-9)/16,y=-(3^(5/2)%i+9)/16],[x=-(3^(5/2)%i+9)/16,y=(3^(5/2)%i-9)/16],[x=(sqrt(3)%i-1)/4,y=-(sqrt(3)%i+1)/2],[x=-(sqrt(3)%i+1)/4,y=(sqrt(3)%i-1)/2]]
    (%i5) length(%);
    (%o5) 10
    Test:
    (%i6) for k thru 10 do
    (subst(sols[k],[p1,p2]),expand(%%),print(%%));
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    [0,0]
    (%o6) done

     

Log in to post a comment.