#3321 algsys regression

[David Billinghurst]

The following system of equations can be solved with version 5.38.1 but later versions return []. Reported by Stefano Milani to mailing list on 2017-07-25. Original report against solve, but the failure is within algsys.

eqs:[k1*x1*x3-j1*x2=0,k2*x2*x3-j2*x1=0,(mu*(x2+x1))/x3-c*x2-b*x1+a*u1=0];
algsys(eqs,[x2,x3,u1]);

(%o6) [[x2 = -sqrt((j2*k1)/(j1*k2))*x1,
        x3 = -(j1*sqrt(j2)*sqrt(k2)-sqrt(j1)*j2*sqrt(k1))
           /(sqrt(j1)*sqrt(k1)*k2-sqrt(j2)*k1*sqrt(k2)),
        u1 = ((b*j1*j2*k2-j2*k1*k2*mu)*x1+sqrt(j1*j2*k1*k2)*(k2*mu-c*j2)*x1)
           /(a*j1*j2*k2)],
       [x2 = sqrt((j2*k1)/(j1*k2))*x1,
        x3 = (j1*sqrt(j2)*sqrt(k2)+sqrt(j1)*j2*sqrt(k1))
           /(sqrt(j1)*sqrt(k1)*k2+sqrt(j2)*k1*sqrt(k2)),
        u1 = -((j2*k1*k2*mu-b*j1*j2*k2)*x1+sqrt(j1*j2*k1*k2)*(k2*mu-c*j2)*x1)
           /(a*j1*j2*k2)]]

git bisect points to

