Elementos de Calculo Numerico, Notas de estudo de Matemática

Baixe Elementos de Calculo Numerico e outras Notas de estudo em PDF para Matemática, somente na Docsity!

ELEMENTOS DE C ´ALCULO NUM´ERICO

Ricardo G. Dur´an, Silvia B. Lassalle y Julio D. Rossi

CAP´ıTULO 1

Punto flotante y redondeo

El objeto de esta secci´on es analizar la representaci´on de los n´umeros en una computadora y la propagaci´on de los errores de redondeo al realizar c´alculos.

Como la cantidad de informaci´on que puede guardarse en una computadora es finita, la m´aquina trabajar´a s´olo con un conjunto finito de n´umeros. A ´estos los llamaremos n´umeros de m´aquina. En consecuencia, toda vez que de nuestros datos o c´alculos surja un n´umero que no pertenece a este conjunto finito, ´este deber´a ser reemplazado por una aproximaci´on (el n´umero de m´aquina m´as cercano). Este reemplazo da lugar a lo que llamamos errores de redondeo.

Al realizar c´alculos estos errores de redondeo se propagan y esto puede llevar a resultados totalmente incorrectos como veremos en algunos ejemplos simples.

En las aplicaciones del c´alculo num´erico es pr´acticamente imposible determinar exactamente la magnitud de los errores de redondeo. Lo que si puede hacerse, y es de fundamental importancia, es identificar las posibles causas de que los errores se propaguen m´as de lo admisible. Esto permite mejorar los algoritmos o determinar que m´etodo es m´as conveniente para resolver un problema. Un claro ejemplo de esto, que veremos m´as adelante, aparece cuando se utiliza el m´etodo de eliminaci´on de Gauss para resolver un sistema de ecuaciones lineales. En este caso, el an´alisis de la propagaci´on de errores permite determinar la forma m´as eficiente de aplicar el m´etodo.

Por otra parte, es fundamental distinguir cuando la propagaci´on excesiva de errores se debe a que el algoritmo utilizado es “malo” o inestable o a que el problema en s´ı mismo est´a “mal condicionado”. En el primer caso se puede (se debe!) tratar de mejorar el m´etodo de resoluci´on mientras que en el segundo caso el problema es m´as esencial. Los ejemplos que presentaremos ilustrar´an estos dos casos.

En lo que sigue supondremos que los n´umeros de m´aquina son los que aparecen en la pantalla. Esto no es exacto pues en realidad la computadora opera internamente con los n´umeros desarro- llados en base 2 y no en base 10. Este abuso de lenguaje es s´olo para mayor claridad (el lector podr´a observar que todo nuestro an´alisis puede repetirse trabajando en base 2).

Observemos primero que un n´umero real cualquiera, x ∈ IR, x > 0, puede escribirse como

x = 0, a 1 a 2 .....ak...... 10 l^ = r × 10 l,

≤ r < 1 (es decir, a 1 6 = 0)

1. PUNTO FLOTANTE Y REDONDEO 3

pr´oximo a x entre x′^ y x′′. A esta forma de aproximar se la conoce por redondeo y es la que usan habitualmente las computadoras.

Veamos que error cometemos al aproximar un n´umero por redondeo. Usaremos la notaci´on

x∗^ = x redondeado

El error absoluto ser´a

|x − x∗| ≤

10 m^

10 l

Mientras que el error relativo se puede acotar de la forma siguiente

|x − x∗| |x|

10 l−m 0 , a 1 a 2 ....am.... 10 l y como

0 , a 1 a 2 ....am.... ≥

se tiene que

|x − x∗| |x|

≤ 5 × 10 −m

Es decir que el error relativo es del orden de 10−m^ si nuestra m´aquina trabaja con m d´ıgitos.

Es importante observar que, si bien el error absoluto que introduce el redondeo depende de la magnitud del n´umero, el relativo, que es el m´as significativo, es independiente de ´esta, est´a con- trolado en t´erminos de la cantidad de d´ıgitos con la que trabaja nuestra computadora (Ejercicio: meditar que tiene que ver esto con la distribuci´on no uniforme de los n´umeros de m´aquina).

Si tenemos

x∗^ = 0, a 1 a 2 ....am 10 l, a 1 6 = 0

decimos que conocemos a x con m d´ıgitos significativos, lo que equivale, seg´un vimos, a conocerlo con un error relativo del orden de 10−m. Observar la importancia de la condici´on a 1 6 = 0: de lo contrario los d´ıgitos dejan de ser significativos.

Observemos que

x∗^ = x(1 + δ)

4 1. PUNTO FLOTANTE Y REDONDEO

con

