Nociones básicas de errores en Cálculo Numérico
Elena E. Álvarez Saiz
Red Educativa Digital Descartes.
Universidad de Cantabria
Título de la obra:
NOCIONES BÁSICAS DE ERRORES
EN CÁLCULO NUMÉRICO
Autores:
Elena E. Alvarez Saiz
Diseño del libro: Juan Guillermo Rivera Berrío
Código JavaScript para el libro: Joel Espinosa Longi, IMATE, UNAM.
Recursos interactivos: DescartesJS
Fuentes: Lato y UbuntuMono
Fórmulas matemáticas: $\KaTeX$
Primera Edición
LICENCIA
Creative Commons Attribution License 4.0.En un curso de cálculo numérico se estudian por lo general métodos o modelos matemáticos que permiten hallar la solución numérica a problemas relacionados entre otros con Ecuaciones diferenciales, Integrales definidas, Sistemas de ecuaciones simultáneas, Raíces de ecuaciones, Interpolación. Como criterio general de trabajo, el objetivo será llegar a la solución correcta (con un error aceptable) de forma eficiente (en un tiempo razonable).
En la aproximación de un fenómeno físico mediante un modelo matemático y su posterior solución numérica, se cometen diferentes clases de errores. Algunos de ellos se debe a errores en los datos iniciales (por ejemplo de medida con algún instrumento), errores de discretización del problema, errores de redondeo debidos a la forma en la que el computador almacena los números, etc. El objetivo de este tema es analizar este último tipo de errores.
El 25 de febrero de 1991, durante la Guerra del Golfo, una batería de misiles Patriot estadounidense en Dharan, Arabia Saudita, no pudo interceptar un misil Scud iraquí entrante. El Scud atacó un cuartel del ejército estadounidense y mató a 28 soldados. Un informe de la entonces Oficina General de Contabilidad titulado "Patriot Missile Defense: Problema de software conducido a una falla del sistema en Dhahran, Arabia Saudita" informó sobre la causa del fallo: un cálculo inexacto del tiempo desde el inicio debido a errores aritméticos del computador.
Específicamente, el tiempo en décimas de segundo, medido por el reloj interno del sistema, se multiplicó por 10 para producir el tiempo en segundos. Este cálculo se realizó utilizando un registro de punto fijo de 24 bits. En particular, el valor 1/10, que tiene una expansión binaria infinita, se cortó a 24 bits después del punto de base. El pequeño error de corte, cuando se multiplica por el gran número que da el tiempo en décimas de segundo, conduce a un error significativo Douglas N. Arnold, “ The Patriot Missile Failure ,” Institute for Mathematics and its Applications,Minneapolis (2000). Fuente http://www.ima.umn.edu/~arnold/disasters/patriot.html.
El Programa Ariane, es un programa emprendido por Europa en 1973 para dotarse de un lanzador que le permitiera un acceso independiente al espacio. El desarrollo del lanzador Ariane se efectúa bajo la dirección de la Agencia Espacial Europea (ESA, por sus siglas en inglés); el Centro Nacional de Estudios Espaciales(CNES) francés fue el contratista principal hasta mayo de 2003, fecha en la que pasó a actuar como tal el consorcio europeo EADS (European Aeronautic Defenceand Space Company). La explotación comercial es gestionada por la sociedad Arianespace, creada en 1980. En total, unas 40 compañías europeas están comprometidas en el desarrollo y construcción del lanzador Ariane.
Todos los lanzamientos se efectúan desde el centro espacial de Kourou, en la Guayana Francesa, que cuenta con varias rampas de lanzamiento, y en el que trabajan de forma permanente varios cientos de personas.
El 4 de junio de 1996, el cohete Ariane 5 Flight 501,de la Agencia Europea del Espacio (ESA) explotó 40 segundos después de su despegue a una altura de 3.7km. tras desviarse de la trayectoria prevista.
Era su primer viaje tras una década de investigación que costó más de 7000 millones de euros. El cohete y su carga estaban valorados en más de 500 millones de euros. La causa del error fue un fallo en el sistema de guiado de la trayectoria provocado 37 segundos después del despegue. Este error se produjo en el software que controlaba el sistema de referencia inercial (SRI).
En concreto, se produjo una excepción debido al intento de convertir un número en punto flotante de 64 bits, relacionado con la velocidad horizontal del cohete respecto de la plataforma de lanzamiento, en un entero con signo de 16 bits. El número mas grande que se puede representarde esta forma es 32767 Douglas N. Arnold, “ The Explosion of the Ariane 5 Rocket ,” Institute for Mathematics and its Applications,Minneapolis (2000). Fuente Fuente http://www.ima.umn.edu/~arnold/disasters/ariane.html.
Más información en: SIAM News, Vol. 29. Number 8, October 1996
El 23 de agosto de 1991, la plataforma petrolífera Sleipner A propiedad de la empresa noruega Statoil situada en el mar del Norte a 82 metros de profundidad se hundió. La causa del error fue un fallo en el modelado numérico de la plataforma utilizando elementos finitos. En concreto se produjo una fuga de agua en una de las paredes de uno de los 24 tanques de aire de 12 metros de diámetro que permitían la flotación de la plataforma de 57000 toneladas de peso que además soportaba a mas de 20 0personas y a equipamiento de extracción con un peso adicional de unas 40000 toneladas. Las bombas de extracción de agua no fueron capaces de evacuar toda el agua. El fallo tuvo un coste económico total de 700 millones de euros.
Para modelar los tanques de la plataforma se utilizó el programa NASTRAN de elementos finitos y una aproximación mediante un modelo elástico lineal. Esta aproximación no era correcta y subestimó en un 47% los esfuerzos que debían soportar las paredes de los tanques. En particular, ciertas paredes fueron diseñadas con un grosor insuficiente. Un análisis con elementos finitos correcto, pero posterior al accidente, demostró que el diseño de la plataforma provocaría fugas en algunas de los tanques cuando ésta estuviese sobre agua a 62 metros de profundidad.La fuga real se produjo cuando ésta estaba sobre 65 metros de agua, lo que ratifica la supuesta causa del fallo. Douglas N. Arnold, “ The sinking of the Sleipner A o®shore platform ,” Institute for Mathematics and itsApplications, Minneapolis (2000). Fuente http://www.ima.umn.edu/~arnold/disasters/sleipner.html.
El objetivo de este capítulo es descrbir como el computador almacena los números y opera con ellos. Esto es fundamental para comprender los errores que comete durante este proceso y cómo afecta al resultado de un algoritmo numérico.
Ver ejemplo 1 a 6 de los apuntes.
Ejemplo. Analiza el siguiente código Matlab .
El resultado debería ser $100$ pero se obtiene un número próximo a $100$: $99.9999999999986e+000$
Ejemplo. Analiza el siguiente código Matlab.
¿Por qué ocurre esto? La razón es que los ordenadores usan un formato (punto flotante binario) que no permite representar de forma precisa números como por ejemplo 0.1, 0.2 o 0.3. Dado que los números están, en general, representados de manera inexacta entonces las operaciones aritméticas también se hacen de manera inexacta.
En este sistema los números se representan utilizando 10 dígitos, es un sistema de numeración posicional en el que las cantidades se representan utilizando como base aritmética las potencias del número diez. Así,
$$x \in \mathbb{R}^{+} \,\,\,\,\, x = \sum\limits_{ - \infty }^n {{a_k} \cdot {{10}^k}} \,\,\,\,\, 0 \le {a_k} \le 9$$Ejemplo. Considera el número $709.3045$
Se puede escribir $$709.3045 = 7 \cdot {10^2} + 9 \cdot {10^0} + 3 \cdot {10^{ - 1}} + 4 \cdot {10^{ - 3}} + 5 \cdot {10^{ - 4}}$$
En general, si consideramos otra base $b$, un número real $x$ se puede escribir de la forma
$$x = \sum\limits_{ - \infty }^n {{a_k} \cdot {b^k}} \,\,\,\,\, 0 \le {a_k} \le b - 1$$ o también $$x = {a_n}{a_{n - 1}}...{a_0}.{a_{ - 1}}{a_{ - 2}}{a_{ - 3}}...$$El Standard de IEEE propone el uso de la base 2 en el almacenamiento y la aritmética sobre el computador. En este sistema de numeración los números se representan utilizando solamente las cifras cero y uno (0 y 1)
$$x = \sum\limits_{ - \infty }^n {{a_k} \cdot {2^k}} \,\,\,\,\, 0 \le {a_k} \le 1$$El procedimiento es sencillo, basta dividir el número en base decimal y sus cocientes entre 2, acumulando los restos que conformarán el número en la nueva base. La lista de ceros y unos leídos de abajo a arriba darán la expresión en binario.
En la escena que se presenta a continuación, se muestra cómo hacer la conversión paso a paso para ciertos ejemplos. Este procedimiento es aplicable a cualquier base. En la escena se puede elegir números enteros y números con parte fraccionaria, observa el procedimiento para este último caso.
Ejemplo. Convierte $0.1_{10}$ a binario.
$$ \begin{aligned}{1 \over {10}} &= {1 \over {{2^4}}} + {1 \over {{2^5}}} + {0 \over {{2^6}}} + {0 \over {{2^7}}} + {1 \over {{2^8}}} + {1 \over {{2^9}}} + {0 \over {{2^6}}} + {0 \over {{2^7}}} + ... \\ &= 0001{\widehat {1001}_2} \end{aligned}$$
Puedes practicar con la siguiente herramienta interactiva introduciendo un número binario sin parte fraccionaria.
Ejemplo. Convierte el número $x=10.111001_2$ a decimal
Se tiene que: $$\begin{aligned} x = &1·2^1 + 0·2^0 + 1·2^{-1}+ 1·2^{-2} + \\ &+1·2^{-3}+ 0·2^{-4} + 0·2^{-5} + 1·2^{-6} \end{aligned}$$
Un número en base 10 puede escribirse como
$$x = {\left( { - 1} \right)^s}× \,{10^e} × \,a.f$$donde
Se denomina:
Ejemplo. Considera el número $x=9028.304$ en base decimal.
Cuando la base es 2, el número se puede escribir
$$x = {\left( { - 1} \right)^s} × \,{2^e} × \,1.f$$Observa que para almacenar $x$ basta almacenar el signo $s$, el exponente $e$ y la mantisa $f$. No es necesario almacenar el 1 (bit oculto).
En el ordenador los números se almacenan en "palabras" de 32 bits (precisión sencilla), 64 bits (doble precisión) o 128 bits (precisión cuádruple). En Matlab la precisión que usaremos para los cálculos será doble precisión. En este caso se utiliza;
Teniendo en cuenta que el exponente se almacena en 11 bits, el valor mayor que se puede almacenar es:
$$\underbrace {0....0}_{11 \,bits}\,\,\,\,\,\,\,\,\text{hasta}\,\,\,\,\,,\,\underbrace {1....1}_{11 \,bits} \equiv 1 + 2 + ... + {2^{10}} = {2^{11}} - 1 = 2047$$Para poder considerar las potencias de 2 negativas se considera la siguiente representación normalizada en 64 bits:
$$x = \left( { - 1} \right)^s × 2^{e - 1023}× \,\,1.f \,\,\,\,\,\, \,\,\,\, 0 < e < 2047 $$ se representa el exponente en exceso a 1023.Ejemplo. Almacenamiento del número $3.125_{10}$
Convirtiendo el exponente en binario, se tendrá $10000000000$. La mantisa ajustada a 52 bits es $1001\underbrace{0...0}_{48 \,\,bits}$.
El número se almacenerá de la forma: $$ \underbrace{0}_{signo} \underbrace{10000000000}_{exponente}\,\,\underbrace{10010...0}_{mantisa, \,48\,\,últimos \,\,bits\,\,0}$$
Herramientas de conversión
Dado que la representación de números reales bajo estos formatos es aproximada, hay dos conceptos importantes en la aritmética en punto flotante:
Mayor número normalizado: realmax en Matlab
$${2^{1023}} × 1.\underbrace {1....1}_{52 \,\,bits} \equiv {2^{1023}}\left( {1 + \sum\limits_{j = 1}^{52} {{1 \over {{2^j}}}} } \right) \approx 1.79 \,× {10^{308}}$$Menor número estrictamente positivo normalizado: realmin en Matlab
$$2^{ - 1022} × 1.00...\underbrace 0_{bit\,\,52} \equiv {2^{ - 1022}} \approx 2.22 \,× {10^{-308}} $$Si consideramos que el bit escondido es cero, podemos obtener números menores que el mínimo número normalizado mediante la siguiente representación: $$x = \left( { - 1} \right)^s × 2^{- 1022} × \,\,0.f$$
Estos números se dicen que están desnormalizados.
Será $${2^{ - 1022}} \cdot 0.\underbrace {00...01}_{52\,bits} = {2^{ - 1022}} \cdot {2^{ - 52}} = {2^{ - 1074}}$$
Todos los números menores que este valor se consideran 0.
Precisión de la máquina: eps en Matlab
El número más pequeño que se puede representar después del 1 utilizando representación binaria en punto flotante con doble precisión es: $$1 = \left( { + 1} \right)× {{\rm{2}}^0}× 1.\underbrace {00...0}_{52}$$ $$1 + \varepsilon = \left( { + 1} \right)× {{\rm{2}}^0}× 1.\underbrace {00...01}_{52}$$
El valor de $\varepsilon$ es la precisión de la máquina:
$$\varepsilon = {2^{ - 52}}$$El número más pequeño $\varepsilon_M$ de forma que $1-\varepsilon_M \ne 1$: $$\varepsilon_M = {2^{ - 53}}$$
Para poder aplicar las operaciones matemáticas básicas a números representados por coma flotante se deben adecuar los números previamente.
Suma y resta: se alinean los bits (se aumenta la mantisa del número de menor exponente) y se suman o restan.
$$123456.7+101.7654$$ $$(1.234567×10^5)+(1.017654×10^2)=$$ $$(1.234567×10^5)+(0.001017654×10^5)=$$ $$(1.234567+0.001017654)×10^5=1.235584654×10^5$$
Multiplicación y división: se multiplican o dividen las mantisas y se suman o restan los exponentes.
El uso del computador para representar valores numéricos y operar con ellos, presenta ciertas limitaciones que no deben olvidarse y que por su importancia pasamos a tratar brevemente.
Si $\widetilde p$ es una aproximación de $p$, se define:
En general, en los métodos numéricos no se conoce el valor verdadero, por lo tanto tampoco el error real. Lo frecuente es tener una estimación del error dada por una cota positiva al tamaño máximo del error.
Cuestión. En 1862 el físico Foucault, utilizando un espejo giratorio, calculó en 298 000 km/s la velocidad de la luz. Aceptando como exacta la velocidad de 299 776 km/s, estimar el error absoluto y el error relativo cometido por Foucault. Por otra parte, la determinación de la constante universal h realizada por Planck en 1913 dio el valor de 6.41×10−27 erg seg. Adoptando el valor de 6.623×10−27 erg seg, estimar el error absoluto y relativo cometido por Planck. ¿Cuál de las dos medidas es más precisa?
Como se ha indicado anteriormente, las limitaciones de almacenamiento de bits generan cálculos aproximados de los números verdaderos, es decir, los números y errores están sujetos al sistema numérico de punto flotante. Si se considera el sistema decimal, un número real puede ser normalizado de la forma:
$$x = a.f_1f_2f_3...f_kf_{k+1}... × 10^e$$ $$1 \le f_i \le 9 \,\,\,\, \text{con}\,\,\,i = 1,2,3,4,...k$$El numero almacenado por el ordenador, que denotaremos por $fl(x)$ (de "floating") se obtiene terminando (recortando) la mantisa de $y$ en $k$ dígitos decimales. Existen dos métodos:
Veamos algunos problemas que se producen al redondear los números almacenados por el computador.
Ejemplo. Calcula con Matlab la siguiente operación $1 - 3*(4/3 - 1)$.
e =
2.2204e-16
Un caso frecuente de cancelación surge cuando se calcula la diferencia de dos números bastante próximos, y el redondeo del ordenador iguala esa diferencia a 0.
Ejemplo. Dados $x=1+2^{-54}$ e $y=1$, calcula con Matlab $x-y$
Nuestro ordenador establecerá que $x-y=0$ cuando en realidad $x - y \approx 5,5 \cdot {10^{ - 17}}$.
>>x=1+2^(-54);y=1;x-y
ans=
0
Ejemplo. Dados $x=2^{60}+2^6$ e $y=2^{60}$, calcula con Matlab $x-y$.
Matlab establecerá que $x-y=0$ cuando en realidad $x - y =64$.
>>x=2^60+2^6;y=2^60;x-y
ans=
0
Ejercicio 5 propuesto del tema 1.
f=@(x) (exp(x)-1).^2+(1./sqrt(1+x.^2)-1).^2;
k=-1:-1:-17;
h=10.^(k);
v=(f(1+h)-f(1))./h;
format long
disp([h' v'])
%h a partir del valor k=-16 es cero
%La función para h pequeño tiende a cero y
%se produce también cancelación en el segundo
%sumando.
Ejemplo. Se sabe que ${e^x} = \sum\limits_{n = 0}^\infty {{{{x^n}} \over {n!}}} $. Calcula una aproximación de $e^{-40}$ con los 100 primeros términos de la serie.
Si $x< 0$ los términos de la serie cambian de signo y si $|x|$ es grande, el valor de los sumandos puede devolver un valor pequeño cuando los sumandos sean muy grandes. En algún momento estaremos con un error de cancelación cuando sumemos "el siguiente" término. la manera de arreglar el problema es usar el hecho de que
${e^{ - x}} = {1 \over {{e^x}}}$. Calcula $e^{-40}$ teniendo esto en cuenta y utilizando la aproximación que darán los 100 primeros sumandos.
x=40; %Inverso del número a calcular
suma=1;
p=1;
for k=1:99
p=p*(x/k);
suma=suma+p;
end
1/suma
x=-40; %Sin considerar el inverso
suma=1;
p=1;
for k=1:99
p=p*(x/k);
suma=suma+p;
end
suma
Ejercicio 7 propuesto del tema 1.
Apartado a
f=@(x) 1./(1+2*x)-(1-x)./(1+x);
%Se realiza la suma para evitar la cancelación
g=@(x) 2*x.^2./((1+2*x).*(1+x));
k=-1:-1:-17;
h=10.^k;
f1=f(h);
g1=g(h);
format long
disp([h' f1' g1'])
Apartado b
f=@(x) (1-cos(x))./x;
%Se multiplica y divide por (1+cos(x))
g=@(x) sin(x).^2./(x.*(1+cos(x)));
k=-1:-1:-17;
h=10.^k;
f1=f(h);
g1=g(h);
format long
disp([h' f1' g1'])
Apartado c
f=@(x) sqrt(x+1./x)-sqrt(x-1./x);
%Se multiplica y divide por el conjugado
g=@(x) 2*x./((sqrt(x+1./x)+sqrt(x-1./x)));
k=-1:-1:-25;
h=10.^k;
f1=f(1+h);
g1=g(1+h);
format long
disp([h' f1' g1'])
Con frecuencia una operación aritmética con dos números válidos da como resultado un número tan grande o tan pequeño que la computadora no puede manejarlo; como consecuencia se produce un overflow o un underflow, respectivamente.
Para este análisis tengamos en cuenta que si
$$x = {\left( { - 1} \right)^s}× {2^{{e_x}}} × 1.f\,\,\,\,\,\,\,\,\,\,\,\,\,\,0 < {e_x} < 1024$$ de acuerdo a las reglas de redondeo se tiene que $$\left| {x - fl\left( x \right)} \right| \le {2^{{e_x}}} \cdot {2^{ - 53}} \Rightarrow {{\left| {x - fl\left( x \right)} \right|} \over {\left| x \right|}} \le {2^{ - 53}}$$Al valor $2^{-53}$ se le llama épsilon de la máquina, que denotaremos por ${\varepsilon _M}$ y nos da una cota superior del error relativo debido a la expresión en formato punto flotante de un número.
Si $\otimes$ es una operáción aritmética el cálculo sigue el siguiente proceso:
$$fl\left( {x \otimes y} \right) = fl\left( {fl\left( x \right) \otimes fl\left( y \right)} \right)$$Errores de redondeo de la suma
Vamos a analizar el error originado a la hora de realizar la operación en punto flotante y representar los números en punto flotante. Denotemos por ${\varepsilon _x}$, el error relativo entre $x$ y su representación en coma flotante $fl(x)$ $${\varepsilon _x} = {{x - fl\left( x \right)} \over x} \,\,\,\,\,\,\,\left| {{\varepsilon _x}} \right| \le {\varepsilon _M}$$ $$fl\left( x \right) = x\left( {1 + {\varepsilon _x}} \right)\,\,\,\,\,\,\,\left| {{\varepsilon _x}} \right| \le {\varepsilon _M}$$
Dados $x$, $y$ se tiene que: $$fl\left( x \right) = x\left( {1 + {\varepsilon _x}} \right)\,\,\,\,\,\,\,\left| {{\varepsilon _x}} \right| \le {\varepsilon _M}$$ $$fl\left( y \right) = y\left( {1 + {\varepsilon _y}} \right)\,\,\,\,\,\,\,\left| {{\varepsilon _y}} \right| \le {\varepsilon _M}$$
y si $z=x+y$ entonces $$fl\left( z \right) = \left( {x\left( {1 + {\varepsilon _x}} \right) + y\left( {1 + {\varepsilon _y}} \right)} \right)\left( {1 + {\varepsilon _z}} \right)$$
Se puede probar que el error relativo cometido al realizar la suma está acotado por: $$E \le {{\left| x \right| + \left| y \right|} \over {\left| {x + y} \right|}}\left( {2 + {\varepsilon _M}} \right){\varepsilon _M}$$
Observamos que el error originado por la operación en sí y por la expresión inexacta de los sumandos será pequeño en cuanto $x+y$ no sea pequeño. Sin embargo, si los números a sumar son aproximadamente iguales en magnitud, pero de signos opuestos, el error relativo es muy grande. Esta situación, conocida como cancelación de dígitos significativos, debe ser evitada tanto como sea posible.
Veamos cómo obtener la expresión anterior.
Dado que $$fl\left( x+y \right) = \left( {x + x{\varepsilon _x} + y + y{\varepsilon _y}} \right)\left( {1 + {\varepsilon _z}} \right)$$ el error relativo $$E = \left| {{{z - fl\left( z \right)} \over z}} \right| = {{\left| {x + y - \left( {x + x{\varepsilon _x} + y + y{\varepsilon _y}} \right)\left( {1 + {\varepsilon _z}} \right)} \right|} \over {\left| {x + y} \right|}} $$Errores de redondeo de la multiplicación y de la división
Análogamente puede demostrarse las acotaciones de errores relativos para las siguientes operaciones:
Octave es una alternativa a Matlab con un aspecto y unos comandos y sintasis similares. Se trata de un lenguaje interpretado de alto nivel. Lo puedes descargar desde su página oficial completamente gratis. Aquí tienes el enlace de descarga.
Pero además, puedes utilizar una versión en línea desde tu móvil en el caso de que no tengas a mano tu ordenador. Esta versión la tienes aquí: Octave online.