为何wxMaxima无法求解明显有解的方程组？

我有一组显然存在解的简单方程组，但wxMaxima调用solve却返回空列表，操作代码如下：

(%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)   []

另一个求解器MERCURY可以正常得到解，结果大致为：d≈8.815，h≈0.185，d1≈8.631。

问题原因 #

wxMaxima的solve函数在处理混合线性/非线性的方程组时，不会自动优先代入所有线性约束来简化非线性方程。原方程组中eq2、eq3、eq4都是线性关系，但eq1是关于d和h的非线性方程，直接用solve求解整个方程组时，可能因非线性项的存在导致求解器无法自动完成约束代入，从而返回空解。

解决方法 #

先手动代入线性约束，简化方程组后再求解：

• 代入d0=9到eq2，得到d1=9-2*h
• 将d0=9和d1=9-2*h代入eq4，得到d=(9 + 9-2h)/2 = 9 - h
• 把d=9-h代入eq1，得到单变量非线性方程：

  eq_simplified:43=%pi/4*(9-h)^2*h;
  

• 再用solve求解这个方程，转换为浮点数解：

  float(solve(eq_simplified,h));
  

  得到h的有效解后，再回代计算d和d1即可。

也可以用subst函数自动完成约束代入：

eqs_substituted:subst(d0=9, [eq1,eq2,eq4]);
eqs_simplified:subst(d1=9-2*h, eqs_substituted);
final_eq:subst(d=9-h, eqs_simplified[1]);
float(solve(final_eq,h));

这样就能得到正确的数值解了。

内容的提问来源于stack exchange，提问作者Michael Hecht