Decimales problemas de lugar con flotantes y decimal.Decimal

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 utilizo 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 utilizo para reorganizar las ecuaciones. b es una matriz de coeficientes e 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 lugar 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 aumenta el número a una potencia muy alta, el error es bastante grande. ¿Hay alguna forma de evitar esto? Decimal clase Decimal pero no funciona bien con poderes (es decir, Decimal(x)**2 ).

¿Algunas ideas?

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 se pueden representar 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 ti? Mientras comienzo con una cadena, en lugar de una flotación, 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 de decimal.Decimal bastante clara, debe verificarlo si aún no lo ha hecho.

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 que se repite, como 1/3 en decimal.

Lea lo que todo científico informático debería 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 lenguaje con racionales exactos

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

En primer lugar, su entrada se puede simplificar mucho. No necesitas leer y analizar un archivo. Simplemente puede declarar sus objetos en notación de Python. Evaluar 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.5999999999999999998z no es inusual. No hay 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 no exactas originales. Esas son ciertas para casi todos los tipos de procesadores de punto flotante que existen.

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

Sí, elevar los números de punto flotante a potencias 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 en su mayoría 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' 

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

En su lugar, si quieres resultados exactos debes usar aritmética racional. Estos representarán números como un par de numerador / denominador, por lo que pueden representar exactamente todos los números racionales. Si solo está utilizando la multiplicación y la división (en lugar de operaciones como las raíces cuadradas que pueden resultar en números irracionales), nunca perderá precisión.

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

Aquí hay una versión un poco ordenada de su código para cargar una matriz que usará gmpy racionals en su lugar:

 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 

No es una respuesta a tu 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 

También vea ¿Qué es un ejemplo simple de error de punto flotante , aquí en SO, que tiene algunas respuestas? El que doy en realidad usa python como lenguaje de ejemplo …