Проблемы с десятичными знаками с плавающей запятой и десятичной дробью.Десятичная дробь

Question

Кажется, я сильно теряю точность при работе с поплавками.

Например, мне нужно решить матрицу:

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]))

Мне нужно решить, используя Гаусса-Зайделя, поэтому мне нужно переставить уравнения для 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).

Есть какие-нибудь идеи?

Answer 1

Я недостаточно знаком с классом Decimal, чтобы помочь вам, но ваша проблема связана с тем фактом, что десятичные дроби часто не могут быть точно представлены в двоичном виде, поэтому то, что вы видите, является наиболее близким приближением; нет способа избежать этой проблемы без использования специального класса (например, Decimal, вероятно).

РЕДАКТИРОВАТЬ: А как насчет десятичного класса не работает для вас должным образом? Пока я начинаю со строки, а не с плавающей запятой, кажется, что полномочия работают нормально.

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

документация по модулю объясняет необходимость использования decimal.Decimal и его использование. довольно ясно, вы должны проверить это, если вы еще не сделали.

Answer 2

IEEE с плавающей запятой является двоичным, а не десятичным.Не существует двоичной дроби фиксированной длины, равной ровно 0,1, или любой кратной ей.Это повторяющаяся дробь, как 1/3 в десятичной системе счисления.

Пожалуйста, прочтите Что должен знать каждый специалист по информатике Об арифметике с плавающей запятой

Другими параметрами, помимо класса Decimal, являются

• используя Common Lisp или Python 2.6 или другой язык с точными рациональными выражениями

• преобразование двойных чисел в близкие рациональные, используя, например, фрап

Answer 3

Во-первых, ваш вклад может быть значительно упрощен. Вам не нужно читать и анализировать файл. Вы можете просто объявить свои объекты в нотации 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'

Answer 4

Я бы предостерег от использования десятичного модуля для подобных задач. Его цель - в большей степени иметь дело с реальными десятичными числами (например, сопоставление с практикой ведения бухгалтерского учета), с конечной точностью, а не с точной математикой точности. Есть числа, которые не могут быть точно представлены в десятичной системе, как в двоичной системе, и выполнение арифметики в десятичной системе также намного медленнее, чем альтернативы.

Вместо этого, если вы хотите получить точные результаты, вы должны использовать рациональную арифметику. Они будут представлять числа в виде пары числитель / деноменатор, поэтому могут точно представлять все рациональные числа. Если вы используете только умножение и деление (а не такие операции, как квадратные корни, которые могут привести к иррациональным числам), вы никогда не потеряете точность.

Как уже упоминали другие, python 2.6 будет иметь встроенный рациональный тип, хотя обратите внимание, что это на самом деле не высокопроизводительная реализация - для скорости вам лучше использовать библиотеки, такие как gmpy . Просто замените вызовы float () на gmpy.mpq (), и ваш код теперь должен давать точные результаты (хотя вы можете отформатировать результаты как плавающие для отображения).

Вот немного приведенная в порядок версия вашего кода для загрузки матрицы, в которой вместо этого будут использоваться gmpy rationals:

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

Answer 5

Это не ответ на ваш вопрос, но связанный с ним:

#!/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

Answer 6

См. также Что является простым примером ошибки с плавающей запятой , здесь на SO, у которого есть несколько ответов Тот, что я даю, на самом деле использует Python в качестве примера языка ...