Aritmetica Finita

Capítulo 2

Aritmética de precisión finita Considere la siguiente operación aritmética: (1 + 10−15

− 1) × 10

30

,

(2.1)

la cual si se evalúa en una calculadora científica se obtiene el valor de cero, mientras que si se usa MATLAB se obtiene el valor 1 valor 1..110223024625157 1015. Sin embargo, nótese que (1 + 10−15 1) 1030 = 10−15 1030 = 1015.

×

− ×

×

Así, tanto la calculadora como MATLAB dan valores erróneos al cálculo (2.1), lo cual se debe a que los cálculos aritméticos realizados por una herramienta computacional, utilizan un número finito de dígitos en sus valores, ya que si se trabajara con una cantidad infinita de dígitos (suponiendo es posible almacenar dicha cantidad), nunca se terminaría de considerar todos los dígitos de dicho número, pues son infinitos. Sin embargo, los número reales, salvo pocas excepciones, cuentan con una cantidad infinita de dígitos, por lo que es importante designar alguna manera para considerar dichos valores. La IEEE (“Institute of Electrical and Electronics Engineer” y se lee i triple e ) define define un formato formato estándar estándar (ANSI/IEEE (ANSI/IEEE standard standard 754-1985, 754-1985, [38]) para representar los números reales en general al considerar una cantidad finita de sus dígitos en alguna base, usualmente base 2. El problema al representar un número con infinitos dígitos por medio de otro número con finitos dígitos, o bien un número con finitos dígitos por un número con una menor cantidad de sus dígitos, es que, por lo general, se introduce un error en la representación de dicho número, el cual es relevante poder medir. Por lo todo lo anterior, en este capítulo se describe la representación de un número en la computadora, denominada representación punto flotante , así como diferentes maneras de medir los errores producidos por dicha representación. Se presenta además algunos ejemplos típicos que producen errores y algunos de ellos con una solución alternativa. Sin embargo, primero se presenta una breve reseña sobre el sistema de numeración binaria, con el objetivo conocer cómo se realizan las operaciones aritméticas en una computadora, de manera más interna.

2.1. 2.1.

Sist Si stem emas as de num umer erac ació ión n

Suponga que al dividir un número entero a por un número entero b, se obtiene un residuo r1 y un cociente q 1 b. Luego, si se divide q 1 por b se obtiene otro residuo r2 , además de otro

≥

105

106

Capítulo 2. Aritmética de precisión finita

cociente q cociente q 2 b. Así, siguiendo siguiendo este proceso proceso hasta que obtener q obtener q n < b, se cumple que a que a se se puede representar de la forma:

≥

+ r2 b + r + r1 , a = q n bn + rn bn−1 + rn−1 bn−2 + . . . + r

(2.2)

y con ello, se denota a en base b de la siguiente manera: a = (a)10 = (q n rn rn−1 . . . r2 r1 )b , 10 (bases mayores que 10), se utilizan donde ri < b, para i = 1, 2, . . . , n. n. En caso de que b > 10 (bases las letras A para el dígito 10, B para el 11 y así sucesivamente. Por ejemplo, para el sistema 10) y hexadecimal (b 16), los dígios se muestran en la tabla 2.1. binario (b (b = 2), decimal (b (b = 10) (b = 16), Sistema Base Dígitos 0, 1 Binario 2 0, 1, . . . , 9 Decimal 10 0, 1, . . . , 9, A , B , . . . , F Hex Hexade adecimal imal 16 Tabla 2.1: Sistemas de numeración más conocidos. Ejemplo 2.1. Considere el número 1511, al cual se desea representar en base 12. Al realizar la división de 1511 entre 12, se obtiene que: 1511 1511 = 125 125 12 + 11, 11,

·

125 y r 1 = 11 = B = B.. Luego, como 125 < 12, 12, se sigue dividiendo de donde q donde q 1 = 125 y



125 = 10 12 + 5, 5,

·

con lo que q 2 = 10 < 12 y r2 = 5. Por lo tanto, se detiene el procedimiento y como = A,, se tiene que q 2 = 10 = A 1511 = (A (A5B )12 . 1511.. Nótese que ( que (A A5B )12 = 10 122 + 5 121 + 11 120 = 1511

·

·

·

Observación: El Observación: El acomodo usual de la división aritmética aritmética brinda gran ayuda cuando se buscan buscan los dígitos que conforman un número en alguna base, como se muestra a continuación: 151 1 150 0

−

11

12 125 120

−

12 10

5

(A5B )12 . A Luego, tomando los número en negrita de derecha a izquierda, se tiene que 1511 = (A base . esta técnica se le conoce como división por base.

106


cociente q cociente q 2 b. Así, siguiendo siguiendo este proceso proceso hasta que obtener q obtener q n < b, se cumple que a que a se se puede representar de la forma:

≥

+ r2 b + r + r1 , a = q n bn + rn bn−1 + rn−1 bn−2 + . . . + r

(2.2)

y con ello, se denota a en base b de la siguiente manera: a = (a)10 = (q n rn rn−1 . . . r2 r1 )b , 10 (bases mayores que 10), se utilizan donde ri < b, para i = 1, 2, . . . , n. n. En caso de que b > 10 (bases las letras A para el dígito 10, B para el 11 y así sucesivamente. Por ejemplo, para el sistema 10) y hexadecimal (b 16), los dígios se muestran en la tabla 2.1. binario (b (b = 2), decimal (b (b = 10) (b = 16), Sistema Base Dígitos 0, 1 Binario 2 0, 1, . . . , 9 Decimal 10 0, 1, . . . , 9, A , B , . . . , F Hex Hexade adecimal imal 16 Tabla 2.1: Sistemas de numeración más conocidos. Ejemplo 2.1. Considere el número 1511, al cual se desea representar en base 12. Al realizar la división de 1511 entre 12, se obtiene que: 1511 1511 = 125 125 12 + 11, 11,

·

125 y r 1 = 11 = B = B.. Luego, como 125 < 12, 12, se sigue dividiendo de donde q donde q 1 = 125 y



125 = 10 12 + 5, 5,

·

con lo que q 2 = 10 < 12 y r2 = 5. Por lo tanto, se detiene el procedimiento y como = A,, se tiene que q 2 = 10 = A 1511 = (A (A5B )12 . 1511.. Nótese que ( que (A A5B )12 = 10 122 + 5 121 + 11 120 = 1511

·

·

·

Observación: El Observación: El acomodo usual de la división aritmética aritmética brinda gran ayuda cuando se buscan buscan los dígitos que conforman un número en alguna base, como se muestra a continuación: 151 1 150 0

−

11

12 125 120

−

12 10

5

(A5B )12 . A Luego, tomando los número en negrita de derecha a izquierda, se tiene que 1511 = (A base . esta técnica se le conoce como división por base.

107

2.1. Sistemas de numeración

Ejemplo 2.2. Muestre que la representación del número 2013 en base 5, corresponde a (31023)5. Solución. Realizando Solución. Realizando las divisiones, se tiene que 2013 2010

−

5 40 2 400

−

3

−

2

5 80 80 0

5 16 15

−

5 3

1

Por lo tanto, 2013 tanto, 2013 = (31023)5 .



En los ejemplos anteriores se ilustra que para convertir un número del sistema decimal a cualquier otra base, basta con realizar la división por base . Mientras que si se desea convertir un número en base b base b a base diez, se puede hacer usando la expresión (2.2), tal y como se muestra en el siguiente ejemplo. Ejemplo 2.3. Los siguientes son ejemplos de cambios a base diez: El número binario (11101) binario (11101)2 corresponde al número 29, pues (11101)2 = 1 24 + 1 23 + 1 22 + 0 21 + 1 20 = 29. 29.

·

·

·

·

·

El número cuaternario (3012) cuaternario (3012)4 corresponde al número 198, pues (3012)4 = 3 43 + 0 42 + 1 41 + 2 40 = 198. 198.

·

·

·

·

El número octal (77) octal (77)8 corresponde al número 63, pues (77)8 = 7 81 + 7 80 = 63. 63.

·

2.1. 2.1.1. 1.

·

Sist Sistem ema a bina binari rio o

El sistema de numeración binaria o sistema binario, consiste en representar los números utilizando solamente las cifras cero y uno (0 y 1). Es el que se utiliza en las computadoras, debido a que trabajan internamente con dos niveles de voltaje (encendido: 1, apagado: 0). En la sección de ejercicios resueltos del capítulo de MATLAB (página 87), se presenta en el problema 1 un fragmento de código que permite representar cualquier número entero en base 10 a la base dos. Con respecto a la aritmética binaria , es importante conocer que la ejecución de cálculos numéricos es esencialmente igual en todos los sistemas de numeración posicional (el (el valor de un dígito

108


depende no sólo de su valor absoluto sino que también de su posición dentro del número). La aritmética binaria se lleva a cabo exactamente igual que la aritmética decimal, excepto que las reglas son más sencillas debido a que son mucho menos las combinaciones posibles de dígitos, lo que la hace ideal para cualquier computadora. Suma y resta binaria A continuación se resumen las cuatro reglas básicas para sumar dígitos en el sistema binario: 0+0 0+1 1+0 1+1

= = = =

0 1 1 10

Suma = Suma = Suma = Suma =

0 1 1 0

Acarreo = Acarreo = Acarreo = Acarreo =

0 0 0 1

Los dígitos que se llevan (acarreo) se manipulan de la misma forma que en la aritmética decimal. Dado que 1 es el dígito más grande en el sistema binario, cualquier suma mayor a 1 requiere que se traslade un dígito a la siguiente posición. Por ejemplo, la suma de (100)2 más (100)2 requiere la suma de dos unos en la tercera posición hacia la izquierda, por lo que un uno se traslada a la siguiente posición. Así, como 1 + 1 = 0 más un 1 que se traslada, la suma de (100)2 y (100)2 es (1000)2 . Ejemplo 2.4. Realice las siguientes sumas 1. (10110111)2 + (110011101)2 2. (1001)2 + (1101)2 + (110)2 + (1011)2 Solución. 10

01 1 01 11 11 01 11 11 1 + 1 1 0 0 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0

+ 1 0

11 1 0 1 0

01 1 1 0 1

01 0 1 1 1

1 1 0 1 1



En el caso de la resta, es necesario establecer un procedimiento para sustraer un dígito grande de uno menor. El único caso de este tipo, es cuando se resta 1 a 0. El residuo es 1 pero es necesario pedir prestado 1 de la siguiente posición a la izquierda, como se muestra a continuación: 0 1 1 0

−0 −0 −1 −1

= = = =

0 1 0 1 (siempre que se pida prestado 1 para obtener 10 − 1 = 1)

109

2.1. Sistemas de numeración

Ejemplo 2.5. Realice las siguientes restas 1. (1001)2

− (101) 2. (10000) − (11) 2

2

2

Solución.

−

10 010 0 1 0 1 0 1 0 1 0 0

−

10 01 01 01 010 0 0 0 1 1 1 1 0 1 

Multiplicación y división binaria La multiplicación en el sistema binario se realiza de la misma forma que la decimal, salvo que es más sencilla, debido a que los únicos productos parciales son: 0 1 0 1

×0 ×0 ×1 ×1

= = = =

0 0 0 1

Ejemplo 2.6. Calcular el valor (1100)2

× (1010) . 2

Solución.

× 1 0 0 + 1 1 0 1 1 1

1 1 0 1 0 0 1

1 0 0 0 0

− 0

0 1 0 0

− 0

0 0 0

− 0



La división binaria también resulta muy sencilla teniendo en consideración que la división por cero no tiene sentido y con ello sólo hay dos posibles divisiones: 0 1

÷1 ÷1

= 0 = 1

110


Ejemplo 2.7. Realizar las divisiones 1. (1010001)2 2. (1110111)2

÷ (11) ÷ (1001) 2

2

Solución.

1

−

0 1 1

−

1 0 1 0 0 1 1 1

−

0

0 1

0 1

1 1 1 1 0 1 1

−

−

0 1 1 1 1 1 0

−

1 1

1 0 1 1

1 0 0 0

−

0 1 1 0 1 1

1 1 1 1 0 0 1 1 1 0 1 1 1 0 1 1 0 0 1 1 0 

Ejercicios (sección 2.1) ♦ Tradicionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. Expresar los siguientes números en base decimal. a ) (10011101)2

c ) (555)7

e ) (ABC )13

b ) (31322)4

d ) (1776)8

f ) (ABC )16

2. Dados los siguientes números en base decimal, encuentre su representación en la base indicada. a ) 199 en base 2

c ) 356281 en base 5

e ) 1000 en base 8

b ) 163 en base 3

d ) 328204 en base 7

f ) 10995 en base 16

