3

我似乎在使用浮点数时失去了很多精度。

例如我需要求解一个矩阵：

4.0x -2.0y 1.0z =11.0
1.0x +5.0y -3.0z =-6.0
2.0x +2.0y +5.0z =7.0

这是我用来从文本文件中导入矩阵的代码：

f = open('gauss.dat')
lines =  f.readlines()
f.close()

j=0
for line in lines:
    bits = string.split(line, ',')
    s=[]
    for i in range(len(bits)):
        if (i!= len(bits)-1):
            s.append(float(bits[i]))
            #print s[i]
    b.append(s)
    y.append(float(bits[len(bits)-1]))

我需要使用 gauss-seidel 求解，所以我需要重新排列 x、y 和 z 的方程：

x=(11+2y-1z)/4
y=(-6-x+3z)/5
z=(7-2x-2y)/7

这是我用来重新排列方程式的代码。b是一个系数矩阵，y是答案向量：

def equations(b,y):
    i=0
    eqn=[]
    row=[]
    while(i<len(b)):
        j=0
        row=[]
        while(j<len(b)):
            if(i==j):
                row.append(y[i]/b[i][i])
            else:
                row.append(-b[i][j]/b[i][i])
            j=j+1
        eqn.append(row)
        i=i+1
    return eqn

然而，我得到的答案并不精确到小数位。

例如，从上面重新排列第二个等式后，我应该得到：

y=-1.2-.2x+.6z

我得到的是：

y=-1.2-0.20000000000000001x+0.59999999999999998z

这似乎不是一个大问题，但是当您将数字提高到非常高的幂时，误差会很大。有没有解决的办法？我尝试了这Decimal门课，但它不适用于权力（即，Decimal(x)**2）。

有任何想法吗？

4

6 回答 6

14

IEEE 浮点数是二进制的，而不是十进制的。没有固定长度的二进制小数正好是 0.1 或其任何倍数。它是一个重复的分数，就像十进制的 1/3。

请阅读每位计算机科学家应了解的浮点运算知识

除 Decimal 类之外的其他选项是

于 2008-11-13T02:11:58.567 回答
12

我不太熟悉 Decimal 类来帮助你，但你的问题是由于十进制分数通常不能准确地用二进制表示，所以你看到的是最接近的近似值；如果不使用特殊类（可能是 Decimal），就无法避免这个问题。

EDIT:十进制类不适合您怎么办？只要我从一个字符串而不是一个浮点数开始，权力似乎就可以正常工作。

>>> import decimal
>>> print(decimal.Decimal("1.2") ** 2)
1.44

模块文档非常清楚地解释了对的需求和用法，decimal.Decimal如果你还没有，你应该检查一下。

于 2008-11-13T02:09:08.667 回答
4

首先，您的输入可以简化很多。您不需要读取和解析文件。你可以只用 Python 表示法声明你的对象。评估文件。

b = [
    [4.0, -2.0,  1.0],
    [1.0, +5.0, -3.0],
    [2.0, +2.0, +5.0],
]
y = [ 11.0, -6.0, 7.0 ]

其次，y=-1.2-0.20000000000000001x+0.59999999999999998z 并不罕见。0.2 或 0.6 没有二进制表示法的精确表示。因此，显示的值是原始不精确表示的十进制近似值。几乎所有类型的浮点处理器都是如此。

您可以尝试 Python 2.6分数模块。有一个较旧的理性包可能会有所帮助。

是的，将浮点数提高到幂会增加错误。因此，您必须确保避免使用浮点数的最右边位置，因为这些位主要是噪声。

显示浮点数时，您必须对它们进行适当的舍入以避免看到噪声位。

>>> a
0.20000000000000001
>>> "%.4f" % (a,)
'0.2000'
于 2008-11-13T02:49:01.113 回答
4

对于这样的任务，我会注意十进制模块。它的目的实际上是更多地处理现实世界的十进制数（例如，匹配人类簿记实践），具有有限的精度，而不是执行精确的精度数学。有些数字不能像二进制那样精确地用十进制表示，而且用十进制执行算术也比其他方法慢得多。

相反，如果你想要精确的结果，你应该使用有理算术。这些将数字表示为分子/分母对，因此可以准确地表示所有有理数。如果您只使用乘法和除法（而不是可能导致无理数的平方根之类的运算），您将永远不会失去精度。

正如其他人所提到的，python 2.6 将有一个内置的有理类型，但请注意，这并不是一个真正的高性能实现——为了速度，你最好使用像gmpy这样的库。只需将您对 float() 的调用替换为 gmpy.mpq()，您的代码现在应该会给出准确的结果（尽管您可能希望将结果格式化为浮点数以进行显示）。

这是您代码的稍微整理的版本，用于加载将使用 gmpy 有理数的矩阵：

def read_matrix(f):
    b,y = [], []
    for line in f:
        bits = line.split(",")
        b.append( map(gmpy.mpq, bits[:-1]) )
        y.append(gmpy.mpq(bits[-1]))
    return b,y
于 2008-11-13T14:30:33.843 回答
2

这不是您问题的答案，而是相关的：

#!/usr/bin/env python
from numpy import abs, dot, loadtxt, max
from numpy.linalg import solve

data = loadtxt('gauss.dat', delimiter=',')
a, b = data[:,:-1], data[:,-1:]
x = solve(a, b) # here you may use any method you like instead of `solve`
print(x)
print(max(abs((dot(a, x) - b) / b))) # check solution

例子：

$ cat gauss.dat
4.0, 2.0, 1.0, 11.0
1.0, 5.0, 3.0, 6.0 
2.0, 2.0, 5.0, 7.0

$ python loadtxt_example.py
[[ 2.4]
 [ 0.6]
 [ 0.2]]
0.0
于 2008-11-25T12:28:41.757 回答
0

另请参阅What is a simple example of floating point error，在 SO 上，其中有一些答案。我给出的那个实际上使用python作为示例语言......

于 2008-11-13T02:49:57.320 回答