浮点数和小数的小数位问题。Decimal
-
08-07-2019 - |
题
我似乎失去了很多浮动的精度。
例如我需要求解一个矩阵:
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
).
有任何想法吗?
解决方案
我不熟悉Decimal类来帮助你,但你的问题是由于十进制分数通常不能用二进制精确表示的事实,所以你看到的是最接近的可能近似值;如果不使用特殊的类(比如Decimal),就没有办法避免这个问题。
编辑:
为什么十进制类不能正常工作?只要我从一个字符串开始,而不是浮动,权力似乎工作正常。
>>> import decimal
>>> print(decimal.Decimal("1.2") ** 2)
1.44
模块文档解释了 decimal.Decimal的必要性和用法
非常清楚,如果你还没有,你应该检查一下。
其他提示
IEEE 浮点数是二进制的,而不是十进制的。不存在恰好为 0.1 或其任意倍数的固定长度二进制小数。它是一个重复分数,如十进制的 1/3。
除了 Decimal 类之外的其他选项还有
使用 Common Lisp 或 Python 2.6 或另一种具有精确理性的语言
使用以下方法将双精度数转换为闭有理数: 弗拉普
首先,您的输入可以简化很多。您不需要读取和解析文件。您可以用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'
对于像这样的任务,我要小心十进制模块。其目的实际上是处理实际十进制数(例如,匹配人工簿记实践),具有有限的精度,而不是执行精确的精确数学。有些数字在十进制中不能完全表示,就像在二进制中一样,并且在十进制中执行算术也比替代算法慢得多。
相反,如果你想要精确的结果,你应该使用理性算术。这些将表示数字作为分子/解构器对,因此可以准确地表示所有有理数。如果你只使用乘法和除法(而不是像平方根这样的操作会导致无理数),你就永远不会失去精确度。
正如其他人所提到的,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
这不是你问题的答案,而是相关的:
#!/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
另见浮点错误的简单示例是什么,这里有SO,有一些答案。我给出的那个实际上使用python作为示例语言...