lugar questões decimais com carros alegóricos e decimal.Decimal

Question

Parece que estou perdendo um monte de precisão com flutuadores.

Por exemplo, eu preciso resolver uma matriz:

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

Este é o código que usei para importar a matriz a partir de um arquivo de texto:

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

Eu preciso resolver usando Gauss-Seidel, então eu preciso reorganizar as equações para x, y, e z:

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

Aqui está o código que eu uso para reorganizar as equações. b é uma matriz de coeficientes e y é o vetor de resposta:

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

No entanto as respostas eu voltar não são precisos para a casa decimal.

Por exemplo, ao reorganizar a segunda equação de cima, devo começar:

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

O que eu recebo é:

y=-1.2-0.20000000000000001x+0.59999999999999998z

Isto pode não parecer um grande problema, mas quando você aumentar o número a uma potência muito alta o erro é muito grande. Existe uma maneira de contornar isso? Tentei a classe Decimal mas ele não funciona bem com poderes (ou seja, Decimal(x)**2).

Todas as idéias?

Answer 1

Eu não sou suficientemente familiarizado com a classe Decimal para ajudá-lo, mas seu problema é devido ao fato de que frações muitas vezes não podem ser precisas representado em binário, então o que você está vendo é a melhor aproximação possível; não há nenhuma maneira de evitar esse problema sem o uso de uma classe especial (como Decimal, provavelmente).

EDIT: E sobre a classe decimal não está funcionando corretamente para você? Enquanto eu começar com uma corda, em vez de um flutuador, poderes parecem funcionar bem.

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

O módulo de documentação explica a necessidade e uso de decimal.Decimal muito claramente, que deveria verificá-la, se você ainda não tiver.

Answer 2

IEEE ponto flutuante é binário, não decimal. Há não é fixo fracção binário comprimento que é exactamente 0,1, ou qualquer dos seus múltiplos. É uma fração repetindo, como 1/3 em decimal.

Por favor, leia O que cada cientista computador deve saber sobre Floating-Point Arithmetic

Outras opções além de uma classe Decimal são

• usando Lisp comum ou Python 2.6 ou outra língua com rationals exatas

• converter as duplas para rationals perto utilizando, por exemplo, frap

Answer 3

Em primeiro lugar, a sua entrada pode ser muito simplificado. Você não precisa de ler e analisar um arquivo. Você pode apenas declarar seus objetos na notação Python. Eval o arquivo.

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 ]

Em segundo lugar, y = -1.2-0.20000000000000001x + 0.59999999999999998z não é incomum. Não há nenhuma representação exata em notação binária para 0,2 ou 0,6. Consequentemente, os valores apresentados são as aproximações decimais do original não representações exatas. Aqueles são verdadeiras para praticamente todo o tipo de processador de ponto flutuante que existe.

Você pode tentar o Python 2.6 frações módulo. Há um pacote racional mais velhos que podem ajudar.

Sim, levantando números de ponto flutuante para poderes aumenta os erros. Consequentemente, você tem que ter a certeza de evitar o uso de mais à direita posições do número de ponto flutuante, uma vez que os bits são na maior parte do ruído.

Ao exibir números de ponto flutuante, você tem que apropriadamente volta deles para evitar ver os bits de ruído.

>>> a
0.20000000000000001
>>> "%.4f" % (a,)
'0.2000'

Answer 4

Eu alertam contra o módulo decimal para tarefas como esta. Sua finalidade é realmente mais lidar com números decimais do mundo real (por exemplo. Combinando as práticas de contabilidade humanos), com precisão finita, não realizar exata matemática precisão. Há números não exatamente representáveis ??em decimal, assim como há em binário, e realizar aritmética em decimal é também muito mais lento do que as alternativas.

Em vez disso, se você quer resultados exatos que você deve usar a aritmética racional. Estes irão representar números como um par numerador / denomentator, isso pode representar exatamente todos os números racionais. Se você estiver usando apenas multiplicação e divisão (em vez de operações como raízes quadradas que podem resultar em números irracionais), você nunca vai perder precisão.

Como já foi mencionado, python 2.6 terá um built-in tipo racional, embora note que este não é realmente um desempenho elevado de implementação - para a velocidade que você está usando melhor bibliotecas como gmpy . Basta substituir as chamadas para float () para gmpy.mpq () e seu código deve agora dar resultados exatos (embora você pode querer formatar os resultados como carros alegóricos para fins de exibição).

Aqui está uma versão ligeiramente arrumado do seu código para carregar uma matriz que vai usar rationals gmpy vez:

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

Não é uma resposta à sua pergunta, mas relacionado:

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

Exemplo:

$ 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

Veja também O que é um simples exemplo de erro de ponto flutuante , aqui no SO, que tem algumas respostas. O que eu realmente dar usa python como língua exemplo ...