3. Dado que el número 177 se representa en base b como (2301)b , halle el valor de b. 4. Realice las siguientes operaciones en binario a ) (111111)2 + (11101111)2 b ) (11001)2 + (11100)2 + (110011)2 c ) (1100011)2 d ) (101101)2

− (1101111) − (1111)

2

2

e ) (111011)2

× (1011) f ) (11100111) × (11) g ) (111001) ÷ (1001) h ) (1011011) ÷ (111) 2

2

2

2 2 2 2

111

2.2. Números reales de precisión finita

2.2. Sea x para i

Números reales de precisión finita

∈ −{0} y β ∈ ≤ t, tales que:

≥ 2, entonces se pueden encontrar dígitos a ∈ {0, 1, . . . , β − 1},

con β

i

t

s

x = ( 1) ( at at−1 . . . a1 a0 . a−1 a−2 a−3 . . . )β = ( 1)

−

−

      parte entera

parte fraccionaria



s

i=

ai β i ,

−∞

0, 1 , según el signo de x. Ahora, suponiendo que se tienen únicamente N posiciones para s de memoria para almacenar x, entonces una de esas posiciones se utilizará para almacenar el signo de x, es decir, el valor de s. Con las restantes N 1 posiciones, suponiendo que k es la cantidad de dígitos en la parte fraccionaria, entonces se tienen N k 1 (t = N k 2) dígitos en la parte entera y haciendo el cambio de variable d i = a i−k−1 se obtiene que:

∈{ }

−

− −

− −

N 1

x

s

≈ (−1) (d

− d − . . . dk+1 .dk dk−1 . . . d1 )β = ( 1)

−

N 1 N 2 N 1

= ( 1)s β −k

−

−

 i=1

s

−



diβ i−k−1 ,

i=1

di β i−1 = ( 1)s (dN −1 . . . d2 d1 )β β −k .

−

   · m

Observación: La parte entera de x se representa de la forma x , mientras que la parte fraccionaria se puede considerar de la forma x x.

−





Definición 2.1.

≥ 2, entonces la representación punto

Sea x un número de t dígitos y β un entero, β flotante de x se denota f l(x) y se define como

f l(x) = ( 1)s m

−

e

× β ,

donde s es el bit de signo, e es llamado exponente y m se denomina mantisa, la cual 0, 1, . . . , β 1 . Además, si se escribe de la forma: consta de t dígitos d i

∈ {

− }

f l(x) = ( 1)s 0.m

−

e+t

× β

,

se dice que está en su forma normalizada y esta representación es única para x.

Ejemplo 2.8. Considere el número x = 13 , el cual puede ser expresado con 4 dígitos en base 10, normalizado, de la forma f l(x) = 0.3333 100 , o bien no normalizado, de las formas: f l(x) = 3.333 10−1 o f l(x) = 33.33 10−2 .

×

×

×

112


Observación: Nótese que la notación punto flotante se parece a la notación científica, la cual consiste de una manera rápida de representar un número utilizando potencias de base diez. Pero en la notación punto flotante la base no tiene porque ser 10, además que generalmente los valores son normalizados. Ejemplo 2.9. Considere x =

