#3461 algsys failure

[David Billinghurst]

Reported to maxima list by Jan-Magnus Økland on 2018-08-25.

I have a nonlinear system of two polynomial equations in two variables
x,y with a parameter t that has a double solution for all t (for t=0 one
circle reduces to a point that lies on the other circle).

eq1:expand((x-(2+3*(40/t)/sqrt(10)))^2+(y-(-5+(40/t)/sqrt(10)))^2-(40/t)^2);
eq2:expand((x-(-2+3*t/sqrt(10)))^2+(y-(7+t/sqrt(10)))^2-t^2);
solve([eq1,eq2],[x,y]); /* yields [] */

but the following gives a solution sub-optimally:

solve([eq1-eq2],[x]);
x0:rhs(%[1]);
solve([subst(x0,x,eq2)],[y]);
y1:rhs(%[1]);
x1:subst(y1,y,x0); /* ‘ratsimp's to quotient of degree 4 polys */
/* the resulting parametric circle plotted */
plot2d([parametric, x1, y1, [t, -100,100]]);

mathematica can do something like the above from Reduce (I actually used WA).

maple finds quotients of degree two polynomials for both x and y from solve:

[x = (2*(12*t*sqrt(10)+t^2-40))/(t^2+40),
 y = (280-5*t^2+8*t*sqrt(10))/(t^2+40)]

Does anyone know what stops solve from finding this solution?

[David Billinghurst]

solve() delegates the solution of nonlinear systems of polynomial equations to algsys().

algsys() uses heuristics, and is sensitive to the ordering of equations and variables. For some problems, preprocessing the equations with poly_reduced_grobner() is helpful. It works here.

(%i1) display2d:false$
(%i2) load(grobner)$
(%i3) eq1:expand((x-(2+3*(40/t)/sqrt(10)))^2+(y-(-5+(40/t)/sqrt(10)))^2-(40/t)^2)$
(%i4) eq2:expand((x-(-2+3*t/sqrt(10)))^2+(y-(7+t/sqrt(10)))^2-t^2)$
(%i5) poly_reduced_grobner([eq1,eq2],[x,y])$
(%i6) s:algsys(%,[x,y]);
(%o6) [[x = (18*t^6+168*sqrt(10)*t^5-7600*t^4-96*10^(5/2)*t^3+304000*t^2
                   +2688*10^(5/2)*t-1152000)
/(9*t^6-24*sqrt(10)*t^5-200*t^4-8000*t^2+384*10^(5/2)*t+576000),
        y = -(45*t^8+348*sqrt(10)*t^7-28240*t^6-1712*10^(3/2)*t^5+2208000*t^4
+8*10^(11/2)*t^3-56192000*t^2-33024*10^(7/2)*t+161280000)
/(9*t^8+84*sqrt(10)*t^7-3440*t^6-144*10^(3/2)*t^5-576*10^(5/2)*t^3
                 +5504000*t^2+5376*10^(7/2)*t-23040000)]]

Some simple numerical checks indicate the solution is reasonable. Perhaps you can simplify the result.

(%i16) float(subst(s[1],[eq1,eq2])),t=3;
(%o16) [-1.4210854715202E-14,-3.552713678800501E-15]
(%i17) float(subst(s[1],[eq1,eq2])),t=2.79;
(%o17) [-2.842170943040401E-14,7.105427357601002E-15]
(%i18) float(subst(s[1],[eq1,eq2])),t=sqrt(39);
(%o18) [-1.77635683940025E-14,0.0]

Thanks for the simple test case.

[David Billinghurst]

The factor sqrt(10) is part of the problem. With the change of variables t = u*sqrt(10)

(%i33) eq3:eq1,t=u*sqrt(10);
(%o33) y^2-(8*y)/u+10*y+x^2-(24*x)/u-4*x+8/u+29
(%i34) eq4:eq2,t=u*sqrt(10);
(%o34) y^2-2*u*y-14*y+x^2-6*u*x+4*x+2*u+53
(%i35) algsys([eq3,eq4],[x,y]);
(%o35) [[x = (2*u^2+24*u-8)/(u^2+4),y = -(5*u^2-8*u-28)/(u^2+4)]]