Problemi di posizione decimale con float e decimale. Decimale

Question

Sembra che stia perdendo molta precisione con i float.

Ad esempio, devo risolvere una matrice:

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

Questo è il codice che uso per importare la matrice da un file di testo:

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

Devo risolvere usando gauss-seidel, quindi ho bisogno di riorganizzare le equazioni per x, ye z:

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

Ecco il codice che uso per riordinare le equazioni. b è una matrice di coefficienti e y è il vettore di risposta:

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

Tuttavia le risposte che ottengo non sono precise al punto decimale.

Ad esempio, dopo aver riorganizzato la seconda equazione dall'alto, dovrei ottenere:

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

Quello che ottengo è:

y=-1.2-0.20000000000000001x+0.59999999999999998z

Questo potrebbe non sembrare un grosso problema, ma quando aumenti il ??numero a una potenza molto elevata l'errore è piuttosto grande. C'è un modo per aggirare questo? Ho provato la classe Decimal ma non funziona bene con i poteri (cioè Decimal (x) ** 2 ).

Qualche idea?

Answer 1

Non ho abbastanza familiarità con la classe Decimale per aiutarti, ma il tuo problema è dovuto al fatto che le frazioni decimali spesso non possono essere rappresentate con precisione in binario, quindi quello che vedi è l'approssimazione più vicina possibile; non c'è modo di evitare questo problema senza usare una classe speciale (come Decimal, probabilmente).

EDIT: Che dire della classe decimale non funziona correttamente per te? Finché inizio con una stringa, piuttosto che con un galleggiante, i poteri sembrano funzionare bene.

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

La documentazione del modulo spiega la necessità e l'uso di decimal.Decimal abbastanza chiaramente, dovresti verificarlo se non l'hai ancora fatto.

Answer 2

IEEE in virgola mobile è binario, non decimale. Non esiste una frazione binaria di lunghezza fissa esattamente pari a 0,1 o un suo multiplo. È una frazione ripetitiva, come 1/3 in decimale.

Si prega di leggere Che cosa dovrebbero sapere tutti gli scienziati informatici sull'aritmetica a virgola mobile

Altre opzioni oltre a una classe decimale sono

• utilizzando Common Lisp o Python 2.6 o un'altra lingua con razionali esatti

• convertendo i doppi in razionali ravvicinati usando, ad esempio, frap

Answer 3

Innanzitutto, il tuo input può essere semplificato molto. Non è necessario leggere e analizzare un file. Puoi semplicemente dichiarare i tuoi oggetti nella notazione Python. Valuta il file.

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 ]

Secondo, y = -1.2-0.20000000000000001x + 0.59999999999999998z non è insolito. Non esiste una rappresentazione esatta in notazione binaria per 0.2 o 0.6. Di conseguenza, i valori visualizzati sono approssimazioni decimali delle rappresentazioni originali non esatte. Questi sono veri per quasi ogni tipo di processore a virgola mobile che esiste.

Puoi provare il modulo razionale che potrebbe aiutare.

Sì, aumentare i numeri in virgola mobile ai poteri aumenta gli errori. Di conseguenza, devi essere sicuro di evitare di usare la posizione più a destra del numero in virgola mobile, poiché quei bit sono per lo più rumore.

Quando si visualizzano numeri in virgola mobile, è necessario arrotondarli in modo appropriato per evitare di vedere i bit di rumore.

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

Answer 4

Vorrei mettere in guardia contro il modulo decimale per attività come questa. Il suo scopo è davvero quello di affrontare i numeri decimali del mondo reale (ad es. Abbinando le pratiche di contabilità umana), con precisione finita, non eseguendo la matematica di precisione esatta. Esistono numeri che non sono esattamente rappresentabili in decimale proprio come in binario, e anche l'esecuzione dell'aritmetica in decimale è molto più lenta delle alternative.

Invece, se vuoi risultati esatti, dovresti usare l'aritmetica razionale. Questi rappresenteranno i numeri come coppia numeratore / denomentatore, quindi possono rappresentare esattamente tutti i numeri razionali. Se stai utilizzando solo la moltiplicazione e la divisione (anziché operazioni come le radici quadrate che possono causare numeri irrazionali), non perderai mai la precisione.

Come altri hanno già detto, python 2.6 avrà un tipo razionale incorporato, sebbene si noti che questa non è davvero un'implementazione ad alte prestazioni - per la velocità è meglio usare librerie come gmpy . Sostituisci semplicemente le tue chiamate a float () in gmpy.mpq () e il tuo codice dovrebbe ora dare risultati esatti (anche se potresti voler formattare i risultati come float per scopi di visualizzazione).

Ecco una versione leggermente riordinata del codice per caricare una matrice che utilizzerà invece razionali 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

Non è una risposta alla tua domanda, ma correlata:

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

Esempio:

$ 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

Vedi anche Che cos'è un semplice esempio di errore in virgola mobile , qui su SO, che ha alcune risposte. Quello che do in realtà usa Python come linguaggio di esempio ...