#4032 A linear system of equations is correctly solved, but its specialization can't be correctly solved.

[Emmanuel Charpentier]

Found in the Maxima version embedded in Sage, currently 5.45.0.

A system of four linear equations in C can be solved (easily) with Maxima, but its specialization by replacing one variable by a constant gives a wrong answer. Original report : an ask.sagemath.org question, where the issue can be traced to Maxima's solver ; discussed in the sage-devel mailing list.

Consider :

/* Original system, where I as a *variable* */
(%i31) Sys:[I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) = U,
I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) = 0,
I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) = 0,
I4/(w*I*C3) + I4*RL - I3/(w*I*C3) = 0]$

(%i32) IVars:[I1, I2, I3, I4]$
/* Solution of the original system */
(%i33) Sol: solve(Sys, IVars)$
/* Solution check : */
(%i34) map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);
(%o34)                     [true, true, true, true]

However :

/* Specializing Sys by setting I=%i
(%i35) Sys0:map(lambda([x], subst(I=%i, x)), Sys)$
/* Solution */
(%i36) Sol0: solve(Sys0, IVars)$
/* Checking */
(%i37) map(lambda([x], is(ratsimp(subst(Sol0, x)))), Sys0);
(%o37)                    [true, false, false, true]
(%i38)

This suggests that something is wrong in the handling of complex constants.

Since this problem can silently lead to erroneous solutions, I tend to consider that it should be given the highest priority.

HTH,

[Robert Dodier]

Thanks for the report. I see it in Maxima 5.46.0. I experimented a little to try to find a simpler example, but didn't find one. If anyone can find a simpler example, that would be really helpful.

Minor point about subst -- it will work on aggregate expressions such as lists all at once, so you can say something like subst(Sol0, Sys0) and subst(I = %i, Sys) instead of mapping over a list.

[Robert Dodier]

Some investigation shows that solve calls the Lisp function TFGELI in src/mat.lisp.

For the record, linsolve also calls TFGELI, and gets the same result as solve.

[Robert Dodier]

Looks like the (undocumented) linsolve_by_lu, which computes the solution by another path (via LU decomposition) finds a correct solution for the equations substituted with I = %i. Just another bit of evidence pointing towards the problem being in TFGELI and not somewhere else in Maxima.

[Robert Dodier]

I think this is a slightly simpler example:
[%i*I1*w+(%i*I2)/(C1*w)-(%i*I1)/(C1*w) = U, %i*I3*w+%i*I2*w+(%i*I3)/(C2*w)-(%i*I2)/(C2*w)-(%i*I2)/(C1*w)+(%i*I1)/(C1*w) = 0, %i*I3*w+%i*I2*w+(%i*I4)/w-(%i*I3)/(C2*w)-(%i*I3)/w+(%i*I2)/(C2*w) = 0, (-(%i*I4)/w)+(%i*I3)/w+I4 = 0];
which I found by substituting C3 = 1, L1 = 1, L2 = 1, L3 = 1, M23 = 1, RL = 1 and I = %i into Sys.

At this point, w, C1, and C2 are required -- substituting 1 for any of those enables solve to get a correct solution.

[Florian Königstein]

Maybe this is a bit more simple (some divisions are replaced by multiplications and replace U by 1):
[%i*I1*w+(%i*I2)/(C1*w)-(%i*I1)/(C1*w) = 1, %i*I3*w+%i*I2*w+(%i*I3)/(C2*w)-(%i*I2)/(C2*w)-(%i*I2)/(C1*w)+(%i*I1)/(C1*w) = 0, %i*I2*w-(%i*I4)*w-(%i*I3)/(C2*w)+2*(%i*I3)*w+(%i*I2)/(C2*w) = 0, %i*I4*w-(%i*I3)*w+I4 = 0]$

What is remarkable: If you replace the summand 2*(%i*I3)*w by (%i*I3)*w (leave away the constant 2), then Maxima finds the correct solution. If you try to replace other divisions by multiplications, e.g. with C1*w or C2*w, it also finds the correct solution.

Sys:[%i*I1*w+(%i*I2)/(C1*w)-(%i*I1)/(C1*w) = U, %i*I3*w+%i*I2*w+(%i*I3)/(C2*w)-(%i*I2)/(C2*w)-(%i*I2)/(C1*w)+(%i*I1)/(C1*w) = 0, %i*I2*w-(%i*I4)*w-(%i*I3)/(C2*w)+2*(%i*I3)*w+(%i*I2)/(C2*w) = 0, ((%i*I4)*w)-(%i*I3)*w+I4 = 0]$

Another more simple version is the following:
[%i*I1*w+(%i*I2)*w/(C1)-(%i*I1)*w/(C1) = 1, %i*I3*w+%i*I2*w+(%i*I3)/(C2*w)-(%i*I2)/(C2*w)-(%i*I2)*w/(C1)+(%i*I1)*w/(C1) = 0, %i*I2*w-(%i*I4)*w-(%i*I3)/(C2*w)+(%i*I2)/(C2*w) = 0, %i*I4*w-(%i*I3)*w+I4 = 0]$

