linsolve strange behavior



Hello, 

 I have a strange behavior of Maxima (app-sci/maxima-5.9.0-r2 
in Gentoo GNU/Linux, should be 5.9.0 plain maxima for practical 
purposes) compiled against cmucl (dev-lisp/cmucl-18e)
and gcl (dev-lisp/gcl-2.5.2).

 For a particular input, I end up needing to solve 
a linear system of 3 eqn.s in 4 unknowns with 4
parameters (DX,coeff_x,DY,coeff_y):

[coe[4]+coe[3]+coe[2]+coe[1] = 1,
     
-coe[2]*DY/(2*coeff_y)-coe[1]*DY/(2*coeff_y)+coe[4]*DY/2+coe[3]*DY/2= 0,

-coe[4]*DX/(2*coeff_x)+coe[3]*DX/(2*coeff_x)-coe[2]*DX/(2*coeff_x)         +coe[1]*DX/(2*coeff_x) = 0]

with unknowns coe[i].

The system is compatible, you can resolve manually, DX, DY 
and coeff_x should be immediately removed and not 
contribute to the solutions. However, Maxima cannot find the
right solution: tested both cmucl and gcl with the 
.max file I include below and on another Maxima
installation on another computer. 

If I simplify manually, as seen below in
list_eq_2, list_eq_3, list_eq_4,list_eq_5,
finally it finds the correct solution.

This is the only special case where my Maxima code
is not working, so I believe there should not
be a visible mistake. It is true I read the
manuals sometime ago, but really I cannot 
understand if there is an error, and which it is.

I apologize in advance if it is an idiotic 
problem, but I am in the middle of it and
possibly I have tunnel vision and cannot
spot the reason of this.

Funny thing is, I asked a friend to try to solve
the same system with Mathematica 4.1 on Windows, and she
cannot solve it, it returns an empty 
solution vector. 

Is there anybody willing to the test the code 
below? Maybe I have something seriously 
screwed in my system and this is affecting
the outcome of the solution?

In any case, thanks in advance.

Best Regards,

 Raimondo Giammanco

########################################################
# file ml.max #
kill (all)$

assume (DX>0)$
assume (DY>0)$
assume (coeff_x>0)$
assume (coeff_y>0)$

array(coe,4)$

list_eq :  [coe[4]+coe[3]+coe[2]+coe[1] = 1,
          
-coe[2]*DY/(2*coeff_y)-coe[1]*DY/(2*coeff_y)+coe[4]*DY/2+coe[3]*DY/2 =
0,
          
-coe[4]*DX/(2*coeff_x)+coe[3]*DX/(2*coeff_x)-coe[2]*DX/(2*coeff_x) +
coe[1]*DX/(2*coeff_x) = 0]$

list_un : [coe[1],coe[2],coe[3],coe[4]]$

sol : linsolve(list_eq,list_un)$

array(test,3)$
for i : 1 thru 3 do
(
  test[i] : lhs(list_eq[i]),
  for j : 1  thru 4 do
  (
    test[i] : subst(rhs(sol[j]),coe[j],test[i])
    ),
  test[i] : rat(test[i]),
  display(test[i])
  )$


list_eq_2 : [a+b+c+d=1, -b*DY/(2*f)-a*DY/(2*f)+d*DY/2+c*DY/2=0, 
-d*DX/(2*e)+c*DX/(2*e)-b*DX/(2*e)+a*DX/(2*e)=0]$
list_un_2 : [a,b,c,d]$

sol_2 : linsolve(list_eq_2,list_un_2)$

array(test_2,3)$
for i : 1 thru 3 do
(
  test_2[i] : lhs(list_eq_2[i]),
  test_2[i] : subst(rhs(sol_2[1]),a,test_2[i]),
  test_2[i] : subst(rhs(sol_2[2]),b,test_2[i]),
  test_2[i] : subst(rhs(sol_2[3]),c,test_2[i]),
  test_2[i] : subst(rhs(sol_2[4]),d,test_2[i]),
  test_2[i] : rat(test_2[i]),
  display(test_2[i])
  )$

list_eq_3 : [a+b+c+d=1, -b*DY/(2*f)-a*DY/(2*f)+d*DY/2+c*DY/2=0, 
-d*1/(2*e)+c*1/(2*e)-b*1/(2*e)+a*1/(2*e)=0]$
list_un_3 : [a,b,c,d]$

sol_3 : linsolve(list_eq_3,list_un_3)$

array(test_3,3)$
for i : 1 thru 3 do
(
  test_3[i] : lhs(list_eq_3[i]),
  test_3[i] : subst(rhs(sol_3[1]),a,test_3[i]),
  test_3[i] : subst(rhs(sol_3[2]),b,test_3[i]),
  test_3[i] : subst(rhs(sol_3[3]),c,test_3[i]),
  test_3[i] : subst(rhs(sol_3[4]),d,test_3[i]),
  test_3[i] : rat(test_3[i]),
  display(test_3[i])
  )$

list_eq_4 : [a+b+c+d=1, -b*1/(2*f)-a*1/(2*f)+d*1/2+c*1/2=0, 
-d*1/(2*e)+c*1/(2*e)-b*1/(2*e)+a*1/(2*e)=0]$
list_un_4 : [a,b,c,d]$

sol_4 : linsolve(list_eq_4,list_un_4)$

array(test_4,3)$
for i : 1 thru 3 do
(
  test_4[i] : lhs(list_eq_4[i]),
  test_4[i] : subst(rhs(sol_4[1]),a,test_4[i]),
  test_4[i] : subst(rhs(sol_4[2]),b,test_4[i]),
  test_4[i] : subst(rhs(sol_4[3]),c,test_4[i]),
  test_4[i] : subst(rhs(sol_4[4]),d,test_4[i]),
  test_4[i] : rat(test_4[i]),
  display(test_4[i])
  )$

list_eq_5 : [a+b+c+d=1, -f*b*1/(2*f)-f*a*1/(2*f)+f*d*1/2+f*c*1/2=0, 
-d*1/(2*e)+c*1/(2*e)-b*1/(2*e)+a*1/(2*e)=0]$
list_un_5 : [a,b,c,d]$

sol_5 : linsolve(list_eq_5,list_un_5)$

array(test_5,3)$
for i : 1 thru 3 do
(
  test_5[i] : lhs(list_eq_5[i]),
  test_5[i] : subst(rhs(sol_5[1]),a,test_5[i]),
  test_5[i] : subst(rhs(sol_5[2]),b,test_5[i]),
  test_5[i] : subst(rhs(sol_5[3]),c,test_5[i]),
  test_5[i] : subst(rhs(sol_5[4]),d,test_5[i]),
  test_5[i] : rat(test_5[i]),
  display(test_5[i])
  )$

list_eq_6 : [a+b+c+d=1, -b-a+f*d+f*c=0,  -d+c-b+a=0]$
list_un_6 : [a,b,c,d]$

sol_6 : linsolve(list_eq_6,list_un_6)$

array(test_6,3)$
for i : 1 thru 3 do
(
  test_6[i] : lhs(list_eq_6[i]),
  test_6[i] : subst(rhs(sol_6[1]),a,test_6[i]),
  test_6[i] : subst(rhs(sol_6[2]),b,test_6[i]),
  test_6[i] : subst(rhs(sol_6[3]),c,test_6[i]),
  test_6[i] : subst(rhs(sol_6[4]),d,test_6[i]),
  test_6[i] : rat(test_6[i]),
  display(test_6[i])
  )$
##########################################################

-- 
Raimondo Giammanco <giamma@vki.ac.be>
von Karman Institute