Problemas con lugares decimales con flotantes y decimales.

Question

Parece que estoy perdiendo mucha precisión con los flotadores.

Por ejemplo, necesito resolver una 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 es el código que uso para importar la matriz desde un archivo 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]))

Necesito resolver usando gauss-seidel, así que necesito reorganizar las ecuaciones para x, y y z:

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

Aquí está el código que uso para reorganizar las ecuaciones. b es una matriz de coeficientes y y es el vector de respuesta:

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

Sin embargo, las respuestas que recibo no son precisas para el decimal.

Por ejemplo, al reorganizar la segunda ecuación desde arriba, debería obtener:

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

Lo que obtengo es:

y=-1.2-0.20000000000000001x+0.59999999999999998z

Esto puede no parecer un gran problema, pero cuando eleva el número a una potencia muy alta, el error es bastante grande. ¿Hay alguna forma de evitar esto? Probé la clase Decimal pero no funciona bien con poderes (es decir, Decimal (x) ** 2 ).

¿Alguna idea?

Answer 1

No estoy lo suficientemente familiarizado con la clase Decimal para ayudarlo, pero su problema se debe al hecho de que las fracciones decimales a menudo no pueden representarse con precisión en binario, por lo que lo que está viendo es la aproximación más cercana posible; no hay forma de evitar este problema sin usar una clase especial (como Decimal, probablemente).

EDIT: ¿Qué pasa con la clase decimal que no funciona correctamente para usted? Mientras empiezo con una cuerda, en lugar de un flotador, los poderes parecen funcionar bien.

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

La documentación del módulo explica la necesidad y el uso de decimal.Decimal con bastante claridad, deberías echarle un vistazo si aún no lo has hecho.

Answer 2

El punto flotante IEEE es binario, no decimal. No hay una fracción binaria de longitud fija que sea exactamente 0.1, o cualquier múltiplo de la misma. Es una fracción repetitiva, como 1/3 en decimal.

Lea Lo que todo informático debe saber sobre la aritmética de punto flotante

Otras opciones además de una clase Decimal son

• usando Common Lisp o Python 2.6 u otro idioma con razones exactas

• convertir los dobles para cerrar racionales usando, por ejemplo, frap

Answer 3

Primero, su entrada puede simplificarse mucho. No necesita leer y analizar un archivo. Simplemente puede declarar sus objetos en notación Python. Evalúe el archivo.

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 ]

Segundo, y = -1.2-0.20000000000000001x + 0.59999999999999998z no es inusual. No hay una representación exacta en notación binaria para 0.2 o 0.6. En consecuencia, los valores mostrados son las aproximaciones decimales de las representaciones originales no exactas. Esos son ciertos para casi todos los tipos de procesadores de punto flotante que hay.

Puede probar el módulo Python 2.6 fracciones . Existe un paquete racional más antiguo que podría ayudar.

Sí, elevar los números de coma flotante a potencia aumenta los errores. En consecuencia, debe asegurarse de evitar el uso de las posiciones más a la derecha del número de punto flotante, ya que esos bits son principalmente ruido.

Al mostrar números de punto flotante, debe redondearlos adecuadamente para evitar ver los bits de ruido.

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

Answer 4

Advertiría contra el módulo decimal para tareas como esta. Su propósito es realmente más tratar con números decimales del mundo real (por ejemplo, hacer coincidir las prácticas de contabilidad humana), con precisión finita, sin realizar matemática de precisión exacta. Hay números que no son exactamente representables en decimal como lo hay en binario, y realizar aritmética en decimal también es mucho más lento que las alternativas.

En cambio, si desea resultados exactos, debe usar aritmética racional. Estos representarán números como un par numerador / denomentador, por lo que pueden representar exactamente todos los números racionales. Si solo usa multiplicación y división (en lugar de operaciones como raíces cuadradas que pueden dar como resultado números irracionales), nunca perderá precisión.

Como otros han mencionado, python 2.6 tendrá un tipo racional incorporado, aunque tenga en cuenta que esta no es realmente una implementación de alto rendimiento, por velocidad, es mejor usar bibliotecas como gmpy . Simplemente reemplace sus llamadas a float () a gmpy.mpq () y su código ahora debería dar resultados exactos (aunque es posible que desee formatear los resultados como flotantes para fines de visualización).

Aquí hay una versión ligeramente ordenada de su código para cargar una matriz que utilizará en su lugar racionales 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

Answer 5

No es una respuesta a su pregunta, sino relacionada:

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

Ejemplo:

$ 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

Consulte también ¿Cuál es un ejemplo simple de error de coma flotante? , aquí en SO, que tiene algunas respuestas. El que doy en realidad usa Python como el lenguaje de ejemplo ...