提高numpy.dot（Python）的精度

Question

我正在尝试模拟一个特定的物理系统。为了传播解决方案，我需要能够乘以行列式= 1的矩阵，这些矩阵描述系统的每个部分。在下面的代码中，T（变量）是一个行列式为1的二维矩阵，i表示区域编号，其余内容不相关。

当我运行具有超过30个区域的系统的此代码时，最终的Msys的行列式不再等于1。我在整个计算过程中检查了Msys的行列式值，并且在前几次迭代中它为1，但然后开始偏离1。我尝试在创建数组T时将dtype = float64放入，以查看是否会提高精度并防止其崩溃，但我没有看到任何改进。

是否有任何方法可以编写代码以避免错误累积，或者可以增加numpy存储的小数位数，使得对于具有100多个区域的系统，误差可以忽略不计？

for i in range(n):                                  
    if i == 0:                                      
        Msys = T(L[i],i,k)
    else:                                           
        Msys = numpy.dot(T(L[i]-L[i-1],i,k), Msys)
return Msys

- Zykx

---

如果可以的话，sympy 具有任意精度。 - jterrace

谢谢！最终使用了Sympy的一部分——Mpmath。 - Zykx

Answer 1

所有的浮点运算都有限制精度，误差会积累。您需要决定多少精度是“足够好”的，或者多少误差积累是“可忽略的”。如果float64对您来说不够精确，请尝试使用float128。您可以像这样查找浮点类型的精度：

In [83]: np.finfo(np.float32).eps
Out[83]: 1.1920929e-07

In [84]: np.finfo(np.float64).eps
Out[84]: 2.2204460492503131e-16

In [85]: np.finfo(np.float128).eps
Out[85]: 1.084202172485504434e-19

这里有更多关于浮点运算的信息：《计算机科学家应该知道的浮点运算》

- Bi Rico

---

在SO上一定有很多关于浮点数的问题。我记得有些例子可以证明两个方程式是相等的。只不过在其中一个方程式中，由于除以一个非常非常小的数字（正如大多数科学计算教科书所述），导致出现了错误。 - r4.

还要记住机器epsilon（计算机可以表示的最小数字）。 - r4.

那些还不够准确，但我找到了另一个名为mpmath的模块，它可以让你选择所需的精度。感谢您的帮助！ - Zykx