Aritmética finita y fuentes del error¶
Contenidos¶
Sistemas decimal y binario¶
En el sistema decimal el número $105.8125$ significa:
Es decir, es una suma cuyos términos son potencias de $10$ y utiliza $10$ dígitos $0,1,\ldots,9$.
Los ordenadores usan el sistema binario porque así los únicos dígitos que se almacenan son el $0$ y el $1$. En sistema binario los números se representan como sumas de potencias de $2$:
Conversión decimal a binario¶
Para convertir la parte entera, dividimos por $2$ de forma repetida. Los restos de estas divisiones son los dígitos en base $2$, de menos a más significativos.
El número $(105)_{10}$ en sistema binario es $(1101001)_2$.
Nota La función de python bin realiza esta conversión.
bin(105)
Escribir una función de2bi_a(x) que transforme la parte entera
de un número decimal x a binario.
La salida de la función será un string
- Empieza con inicializando el string
s = '' - Añade elementos a la izquierda del string con
s = (str(d)) + sdondedes0o1
Aplicarla al número x = 105.8125 y al x = 310.671875
NOTA: pueden ser útiles las siguientes funciones
np.fix(x): redondeaxhacia0.a//bcalcula el cociente de la división entera.a%bcalcula el resto de la división entera.
La estructura del programa en el fichero será
import numpy as np
def de2bi_a(x):
... # <- escribe aquí tu código
def main():
x = 105.8125
binary1 = de2bi_a(x)
print('La parte entera de ', x, ' en binario es ', binary1)
if __name__ == "__main__":
main()
Así podrás importar la función debi_a del fichero Ejercicio1a.py sin que se ejecute otra vez el código contenido en main
%run Ejercicio1a
Para convertir a binario la parte decimal, multiplicamos por $2$, sustraemos la parte entera, que será nuestro dígito binario y repetimos el proceso.
Escribir una función de2bi_b(x) que transforme la parte
fraccionaria de un número en base $10$ a base $2$.
La salida se almacenará en un string
- Empieza con inicializando el string
s = '' - Añade elementos a la derecha del string con
s = s + (str(d))dondedes0o1
Limitar el número de dígitos fraccionarios en binario a 23 (truncamiento).
Aplicarla a los números x = 105.8125, x = 310.671875 y x = 14.1.
Usando las dos funciones de2bi_a(x) y de2bi_b(x) escribir el número x en binario para los dos ejemplos anteriores
La estructura del programa en el fichero será
import numpy as np
import Ejercicio1a as ej1a
def de2bi_b(x):
... # <-- escribe aquí tu código
def main():
x = 105.8125
binary1 = ej1a.de2bi_a(x)
binary2 = de2bi_b(x)
print('La parte fraccionaria de ', x, ' en binario es ', '0.'+binary2)
print('El número ', x, 'en binario es ', binary1+'.'+binary2)
if __name__ == "__main__":
main()
%run Ejercicio1b
Conversión binario a decimal¶
De acuerdo con $(1)$ para convertir $(1101001.1101)_2$ a base $10$
1*(2**6)+1*(2**5)+0*(2**4)+1*(2**3)+0*(2**2)+0*(2**1)+1*(2**0)+1*(2**-1)+1*(2**-2)+0*(2**-3)+1*(2**-4)
Representación de enteros¶
Si tenemos $m$ dígitos o bits de memoria podemos almacenar $2^m$ números diferentes en binario.
Enteros positivos¶
Si sólo tomamos enteros positivos podríamos representar enteros con valores entre
Por ejemplo, para $m=3$ bits podríamos representar los enteros del $0$ al $7$.
Enteros con signo¶
Si queremos representar enteros con signo reservaremos el primer bit para el signo y entonces tendremos espacio para $m-1$ cifras significativas. Es decir $2^{m-1}-1$ números positivos, otros tantos negativos y dos ceros, uno con signo positivo y otro negativo. El rango de números representados sería de $[-2^{m-1}+1, 2^{m-1}-1]$. Si $m=3$ podremos representar números de $-3$ a $3$.
El primer bit es el signo, $0$ para los positivos y $1$ para los negativos.
Ejemplo
Calcular cómo se almacenará $(-80)_{10}$ si usamos representación binaria con una precisión de 8-bits con signo.
Solución: como es un número negativo, la representación es la misma que en el caso anterior pero el primer bit es el signo y como es negativo, es $1$. Es decir, $(-80)_{10}=(11010000)_2$
Representación al complemento de dos¶
Una forma de evitar la doble representación del cero es usar la representación al complemento de dos: los números negativos se construyen tomando los dígitos de los positivos y cambiando los $0$ a $1$ y los $1$ a $0$ y sumándole $1$. De esta manera la suma de un número y su opuesto es siempre $0$.
Ejemplo
Para representar $(-2)_{10}$ empezamos convirtiendo $(2)_{10}$ a binario $(010)_2$
- invertimos sus dígitos $\rightarrow (101)_2$
- le sumamos $001$ (tener en cuenta que en binario $0+0=0 \quad 0+1=1 \quad 1+1=10$) $\rightarrow (110)_2$
Si sumamos $(010)_2+(110)_2=(000)_2$
En este caso, el primer bit de los negativos es $1$ y el de los positivos es $0$. Hemos representado enteros en el rango de $[-2^{m-1},2^{m-1}-1]$.
Ejemplo
Calcular cómo se almacenará $(-80)_{10}$ si usamos representación binaria representación al complemento de dos con una precisión de 8-bits.
Solución: como $(80)_{10}=(01010000)_2$
- en el número anterior cambiamos los $0$ por $1$ y viceversa $\rightarrow 10101111 $
- y sumándole $1$ ya tenemos la representación $(-80)_{10}=(1011000)_2$
Representación sesgada¶
Otra forma es la representación sesgada. Los números negativos se representan con valores positivos consecutivos, empezando por el menor negativo, y los positivos toman los valores siguientes. La representación de los números la obtenemos añadiendo el sesgo $2^{m-1}$ al número $x_r=x+2^{m-1}\in[0,2^m-1]$.
De nuevo vemos que hemos representado enteros en el rango de $[-2^{m-1},2^{m-1}-1]$.
¿De cual de las formas descritas se almacenan los enteros? La mayor parte de las máquinas almacenan los enteros con la representación del complemento al dos y los exponentes de los números en punto flotante, que son enteros, con representación sesgada. ¿Por qué? Hay varios motivos. Uno de ellos es que los números con representación sesgada son más fáciles de comparar. Los códigos comparan frecuentemente números reales, y se empieza comparando el exponente y sólo si son iguales se compara la mantisa.
Representación de enteros en Python¶
Las variables numéricas en Python pueden ser integers, floating point y complex numbers. Las variables booleanas son un subtipo de los enteros. Los enteros no tienen precisión límite en Python 3.
Por ejemplo, en representación sesgada el mayor entero viene dado por $2^{m-1}-1$. Si $m=64$ es el número de bits el máximo entero sería
m = 64
print(2**(m-1)-1)
print(type(2**(m-1)-1))
Si $m=128$ fuera el número de bits el máximo entero sería
m = 128
print(2**(m-1)-1)
print(type(2**(m-1)-1))
Para $m=256$
m = 256
print(2**(m-1)-1)
print(type(2**(m-1)-1))
Representación de los números reales¶
En el caso de los números reales se utiliza una representación en punto flotante en base $\beta=2$.
donde:
- $s$ es el signo. Se almacena un $1$ si el número es negativo y $0$ si es positivo.
- $m$ es la mantisa que si está normalizada tiene un valor $1\leq m < \beta$ y el primer dígito ha de ser distinto de $0$. Este dígito en base dos solo puede ser $1$ y por lo tanto, en una mantisa normalizada siempre es $a_0=1$, no es necesario almacenarlo y ganamos un bit. Esto se llama técnica del bit escondido.
- $e$ es el exponente que es un entero con signo y utiliza representación binaria sesgada.
Los números se almacenan en "palabras" de 32 bits (precisión sencilla), 64
bits (precisión doble), o incluso 128 bits (precisión cuádruple). En casi todas las máquinas, la
precisión en Python para el tipo float es la precisión doble.
En el caso de precisión doble los bits se reparten:
- 1 bit para el signo.
- 11 bits para el exponente.
- 52 bits para la mantisa.
%run Doble_precision.py
Si tenemos $11$ bits para el exponente con signo, significa que tenemos espacio para $2^{11}=2048$ números en formato binario. Si empezamos en el cero, el $ 0 \le E \le 2047 $.
El primer valor $00000000000$ lo reservamos para cero y números
desnormalizados y el último $11111111111$ para Inf y NaN. Así que el rango de valores que puede tomar el exponente queda reducido a $ 1 \le E \le 2046 $.
Y como el sesgo es $2^{m-1}-1=2^{10}-1=1024-1=1023$ estos valores representan a los exponentes $ -1022 \le E - 1023 \le 1023 $.
Por lo tanto el exponente máximo es $e_{max}=1023$, y el mínimo $e_{min}=-1022$. Si dividimos por el valor mínimo:
y no se supera el valor máximo (overflow).
Escribir un programa que calcule la representación de los números los números 105.8125, 120.875, 7.1 y -1.41 en precisión simple, según la norma IEEE 754, utilizando las funciones del ejercicio 1. Considerar únicamente el caso en el que la parte entera es mayor que cero. Redondear por truncamiento.
Nota: Recordar que en precisión simple se utilizan 32 bits, de los cuales, uno es para el signo, 8 para el exponente y el resto para la mantisa. El sesgo del exponente es 127.
%run Ejercicio2
Las características del tipo float en nuestra máquina nos las da la orden sys.float_info
import sys
sys.float_info
Vemos que el número de dígitos de la mantisa, mant_dig es 53 (52 bits pero 53 dígitos debido al bit escondido). Es decir, el float que usa Python por defecto es el de 64 bits.
Pero si observamos esta salida dice
max_exp=1024 y min_exp=-1021
que no concuerda con lo dicho arriba.
Pero si ejecutamos
help(sys.float_info)
entre otras cosas obtenemos
| max_exp | DBL_MAX_EXP -- maximum int e such that radix**(e-1) is representable
| min_exp | DBL_MIN_EXP -- minimum int e such that radix**(e-1) is a normalized float
Como habíamos dicho, el primer $1$ no hace falta guardarlo y nos quedan 52 bits donde almacenamos los unos de $0.11\ldots 11$. Por lo tanto, en decimal este número será:
Hacer un programa Ejercicio4a.py que calcule esta suma. Como norma general, conviene sumar los elementos de la serie de menores a mayores (más adelante veremos por qué). Su valor tiene que coincidir con el valor obtenido con la orden python sys.float_info.max. Llamar sal al valor de salida.
import sys
%run Ejercicio3
print ('%.16e' % sal) # forma exponencial con 16 decimales
sys.float_info.max
Comprobamos que, efectivamente, son iguales
print(sal == sys.float_info.max)
El valor real positivo mínimo representable en punto flotante de forma normalizada es
sal = 2**-1022
print('%.16e' % sal)
Su valor tiene que coincidir con orden el valor de salida de python sys.float_info.min.
sys.float_info.min
Mayor entero representable exactamente de forma normalizada tal que todos los enteros positivos menores se pueden almacenar de forma exacta en forma normalizada
El mayor entero consecutivo, utilizando el formato en punto flotante, del que podemos almacenar todos los dígitos significativos y que, por tanto, podemos almacenar de forma exacta en binario podría ser:
Es decir, tenemos $53$ cifras de precisión en binario. El entero siguiente:
también podemos almacenarlo de forma exacta. Pero el siguiente a este ya no porque necesitaríamos un bit más a la derecha en la mantisa. Este entero es
print(2**53)
Todos los enteros menores se pueden almacenar de forma exacta. Si contamos, vemos que tiene 16 dígitos. Por lo tanto, podemos almacenar todos los enteros con 15 dígitos decimales y casi todos los enteros con 16 dígitos de forma exacta.
Pero, podría decirse, el número
es más grande y también es entero. Sin embargo sus $10$ últimas cifras serían forzosamente $0$ porque no las podemos almacenar y desde el punto de vista de la precisión no cuentan porque no les podemos dar el valor que queramos.
Valores desnormalizados
¿Qué sucede con los valores de los exponentes fuera del rango $[-1022,1023]$?
Si el exponente es menor que $-1022$ , el número es desnormalizado o cero y estamos usando para los bits correspondientes al exponente el valor especial $00000000000$. Entonces ya no se asume que el bit escondido es $1$ sino $0$.
x = 0.5*(2**-1023)
print(x) # desnormalizado
print(1/x) # su inverso
print(2**-1080 ) # Underflow
El menor valor desnormalizado es
es decir, $2^{-1022-52}$
num_min = 2.**(-1022-52)
print(num_min)
Si el valor es menor, devuelve es Underflow y devuelve cero.
2.**(-1022-53) # Underflow
Si el exponente es mayor que $1023$ devuelve OverflowError.
2.**1024 # Overflow
Precisión de la máquina¶
Vamos a calcular el número más pequeño que se puede representar después de $1$ utilizando la representación binaria en punto flotante con doble precisión:
El número real $1$ en representación en punto flotante de doble precisión normalizada (con $52$ ceros en la mantisa) es
y no podemos representar, de forma exacta en punto flotante, un número real entre los dos. La diferencia entre $1$ y $1+\epsilon$ es
que pasado a decimal es
1*2.**(-52)
o lo que es lo mismo
sys.float_info.epsilon
Este valor, el menor número $\epsilon$ tal que $1+\epsilon>1$ se llama precisión de la máquina y es una cota del error relativo que cometemos al almacenar cuanlquier número real en punto flotante, en este caso, en doble precisión.
Como, en doble precisión, $\epsilon\approx 2.22\cdot 10^{-16}$ se puede decir que corresponde, aproximadamente, a $16$ cifras decimales.
Si queremos el número más pequeño comparable con $x$ hay que sumarle un valor eps
x = 1.
eps = np.spacing(x)
print(eps)
x = 10.
eps = np.spacing(x)
print(eps)
x = 100.
eps = np.spacing(x)
print(eps)
Y vemos que la diferencia entre dos números consecutivos representables de forma exacta en punto flotante, aumenta conforme nos alejamos del cero. Por ello, la densidad de números representados en punto flotante no es uniforme. Es mayor cerca del cero y menor a medida que nos alejamos.
Resultados especiales¶
Algunas operaciones dan resultados especiales
x = 1.e200
y = x*x
print(y)
y/y
nan significa "not a number". Su aparición significa, en general,
que algo ha ido mal con el programa y/o que se ha realizado una operación
no válida.
Overflow. Ocurre cuando el resultado de la operación realizada es finita pero el resultado supera el rango de reales que pueden ser representados.
x = 1.e300
x**2
Aproximación y error de redondeo¶
Redondeo
A menudo, para un número real $x$, no exite una representación exacta en punto flotante y cae entre dos números racionales consecutivos $x^- \lt x \lt x^+$. Como representación de $x$ elegiremos uno de los dos dependiendo del método de redondeo usado. IEEE propone 5 sistemas de redondeo:
- Hacia arriba.
- Hacia abajo.
- Hacia cero ("truncamiento").
- Hacia el más alejado de cero más cercano.
- Hacia el par más cercano.
Este último método de redondeo es el más utilizado.
Pérdida de dígitos¶
Supongamos que estamos usando un sistema de representación decimal que almacena 3 cifras significativas y queremos sumar los números $a=123$ y $b=1.25$. El resultado es $s=124.25$ pero como solo podemos almacenar tres cifras redondeamos, por ejemplo, a $s_r=124$. Hemos perdido cifras y estamos cometiendo un error. En general, cualquier operación aritmética conlleva errores de redondeo.
Ejemplo
Si realizamos la operación $3\times 0.1 \times 100$:
3*0.1*100
El resultado no es $30$ debido al error de redondeo. Estos se pueden propagar y afectar hasta invalidar los resultados. En el caso anterior el error se debe a que, como vimos al principio, el número $x=0.1$ en base dos es periódico y por ello su representación en el ordenador no es exacta.
Ejemplo
Si calculamos $$\sum_{k=1}^{10000} 0.1$$
suma = 0.
for i in range(10000):
suma += 0.1
print ('%.16f'% suma) # En forma decimal con 16 decimales
y el error absoluto es
abs(suma-1000.)
Ejemplo
Si sumamos o restamos dos números muy diferentes en magnitud se pierde exactitud debido a los errores de redondeo. Por ejemplo
a = 1.e+9
epsilon = 1.e-8
suma = a
for i in range(10000):
suma += epsilon
print('%.16e'% suma) # En forma exponencial con 16 decimales
Sin embargo el resultado debería haber sido
suma = a + 10000*epsilon
print('%.16e'% suma)
Ejemplo
Vamos a sumar los $N$ primeros términos de la serie armónica $$\sum_{n=1}^N \frac{1}{n}$$
El resultado teórico para 1000 términos sería
from scipy import special
N = 10000
suma = special.polygamma(0,N+1)-special.polygamma(0,1)
print('%.20f'% suma)
Empezando la suma por el término $n=1$
suma1 = 0.;
for n in range(1,N+1):
suma1 += 1./n
print('%.20f'% suma1)
error1 = abs(suma1 - suma)
print(error1)
Empezando por el término final
suma2 = 0.;
for n in range(N,0,-1):
suma2 += 1./n
print('%.20f'% suma2)
La diferencia entre ambos resultados es
error2 = abs(suma2 - suma)
print(error2)
Es mayor el error en el caso primero porque los términos de la sucesión son cada vez más pequeños mientras que la suma es cada vez más grande. Por lo tanto, conforme avanzamos en la suma se producirá pérdida de dígitos por sumar números de magnitud muy distinta.
Error de cancelación¶
Al restar números muy parecidos usando aritmética finita donde, obviamente usamos la misma precisión, se produce el error de cancelación. Un ejemplo, suponiendo base decimal y precisión 7, si restamos dos números muy parecidos:
pasamos de tener precisión 7 a tener precisión 1 porque sólo nos queda una cifra significativa. Podemos añadir las demás cifras como ceros, pero realmente desconocemos el valor que tienen, porque tanto el minuendo como el sustraendo se supone que son aproximaciones de números almacenados o resultados de operaciones cuyo número de cifras reales pueden ser muchas más.
Ejemplo¶
Si usamos la fórmula
$$\sqrt{x^2+\epsilon^2}-x \quad \mathrm{con} \quad x \gg \epsilon$$tendremos un error de cancelación. Es decir, si usamos esta fórmula, tendremos un error grande
import numpy as np
x = 1000000.
ep = 0.01
a = np.sqrt(x**2 + ep**2) - x
print(a)
En este caso el error de calcelación se puede evitar utilizando una expresión matemática equivalente
b = ep**2 / (np.sqrt(x**2 + ep**2) + x)
print(b)
O bien
b = ep**2 / (2*x)
print(b)
El error relativo cometido al calcularlo con la primera expresión es muy grande
Er = abs((a - b)/b)
print(Er)
Además, este ejemplo muestra que expresiones matemáticamente equivalentes no son, necesariamente, equivalentes desde el punto de vista computacional. Y viceversa, que expresiones matemáticas que no son iguales, pueden tener el mismo resultado en el ordenador y ser, por tanto, computacionalmente equivalentes.
Ejemplo¶
Al resolver la ecuación de segundo grado
utilizando la fórmula
como
$$a=1,\,b=10^{8},\:c=1\qquad b^{2}\gg 4ac$$Y es de esperar que tengamos un error de cancelación. La solución que obtenemos es
a = 1.; b = 10.**8; c = 1.
x11 = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
print(x11)
Pero al sustituir en el polinomio para calcular el residuo
residuo11 = a*x11**2 + b*x11 + c
print(residuo11)
es muy grande comparado con la solución. No es una buena solución porque deberíamos haber obtenido
$$ax_{1}^{2}+bx_{1}+c\approx0$$y tenemos
$$ax_{1}^{2}+bx_{1}+c\approx0.255$$Reescribimos la fórmula
Multiplicamos y dividimos por la diferencia
como $(x+y)(x-y)=x^2-y^2$
Y la fórmula siguiente es equivalente a (1) y no tenemos el problema del error de cancelación
Como $b^{2}\gg 4ac$ podemos aproximar con $\sqrt{b^2-4ac}\approx \sqrt{b^2} \approx b$ y
x12 = -c/b
residuo12 = a*x12**2 + b*x12 + c
print(x12)
print(residuo12)
Ahora el residuo es mucho más pequeño que la solución
utilizando la fórmula
como
$$a=1,\,b=-10^{8},\:c=1\qquad b^{2}\gg 4ac$$como $b$ es negativo, $-b$ es positivo y $-\sqrt{b^2-4ac}$ es negativo y con un valor próximo a $-b$ se producirá un error de cancelación. La solución que obtenemos es
a = 1.; b = -10.**8; c = 1
x21 = (- b - np.sqrt(b**2 - 4*a*c)) / (2*a)
print(x21)
Pero al sustituir en el polinomio
residuo21 = a*x21**2 + b*x21 + c
print(residuo21)
es relativamente grande.
Reescribimos la fórmula
Teniendo en cuenta que si $b$ es negativo $-b$ es positivo y que $-\sqrt{b^2-4ac}$ es siempre negativo, la fórmula $(-b)-\sqrt{b^2-4ac}$ con $b^{2}\gg4ac$ dará un error de cancelación.
Multiplicamos y dividimos por la suma
Como $(x+y)(x-y)=x^2-y^2$
Y la fórmula siguiente es equivalente matemáticamente a la fórmula (2) pero no tiene el problema del error de cancelación (el denominador es la suma de dos números positivos porque $b$ es negativo)
Podemos simplificar la fórmula teniendo en cuenta que como $b^{2}\gg4ac$ podemos decir que $\sqrt{b^2-4ac}\approx \sqrt{b^2} = |b|=-b$
x22 = -c/b
residuo22 = a*x22**2 + b*x22 + c
print(x22)
print(residuo22)
Fuentes de error¶
Hay dos tipos de métodos numéricos:
- Métodos directos. Llegamos a la solución en un número finito de pasos. Por ejemplo, la solución de una ecuación de segundo grado que acabamos de ver.
- Métodos iterativos. Generamos una secuencia de soluciones aproximadas que tienden a la solución exacta.
Dependiendo de la fuente podemos tener:
- Errores de redondeo. Afecta a todos los métodos si usamos aritmética finita.
- Errores del método. Solo afecta a los métodos iterativos, no a los directos.
- Errores en los datos. Afectan a todos los métodos.
Error de truncamiento¶
Es un error del método. Es el error que obtenemos cuando aplicamos un método iterativo. En lugar de generar infinitos términos de la sucesión que tiende a la solución exacta, calculamos solo un número finito de elementos de la sucesión y asumimos que la solución tiene un error pero que es "suficientemente buena."
Ejemplo
Tenemos que
Si calculamos un número infinito de términos (y usamos precisión infinita) tendremos la solución exacta. Pero si solo calculamos un número finito de términos tendremos un error de truncamiento.
Aproximemos el valor de $e$, usando $x=1$, con un número finito de términos de la suma, por ejemplo, $50$
suma = 0.
x = 1.
fact = 1.
for n in range(100):
suma += x**n / fact
fact *= n+1
if n in np.array([5,10,15,20,40,60,80,100])-1:
error = abs(np.exp(1.) - suma)
print("Número de términos %i" % (n+1))
print("Error %.16e\n" % error)
Vemos que primero el error disminuye y luego se estanca. ¿Por qué? Si estuviéramos usando aritmética infinita el error disminuiría con las iteraciones. Pero como estamos usando aritmética finita, el error, en este caso, está limitado por el espaciado entre números representables en la zona de la solución $e$ que es
Ea = np.spacing(np.exp(1.))
print (Ea)
Calculemos el menor número de términos necesarios para obtener el número $e$ con el menor error posible.
Ea = np.spacing(np.exp(1.))
n = 0
Sum = 1.
fact = 1.
itermax = 100
error = np.abs(np.exp(1.) - Sum)
while (error > Ea and n < itermax):
n += 1
Sum += 1. / fact
fact *= n+1
error = np.abs(np.exp(1.) - Sum)
print("Número de iteraciones %i" % n)
print("Error %.16e" % error)
print("Espaciado alrededor de e %.16e" % Ea)
Y aunque sabemos que con $17$ iteraciones todavía cometemos error, este ya no es apreciable por la máquina y aunque realicemos más iteraciones este error computado no disminuye.
calcular el número de términos necesarios para calcular $\log 1.5$ con el menor error posible.
%run Ejercicio4
print("Número de iteraciones %i" % n)
print("Error %.16e" % error)
print("Espaciado alrededor de ln(1.5) %.16e" % Ea)
(a)
(b)
(c)
%run Ejercicio5