commit 640ca759d4b7dc54de9fd0d15149bd4471867546
Author: David Billinghurst <dbmaxima@gmail.com>
Date:   Tue Oct 4 14:14:59 2016 +1100

    SF [bug:#1266][bug:#3208] Don't use $radcan when simplifying in algsys.

    Inappropriate use of $radcan when simplifying expressions caused a
    number of failures in algsys.  This commit follows on from recent
    work that removed one call to radcan and introduced use of sqrtdenest.
    This addresses a number of long-standing solve/algsys failures.

    The recent changes to ebaksubst simpified some expressions to zero
    and allowed the solution to continue until it crashed in a call to
    radcan from pfreeofmainvarsp.

    * src/algsys.lisp:
      New function simplify-after-subst that uses sqrtdenest, $rectform
      and $ratsimp.  This replaces uses of radcan to simplify expressions.
      Call simplify-after-subst from ebaksubst1, pfreeofmainvarsp and
      pfreeofmainvarsp.

[David Billinghurst]

Reduce the problem. Rename constants and variables, then remove constants while problem remains unsolvable with maxima 5.39.0

(%i31) neqs:eqs,[k1=d,x1=e,j1=f,k2=g,j2=h,mu=m,x2=x,x3=y,u1=z];
(%o31) [d*e*y-f*x = 0,g*x*y-e*h = 0,a*z+(m*(x+e))/y-c*x-b*e = 0]
(%i32) subst([a=1,b=1,c=1,d=1,e=1,f=1,m=1],neqs);
(%o32) [y-x = 0,g*x*y-h = 0,z+(x+1)/y-x-1 = 0];
(%i33) algsys(%,[x,y,z]);
(%o33) []

This is solvable in 5.38.1, with solution

(%o18) [[x = sqrt(h/g),y = (h+sqrt(g)*sqrt(h))/(sqrt(g)*sqrt(h)+g),
         z = sqrt(1/(g*h))*h-g*sqrt(1/(g*h))],
        [x = -sqrt(h/g),y = -(h-sqrt(g)*sqrt(h))/(sqrt(g)*sqrt(h)-g),
         z = g*sqrt(1/(g*h))-sqrt(1/(g*h))*h]]

[David Billinghurst]

Tracing this with

 display2d:false$
(?trace(?algsys),
?trace(?algsys0),
?trace(?algsys1),
?trace(?condensesolnl),
?trace(?distrep),
?trace(?lofactors),
?trace(?findleastvar),
?trace(?presultant),
?trace(?bakalevel),
?trace(?bakalevel1),
?trace(?mergesoln),
?trace(?callsolve),
?trace(?callsolve1),
?trace(?callsolve2),
?trace(?baksubst),
?trace(?ebaksubst),
?trace(?ebaksubst1),
?trace(?$resultant),
?trace(?solve),
?trace(?simplify\-after\-subst)
)$

eq1:y-x = 0;
eq2:g*x*y-h = 0;
eq3:z+(x+1)/y-x-1 = 0;

algsys([eq1,eq2,eq3],[x,y,z]);

Looking quickly through the output I found an expression that is simplified to zero by $radcan but not by lisp function simplify-after-subst. If we restrict g and h to positive values then it can be simplified

ex:h*(1-g*sqrt(1/(g*h))*sqrt(h/g))+g^2*sqrt(1/(g*h))*sqrt(h/g)-g$

radcan(ex);
0

?simplify\-after\-subst(ex);
sqrt(h/g)*(g^2*sqrt(1/(g*h))-g*sqrt(1/(g*h))*h)+h-g

assume(g>0,h>0)$

?simplify\-after\-subst(ex);
0

If we assume all the constants in the original problem are positive then current maxima gives solutions that seem equivalent to 5.38.1.

I haven't tried to find the minimum assumptions necessary.

(%i1) eqs:[k1*x1*x3-j1*x2=0,k2*x2*x3-j2*x1=0,(mu*(x2+x1))/x3-c*x2-b*x1+a*u1=0];

(%o1) [k1*x1*x3-j1*x2 = 0,k2*x2*x3-j2*x1 = 0,(mu*(x2+x1))/x3-c*x2-b*x1+a*u1
                                               = 0]
(%i2) assume(k1>0,x1>0,j1>0,k2>0,j2>0,mu>0,c>0,b>0,a>0);

(%o2) [k1 > 0,x1 > 0,j1 > 0,k2 > 0,j2 > 0,mu > 0,c > 0,b > 0,a > 0]
(%i3) algsys(eqs,[x2,x3,u1]);

(%o3) [[x2 = -(sqrt(j2)*sqrt(k1)*x1)/(sqrt(j1)*sqrt(k2)),
        x3 = (j1*j2*k2-sqrt(j1)*j2^(3/2)*sqrt(k1)*sqrt(k2))
           /(j2*k1*k2-sqrt(j1)*sqrt(j2)*sqrt(k1)*k2^(3/2)),
        u1 = -(((-sqrt(j1)*sqrt(j2)*sqrt(k1)*k2^(3/2)*mu)
           +j2*k1*k2*mu-b*j1*j2*k2+c*sqrt(j1)*j2^(3/2)*sqrt(k1)*sqrt(k2))
           *x1)
           /(a*j1*j2*k2)],
       [x2 = (sqrt(j2)*sqrt(k1)*x1)/(sqrt(j1)*sqrt(k2)),
        x3 = (j1*j2*k2+sqrt(j1)*j2^(3/2)*sqrt(k1)*sqrt(k2))
           /(sqrt(j1)*sqrt(j2)*sqrt(k1)*k2^(3/2)+j2*k1*k2),
        u1 = -((sqrt(j1)*sqrt(j2)*sqrt(k1)*k2^(3/2)*mu
           +j2*k1*k2*mu-b*j1*j2*k2-c*sqrt(j1)*j2^(3/2)*sqrt(k1)*sqrt(k2))
           *x1)
           /(a*j1*j2*k2)]]

[David Billinghurst]

This isn't a bug. The solution from 5.38.1 is wrong for certain choices of constants.

Using the reduced test case above, then taking g=1 and h=-1, we see that the "solution" doesn't satisfy the original equations.

(%i1) display2d:false$

(%i2) eqs:[y-x = 0,g*x*y-h = 0,z+(x+1)/y-x-1 = 0]$

(%i3) s:algsys(eqs,[x,y,z]);

(%o3) [[x = sqrt(h/g),y = (h+sqrt(g)*sqrt(h))/(sqrt(g)*sqrt(h)+g),
        z = sqrt(1/(g*h))*h-g*sqrt(1/(g*h))],
       [x = -sqrt(h/g),y = -(h-sqrt(g)*sqrt(h))/(sqrt(g)*sqrt(h)-g),
        z = g*sqrt(1/(g*h))-sqrt(1/(g*h))*h]]

(%i5) eqs,s;

(%o5) [(h+sqrt(g)*sqrt(h))/(sqrt(g)*sqrt(h)+g)-sqrt(h/g) = 0,
       (g*sqrt(h/g)*(h+sqrt(g)*sqrt(h)))/(sqrt(g)*sqrt(h)+g)-h = 0,
       ((sqrt(g)*sqrt(h)+g)*(sqrt(h/g)+1))/(h+sqrt(g)*sqrt(h))
        -sqrt(h/g)+sqrt(1/(g*h))*h-g*sqrt(1/(g*h))-1
         = 0]

(%i6) %,g=1,h=-1;

(%o6) [(%i-1)/(%i+1)-%i = 0,((%i-1)*%i)/(%i+1)+1 = 0,(%i+1)^2/(%i-1)-3*%i-1
                                                       = 0]
(%i7) ratsimp(%);

(%o7) [0 = 0,0 = 0,(4*%i+4)/(%i-1) = 0]

[David Billinghurst]

Not a bug. Test cases added to rtest_algsys.mac as commit [cd0fd4].

[David Billinghurst]

Not a bug. It is unfortunate that algsys returns [] with no message when a solution requires constraints on a constant.