−4.75, represente este valor en base 2 con todos los dígitos posibles.

Solución. Nótese que 4 = (100)2 , mientras que: 0.75

× 2 = 1.5, Genera 1 y sigue con los restantes, 0.5 × 2 = 1.0, Genera 1 y termina, con lo que 0.75 = (0.11) . Por lo tanto x = −(100.11) y 2

2

f l(x) =

2

−10011 × 2− .



Observación: Del ejemplo anterior, se puede apreciar que también es posible considerar f l(x) = 0.10011 23 en forma normalizada. Es decir, que se pueden considerar diferentes representaciones con la misma cantidad de dígitos al cambiar el exponente, sin embargo, lo usual es considerar el exponente en un intervalo finito de valores permitidos: L e U (generalmente L < 0 y U > 0). De esta forma, las N posiciones de memoria serán distribuidas de la siguiente manera: una posición para el signo, t posiciones para la mantisa y las restantes N t 1 posiciones son las cifras del exponente (entre L y U ).

−

×

≤ ≤

− −

Definición 2.2. Sean β, t con β 2 y L, U precisión finita se define por:

∈

≥

(β,t,L,U ) =

∈

{0} ∪

con L



x

≤ U . El conjunto de números reales de t

s

∈{

e

∈ − {0} / x = (−1) β

0.1 , di 0, 1, . . . , β 1 con d1 = 0 y e donde s número t se le conoce como la precisión.

∈ { }

 

− }



di β −i ,

i=1

∈ {L, L + 1, . . . , U − 1, U }. Al

Con respecto a las N posiciones de memoria, la IEEE admite dos valores N = 32 denominada precisión simple o N = 64 denominada precisión doble . En la tabla 2.2 se muestra la división de estas posiciones. Precisión simple Signo Mantisa Exponente 1 23 8

Precisión doble Signo Mantisa Exponente 1 52 11

Tabla 2.2: Estándar IEEE para punto flotantes.

113


Nótese que al ser base β = 2, entonces la precisión simple provee posiciones para aproxima1.2 10−7 , mientras que la precisión doble damente siete dígitos decimales, ya que 2−23 2.2 10−16 . La precisión ofrece aproximadamente dieciséis dígitos decimales, ya que 2−52 simple se presenta en la mayoría de calculadoras científicas, mientras que la doble esn la de las computadoras.

≈

×

≈

×

Ejemplo 2.10. Considere x = 1.7 y encuentre f l(x)

∈

(2, 8, 0, 1).

Solución. Se tiene que para la parte entera 1 = (1)2 , mientras que para la parte fraccionaria:

0.7 0.4 0.8

×2 ×2 ×2

= 1.4 = 0.8 = 1.6

→ → →

0.6 0.2 0.4

1 0 1

×2 ×2 ×2

= 1.2 = 0.4 = 0.8

→ → →

1 0 0

Por lo tanto, 0.7 = (0.10110)2 y como se requieren 8 dígitos de precisión, entonces se tiene que  f l(x) = 1.1011001 20 = 0.11011001 21 .

×

×

Ejemplo 2.11. Indique todos los números del conjunto

Solución. Los elementos de 2

±0.111 × 2 ±0.111 × 2 ±0.111 × 2 ±0.111 × 2− 1

0

1

7 2 7 = 4 7 = 8 7 = 16 =

(2, 3, 1, 2).

−

(2, 3, 1, 2) son 0 y los siguientes valores:

−

2

=3

1

=

±0.110 × 2 ±0.110 × 2 ±0.110 × 2 ±0.110 × 2− 0

1

3 2 3 = 4 3 = 8

2

±0.101 × 2 ±0.101 × 2 ±0.101 × 2 ±0.101 × 2− 1

0

1

5 2 5 = 4 5 = 8 5 = 16 =

2

=2

1

=1

0

=

±0.100 × 2 ±0.100 × 2 ±0.100 × 2 ±0.100 × 2−

1

1 2 1 = 4

Es decir, card( ) = 33. Además, en la figura 2.1 se pueden apreciar estos valores en la recta numérica, donde es claro que hay mayor acumulación de los mismos cerca del cero. 

Figura 2.1: Distribución de los elementos de

(2, 3, 1, 2).

−

114


Teorema 2.1. Para todo x

∈

(β,t,L,U ), con x = 0, se cumple que



0 < xm´ın

≤ |x| ≤

xmáx ,

donde: xm´ın = β L−1

− 

xmáx = β U 1

y

β −t .

Demostración. Nótese que el menor valor se puede obtener al tomar e = L, d1 = 1 y di = 0 para i = 2, 3, . . . , t, es decir, xm´ın = β −1 β L = β L−1 .

·

Luego, el mayor valor ocurre cuando e = U y d i = β t

U

xmáx = β



(β

= β

·

β

t 1

t

i=1

U

− 1 para i = 1, 2, . . . , t, con lo que:

−

1)β −i = β U (β t 1

−      1

β

−

1 β

i=0

− 1)

i

U

= β



β −i = β U (β

i=1

− 1)

t

− 1) · −β 1 · β (β (β − 1)β β

t

suma geométrica

  −  −

i=0

1 β

= β U 1

i+1

,

β −t . 

Ejemplo 2.12. De los resultados del ejemplo 2.11 se puede apreciar que: 1 xm´ın = β L−1 = 2−1−1 = , 4

−  − 

xmáx = β U 1

β −t

= 22 1

2−3

=

7 . 2

Teorema 2.2. Considere (β,t,L,U ) para β, t con β 2 y L, U con L U , entonces se cumple que: card ( ) = 2(β 1)β t−1 (U L + 1) + 1.

∈

≥

−

∈

≤

−

Demostración. (ejercicio) Observación: De los resultados obtenidos en el ejemplo 2.10, se puede apreciar que: 0.7

≈

(1.1011001)2 = 20 + 2−1 + 2−3 + 2−4 + 2−7 = 1.6953125.

Así, surge la interrogante de cómo medir el error en la aproximación de f l(x) y por tal razón a continuación se consideran los estimados usuales para medir el error.

115


Definición 2.3 (Errores). Existen dos maneras de medir el error al representar x por f l(x), los cuales son: Error absoluto: x

| − f l(x)|, |x − f l(x)| , para x = 0. Error relativo: |x| El error relativo tiene más sentido que el error absoluto a la hora de comparar representaciones. Debido a que este considera que tan lejos están los valores de cero, ya que si estan lejos de cero, hay menos posibles valores para usar como representaciones. Esto se ilustra en el siguiente ejemplo, donde el valor x 2 esta más cerca de cero que x 1 , por lo que se presenta una diferencia en los errores relativos, mientras que no en los errores absolutos. Ejemplo 2.13 (Error relativo contra error absoluto). Considere x1 = 1.31 con f l(x1 ) = 1.30 y x2 = 0.12 con f l(x2 ) = 0.11. Los errores absolutos en ambos casos son iguales:

|x − f l(x )| 1

1

= x2

| − f l(x )| 2

= 0.01.

Por otro lado, el error relativo en el primer caso es:

|x − f l(x )| |x | 1

1

= 0.0076335,

1

mientras que para el segundo valor corresponde a:

|x − f l(x )| |x | 2

2

= 0.08333333.

2

Así, los errores relativos muestran que f l(x1 ) está más cerca de x1 que f l(x2 ) de x2 , mientras que los errores absolutos dados no indican esto. Definición 2.4. Se dice que f l(x) aproxima a x con q dígitos significativos, si q es el mayor entero no negativo para el cual el error relativo cumple que:

|x − f l(x)| |x|

< 5

q

× 10− .

Además, nótese que se tiene que: q =

−  | − |   log

x

f l(x) x

||

+

1 . 2

116


Ejemplo 2.14. Según el ejemplo 2.13, se tiene que f l(x1 ) aproxima a x1 con dos dígitos significativos, mientras que f l(x2 ) aproxima a x 2 con un único dígito significativo. Observación: El término dígito significativo suele confundirse con el de dígito correcto, sin embargo, no son lo mismo ya que se dice que una aproximación tiene p dígitos correctos , si el error relativo es menor a 12 10− p , es decir, si tiene p + 1 dígitos significativos.

×


1. Represente los siguientes números en su notación punto flotante normalizada, aplicando un corte a los tres dígitos decimales. a )

