search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

algsys/solve cannot find solutions #1244

Open rtoy opened 4 months ago

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:55:52 Created by *anonymous on 2007-02-04 21:07:23 Original: https://sourceforge.net/p/maxima/bugs/1106

(%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?

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:55:53 Created by willisbl on 2007-02-05 12:45:32 Original: https://sourceforge.net/p/maxima/bugs/1106/#462c

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]]

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:55:57 Created by robert_dodier on 2007-02-08 01:50:20 Original: https://sourceforge.net/p/maxima/bugs/1106/#6a2b

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:00 Created by billingd on 2016-10-02 06:55:11 Original: https://sourceforge.net/p/maxima/bugs/1106/#4222

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:04 Created by billingd on 2016-10-02 06:55:11 Original: https://sourceforge.net/p/maxima/bugs/1106/#b1c2

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.
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:08 Created by billingd on 2016-10-04 03:34:02 Original: https://sourceforge.net/p/maxima/bugs/1106/#4f4b

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]]
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:11 Created by billingd on 2016-10-07 06:53:13 Original: https://sourceforge.net/p/maxima/bugs/1106/#1a42

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]

Attachments:

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:15 Created by billingd on 2016-10-07 07:01:09 Original: https://sourceforge.net/p/maxima/bugs/1106/#09a7

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 */
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-04 06:56:18 Created by aleksasd on 2016-10-07 09:23:50 Original: https://sourceforge.net/p/maxima/bugs/1106/#9521

(%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+7073x^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