|δ| ≤ ε = 5 × 10 −m

Este ε es el error de redondeo unitario, es decir, que el error que se comete al aproximar x por x∗^ es xδ con |δ| ≤ ε, de esta forma el error |x − x∗| ser´a menor o igual que ε|x|.

Este valor, ε, depende de la m´aquina y se lo llama el ε de la m´aquina. Recordemos que, seg´un lo comentado antes, el verdadero ε de m´aquina no es exactamente ´este debido a que la computadora trabaja en base 2. Pero desde el punto de vista pr´actico s´olo interesa el orden de ε (y ´este s´ı es correcto!).

Ejercicio: calcular el ε exacto suponiendo que la m´aquina trabaja en base 2 con k d´ıgitos.

Otra forma de interpretar el significado del ε de la m´aquina es la siguiente: ε nos dice cu´al es el menor n´umero tal que, en la m´aquina,

1 + ε 6 = 1

O sea, si se le suma a 1 algo menor que ε, “desaparece” debido al redondeo. En efecto, seg´un la m´aquina tendremos

1 + 4 × 10 −m^ = 1, 0 ....04 = 0, 10 .... 00 × 10 = 1

en cambio,

1 + ε = 1, 0 ....06 = 0, 10 .... 01 × 10 6 = 1

Si sumamos exactamente ε el resultado depender´a de como redondee la m´aquina en el caso en que x equidista de dos n´umeros de m´aquina, es decir, cuando la cifra m + 1 es justo 5.

M´as en general, el orden del menor n´umero que sumado a un x da, en la m´aquina, un resultado distinto de x es ε|x|.

Es fundamental observar que ε no es el menor n´umero que se puede representar en la m´aquina (y est´a muy lejos de ´este!). Este ´ultimo depende de M 1 y no de m.

Veamos ahora algunos de los problemas que surgen debido al redondeo.

Empecemos por sumar dos n´umeros. Como vimos, si sumamos dos n´umeros de ordenes muy distintos, el m´as chico puede “desaparecer”. Por ejemplo, si m = 5 y

x = 78473000; y = 24

tenemos

6 1. PUNTO FLOTANTE Y REDONDEO

y =

x^2 + 1 − 1

para valores de x peque˜nos. Si lo hacemos directamente, estamos en la situaci´on anterior. En cambio podemos proceder de la siguiente manera:

y =

x^2 + 1 − 1 =

x^2 + 1 − 1)(

x^2 + 1 + 1) (

x^2 + 1 + 1)

x^2 (

x^2 + 1 + 1)

y utilizar la ´ultima expresi´on para calcular y. Si los c´alculos fueran exactos ambas f´ormulas dar´ıan el mismo resultado pero, debido al redondeo, dan resultados distintos. Por ejemplo, trabajando con 5 d´ıgitos, si x = 0, 0312 obtenemos con la primera f´ormula y = 0, 0004 (un solo d´ıgito significativo si bien conociamos x exactamente). Mientras que con la segunda, y = 0 , 00048662 (que tiene cuatro d´ıgitos significativos correctos).

El mismo problema surge al calcular y = x − sin x para x peque˜no. En este caso se puede usar el desarrollo de Taylor,

y = x − sin x =

x^3 3!

x^5 5!

x^7 7!

y calcular y sumando los primeros t´erminos de la serie.

Otro caso en el que la cancelaci´on de d´ıgitos significativos puede presentarse es en la resoluci´on de una ecuaci´on cuadr´atica

ax^2 + bx + c = 0

si utilizamos la f´ormula habitual

x 1 =

−b +

b^2 − 4 ac 2 a

y x 2 =

−b −

b^2 − 4 ac 2 a

Consideremos por ejemplo,

x^2 − 105 x + 1 = 0

Los primeros d´ıgitos exactos de las ra´ıces son

x 1 = 99999. 99999

y x 2 = 0. 000010000000001

Usando la f´ormula para x 2 tenemos

1. PUNTO FLOTANTE Y REDONDEO 7

x 2 =

Si trabajamos con ocho d´ıgitos el 4 desaparece y x 2 resulta igual a cero. Otra vez hemos restado dos n´umeros parecidos!

Esto se puede solucionar calculando primero x 1 y luego obtener x 2 usando que x 2 = (^) axc 1.

En general, la forma correcta de encontrar las ra´ıces en el caso en que ac sea despreciable respecto de b, es calculando primero

x 1 =

−b − sign(b)

b^2 − 4 ac 2 a y luego la otra raiz como hicimos en el ejemplo. De esta forma se evita la p´erdida de d´ıgitos significativos.

