Search code examples
solvermaxima

Why no solution in wxMaxima


I have this simple set of equations which obviously has a solution. Why doesn't wxMaxima find the solution?

(%i24)  kill(all);
(%o0)   done
(%i1)   eq1:43=%pi/4*d*d*h;
(eq1)   43=(%pi*d^2*h)/4
(%i2)   eq2:d1=d0-2*h;
(eq2)   d1=d0-2*h
(%i3)   eq3:d0=9;
(eq3)   d0=9
(%i4)   eq4:d=(d0+d1)/2;
(eq4)   d=(d1+d0)/2
(%i8)   float(solve([eq1,eq2,eq3,eq4],[d,d1,h]));
(%o8)   []

The same in my favourite equation system solver MERCURY enter image description here


Solution

  • This looks like a bug in solve. I'm not sure what the problem is, but solve can find the solutions with just equations eq2, eq3, and eq4:

    (%i90) solve([eq2,eq3,eq4],[d,d1,h,d0]);
    (%o90)         [[d = %r3, d1 = 2 %r3 - 9, h = 9 - %r3, d0 = 9]]
    

    Then substitute this into eq1:

    (%i91) subst(%o90, eq1);
                                                      2
                                     %pi (9 - %r3) %r3
    (%o91)                      43 = ------------------
                                             4
    

    Finally, solve can produce roots of this cubic equation. It's really messy, and it seems you want a numerical solution, so:

    (%i93) allroots(float(%o91));
    (%o93) [%r3 = 3.027752659221862, %r3 = - 2.209973166894623, 
                                                            %r3 = 8.18222050767276]
    

    We can substitute these 3 roots back into the solutions found in %o90 to get the desired results:

    (%i94) map(lambda([r], subst(r,%o90)), %o93);
    
    (%o94) [[[d = 3.027752659221862, d1 = - 2.944494681556276, 
    h = 5.972247340778138, d0 = 9]], [[d = - 2.209973166894623, 
    d1 = - 13.419946333789246, h = 11.209973166894624, d0 = 9]], 
    [[d = 8.18222050767276, d1 = 7.36444101534552, h = 0.8177794923272401, 
    d0 = 9]]]
    

    Maxima should be able to figure this out for itself, but it can't. The last solution [d = 8.18222050767276, d1 = 7.36444101534552, h = 0.8177794923272401, d0 = 9] is the answer you get with the other software.

    This issue is being tracked as Maxima issue 4247. And curiously, if you replace %pi by, say, pp. solve can find a solution. It's quite messy since it involves the solution of a cubic, but converting the answer to floating point (and rectform) produces the same answer as above, except there are some spurious imaginary parts that are very close to zero.