0.6523689

c )

b)

−52324.2365

d )

0.0000020035

−1.234000569

2. Represente los siguientes números en su notación punto flotante correspondiente, para que cumplan estar en el conjunto (2, 10, 1, 6).

−

a )

27.25

c )

2.625

b)

−0.40625

d )

−49.8125

3. ¿Cuál es el error relativo y absoluto al aproximar: a ) π por

22 ? 7

b)

1 por 0.333? 3

c )

1 por 0.166? 6

¿Cuántos dígitos significativos hay en cada representación? ♦ De

programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. Utilizando los dos tipos de errores (absoluto y relativo), analice la calidad de la aproximación de Stirling

√



n n S n = 2πn , e al factorial de un número n. Para ello, elabore una tabla de errores (puede ser con ayuda de MATLAB) para diferentes valores de n (por ejemplo n = 1, . . . , 15). ¿Qué se puede concluir acerca del uso de estos dos tipos de errores? 5. ¿Cuál será el resultado al ejecutar el siguiente fragmento de código en MATLAB? x = 2^(30); x + 2^(-23) == x x - 2^(-23) == x

Comente los resultados.

117

2.3. Truncamiento y redondeo

6. Utilice MATLAB para aproximar el valor de e −5 utilizando una serie de Taylor de orden 9 de las siguientes dos formas distintas: 9

e−5

( 5)n n!

− ≈ n=0

e−5 =

y

1 e5

≈

1 9



n=0

.

(5)n n!

¿Cuál es la mejor aproximación? (para ello puede usar la instrucción exp(-5) de MATLAB). ¿Qué se puede concluir de las dos estrategias utilizadas para aproximar e −5 ?

2.3.

Truncamiento y redondeo

Es posible que durante los cálculos entre valores punto flotantes se obtenga un nuevo valor que no es representable en el conjunto de números reales de precisión finita. Por ejemplo, considere (10, 3, 2, 3) definidos como: a, b

∈

−

3

a = 0.611

3

× 10 y b = 0.525 × 10 , ∈ (10, 3, −2, 3) ya que consta de 4 y no de 3 dígitos y nótese que c = a + b = 1.136 × 10  3

en la mantisa. Por tal razón, uno de los dígitos de c, el último, debe ser removido para que siga perteneciendo al conjunto de precisión finita. Para remover este último dígito se pueden considerar dos métodos de aproximación: el truncamiento y el redondeo. Para definirlos, considere x = ( 1)s 0.d1 . . . dtdt+1 . . . ,

−

con lo que, el método de truncamiento consiste en eliminar todos los dígitos desde dt+1 incluyendo el mismo, mientras que el método de redondeo consiste no sólo en eliminar todos los dígitos desde dt+1 , sino que el dígito d t puede ser aumentado en uno, dependiendo del valor de dt+1 . Así, se tiene que: Por truncamiento: f l(x) = ( 1)s 0.d1 d2 . . . dt−1 dt Por redondeo: f l(x) =

Ejemplo 2.15. Considere

 

− (−1) 0.d d . . . d − d (−1) (0.d d . . . d − d + β − ) s

1 2

s

t 1 t

1 2

t 1 t

t

si dt+1 < si dt+1

(10, 3, 2, 2) y el valor

−

√

x = 2 2 = 0.282842712474 . . . entonces, se tiene que: Por truncamiento: f l(x) = 0.282

× 10

Por redondeo: f l(x) = 0.283

1

× 10

1

× 10

1

≥

β 2 β 2

118


Ejemplo 2.16.

−∞, +∞) y x = 3.141596, entonces por redondeo se tiene que: Para t = 1: f l(x) = 0.3 × 10 Para t = 4: f l(x) = 0.3142 × 10 Para t = 2: f l(x) = 0.31 × 10 Para t = 5: f l(x) = 0.31416 × 10 Para t = 3: f l(x) = 0.314 × 10 Para t = 6: f l(x) = 0.314160 × 10

Considere

(10, t,

1

1

1

1

1

1

Debido a que se cuenta con dos métodos de aproximación, surge la interrogante de ¿qué es mejor usar, entre truncamiento y redondeo? . El siguiente resultado responde a esta pregunta, al presentar una cota para el error relativo en ambos métodos de aproximación, el cual muestra que las aproximaciones por redondeo son mejores que las aproximaciones por truncamiento. Teorema 2.3 (Cota superior para el error relativo). Sea f l(x) la representación punto flotante de un número real x, entonces

|x − f l(x)| ≤ |x|



µ =

1 1−t β 2

por redondeo,

β 1−t

por truncamiento.

(2.3)

Demostración. Considere el número x = ( 1)s (0.d1 d2 . . . dt . . .)

−

e

× β ,

1. En el caso de redondeo, puede pasar una de las siguientes dos posibilidades: f l(x) = x1 = ( 1)s (0.d1 d2 . . . dt )

− (−1)

f l(x) = x2 = Es claro que x x2 , con lo que:

s



e

× β ,

(0.d1 d2 . . . dt ) + β −t

×

β e .

∈ ]x , x [ y sin pérdida de generalidad suponga que x está más cerca de 1

2

|x − x | ≤ 12 |x − x | = 12 β − . Luego, como 0.1 × β = β − ≤ |x|, entonces el error relativo, viene dado por: |x − f l(x)| = |x − x | ≤ β − = 1 β − . |x| |x| 2 β − 2

e

1

e t

2

e 1

2

1 2

e t

1 t

e 1

2. En el caso de truncamiento, se tiene que f l(x) = ( 1)s (0.d1 d2 . . . dt )

−

e

× β ,

y con ello:

|x − f l(x)|

= (0.0 . . . 0dt+1 . . .)

e

× β

= (0.dt+1 . . .)

e t

× β − ≤

(β

e t

− 1)(0.11 . . .) × β − .

119

2.3. Truncamiento y redondeo

Luego, se tiene que (β

− 1)(0.11 . . .)

= (β

e t

con lo que x

| − f l(x)| ≤ β −

− 1)

+∞

  1 β

k =0

y como β e−1

k +1

=

β

−1 β

+∞

  k =0

1 β

k

= 1,

≤ |x|, entonces

|x − f l(x)| ≤ |x|

β e−t = β 1−t. − 1 e β



Observación: Debido a que el resultado anterior muestra que la cota de error en redondeo siempre es menor que la cota de error en truncamiento, entonces en caso de no referirse a un método de aproximación específico, se conciderará el de redondeo. Ejemplo 2.17. Considere x = 0.2346 en

(10, 3,

−∞, +∞), entonces:

Por redondeo: f l(x) = 0.235, con lo que el error relativo es 0.001705 y nótese que 0.001705 < 12 10−2 = 0.005.

·

Por truncamiento: f l(x) = 0.234 y el error relativo es 0.0025575 < 10 −2 = 0.01. Nótese además que en ambas representaciones se tienen tres dígitos significativos. Definición 2.5. El número µ en (2.3) se denomina precisión de máquina o épsilon de máquina (en este caso µ se representa como eps) y corresponde al menor punto flotante positivo tal que: f l(1 + µ) > 1. Ejemplo 2.18. Escriba en MATLAB una rutina para aproximar el valor numérico del número real positivo eps más grande tal que 1 + eps = 1. Solución. Considere la siguiente rutina: eps = 1.0; while (1.0 + eps) ~= 1.0 eps = eps / 2.0; end

disp([' eps = ' num2str(eps, 16)]); % Mostrar eps con 16 dígts disp([' 1 + e p s = ' num2str(1.0 + eps, 16)]); % Mostrar 1 + eps

120


Al ejecutarse produce el resultado: eps = 1.110223024625157e-016 1 + eps = 1




1. Para cada uno de los valores que se presentan a continuación, determine su representación en (10, 5, , + ) por redondeo y por truncamiento. Calcule además los errores relativos respectivos.

−∞ ∞

a )