Un problema fundamental del c´alculo num´erico es la resoluci´on de sistemas de ecuaciones lineales. Veamos como los errores de redondeo pueden afectar la soluci´on a´un en problemas de dos ecuaciones con dos inc´ognitas. Tratemos de resolver el siguiente sistema utilizando el m´etodo de eliminaci´on de Gauss,

( ε 1 1 1

x y

Pongamos, a modo de ejemplo, ε = 10−^6 y supongamos que la m´aquina trabaja con cinco d´ıgitos.

Multiplicando la primera fila por (^1) ε = 10^6 y rest´andosela a la segunda obtenemos

( 10 −^6 0 1 − 106

x y

pero, con el redondeo a cinco cifras nos queda

( 10 −^6 0 − 106

x y

(perdimos la informaci´on de la segunda ecuaci´on!).

Mientras que el resultado exacto es

( 10 −^6 0 − 999999

x y

Hasta aqui el error no parece grave. Pero veamos: si utilizamos la matriz obtenida con la m´aquina y despejamos de la segunda ecuaci´on obtenemos la soluci´on y′^ = 1 y luego, reemplazando en la primera, x′^ = 0.

1. PUNTO FLOTANTE Y REDONDEO 9

Hay casos de naturaleza esencialmente diferente en los cuales el problema que se quiere resolver est´a “mal condicionado”. Esto significa que peque˜nos cambios en los datos implican grandes cambios en la soluci´on. Esto hace que los errores de redondeo puedan amplificarse mucho independientemente del m´etodo que usemos para resolverlo.

Veamos un ejemplo de esto. Supongamos que nuestra m´aquina trabaja con 3 d´ıgitos y trunca. Resolvamos el sistema

(

1. 780 0. 563
2. 913 0. 659

x y

La soluci´on exacta es x = 1, y = −1.

Teniendo en cuenta lo visto antes, intercambiamos filas antes de hacer la eliminaci´on de Gauss. Obtenemos

(

1. 913 0. 659 0 0. 001

x y

y en consecuencia y′^ = 1 y x′^ = − 0 .443 que no tiene nada que ver con la soluci´on exacta. En particular, el error es mayor que la soluci´on misma!

Lo que sucede en este ejemplo es que la matriz est´a “mal condicionada” (m´as adelante precisare- mos lo que esto significa) y habr´a problemas independientemente del algoritmo que utilicemos.

Otro ejemplo de problema “mal condicionado” es el siguiente. Las ra´ıces de

(x − 2)^2 = 10−^6

son x 1 = 2 + 10−^3 x 2 = 2 − 10 −^3

en cambio, las ra´ıces de

(x − 2)^2 = 0

son x 1 = x 2 = 2.

Este ejemplo trivial muestra que un peque˜no cambio en un coeficiente de la ecuaci´on polinomial puede dar lugar a un cambio de otro orden en las ra´ıces. En este caso, un cambio de 10−^6 en el t´ermino independiente origina un cambio de 10−^3 en las ra´ıces.

Un ejemplo m´as interesante es el estudiado por Wilkinson en 1963. Se trata de calcular las ra´ıces de

p(x) = (x − 1)(x − 2)......(x − 19)(x − 20) = x^20 − 210 x^19 + ....

10 1. PUNTO FLOTANTE Y REDONDEO

Wilkinson demostr´o que al cambiar el coeficiente −210 por − 210 − 2 −^23 las ra´ıces 16 y 17 se transforman en el par complejo

16. 73 ... + i 2. 812 ... 16. 73 ... − i 2. 812 ...

Para finalizar, veamos otro ejemplo de algoritmo inestable. El problema consiste en calcular

En =

0

xnex−^1 dx n = 1, 2 , ...

Integrando por partes se obtiene

En =

0

xnex−^1 dx = xnex−^1 |^10 −

0

nxn−^1 ex−^1 dx

es decir

En = 1 − nEn− 1

y es f´acil ver que

E 1 = 1/e

con lo que tenemos definida la sucesi´on En en forma recursiva.

Usando esta relaci´on recursiva para calcular con una m´aquina que trabaje con seis d´ıgitos se obtiene,

E 1 ∼ 0. 367879

E 2 ∼ 0. 264242

E 3 ∼ 0. 207274

E 9 ∼ − 0. 0684800

cuando en realidad

E 9 ∼ 0. 0916

con lo que el resultado computacional es p´esimo.

En este caso lo que sucede es que el error de En− 1 se amplifica multiplic´andose por n. Entonces en nueve pasos se multiplica por 9!, obteni´endose un error de

9! × error inicial = 9! × 4. 412 × 10 −^7 ∼ 0. 1601

que resulta mayor que el verdadero E 9.

CAP´ıTULO 2

Normas y condicionamiento de una matriz.

Consideramos el sistema de n ecuaciones con n inc´ognitas

Ax = b

con A ∈ IRn×n, x ∈ IRn^ y b ∈ IRn^ y nos preguntamos cu´anto afectar´a a la soluci´on un error en el dato b. Para poder dar una respuesta debemos decidir primero c´omo medir el error. Es decir, necesitamos dar alguna forma de medir vectores de IRn. Una forma posible es utilizar la longitud o norma eucl´ıdea del vector, o sea,

‖x‖ 2 =

x^21 + ..... + x^2 n

Pero ´esta no es la ´unica medida razonable y en muchos casos es conveniente trabajar con otras. Por ejemplo, podemos decir que un vector es “chico” si lo son todas sus componentes y tomar entonces como medida de x la siguiente, llamada “norma infinito”,

‖x‖∞ = max 1 ≤i≤n

|xi|

Otra elecci´on natural es la “norma uno”,

‖x‖ 1 = |x 1 | + ..... + |xn|

o m´as en general la “norma p”,

‖x‖p = (|x 1 |p^ + ..... + |xn|p)

1 p

con 1 ≤ p < ∞.

Todas estas formas de medir resultan equivalentes en el sentido de que, si x es “chico” en una de las normas entonces lo es en cualquier otra, puesto que una norma mayora a la otra salvo una constante que depende s´olo de n. Por ejemplo, utilizando la desigualdad de Schwartz se obtiene

‖x‖ 1 ≤

n‖x‖ 2

y por otra parte, es f´acil ver que,

‖x‖ 2 ≤ ‖x‖ 1

2. NORMAS Y CONDICIONAMIENTO 15

‖ym‖ 2 ≥ m‖ym‖

y llamando xm = ym/‖ym‖ 2 obtenemos ‖xm‖ 2 = 1 y

‖xm‖ ≤

m

pero, toda sucesi´on acotada en la norma eucl´ıdea tiene una subsucesi´on convergente. Entonces existe una subsucesi´on de (xm), (x′ m) tal que

‖x′ m − x‖ 2 → 0

para cierto x ∈ IRn. Pero entonces por (2.1), tambi´en vale que

‖x′ m − x‖ → 0

Por otra parte, por (2.3) tenemos que ‖x′ m‖ → 0 y en consecuencia, por unicidad del l´ımite, resulta x = 0. Pero observemos finalmente que, por la desigualdad triangular,

∣∣ ∣‖x′ m‖ 2 − ‖x‖ 2

∣ ≤ ‖x′ m −^ x‖ 2

y entonces se llega a la contradicci´on 1 = ‖x′ m‖ 2 → ‖x‖ 2 = 0, finalizando la demostraci´on.

Ahora s´ı, estamos en condiciones de abordar el problema de c´omo afecta el error en los datos a la soluci´on de un sistema lineal cuya matriz A es inversible. Si se reemplaza el dato b por b + ∆b, la soluci´on x del sistema ser´a modificada de tal forma que tendremos

A(x + ∆x) = (b + ∆b)

o equivalentemente, A∆x = ∆b

y nos preguntamos que relaci´on hay entre

‖∆x‖ ‖x‖

y

‖∆b‖ ‖b‖ Veamos primero el siguiente ejemplo simple,

(

4. 1 2. 8
5. 7 6. 6

x 1 x 2

La soluci´on es

x =

16 2. NORMAS Y CONDICIONAMIENTO

Observemos que

‖b‖ 1 = 13. 8 ‖x‖ 1 = 1

Si modificamos b poniendo

b′^ = b + ∆b =

entonces la soluci´on es

x′^ = x + ∆x =

y se obtiene en consecuencia

‖∆b‖ 1 = 0. 01 ‖∆x‖ 1 = 1. 63

con lo que el error relativo se amplific´o mucho, en efecto,

‖∆b‖ 1 ‖b‖ 1

= 0. 7246 × 10 −^3

mientras que

‖∆x‖ 1 ‖x‖ 1

Nuestro objetivo es tratar de entender a qu´e se debe este comportamiento y poder predecir, dada una matriz A, cu´al ser´a el factor de amplificaci´on del error relativo o, al menos, dar una cota de ´este en t´erminos de A.

Analicemos primero el caso de una matriz diagonal.

( λ 1 0 0 λ 2

x 1 x 2

b 1 b 2

La soluci´on es (si λ 1 , λ 2 6 = 0)

x 1 =

b 1 λ 1

x 2 =

b 2 λ 2 Por ejemplo, si

A =

entonces,

x 1 =

b 1 1000

x 2 = 100b 2