wolfram-mathematica - 这是数学中 NSolve 的错误吗？

Question

人们会期望并希望，如果您要求 Mathematica 求多项式的根，无论您是用符号来做，它都应该给出相同的(近似)答案，然后找到这些精确的数值近似值答案，或者你是否用数字来做。这是一个严重失败的示例(在 Mathematica 7 中，在 OS X 上运行):

poly = -112 + 1/q^28 + 1/q^26 - 1/q^24 - 6/q^22 - 14/q^20 - 25/q^18 -
  38/q^16 - 52/q^14 - 67/q^12 - 81/q^10 - 93/q^8 - 102/q^6 - 108/
  q^4 - 111/q^2 - 111 q^2 - 108 q^4 - 102 q^6 - 93 q^8 - 81 q^10 -
  67 q^12 - 52 q^14 - 38 q^16 - 25 q^18 - 14 q^20 - 6 q^22 - q^24 +
  q^26 + q^28;

Total[q^4 /. NSolve[poly == 0, q]] - Total[q^4 /. N[Solve[poly == 0, q]]]

(注意:这实际上是一个 Laurent 多项式，如果您乘以 q 的大幂，问题就会消失。)

这里的最后一行只是证明找到的解决方案非常不同；事实上，这是我们在解决问题时试图计算的数量。

如果您仔细观察 NSolve[poly == 0, q] 和 N[Solve[poly == 0, q] 的输出，您会发现您将看到 NSolve 只给出 54 根而不是预期的 56。这并不是说它只是错过了重复的根或其他任何东西；它缺少两个最大的幅度根(大约 +/- 1.59)

Is this a bug in Mathematica? Does anyone have an explanation for why this is happening?

Answer 1

很遗憾，您的期望是没有根据的。

正如您所说，Solve[] 给出了一个精确的解决方案，而 N[] 引入了一个小错误，但只有一次，当您评估它时。另一方面，NSolve[] 从一开始就使用数值近似值，因此会累积显着的舍入误差。

您还受到默认计算精度的限制，这可能导致数值方法完全失败，例如，缺少根(参见 Wilkinson's polynomial )。您可以通过告诉 NSolve[] 使用更高的精度来抵消这种情况，如下所示:

In[1]  := Total[q^4 /. NSolve[poly == 0, q, WorkingPrecision -> 50]] - 
          Total[q^4 /. N[Solve[poly == 0, q]]]
Out[1] := 0. - 3.66374*10^-15 I

在使用数值方法时，牢记错误始终很重要。因为对于大量的 numerical analysis 都是如此。从求解长多项式到对角化大矩阵到奇怪函数的积分等问题。没有一种正确的方法，需要告诉 Mathematica，例如，提高 WorkingPrecision，或应用不同的数值技术。