32.419

b)

−19.0686

c ) 712342

× 10−

4

0.999999

e )

d ) 35364500001

f )

−0.000000001

2. Calcule los siguientes valores en la calculadora. Luego, represente cada resultado final en el sistema (10, 4, 0, 1). Utilice ambos métodos de aproximación y determine su respectivo error relativo. a ) π +

5 b) 12 c )

3.

√

2

d )

   e +

5

1 3 √ + e e 2

√

 √

3 11

1 π + 0.725 e ) + 3 11 1 f ) 1 + 1+ 1 1

0.6 + 1 ln 2

π

π

− e + 1.267



1+ 5 9

a ) Utilice (2.3) para probar que

x

− µ · |x| ≤

f l(x)

≤

·| |

x + µ x .

b ) Para cada uno de los valores que se presentan a continuación, determine el intervalo

de valores representados por redondeo al número de dígitos indicados en cada caso: 1)

1.325 a 4 dígitos

3)

2)

0.00123 a 3 dígitos

4)

−15.002 a 5 dígitos 0.1 a 1 dígito

c ) Para cada uno de los valores que se presentan a continuación, determine el intervalo

de valores representados por truncamiento al número de dígitos indicados en cada caso: 1)

−0.00010001 a 5 dígitos

2) 15012.06 a 7 dígitos

3)

0.2365 a 4 dígitos

4)

0.01 a 1 dígito

121

2.4. Operaciones aritméticas en precisión finita

4. Para n

∈

considere la integral



1

I n =

tn et−1 dt.

0

a ) Verifique que se cumple la siguiente relación de recurrencia:

I n = 1

− nI − , para todo n > 1. n 1

b ) Utilizando esta relación de recurrencia. Calcule I n para n = 1, . . . , 8 en el sistema

de precisión finita (10, 4, 10, 10) utilice ambos métodos de aproximación. Muestre los resultados en una tabla de la siguiente forma:

−

n I n por truncamiento 1 .. .

I n por redondeo

.. .

.. .

8 c ) Determine los errores relativos para cada valor I n , con n = 1, . . . , 8 y comente sus

resultados. ♦ De

programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5. Implemente en MATLAB una rutina que reciba un número real x y un entero positivo t y retorne el valor de f l(x) en (10, t, , + ), utilizando truncamiento.

−∞ ∞

♦ Avanzados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

∈ − {0} se cumple que |δ | ≤ µ.

6. Muestre que para cualquiera de las representaciones de x f l(x) = x(1 + δ ),

2.4.

con

Operaciones aritméticas en precisión finita

Como ya se ha mencionado, es posible que durante los cálculos entre valores punto flotantes se obtenga un nuevo valor que no se representa en el conjunto de números reales de precisión finita. Por ejemplo, considere el conjunto (10, 3, 2, 3) con todos sus elementos normalizados, de entre los cuales dos de ellos corresponden a:

−

a = 0.111

3

× 10

y

b = 0.120

3

× 10 .

(10, 3, 2, 3) ya que el exponente de c excede a U = 3 Nótese que c = ab = 0.133 105 lo cual se conoce como desbordamiento o “overflow”. Similarmente, si a = 0.1 10 −1 y b = 0.2 10−1 , entonces c = ab = 0.2 10−3 tiene exponente menor a L = 2, por lo que c no esta en el conjunto de precisión finita debido a un subdesbordamiento o “underflow”.

×

×

∈

−

×

−

×

Para tratar de evitar el desbordamiento y el subdesbordamiento de las operaciones básicas: suma, resta, multiplicación y división, lo que se hace es considerar las representaciones después de efectuar cada operación, tal y como se muestra en la tabla 2.3.

122


Operación Representación x + y

f l (f l(x) + f l(y))

x

f l (f l(x)

−y x·y x÷y

− f l(y)) f l (f l(x) · f l(y)) f l (f l(x) ÷ f l(y))

Tabla 2.3: Representaciones de las operaciones aritméticas. Ejemplo 2.19. Considere la operación 1 + ε en

(10, 2, 3, 3) con ε = 10−3 . Entonces se tiene que:

−

1

f l(1) = 0.10

× 10 , f l(ε) = 0.10 × 10− , 2

f l(1) + f l(ε) = 1.001, f l (f l(1) + f l(ε)) = f l(1.001) = 0.10

1

× 10

= 1.

Por lo tanto, se tiene que 1 + ε = 1. Ejemplo 2.20. Considere el vector x = (1, 0.01, 0.01, . . . , 0.01)T 1+10 y determine la norma Euclidea x 2 en (10, 2, 3, 3), además del error relativo de la representación con redondeo.

 

−

Solución. Considere

4

∈

(10, 2, 3, 3) y nótese que

−

1

f l(1) = 0.10

× 10 , f l(0.01) = 0.10 × 10− . √ Luego, como x = 1 + 10 · 0.01 , entonces se tiene que: f l(1 ) = f l (f l(1) · f l(1)) = f l(1) = 0.10 × 10 , f l (0.01 ) = f l (f l(0.01) · f l(0.01)) = f l(0.0001) = 0.10 × 10− , f l (10 ) = 0.99 × 10 (¡desbordamiento!). f l (10 · 0.01 ) = f l (f l(10 ) · f l(0.01 )) = f l(0.099) = 0.99 × 10− . f l (1 + 10 · 0.01 ) = f l (f l(1 ) + f l (10 · 0.01 )) = f l(1.099) = 0.11 × 10 . √ f l 1 + 10 · 0.01 = f l f l (1 + 10 · 0.01 ) = f l(1.0488) = 0.10 × 10 . 1

2

2

4

2

2

1

2

3

4

3

4

2

2



4

2

4

2

4

2

2

2

 

1

4

2

2

4

1

2



1

123


Por lo tanto, f l ( x 2 ) = 1 en (10, 2, 3, 3), mientras que el valor exacto corresponde a x 2 = 2. Finalmente, el error relativo corresponde a

  √



−

  −       x

2

f l ( x 2)

x

=

2

√ 2 − 1 √ 2 ≈

0.2928932188 < 5

1

× 10− .

Nótese que la representación f l ( x 2 ) sólo tiene un dígito significativo.





Observación: Del ejemplo anterior se puede apreciar que en las operaciones aritméticas con números punto flotantes se cometen errores, debido a que la representación de un número no siempre es exacta. Así, al aplicar una operación y luego encontrar su representación se comete un nuevo error de aproximación, es decir, se da una propagación del error , donde al propagarse el error, este se incrementa ocasionando que los resultados no reflejen el esperado, tal y como se muestra en el siguiente fragmento de código de MATLAB: clc, format long e H = hilb(9); % Matriz de Hilbert de 9x9 B = inv(H); % Inversa de H I = H*B % Matriz identidad

en el cual se calcula la inversa B de la matriz H de Hilbert y luego se multiplican ambas matrices. El resultado esperado debería corresponder a la matriz identidad, sin embargo, esto no ocurre debido a que el resultado es el siguiente: I = Columns 1 through 3 1.000000000018190e+000 2.810847945511341e-007 3.333334461785853e-007 3.374734660610557e-007 3.271597961429507e-007 3.122804628219456e-007 2.963952283607796e-007 2.809101715683937e-007 2.663236955413595e-007

