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

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

Мне нужно решить, используя 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).

Есть идеи?

Answer 1

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

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

Другие параметры, кроме десятичного класса,

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

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

Answer 2

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

EDIT: А как насчет десятичного класса не работает для вас должным образом? Пока я начинаю со строки, а не с плавающей запятой, полномочия, кажется, работают нормально.

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

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

Answer 3

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

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

Как уже упоминали другие, 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 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'

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 в качестве примера языка ...