and more simple:
[%i*I1*w+(%i*I2)*w/(C1)-(%i*I1)*w/(C1) = 1, %i*I3*w+%i*I2*w+(%i*I3)/(C2*w)-(%i*I2)/(C2*w)-%i*(I1)*w/(C1)= 0, %i*I2*w-(%i*I4)*w-(%i*I3)/(C2*w)+(%i*I2)/(C2*w) = 0, %i*I4*w-(%i*I3)*w+I4 = 0]$

and more simple:
[-(%i*I1)*w/(C1) = 1, %i*I3*w+%i*I2*w-%i*(I1)*w/(C1)= 0, %i*I2*w-(%i*I3)/(C2*w) = 0, %i*I4*w+I4 = 0]
In this version the solution for I4 is 0 (assuming w is real).

This is more simple:
[-(%i*I1)*w/(C1) = 1, %i*I3*w+%i*I2*w-%i*(I1)*w/(C1)= 0, %i*I2*w-(%i*I3)/(C2*w) = 0, %i*I4*w+I4 = 0]$

From the first equation you can see that I1 must be C1 / (-%i*w), but Maxima finds the solution I1 = -(C1*w^2-2*%i*C1*w-C1)/(w^2-%i*w) .

This is more simple:
[-(%i*I1)*w/(C1) = 1, %i*I2*w-%i*(I1)*w/(C1)= 0, (%i*I3)/(C2*w) = 0, %i*I4*w+I4 = 0]$

I3 and I4 are correctly found to be zero, but I1 and I2 are incorrect.

And more simple:
-(%i*I1)*w/(C1) = 1, %i*I2*w-%i*(I1)*w/(C1)= 0, I3 = 0, %i*I4*w+I4 = 0]$
Again, I1 and I2 are incorrect.

And more simple by replacing mult with div: I list the whole cell again:

Sys:[-(%i*I1)*w*C1 = 1, %i*I2*w-%i*(I1)*w*C1= 0, I3 = 0, %i*I4*w+I4 = 0]$
IVars:[I1, I2, I3, I4]$
Sol: solve(Sys, IVars);
/* Solution check : */
map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);

The output is:
[[I1=-(w^2-2*%i*w-1)/(%i*C1*w^2+C1*w),I2=-(w-%i)/(%i*w^2+w),I4=0]] [false,false,true,true]

And you can drop the equation for I3=0 and do not solve for I3:

Sys:[-(%i*I1)*w*C1 = 1, %i*I2*w-%i*(I1)*w*C1= 0, %i*I4*w+I4 = 0]$
IVars:[I1, I2, I4]$
Sol: solve(Sys, IVars);
/* Solution check : */
map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);

Output:

[[I1=-(w^2-2*%i*w-1)/(%i*C1*w^2+C1*w),I2=-(w-%i)/(%i*w^2+w),I4=0]]
[false,false,true]

Changing some minus to plus and setting C2=1 makes the solution for I1 correct, but I2 is wrong:

Sys:[%i*I1*w*C1 = 1, %i*I1*w+%i*w*I2= 0, %i*I4*w+I4 = 0]$
IVars:[I1, I2, I4]$
Sol: solve(Sys, IVars);
/* Solution check : */
map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);

Output:
[[I1=(w-%i)/(%i*C1*w^2+C1*w),I2=-(w^2-2*%i*w-1)/(%i*C1*w^2+C1*w),I4=0]] [true,false,true]

After replacing one %i by 1 and I1 by I5 that is not solved for (only in the second equation):

Sys:[%i*I1*w*C1 = 1, I2*w = I5*w, %i*I4*w+I4 = 0]$
IVars:[I1, I2, I4]$
Sol: solve(Sys, IVars);
/* Solution check : */
map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);

Output:

[[I1=-(%i*w^2+2*w-%i)/(C1*w^2-%i*C1*w),I2=I5,I4=0]]
[false,true,true]

Remarkable in this system is that initially all 3 equations are independent, so that their solution could easily be found without any Gauss-elimination steps.

I have also tried to make some assumptions about the parameters (not being solved for), but the solution stays wrong:
assume(w>0, C1>0, I5>0);

Also remarkable is: Replacing w by w2 in the second equation keeps the solution wrong. But replacing w by w3 in the third equation makes it correct.

And slightly simpler:

assume(a>0, b>0, c>0, d>0, f>0)$
Sys:[I1*a*b = %i, I2*c = d*f, I3*(%i*a + 1) = 0]$
IVars:[I1, I2, I3]$
Sol: solve(Sys, IVars);
/* Solution check : */
map(lambda([x], is(ratsimp(subst(Sol, x)))), Sys);

Output:

[[I1=-(%i*a^2+2*a-%i)/((%i*a^2+a)*b),I2=(d*f)/c,I3=0]]  
[false,true,true]

[Robert Dodier]

Great work, Florian! Your final example can be simplified just a little bit more. It appears that the assume isn't needed -- the result from solve is incorrect even when assume is not present. (That's what I would expect since solve ignores assume.)