6.984919309616089e-010 1.000000002793968e+000 -6.984919309616089e-009 -1.396983861923218e-008 -1.373700797557831e-008 -1.490116119384766e-008 -1.769512891769409e-008 -1.676380634307861e-008 -1.618172973394394e-008

-2.235174179077148e-008 -7.450580596923828e-009 9.999999888241291e-001 2.980232238769531e-008 -3.352761268615723e-008 -3.352761268615723e-008 -7.450580596923828e-009 3.725290298461914e-009 -7.450580596923828e-009

4.768371582031250e-007 2.384185791015625e-007 0 4.768371582031250e-007 1.000000476837158e+000 5.960464477539063e-007 2.980232238769531e-007 0

0 -2.384185791015625e-007 1.192092895507813e-006 0 -2.384185791015625e-007 1.000000953674316e+000 2.384185791015625e-007 -1.192092895507813e-006

Columns 4 through 6 0 -2.980232238769531e-008 2.980232238769531e-008 1.000000119209290e+000 5.960464477539063e-008 8.940696716308594e-008 -8.940696716308594e-008 0

124

Capítulo 2. Aritmética de precisión finita 5.960464477539063e-008

5.960464477539063e-008

3.576278686523438e-007

4.768371582031250e-007 1.668930053710938e-006 5.960464477539063e-007 3.576278686523438e-007 1.192092895507813e-007 -1.192092895507813e-007 -1.192092895507813e-007 1.000001311302185e+000 5.960464477539063e-007

5.960464477539063e-008 1.192092895507813e-007 -2.980232238769531e-007 1.192092895507813e-007 5.960464477539063e-008 1.192092895507813e-007 -8.940696716308594e-008 0 1.000000029802322e+000

Columns 7 through 9 1.192092895507813e-006 0 7.152557373046875e-007 2.861022949218750e-006 1.430511474609375e-006 1.192092895507813e-006 1.000001192092896e+000 -2.384185791015625e-007 -4.768371582031250e-007

Conocer estos errores introducidos por la aritmética de precisión finita es muy importante, sobre todo en aplicaciones en las que se realizan muchas operaciones, ya que conforme se realizan las operaciones el error se propaga, se acumula y en algunos casos puede llegar a crecer exponencialmente, provocando grandes errores en el resultado final. Dos ejemplos de esto son: el fallo de un misil Patriot: murieron 28 soldados norteamericanos a causa de los errores numéricos por utilizar truncamiento en lugar de redondeo en el sistema que calcula el momento exacto en que debe ser lanzado el misil y la explosión del cohete Ariane 5: explotó 40 segundos después de su despegue, tomando un costo de 75000 millones de euros, a causa de un fallo en el sistema de guiado de la trayectoria, provocado al intentar convertir un número en punto flotante. En algunos casos el error puede ser identificado y corregido, por ejemplo MATLAB ofrece la instrucción invhilb() para aproximar de una mejor manera la inversa de una matriz de Hilbert. Además, el ejemplo siguiente muestra como evitar otro error usual. Ejemplo 2.21 (Evitar errores en los cálculos). Determine en

(10, 3, 5, 5) las soluciones de la ecuación x2 + 28x + 1 = 0.

−

Solución. Las soluciones de la ecuación ax2 + bx + c = 0 están dadas por las fórmulas:

−b + √ b − 4ac 2

x1 =

2a

−b − √ b − 4ac . 2

x2 =

y

2a

Y nótese que fl

√ −  b2

4ac

= fl

√

0.784

×

103

= f l(27.92848009)

101

 √

− 0.400 × = 0.279 × 10 .

= fl

0.780

2

Así utilizando truncamiento, los valores de las soluciones son dadas por: f l(x1 ) = f l

−

0.280

× 10

2

+ 0.279 0.2000 101

×

× 10

2



=

1

−0.500 × 10− ,

×

103



,

125


−

2

0.280

2



× 10 − 0.279 × 10 = −0.279 × 10 . f l(x ) = f l 0.2000 × 10 Sin embargo, los valores exactos son: x = −0.03575995623 . . . y x = −27.96424004 . . . con lo 2

1

2

1

2

que los errores relativos vienen dados por:

|x − f l(x )| ≈ |x | 1

1

0.3982120022

y

1

|x − f l(x )| ≈ |x | 2

2

0.0022972210.

2

Así, x2 está bien representado con tres dígitos significativos, pero x1 sólo tiene un dígito significativo. La razón es debido a que el numerador de la fórmula de x1 da un valor pequeño, provocado por una resta . Esta resta puede ser eliminada al racionalizar el numerador como se muestra: 2c b + b2 4ac b + b2 4ac = x1 = . 2a b + b2 4ac b + b2 4ac Con esta nueva fórmula se tiene que

− √ −

f l(x1 ) = f l



√ − √ −

·

0.200 101 0.100 101 0.280 102 + 0.279 102

× ×

·

× ×

√ − −



=

1

−0.357 × 10− ,

la cual tiene un error relativo de 0.0016766304. Esta nueva aproximación para x1 es mucho  mejor debido a que tiene tres dígitos significativos. Observación: Las operaciones de suma y multiplicación flotantes son conmutativas, pero no son ni asociativas ni distributivas. Es decir, por ejemplo, aunque x + y + z = (x + y) + z = x + (y + z ), se puede llegar a obtener que f l(f l(x + y) + z ) = f l(x + f l(y + z )),



y por lo tanto el orden con el que se realizan las operaciones es muy importante tanto por cómo se afecta el valor, a cómo se acumulan los errores. Un orden adecuado puede llegar a reducir los errores en el resultado final, mientras que un orden inadecuado puede incrementar éstos, tal y como se puede apreciar en el ejercicio 7 de la página 126.


1. Evalúe utilizando truncamiento el polinomio p(x) = x 3 5x2 + 6x + 0.55 en x = 2.73, con (10, 3, p(x) , + ). Además, calcule el error relativo.

∈

−∞ ∞

−

2. El polinomio p(x) del ejercicio anterior se puede reescribir como p(x) = ((x 5)x + 6) x+ 0.55. Así, con esta nueva formulación repita el ejercicio anterior y compare ambos resultados.

−

3. Construya ejemplos que muestren que la ley asociativa y la ley distributiva de la suma y la multiplicación no se cumplen para números punto flotantes.

126


4. Considere (10, t,

−∞, +∞) y calcule f l A =

 

 

AT A para t = 4 y

1 1 10−4 0 0 10−4

 

.

Repita los cálculos para t = 9 y compare los resultados. 5. Considere (10, 4,

−∞, +∞) y

 −  ÷

1 0.1666 6 ¿cuántos dígitos significativos se obtuvieron? a =

0.1666,

6. Considere los puntos (x1 , y1 ) y (x2 , y2 ) con y1 = y2 . Para la recta que pasa por ellos, se tienen dos formas de calcular la abscisa de cero, dadas por: (x2 x1 )y1 x1 y2 x2 y1 y x = x1 x = . y2 y1 y2 y1 a ) Demuestre que ambas fórmulas son algebraicamente correctas.



− −

−

− −

b ) Utilice los valores x 1 = 1.31, y 1 = 3.24, x 2 = 1.93 y y 2 = 4.76 en

(10, 3, , + ), para determinar la abscisa de cero por redondeo con ambas fórmulas. ¿Cuál método es mejor y por qué?

♦ De

7.

−∞ ∞

programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

a ) Implemente una función en MATLAB [nrm] = L2normaIzq(x) la cual recibe un vector x y devuelve x 2 , pero esta es calculada de la siguiente manera:

     · · · n

x

2

=

x2i =

(

((x21 + x22 ) + x23 ) + . . . + x2n ).

i=1

b ) Implemente una función en MATLAB [nrm] = L2normaDer(x) la cual recibe un vector x y devuelve

    n

x

2

=

x2i =

i=1



x21 + . . . + x2n−2 + (x2n−1 + x2n )

 ·· · 

.

c ) Muestre mediante un ejemplo numérico que estas funciones no retornan el mismo

