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

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

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

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

Любые идеи?


person darudude    schedule 13.11.2008    source источник


Ответы (6)


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

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

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

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

person Jeremy    schedule 13.11.2008
comment
Проблема, с которой я столкнулся ранее с мощностью, заключалась в том, что я пытался поднять ее до степени 0,5. Для этого мне пришлось написать decimal.Decimal (1.2) ** decimal.Decimal (0,5) - person darudude; 13.11.2008
comment
Вы упустили цитату. Может, в этом и была проблема. Если вы нашли его слишком подробным, вы всегда можете D = decimal.Decimal; D (1,2) ** D (0,5) - person recursive; 02.01.2009

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

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

Другие варианты помимо класса Decimal:

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

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

person Doug Currie    schedule 13.11.2008
comment
@Doug: Я добавил ссылку на Python 2.6 и его модуль фракций. - person tzot; 13.11.2008
comment
Не нужно менять язык - просто используйте библиотеку рациональной арифметики, например gmpy. - person Brian; 13.11.2008
comment
Что значит кратное? Есть двоичная дробь фиксированной длины, равная точно 0,5. - person Nosredna; 02.06.2009
comment
Да, но не существует двоичной дроби фиксированной длины, равной 0,1 (одной десятой) или кратной 0,1. - person Doug Currie; 02.06.2009
comment
@Doug Currie 0,5 кратно 0,1, поэтому последняя часть неверна;) - person Voo; 08.09.2011

Во-первых, ваш ввод можно значительно упростить. Вам не нужно читать и разбирать файл. Вы можете просто объявить свои объекты в нотации 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'
person S.Lott    schedule 13.11.2008
comment
Мне нужно ввести данные из файла, потому что моя программа должна быть достаточно универсальной, чтобы обрабатывать матрицу NxN. Также не% .4f формат форматирования строки? это не кажется очень чистым - person darudude; 13.11.2008
comment
Обозначения Python будут описывать матрицу ЛЮБОГО размера. Зачем писать собственный парсер? Зачем искать и конвертировать числа с плавающей запятой? Просто используйте нотацию Python для своей матрицы. - person S.Lott; 13.11.2008
comment
Все числа преобразуются в строки для печати. Ваш источник находится в десятичной системе счисления. (0,2). Python работает в двоичном приближении к этому. Когда вы запрашиваете вывод - любой вывод - он конвертируется обратно в строку. Всегда. Запрашивайте в этой строке значимые цифры, а не ВСЕ цифры. - person S.Lott; 13.11.2008

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

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

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

Вот слегка урезанная версия вашего кода для загрузки матрицы, которая вместо этого будет использовать рациональные числа 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
person Brian    schedule 13.11.2008

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

#!/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
person jfs    schedule 25.11.2008

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

person Matthew Schinckel    schedule 13.11.2008