valor para todos los vectores dados. ♦ Avanzados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8.

a ) Calcule la norma Euclidea del vector x = (1, 101 , 102 ) en aritmética exacta. b ) Realice el mismo cálculo que del inciso anterior pero en aritmética de precisión finita

utilizando el sistema (10, 3, 4, 4). El método de aproximación es por truncamiento. Explique su resultado.

−

c ) Encuentre un método para evitar el problema de cálculo que se presentó en el inciso n

ax xi . anterior. Sugerencia: considere m = m´ i=1

| |

127

2.5. Errores debido a la cancelación y cálculos recursivos

2.5.

Errores debido a la cancelación y cálculos recursivos

Es conocido que la gran cantidad de cálculos con números punto flotantes acumulan el error considerablemente. Sin embargo, los errores de redondeo pueden ser desastrosos, incluso en un sólo cálculo. Por ejemplo, considere la resta de los números: a = 0.54617

b = 0.54601,

y

cuyo valor exacto es: d = a Ahora, considere la misma operación en se tiene que:

−b

= 0.00016.

(10, 4,

−∞, +∞) por medio de redondeo, con lo que

f l(a) = 0.5462 (cuatro dígitos significativos), f l(b) = 0.5460 (cuatro dígitos significativos), f l(d) = f l(a)

− f l(b) = 0.0002.

Pero ¿qué tan buena es esta aproximación? Para responder a esto considere el error relativo dado por d f l(d) = 0.25 (el cual se considera bastante grande). d

| − | ||

La razón de obtener un error “tan grande” es debido a que en aritmética de cuatro dígitos los números 0.5462 y 0.5460 son de casi del mismo tamaño (valores casi idénticos salvo en los últimos dígitos). Así, cuando se resta el segundo número al primero, los dígitos más significativos se cancelaron y el dígito menos significativo fue dejado en la respuesta . A este fenómeno se le denomina cancelación catastrófica y ocurre cuando dos números de aproximadamente el mismo tamaño se restan.

2.5.1.

Evitar la cancelación

Como se muestra en el ejemplo 2.21 la solución de la ecuación cuadrática ax 2 + bx + 1 = 0, con a = 0, dada por b + b2 4ac b b2 4ac x1 = y x2 = , 2a 2a



− √ −

− − √ − √ puede producir una cancelación catastrófica cuando −b y b − 4ac tienen aproximadamente 2

el mismo tamaño, con respecto a las operaciones aritméticas. En el mismo ejemplo se puede apreciar que la forma de evitar la cancelación catastrófica es utilizando las fórmulas equivalentes: x1 =

−

√ − 4ac

b + sgn(b) b2 2a

y

x2 =

c , ax1

con x1 = 0,



(2.4)

128


donde sgn(x) es la función signo definida por: sgn(x) =

 −

1 si x > 0, 0 si x = 0, 1 si x < 0.

En el caso de que x1 = 0, entonces se tiene que x2 =

−

b a

.

Ejemplo 2.22. Resuelva x2

5

− 10 x + 1 = 0 en

(10, 8, 50, 50).

−

Solución. Utilizando la fórmula original se tiene que: x1 = x2 =

√ 10 − 4 2 √ − 10 − 4

105 +

10

105

10

= =

2

105 + 105 = 105 (respuesta correcta), 2 105

− 10 2

5

= 0

(respuesta incorrecta).

Ahora, utilizando las fórmulas (2.4) se tiene que: x1 = 100000.00, x2 =

1.0000000 = 0.000010000 (respuesta correcta). 100000.00



Ejemplo 2.23. Evalúe la función f (x) = e x

− x − 1 en x = 0.01, utilizando el sistema

(10, 5,

−∞, +∞).

Solución. Al evaluar f (0.01) directamente en la expresión se tiene que: f (0.0 1 ) = 1.0101

− (0.01) − 1

= 0.00001.

Luego, la respuesta correcta es 0.000050167, con lo que el error relativo es dado por: 0.0001 0.000050167 = 0.99 0.00005016

−

0

× 10 ,

e indica que no se tienen dígitos significativos. Afortunadamente, una manera de evitar la cancelación es utilizando la serie convergente para e x dada por x

e

=

+∞

 n=0

xn x2 x 3 = 1 + x + + + . . . , 2 3! n!

y con ella se tiene que x

e

−x−



−

x 2 x 3 1 = 1 + x + + + . . . 2 3! x2 x 3 x 4 = + + + . . . 2 3! 4!

x

− 1,

129


Ahora, para x = 0.01, esta fórmula da: x

e

−x−

(0.01)2 (0.01)3 (0.01)4 1 = + + + . . . , 2 3! 4! = 0.00005 + 0.000000166666 + 0.00000000004166 + . . . , = 0.000050167 (cinco dígitos significativos). 

2.5.2.

Cálculos recursivos

Los cálculos recursivos son aquellos que se realizan de modo que el cálculo de un paso depende de los resultados de los pasos anteriores. En tales casos, incluso si el error en el primer paso es insignificante, debido a la acumulación y magnificación del error a cada paso, el error final puede ser bastante grande, obteniendo una respuesta completamente errónea. Ejemplo 2.24. Determine el valor de la siguiente integral para varios valores de n.



1

I n =

tn et−1 dt.

0

Solución. Integrando por partes se tiene que:



1

I n =

n t 1

t e − dt = (tn et−1 )

0

 − 

1

1 0

n 1 t 1

nt − e − dt = 1

0

 −

1

n

tn−1 et−1 dt,

0

es decir, se tiene la fórmula por recurrencia: I n = 1

− nI − , para n = 2, 3, . . . n 1

Así, si I 1 es conocido, entonces para diferentes valores de n, I n puede ser calculado utilizando la fórmula recursiva anterior. Luego, para (10, 6, , + ) y el valor inicial I 1 = 0.367879, el cual se obtiene al representar e−1 en el sistema, se pueden determinar los valores:

−∞ ∞

I 1 =

0.367879,

I 2 =

0.264242,

I 3 =

0.207274,

I 4 = .. .

0.170904, .. .

I 9 =

−0.068480 (valor incorrecto).

La integral I n es positiva sobre [0, 1], debido a esto el resultado de I 9 es completamente incorrecto por ser negativo. La razón de esto, es que el error al calcular I 2 fue 2 veces el error al calcular I 1 y el error al calcular I 3 fue 3 veces el error en I 2 (por lo tanto, el error en este paso fue

−

−

130


exactamente 6 veces el error en I 1 ). Así, el error al calcular I 9 fue ( 2)( 3)( 4) . . . ( 9) = 9! veces el error en I 1 , el cual se debió a redondear e−1 usando seis dígitos significativos y corresponde a 4.412 10−7 . Sin embargo, este pequeño error multiplicado por 9! corresponde a 0.1601 el cual es bastante grande, ya que no tiene dígitos significativos.

− − − · ·−

×

Para corregir este problema, mejorando el resultado, basta con reordenar la recurrencia para disminuir el incremento del error en cada paso, de la forma: I n−1 =

1

− I , para n = . . . , 3, 2. n

n Con esta nueva recurrencia el error en cada paso será reducido por un factor de n−1 . Así, iniciando con un valor grande de n (por ejemplo, n = 20) y recorriendo la recursión hacia atrás, se calcula I 9 con seis dígitos de precisión. Para obtener el valor inicial, nótese que:



1

I n =

n t 1

t e − dt

0

 ≤

1

tn dt =

0

1 . n + 1

1 = 0 y con ello se obtiene que Así, para n = 20 se tiene que I 20 21 , por lo que considere I 20  I 9 = 0.0916123, el cual es el valor exacto con seis dígitos de precisión.

≤

Observación: La razón de obtener toda la precisión en la solución del ejemplo anterior, fue 1 1 debido a que el error en I 20 fue de a lo más 21 , el cual fue multiplicado por 20 al calcular I 20 , 1 1 dando un error de a lo más 20 21 = 0.0024 al calcular I 19 y así sucesivamente.

·

2.5.3.

Errores en la ortonormalización de Gram-Schmidt

El proceso de ortogonalización de Gram-Schmidt B.4, consiste en a partir de una base cualquiera n de vectores linealmente independientes x1 , x2 , . . . , xm , con xi , para así obtener una n nueva base q1 , q2 , . . . , qm ortonormal de .

{

{

}

}

∈

En la sección 1.7.5 se muestra el algoritmo del proceso de Gram-Schmidt, el cual se resume a continuación: Algoritmo 2.1: GramSchmidt. n Recibe: x1 , x2 , . . . , xm linealmente independientes. n× m Retorna: Q = ( q 1 , q2 , . . . , qm ) ortonormal y R

∈

1 2 3

4 5

Para j 1 hasta m hacer Para i 1 hasta j 1 hacer rij x j , qi j −1 rij qi q j x j

←

← ← 

← − r ← q  q q ← r j

6



 i=1

j 2

jj j

jj

−

·

∈

∈

m m

×

triangular superior.

131


Nótese que los vectores qi y las constantes rij forman una factorización de la matriz X = n×m ( x 1 , x2 , . . . , xm ) , de la forma X = QR, cuando rank(X) = m. La matriz Q es T una matriz ortogonal (Q Q = I m) y R es una matriz triangular superior. A esta factorización se le conoce como la factorización QR . Los principales usos de esta factorización es en la resolución de sistemas rectangulares y el cálculo de valores propios.

∈

El método de Gram-Schmidt es muy útil en la práctica, pero tiene el inconveniente de ser un método numéricamente inestable , lo que significa que tiene serias dificultades numéricas, entre ellas la gran variedad de cálculos en valores punto flotantes. Durante los cálculo de los qi , la cancelación se puede llevar acabo y en consecuencia el resultado final para Q no corresponde a una matriz ortogonal , tal y como se puede apreciar en el siguiente ejemplo.

Ejemplo 2.25. Sea ε > 0 tal que f l(1 + ε 2 ) = 1 y X = ( x 1 , x2 , x3 ), con x1 = ( 1, ε, 0, 0 )T , T T x2 = ( 1, 0, ε, 0 ) y x3 = ( 1, 0, 0, ε ) . Utilice el proceso de Gram-Schmidt para obtener la ortonormalización de los vectores x1 , x2 y x3 . Solución. Al proceder con el método de Gram-Schmidt usual, se obtiene:

x  = √ 1 1 2

2

+ ε2 + 02 + 02 = 1, por lo que: T

q1 = x1 = ( 1, ε, 0, 0 ) .

ˆ 2 = x 2 q

− x , q  q = ( 0, −ε, ε, 0 ) , qˆ  = √ 0 + ε + ε + 0 = ε√ 2, por lo que: ˆ −1 1 q = 0, √ , √ , 0 q = . qˆ  2 2 −x , q  q −x , q  q = ( 0, −ε, 0, ε ) , qˆ  = √ 0 + ε + 0 + ε = ε√ 2, 2

1

T

1

2 2

ˆ 3 = x 3 q

3

1

1

3

2

q3 =



=

qˆ 

3 2

2

T

T

por lo que:

2



2

ˆ3 q

2



2

2

2

2 2

0,

2

3 2

−1 , 0, √ 1 √ 2 2



2

2

2

T

.

Luego, se tiene que

T

Q Q =

 

1 0 0

ε 1 √ 2 1 √

− −

2

0 1 √

2

0

0 0 1 √

2

   · 

1 ε 0 0

0 1 √ 2 1 √

−

2

0

0 1 √

−

2

0 √ 1

2

    −  − =

1 √ ε 2 √ ε 2

ε

ε

− √ − √ 2

1

1 2

1 2

1

2

 

,

0, QT Q no es la identidad, por lo que Q no es ortogonal, por lo que nótese que cuando ε produce una factorización incorrecta. La razón de esto es la cancelación catastrófica ocasionada porque los vectores columnas de X están muy juntos. 

→

132


Con el fin de obtener mejores resultados numéricos, el algoritmo 2.1 puede ser modificado al cambiar algunos de sus cálculos, obteniendo el siguiente algoritmo conocido como GramSchmidt Modificado (MGS), el cual determina la factorización QR de X, donde en el k-ésimo paso se calcula la k-ésima columna de Q y la k-ésima fila de R (en el algoritmo 2.1 se calcula la k-ésima columna de Q y la k-ésima columna de R). Algoritmo 2.2: GramSchmidtModificado. n Recibe: x1 , x2 , . . . , xm linealmente independientes. n× m Retorna: Q = ( q 1 , q2 , . . . , qm ) ortonormal y R

∈

1 2 3 4 5

6 7

∈

∈

m m

×

triangular superior.

Para j

← 1 hasta m hacer q ← x Para i ← 1 hasta j − 1 hacer r ← q , q  q ← q − r · q r ← q  q ← j

j

ij

j

j

i

j

ij

i

j 2

jj

qj rjj

j

Así, si se resuelve el ejemplo 2.25 usando el método MGS, se obtiene que:

T

Q Q =

y cuando ε

 

1 0

ε √ 1

0 √ 1

− − √ − √ 2 1 6

0

T

→ 0, Q

2 1 6

    ·  0 0

2 3

1 ε 0 0

Q es la matriz identidad.

0 1 √ 2 1 √

−

2

0

  −   −  −   − 0 1 √ 6 1 √ 6 2 3

=

1 √ ε 2 ε √ 6

−

√ ε

1 0

2

−

√ ε

6

0 1

 

,


1. ¿Cuál problema se pueden notar al resolver las siguientes ecuaciones cuadráticas: a ) x2

6

− 10 x + 1 = 0, b ) 10− x − 10 x + 10 10 2

10

10

= 0,

usando la conocida fórmula

−b ± √ b − 4ac ? 2

x =

2a ¿Qué estrategia se puede seguir para corregir el problema? Resuelva las ecuaciones con la nueva estrategia sugerida.z

133

2.5. Errores debido a la cancelación y cálculos recursivos ♦ De

2.

programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

a ) Explique que hace el siguiente fragmento de código de MATLAB:

A = 1/3 * [-20e03 - 1e20, 2e20 + 10e03; -10e03 - 2e20, 4e20 + 5e03]; x = [1; 1]; b = A * x; A \ b b ) Ejecute el fragmento anterior en MATLAB y explique qué anormalidad encuentra

en la respuesta. Comente significativamente porqué se produce esta anormalidad. ♦ Avanzados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. Muestre cómo reescribir las siguientes operaciones para evitar la pérdida de dígitos significativos. Realice un experimento numérico en cada uno para darle soporte a su respuesta. a ) ex

− x − 1 para valores negativos de x. √ b ) x + 1 − x para valores grandes de x. c ) x − sen(x) para valores de x cercanos a cero. d ) 1 − cos(x) para valores de x cercanos a cero. e −1 e ) para |x| ≪ 1 (|x| mucho menor a 1). x 1 − cos(x) f ) para valores pequeños de x. 4

2

x

x2 4. Muestre que la integral



1

ti yi = dt, 0 t + 5 puede ser calculada utilizando la fórmula por recurrencia: 1 5yi−1 . yi = i Determine y1 , y2 , . . . , y10 usando esta fórmula, tomando

−

y0 = ln(t + 5)



1 0

= ln6

− ln5 = ln(1.2).

¿Qué anormalidades se pueden observar en los cálculos? Explique que está sucediendo y luego reordene la recurrencia para que los valores de yi pueden ser determinados con mayor precisión. 5. Encuentre la factorización QR de la matriz A =

 

9 12 0

−6 −8

20

 

,

y utilícela para encontrar la solución del sistema Ax = b , donde b = (300, 600, 900)T .

Aritmetica Finita

Recommend Documents