Apuntes Programación 2 Escom Ipn

´ INSTITUTO POLITECNICO NACIONAL Escuela Superior de C´ omputo

Apuntes de Programaci´ on II

autor: M en C. Benjam´ın Luna Benoso

M´exico, D. F.

Diciembre 2008

ii

Introducci´ on Este libro est´a estructurado tomando como referencia el temario para el curso de Programaci´on II impartido en la Escuela Superior de C´omputo del Instituto Polit´ecnico Nacional. Los temas analizados difieren s´olo un poco en el nombre y en el orden en que son presentados en el temario correspondiente. El estudio del an´alisis num´erico se basa en buscar alternativas para resolver problemas cuya soluci´on es dif´ıcil usando m´etodos directos, para ello, el an´alisis num´erico se basa en la alternativa de m´etodos num´ericos. En varias ´areas de la ingenier´ıa y las ciencias, se busca resolver problemas utilizando herramientas del ´algebra lineal o de las ecuaciones diferenciales; estos problemas, por ejemplo, suponen que se trabajan sobre funciones continuas, que para efectos pr´acticos de ejemplificaci´on, se obtiene lo que se desea, sin embargo, en la vida real, varios problemas de los que se plantean, no se conocen muy bien su naturaleza, generalmente se conocen o se obtienen mediciones por medio de aparatos y estos son registrados. Por ejemplo, un f´ısico experimental, en un laboratorio registra datos de un experimento el cual repite varias veces, registra sus datos y ahora lo que le interesa conocer, es una funci´on continua de tal manera, que si se restringe a los puntos que ´el registro, exactamente se comporte como en el experimento realizado. De est´a manera, pasamos de los problemas continuos a los totalmente discretos, que es como se obtienen mediciones y a partir de estas mediciones discretas, pasamos nuevamente al espacio continuo, de ah´ı la importancia del an´alisis num´erico. Otra problem´atica interesante donde se puede notar la importancia del an´alisis num´erico es la siguiente: en un curso de C´alculo, se hace notar √ que la recta esta compuesto por una umeros irracionales) no tienen una infinidad de n´ umeros, y que n´ umeros como π, e y 2 (n´ representaci´on finita en su notaci´on decimal, esto significa que a ninguna persona le bastar´ıa su vida completa para escribir en notaci´on decimal alguno de estos n´ umeros, entonces, como lo hace una computadora. La parte que nos interesa de una computadora, es aquella donde almacena la informaci´on para hacer c´alculos, para representar n´ umeros, sumas, etc., la parte que nos interesa es la memoria. Si finalmente una computadora est´a compuesta en su totalidad de circuitos, con una cantidad finita de estos, entonces como le hace para almacenar un n´ umero irracional como los presentados si su representaci´on es infinita y la computadora s´olo puede guardar n´ umeros binarios (ceros y unos). En realidad una computadora, no es capaz de almacenar todos los n´ umeros reales tal y como se conocen en su abstracci´on. Si una computadora requiere almacenar un n´ umero irracional como π, este lo hace utilizando alg´ un m´etodo de truncamiento que se ver´a posteriormente y en realidad la computadora no estar´a almacenando el n´ umero irracional π, si no mas bien, una aproximaci´on a dicho valor.

iii

Si una computadora esta realizando c´alculos num´ericos donde ha utilizado valores aproximados a n´ umeros reales, esto podr´ıa generar una serie de errores que vendr´ıan generandose en cada operaci´on. Obteniendo finalmente una aproximaci´on al valor real de la operaci´on que se deseaba conocer. Una problem´atica importante en este aspecto es el error que se va obteniendo y saber decidir si es aceptable o no. Esto es algo que tambi´en es tratado en el an´alisi num´erico. Los ejemplos anteriores, muestran la importancia que tiene el an´alisis num´erico dentro de las areas de la ciencias y la ingenier´ıa.

iv

Contenido Introducci´ on 1 Aritm´ etica de Punto Flotante 1.1 N´ umeros de Punto Flotante . . . . . . . . 1.1.1 Sistemas de N´ umeraci´on . . . . . . 1.1.2 Representaci´on Binaria de N´ umeros 1.1.3 Norma IEEE 754 . . . . . . . . . . 1.2 Errores . . . . . . . . . . . . . . . . . . . . 1.2.1 Error Relativo y Error Absoluto . . 1.2.2 Aritm´etica de Punto Flotante . . . 1.2.3 Errores de Representaci´on . . . . . 1.2.4 Reducci´on de Errores . . . . . . . . 1.2.5 Propagaci´on de Errores . . . . . . . 1.2.6 N´ umero de Condicion . . . . . . . .

ii

. . . . . . . . . . . . Negativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 Algoritmos para dar Soluci´ on a la Ecuaci´ on de Una 2.1 Algoritmos Para B´ usqueda de Ra´ıces . . . . . . . . . 2.1.1 Algoritmos . . . . . . . . . . . . . . . . . . . . 2.1.2 M´etodo de la Bisecci´on . . . . . . . . . . . . . 2.1.3 Iteraci´on de Punto Fijo . . . . . . . . . . . . . 2.1.4 M´etodo de Newton-Raphson y de la Secante . 2.2 Ra´ıces de Polinomios . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

1 1 1 7 9 14 14 16 16 20 23 24

Variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

25 25 25 26 28 30 33

. . . . . . .

37 37 38 40 41 44 45 47

. . . . . . . . . . .

. . . . . . . . . . .

3 Algoritmos de Interpolaci´ on 3.1 Algoritmos de Interpolaci´on . . . . . . . . . . . . . . . . 3.1.1 Polinomio de Lagrange . . . . . . . . . . . . . . . 3.1.2 M´etodo de Neville . . . . . . . . . . . . . . . . . 3.1.3 Diferencias Divididas de Newton . . . . . . . . . 3.2 Algoritmos de Ajuste de Curvas . . . . . . . . . . . . . . 3.2.1 Ajuste por M´ınimos Cuadrados a una Recta . . . 3.2.2 Ajuste por M´ınimos Cuadrados a un Polinomio de

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . Grado

. . . . . . n

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

CONTENIDO

v

4 Algoritmos para Derivaci´ on e Integraci´ on Num´ erica 4.1 Algoritmos para Derivaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Algoritmos de Integraci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49 49 51

5 Algoritmos Num´ ericos Utilizados en Algebra Lineal 5.1 Eliminaci´on Gaussiana con Sustituci´on hacia Atr´as . 5.2 Algoritmos Basados en T´ecnicas Iterativas . . . . . . 5.2.1 Normas de Vectores y Matrices . . . . . . . . 5.2.2 Algoritmo de Jacobi . . . . . . . . . . . . . . 5.2.3 Algoritmo de Gauss-Seidel . . . . . . . . . . .

55 55 58 58 59 61

6 Algoritmos para dar Soluci´ on Num´ erica a narias 6.1 Algoritmos por M´etodos de Taylor . . . . 6.1.1 M´etodo de Euler . . . . . . . . . . 6.1.2 M´etodo de Taylor de Orden N . . . 6.2 Algoritmo de Runge-Kutta . . . . . . . . . 6.2.1 Algoritmo de Runge-Kutta . . . . . Bibliograf´ıa

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

Ecuaciones Diferenciales Ordi. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

63 63 63 64 65 65 67

1

Cap´ıtulo 1 Aritm´ etica de Punto Flotante 1.1 1.1.1

N´ umeros de Punto Flotante Sistemas de N´ umeraci´ on

A diario utilizamos n´ umeros para el conteo de cosas; contamos objetos, personas, dinero, etc., estamos muy familiarizados con los n´ umeros naturales 1, 6, 98.32, 23784.12, 20993882 por mencionar algunos de una infinidad. Sin embargo, que significan exactamente. El n´ umero 97 por ejemplo significa nueve veces diez m´as siete unidades, esto es: 97 = (9 × 10) + 7 haciendo una analog´ıa a lo anterior, 7583 significa 7583 = (7 × 103 ) + (5 × 102 ) + (8 × 101 ) + (3 × 100 ) y 362.15 significa 362.15 = (3 × 102 ) + (6 × 101 ) + (2 × 100 ) + (1 × 10−1 ) + (5 × 10−2 ), en general, si z = xn xn−1 . . . x3 x2 x1 x0 .x−1 x−2 . . ., se tiene: X z= xi × 10i (1.1.1) i

A los n´ umeros representables en el modo anterior, se le conoce como sistema decimal. En general, para cada n´ umero natural b mayor que uno, se tiene un sistema de numeraci´on cuya base es b. Sin embargo, en la pr´actica, a parte del sistema decimal, suelen utilizarse los sistemas binario, octal y hexadecimal. El sistema binario o base 2 se representa por medio de dos d´ıgitos: 0 y 1. Para evitar confusiones, 02 y 12 representar´an el cero y uno en base binaria, y en general, el n´ umero zb representar´a el n´ umero z en base b. Se tiene:

1.1. N´ umeros de Punto Flotante

2

02 = 010 12 = 110 De manera an´aloga a lo que se hizo con la base decimal, se tiene lo siguiente: 1012 = (1 × 22 ) + (0 × 21 ) + (1 × 20 ) = 510 1101012 = (1 × 25 ) + (1 × 24 ) + (0 × 23 ) + (1 × 22 ) + (0 × 21 ) + (1 × 20 ) = 5310 110.012 = (1 × 22 ) + (1 × 21 ) + (0 × 20 ) + (0 × 2−1 ) + (1 × 2−2 ) = 6.2510 , en general, si z2 = xm . . . x2 x1 x0 .x−1 x−2 . . . x−n representa un n´ umero en base binaria, entonces, su representaci´on en base decimal queda como sigue: i=−n X

(xi × 2i ))10

z2 = (

(1.1.2)

i=m

La ecuaci´on (1.1.1) muestra la forma para convertir un n´ umero de base binaria a base decimal. Ahora, analicemos el caso contrario, es decir, la conversi´on de un n´ umero de base decimal a base binaria. Para convertir el n´ umero entero positivo N = N10 a base binaria se hace lo siguiente: Paso 1: Dividimos N entre 2, de aqu´ı obtenemos un cociente Q1 y un residuo R1 , con R1 = 0 o 1. Entonces se tiene: N = 2 × Q1 + R1 Paso 2: Enseguida dividimos el cociente Q1 por 2, obteniendo un cociente Q2 y un residuo R2 , con R2 = 0 o 1. Se tiene: Q1 = 2 × Q2 + R2 Paso i: Se sigue el proceso anterior, de tal manera que en el paso i se tiene: Qi−1 = 2 × Qi + Ri Paso n+1: Se sigue el proceso anterior hasta que el residuo del cociente en el proceso continuo de hacer divisiones sea cero, esto es, hasta que en el proceso se tenga un cociente igual a cero Qn = 2 × Qn+1 + Rn+1 = 2 × (0) + Rn+1 = Rn+1 para alg´ un n + 1, de lo anterior, se tiene:

1.1. N´ umeros de Punto Flotante

3

N = 2 × Q1 + R1 = 2 × (2 × Q2 + R2 ) + R1 = 22 × Q2 + 2 × R2 + R1 = 22 × (2 × Q3 + R3 ) + 21 × R2 + 20 × R1 = 23 × Q3 + 22 × R3 + 21 × R2 + 20 × R1 ··· = 2i × Qi + · · · + 22 × R3 + 21 × R2 + 20 R1 ··· = 2n × Qn + · · · + 22 × R3 + 21 × R2 + 20 R1 = 2n × Rn+1 + · · · + 22 × R3 + 21 × R2 + 20 R1 , en conclusi´on, Paso n+2: N10 = (Rn+1 · · · R3 R2 R1 )2 Ejemplo 1.1.1 Siguiendo la metodolog´ıa anterior, para convertir el n´ umero 2710 a base binaria se tiene lo siguiente: 27 = 2 × 13 + 1, Q1 = 13 y R1 = 1, 13 = 2 × 6 + 1, Q2 = 6 y R2 = 1, 6 = 2 × 3 + 0, Q3 = 2 y R3 = 0, 3 = 2 × 1 + 1, Q4 = 1 y R4 = 1, 1 = 2 × 0 + 1, Q5 = 0 y R5 = 1, luego, 2710 = (R5 R4 R3 R2 R1 )2 = 110112 X Para convertir un n´ umero decimal fraccionario a base binaria, se hacen repetidas multiplicaciones por 2, como se ilustra en ejemplo 1.1.2. En cada caso se multiplica por 2 la parte fraccionaria del n´ umero decimal. El d´ıgito del producto a la izquierda del punto decimal ser´a 1 o 0 y contribuye a la representaci´on binaria, comenzando con el bit m´as significativo 1 . La parte fraccionaria del producto se utiliza como en el siguiente paso. Para mostrar lo anterior, consideremos el siguiente ejemplo. Ejemplo 1.1.2 Para convertir el n´ umero 0.3910 a base binaria, se tiene el siguiente procedimiento: 0.39 × 2 = 0.78, parte 0.78 × 2 = 1.56, parte 0.56 × 2 = 1.12, parte 0.12 × 2 = 0.24, parte 0.24 × 2 = 0.48, parte 0.48 × 2 = 0.96, parte 0.96 × 2 = 1.92, parte ... 1

entera entera entera entera entera entera entera

bit que se encuentra a la izquierda de la representaci´on binaria

= = = = = = =

0, 1, 1, 0, 0, 0, 1,

1.1. N´ umeros de Punto Flotante

4

luego 0.3910 = 0.01100012 (aproximado) X Ejemplo 1.1.3 Para convertir el n´ umero 14.2510 a base binaria primeramente, se convierte la parte entera a base binaria quedando de la siguiente manera: 14 = 2 × 7 + 0, 7 = 2 × 3 + 1, 3 = 2 × 1 + 1, 1 = 2 × 0 + 1,

Q1 Q2 Q3 Q4

=7 =3 =1 =0

y y y y

R1 R2 R3 R4

= 0, = 1, = 1, = 1,

por tanto 1410 = (R4 R3 R2 R1 )2 = 11102 Convertida la parte entera a base binaria, procedemos a convertir la parte fraccionaria a base binaria obteni´endose lo siguiente: 0.25 × 2 = 0.5, parte entera = 0 0.5 × 2 = 1.0, parte entera = 1 luego 0.2510 = 0.102 . Por tanto 14.2510 = 1410 + 0.2510 = 11102 + 0.102 = 1110.102 X Los ejemplos anteriores muestran la conversi´on de un n´ umero entre la base decimal y la base binaria. Sin embargo el procedimiento es el mismo para convertir un n´ umero de base 10 a base b. Generalicemos entonces el concepto anterior. Supongamos que se tiene el n´ umero (xm . . . x1 x0 .x−1 x−2 x−3 . . . x−n )10 y se desea conocer su representaci´on en base b. Para ello, primero observemos que la base b est´a compuesta por b s´ımbolos, siendo estos: 0, 1, . . . , b − 1. Paso 1. Se considera la parte entera, esto es P Entera = (xm . . . x1 x0 )10 y se divide entre b, de aqu´ı obtenemos un cociente Q1 y un residuo R1 , con R1 ∈ {0, 1, . . . , b − 1}. Entonces se tiene: (xm . . . x1 x0 )a = b × Q1 + R1 Paso 2: Enseguida dividimos el cociente Q1 por b, obteni´endose un cociente Q2 y un residuo R2 , con R2 ∈ {0, 1, . . . , b − 1}. Se tiene Q1 = b × Q2 + R2 Paso i: Se sigue el proceso anterior, de tal manera que en el paso i se tiene: Qi−1 = b × Qi + Ri Paso n+1: Se sigue el proceso anterior, hasta que el residuo del cociente en el proceso continuo de hacer divisiones sea cero, esto es hasta que en el proceso se tenga un cociente igual a cero

1.1. N´ umeros de Punto Flotante

5

Qn = b × Qn+1 + Rn+1 = b × (0) + Rn+1 = Rn+1 para alg´ un n + 1, de lo anterior, se tiene: P Entera = b × Q1 + R1 = b × (b × Q2 + R2 ) + R1 = b2 × Q2 + b × R2 + R1 = b2 × (b × Q3 + R3 ) + b1 × R2 + b0 × R1 = b3 × Q3 + b2 × R3 + b1 × R2 + b0 × R1 ··· i 2 = b × Qi + · · · + b × R3 + b1 × R2 + b0 R1 ··· n 2 = b × Qn + · · · + b × R3 + b1 × R2 + b0 R1 = bn × Rn+1 + · · · + b2 × R3 + b1 × R2 + b0 R1 , en conclusi´on, Paso n+2: P Entera = (Rn+1 · · · R3 R2 R1 )b Enseguida se considera la parte fraccionaria, esto es: P F racc = (0.x−1 x−2 x−3 . . . x−n )10 , Paso n+3: Se multiplica por b obteniendo una parte entera p1 y una parte fraccionaria f1 , esto es: P F racc × b = p1 .f1 Paso n+4: Se multiplica 0.f1 por b obteniendo una parte entera p2 y una parte fraccionaria f2 , esto es: 0.f1 × b = p2 .f2 Paso n+5: Se multiplica 0.f2 por b obteniendo una parte entera p3 y una parte fraccionaria f3 , esto es: 0.f2 × b = p3 .f3 , se realiza el proceso anterior repetidas veces tantas como sea necesario, quedando la parte fraccionaria de la siguiente manera P F racc = (0.p1 p2 p3 . . . pk )b . Finalmente se junta la parte entera y la parte fraccionaria, quedando la representaci´on del n´ umero de base decimal a binaria de la siguiente manera:

1.1. N´ umeros de Punto Flotante

6

(xm . . . x1 x0 .x−1 x−2 x−3 . . . x−n )10 = (Rn+1 R3 R2 R1 .p1 p2 p3 . . . pk )b Ahora analicemos el caso rec´ıproco, supongamos que se tiene el n´ umero Mb = (xm . . . x1 x0 .x−1 x−2 . . . x−n )b en base b y se desea convertir a base decimal. Esta conversi´on es directa obteni´endose de la siguiente manera: i=−n X

(xi × bi ))10

Mb = (

i=m

= (xm × bm + · · · + x1 × b1 + x0 × b0 + x−1 × b−1 + x−2 × b−2 + · · · + x−n × b−n )10 Los procedimientos anteriores son una generalizaci´on para convertir un n´ umero de su representaci´on de base decimal a base b y viceversa, sin embargo, generalmente se hacen uso u ´nicamente de la base binaria, la base octal y la hexadecimal. La base binaria consta de los d´ıgitos 0 y 1; la base octal de los d´ıgitos 0, 1, 2 , 3, 4, 5, 6, y 7; y la base hexadecimal hace uso de los s´ımbolos 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E y F . Ejemplo 1.1.4 Convertir el n´ umero 3510 a base binaria, octal y hexadecimal. Se tiene: 35 = 2 × 17 + 1, Q1 = 17 y R1 = 1, 17 = 2 × 8 + 1, Q2 = 8 y R2 = 1, 8 = 2 × 4 + 0, Q3 = 4 y R3 = 0, 4 = 2 × 2 + 0, Q4 = 2 y R4 = 0, 2 = 2 × 1 + 0, Q5 = 1 y R5 = 0, 1 = 2 × 0 + 1, Q6 = 0 y R6 = 1, luego 3510 = (R6 R5 R4 R3 R2 R1 )2 = 1000112 . Por otro lado, se tiene: 35 = 8 + 3, Q01 = 4 y R10 = 3, 4 = 8 + 4, Q02 = 0 y R20 = 4, luego 3510 = (R20 R10 )8 = 438 . Se tiene adem´as: 35 = 16 × 2 + 3, Q001 = 2 y R100 = 3, 2 = 16 × 0 + 2, Q002 = 0 y R200 = 2,

1.1. N´ umeros de Punto Flotante

7

luego 3510 = (R200 R100 )16 = 2316 . X Del problema anterior, se tiene: 3510 = 1000112 = 438 = 2316 . El procedimiento para encontrar el n´ umero decimal a octal y hexadecimal es el mostrado en la secci´on, sin embargo, observemos que: 4

3

8 8 z}|{ z}|{ 3510 = 1000112 = ( 100 011 )2 = 438 y 216 316 z}|{ z}|{ 3510 = 1000112 = (0010 0011)2 = 238

Ejemplo 1.1.5 Para convertir el n´ umero 2378 a base decimal se hace lo siguiente: 2378 = 2 × 82 + 3 × 81 + 7 × 80 = 15910 , por tanto 2378 = 15910 . X Ejemplo 1.1.6 Para convertir el n´ umero A2D16 a base decimal se hace lo siguiente: A2D16 = A × 162 + 2 × 161 + D×0 = 10 × 162 + 2 × 161 + 13×0 = 260510 , por lo tanto A2D16 = 260510 . X Cuando no haya confusi´on, se omitir´a colocar el sub´ındice para representar un n´ umero en base decimal.

1.1.2

Representaci´ on Binaria de N´ umeros Negativos

La secci´on anterior muestra el procedimiento para representar un n´ umero de base decimal a base b y viceversa. Sin embargo, u ´nicamente se muestra dicha conversi´on para n´ umeros positivos. La problem´atica a resolver a continuaci´on es la representaci´on binaria de enteros negativos. Para ello, se hacen uso principalmente de tres representaciones: • Signo-magnitud • Complemento a uno • Complemento a dos

1.1. N´ umeros de Punto Flotante

8

Representaci´ on Signo-Magnitud Las cantidades se representan de izquierda a derecha (el bit m´as significativo a la izquierda y el menos significativo a la derecha) reservando el bit m´as significativo para representar el signo, siendo 0 para valores positivos y 1 para valores negativos. Por ejemplo, en una palabra de 8 bits, la representaci´on de 61 y -61 son: +61 = 0011 1101SM −61 = 1011 1101SM Sin embargo, observemos que 0000 0000SM = +0 = −0 = 1000 0000SM ; esto es, el cero tiene dos representaciones. Si consideramos una palabra de 8 bits para representar un n´ umero en signo-magnitud, un Byte puede representar n´ umeros en el rango: −127 ≤ n ≤ 127, siendo 127 = 0111 1111SM el entero positivo m´as grande y −127 = 1000 0001SM el entero negativo m´as grande.

Representaci´ on Complemento a Uno Para enteros positivos el n´ umero se representa como en la representaci´on signo-magnitud con el bit m´as significativo 0. Para enteros negativos, se hace uso del complemento a uno, esto es, primero se considera el entero positivo y despu´es se cambian los bits 1 por 0 y los bits 0 por 1. En una palabra de 8 bits, los n´ umeros +61 y −61 quedan representados como: +61 = 0011 1101C1 −61 = 1100 0010C1 Observemos que 0000 0000C1 = +0 = −0 = 1111 1111C1 ; esto es, el 0 tiene dos representaciones nuevamente.

1.1. N´ umeros de Punto Flotante

9

Representaci´ on Complemento a Dos Los enteros positivos se representan de igual manera a la representaci´on signo-magnitud y complemento uno, con el bit m´as significativo igual a 0. Para enteros negativos se hace uso del complemento a dos, esto es, primero se considera su representaci´on positiva, se aplica el complemento a uno y despu´es se le suma un bit. El sistema sigue reservando el bit m´as significativo para el signo, siendo 0 para entero positivos y 1 para negativos. En una palabra de 8 bits, el n´ umero +61 queda representado como: +61 = 0011 1101C2 y para representar el n´ umero −61 se realiza lo siguiente: 1. Se considera su representaci´on positiva, esto es: +61 = 0011 1101 2. Se aplica complementa a uno: 1100 0010 1100 0010 1 3. Se suma un bit: + 1100 0011 por lo tanto −61 = 1100 0011C2 . Observemos que +0 = 0000 0000C2 = −0. Observaci´ on: Algunos compiladores, como el Borland C++ utilizan tanto para palabras de 16 y 32 bits la representaci´on complemento a dos para representar de manera interna los n´ umeros enteros.

1.1.3

Norma IEEE 754

En el a´rea de ciencias e ingenier´ıa es usual utilizar n´ umeros muy grandes y n´ umeros muy peque˜ nos. La norma IEEE 754 fue desarrollada en 1985 por el Insitute for Electrical and Electronic Engineers, IEEE (Instituto para Ingenieros El´ectricos y Electr´onicos) y sirve para representar en un computador los n´ umeros de punto flotante. Esta norma fue desarrollada con la finalidad para facilitar la portabilidad de programas num´ericos. Para facilitar su comprensi´on del formato IEEE 754, primero se considera un n´ umero de punto flotante en su notaci´on cient´ıfica, esto es, en la que un n´ umero es representado mediante dos cantidades, la mantisa y la caracter´ıstica o exponente. De esta manera, cualquier n´ umero puede ser representado en la forma m ∗ be , siendo m la mantisa y e la caracter´ıstica.

1.1. N´ umeros de Punto Flotante

10

Por ejemplo, el n´ umero 523.789 es lo mismo que los numeros 0.523789 ∗103 , 5237.9 ∗ 10−1 , 5.23789 ∗ 102 o que 523789 ∗ 10−3 . Con la notaci´on cient´ıfica se pueden tener varias formas de representar un n´ umero. Sin embargo la representaci´on de la que generalmente hacen uso los computadores es la notaci´ on cient´ıfica normalizada, que es aquella donde el primer d´ıgito fraccionario es diferente de cero, por ejemplo para el n´ umero 523.789, su notaci´on cient´ıfica normalizada es 0.523789 ∗ 103 . Si la base es 2, significa que el primer bit fraccionario despu´es del punto es 1, la norma IEEE 754 aprovecha este hecho de tal manera que al almacenar un n´ umero real, se emplea la siguiente notaci´on cient´ıfica: x = ±(1.b1 b2 b3 . . .) × 2m con bi = 0 o 1 para cada i

(1.1.3)

La norma IEEE 754 define dos formatos b´asicos para almacenar n´ umeros de punto flotante: un formato simple de 32 bits y un formato doble de 64 bits. La base ´ımplicita es 2. Formato Simple Doble

bits para Signo bits para la Caracter´ıstica 1 8 1 11

bits para Mantisa 23 52

Tabla 1.1.1: Formato IEEE

El bit m´as significativo contiene el signo del n´ umero, siendo 0 para valores positivos y 1 para valores negativos. El valor de la caracter´ıstica se almacena en el siguiente bloque de bits, siendo 8 para el formato simple y 11 para el formato doble. La representaci´on de la caracter´ıstica se conoce como el termino sesgado. Un valor fijo llamado sesgo se resta de este campo para conseguir el valor de la caracter´ıstica verdadera. El sesgo tiene una valor de 2k−1 − 1 con k el n´ umero de bits en el exponente binario. As´ı, para el formato simple, el sesgo tiene un valor de 28−1 − 1 = 127 de tal manera que los valores verdaderos var´ıan entre -127 y 128; para el formato doble, el sesgo tiene un valor de 211−1 − 1 = 1023 de tal manera que los valores verdaderos var´ıan entre -1023 y 1024. En los u ´ltimos 23 bits para el formato simple y 52 para el formato doble se almacena la mantisa. Para almacenar un n´ umero de punto flotante se siguen los siguientes pasos: 1. El n´ umero binario se escribe como en la notaci´on cient´ıfica 1.1.3 2. El signo se almacena en el bit m´as significativo: 0 para valores positivos y 1 para negativos. 3. Se suma el valor del sesgo al exponente original y se almacena en el bloque de bits correspondiente a la caracter´ıstica.

1.1. N´ umeros de Punto Flotante

11

4. Se trunca la mantisa al tama˜ no de la palabra para la mantisa. Ejemplo 1.1.7 Para convertir el n´ umero 20 a formato simple IEEE 754 se siguen la serie de pasos siguientes: 1. Se convierte en n´ umero a base binaria y se escribe en la forma cient´ıfica 1.1.3; esto es: 20 = 101002 = 1.0100100 2 . 2. Puesto que el n´ umero a almacenar es positivo, entonces el bit m´as significativo tendr´a un valor de 0. 3. Se suma el sesgo (127) al valor del exponente (1002 = 4) y se almacena en el siguiente bloque de 8 bits; esto es: 127 + 4 = 131 = 100000112 , as´ı, en el bloque de la caracter´ıstica se almacenan los bits 10000011. 4. En el bloque de la mantisa se almacenan los bits 01000000000000000000000 finalmente el n´ umero 20 en formato simple IEEE 754 es 0 10000011 01000000000000000000000. X Ejemplo 1.1.8 Convertir el n´ umero -125.32 al formato simple IEEE 754. Se tiene: 1. Se convierte el n´ umero a base binaria y se escribe en la forma cient´ıfica 1.1.3; esto es: −125.32 = −1111101.0101000111101012 = −1.11110101010001111010111110 2 . 2. Puesto que el n´ umero a almacenar es negativo, entonces el bit m´as significativo tendr´a un valor de 1. 3. Se suma el sesgo (127) al valor del exponente (1102 = 6) y se almacena en el siguiente bloque de 8 bits; esto es: 127 + 6 = 133 = 100001012 , as´ı, en el bloque de la caracter´ıstica se almacenan los bits 10000101. 4. En el bloque de la mantisa se almacenan los bits 11110101010001111010111 finalmente el n´ umero -125.32 en formato simple IEEE 754 es 1 10000101 11110101010001111010111. X Para convertir un n´ umero de formato IEEE 754 a su representaci´on binaria se hace la siguiente operaci´on: (−1)bit mas signif icativo × 2caracteristica × 2−sesgo × 1.mantisa

(1.1.4)

1.1. N´ umeros de Punto Flotante

12

Ejemplo 1.1.9 Para convertir el n´ umero 0 10000011 01000000000000000000000 de formato simple IEEE 754 a su representaci´on binaria se hace la operaci´on 1.1.4, de donde se tiene lo siguiente: sesgo=127

z }| { (−1)0 × 210000011 × 2− 1111111 × 1.01000000000000000000000 = 2100 × 1.01000000000000000000000 = 1.01000000000000000000000 × 2100 = 10100 = 2010 obteni´endose exactamente el n´ umero exacto 20 que as´ı era de esperar ya que este ejercicio muestra la parte inversa al ejercicio 1.1.7 X La norma IEEE 754 reserva ciertos valores para casos especiales como lo son ±0 y ±∞ y N aN 2 (Not a Number ). Es necesario considerar este tipo de valores para prevenir errores que se puedan presentar. Las siguientes tablas muestran toda la gamma de interpretaci´on de los n´ umeros en el formato IEEE 754 para el formato simple y para el formato doble. Precisi´on Simple Signo Exponente Sesgado Cero Positivo 0 0 Cero Negativo 1 0 M´as Infinito 0 255 (todos unos) Menos Infinito 1 255 (todos unos) NaN Silencioso 0o1 255 (todos unos) NaN Indicador 0o1 255 (todos unos) Positivo Normalizado 6= 0 0 0 < e < 255 Negativo Normalizado 6= 0 1 0 < e < 255 Positivo Denormalizado 0 0 Negativo Denormalizado 1 0

Parte Fraccionaria 0 0 0 0 6= 0 6= 0 f f f 6= 0 f 6= 0

Tabla 1.1.2: Interpretaci´on de n´ umeros en formato simple IEEE 754 2

Para una mayor referencia a los NaN pueden ver [5]

Valor 0 -0 ∞ −∞ NaN NaN (1.f ) ∗ 2e−127 −(1.f ) ∗ 2e−127 (0.f ) ∗ 2e−126 −(0.f ) ∗ 2e−126

1.1. N´ umeros de Punto Flotante Precisi´on Doble Signo Exponente Sesgado Cero Positivo 0 0 Cero Negativo 1 0 M´as Infinito 0 2047 (todos unos) Menos Infinito 1 2047 (todos unos) NaN Silencioso 0 o 1 2047 (todos unos) NaN Indicador 0 o 1 2047 (todos unos) Positivo Normalizado 6= 0 0 0 < e < 2047 Negativo Normalizado 6= 0 1 0 < e < 2047 Positivo Denormalizado 0 0 Negativo Denormalizado 1 0

13

Parte Fraccionaria 0 0 0 0 6= 0 6= 0 f f f 6= 0 f 6= 0

Valor 0 -0 ∞ −∞ NaN NaN (1.f ) ∗ 2e−1023 −(1.f ) ∗ 2e−1023 (0.f ) ∗ 2e−1022 −(0.f ) ∗ 2e−1022

Tabla 1.1.3: Interpretaci´on de n´ umeros en formato doble IEEE 754

En la norma est´andar IEEE 754, para el caso del formato simple, los exponentes est´an restringidos por la desigualdad: −126 < e − 127 < 127 y la mantisa est´a restringida por la desigualdad: 1 ≤ (1.f )2 ≤ (1.11111111111111111111111)2 = 2 − 2−23 por lo tanto el n´ umero m´as grande representado en el formato simple es: (2 − 2−23 ) ∗ 2127 ≈ 3.4 × 1038 y el n´ umero m´as peque˜ no representable en el formato simple es: 2−126 ≈ 1.2 × 10−38 Ejercicio Encontrar el n´ umero m´as grande y m´as peque˜ no representable en el formato doble IEEE 754.

1.2. Errores

1.2

14

Errores

A lo largo de esta obra, se emplear´an t´ecinas num´ericas que ayudaran a resolver problemas matem´aticos. Mediante estas t´ecnicas se tendr´a una aproximaci´on a la soluci´on real, sin embargo, el hecho de que no sea la soluci´on real, nos genera un margen de error. Este error en la pr´actica debe de ser considerado, pues nos puede acarrear situaciones inesperadas. Los computadores son un ejemplo claro en la generaci´on de errores. La recta real esta compuesta de una infinidad de n´ umeros, sin embargo un computador tiene una cantidad finita de circuiter´ıa, y por ende, no es capaz de almacenar toda la infinidad de n´ umeros reales. Por citar un ejemplo, consideremos el n´ umero irracional π = 3.141592654 . . .. Un n´ umero irracional tiene la caracter´ıstica de que no puede ser expresado en forma fraccionaria de una cantidad finita de d´ıgitos con cierto periodo. Para que un computador maneje el valor de π, debe de hacer uso de alg´ un m´etodo que le permita registrar una cierta cantidad finita de la parte fraccionaria del n´ umero y almacenarla. Esto, por supuesto nos generar´a un cierto error, ya que no se estar´a trabajando con el valor de la cantidad exacta, si no mas bien con una aproximaci´on. De est´a manera, en esta secci´on se dedicar´a a analizar la parte correspondiente a los errores que se producen al hacer uso de una cantidad aproximada a su respectivo valor real.

1.2.1

Error Relativo y Error Absoluto

En lo que sigue, suponemos un computador imaginario que trabaja con el sistema decimal con una palabra de k d´ıgitos en la mantisa. Adem´as suponemos que los n´ umeros est´an representados en la forma normalizada de punto flotante: mantisa

z }| { ±0. d1 d2 · · · dk ×10n , con 1 ≤ d1 ≤ 9 y 0 ≤ di ≤ 9, i = 2, . . . k

(1.2.1)

A los n´ umeros de la ecuaci´on 1.2.1 se les llama n´ umeros de m´aquina decimales con kd´ıgitos. En general, cualquier n´ umero real positivo y puede ser normalizado de la forma: 0.d1 d2 · · · dk dk+1 · · · × 10n , con 1 ≤ d1 ≤ 9 y 0 ≤ di ≤ 9, i = 2, . . . k

(1.2.2)

Si nuestro computador imaginario de k d´ıgitos en la mantisa requiere almacenar un n´ umero cuya forma es como el de la ecuaci´on 1.2.2, entonces requerir´a trabajar u ´nicamente con k d´ıgitos de la mantisa, eliminando el resto, a esto se le conoce como forma de punto flotante. Denotamos mediante f l(y) a la forma de punto flotante de y. El valor de f l(y) se obtiene

1.2. Errores

15

terminando el valor de la mantisa de y en k cifras decimales. Para ello hay dos formas de realizarlo que se muestran en las siguientes definiciones. Definici´ on 1.2.1 Si y es un n´ umero real expresado de la forma como en la ecuaci´on 1.2.2, entonces f lc (y) = 0.d1 d2 · · · dk , este m´etodo es llamado m´etodo de corte o truncamiento a k-d´ıgitos. Definici´ on 1.2.2 Si y es un n´ umero real expresado de la forma como en la ecuaci´on 1.2.2, entonces

f lr (y) =

0.d1 d2 · · · d0k ,

con

d0k

 =

dk si dk+1 < 5 dk + 1 si dk+1 ≥ 5

este m´etodo es llamado m´etodo de redondeo a k-d´ıgitos. √ Ejemplo 1.2.1 3 = 1.732050808 · · · , normalizado es igual a 0.1732050808 × 101 , por √ tanto √ el m´etodo de corte y redondeo a 5-d´ıgitos son respectivamente f lc ( 3) = 0.17320 y f lr ( 3) = 0.17321. X Definici´ on 1.2.3 El error resultante al reemplazar un n´ umero con su forma de punto flotante recibe el nombre de error de redondeo (independientemente del m´etodo que se aplique). La siguiente definici´on describe m´etodos para medir los errores de aproximaci´on. Definici´ on 1.2.4 Si p? es una aproximaci´on de p, el error absoluto es EA = |p − p? |. Si |p−p? | p 6= 0 el error relativo es ER = EA = , y el error relativo porcentual o simplemente |p| p error porcentual es ERP = ER × 100 (%). Ejemplo 1.2.2 Si p = 0.17320 × 101 y p? = 0.17321 × 101 , entonces EA = |0.17320 × 101 − −3 EA 0.17321×101 | = 0.1×10−3 , ER = |0.17320| = 0.1×10 = 0.57736×10−3 y ERP = ER×100 = 0.17320 0.057%. X Definici´ on 1.2.5 Se dice que el n´ umero p? se aproxima a p con t-d´ıgitos o cifras significativas, si t es el entero positivo m´as grande para el cual ER =

|p−p? | |p|

< 5 × 10−t

Ejemplo 1.2.3 Si p = 0.51 × 10−1 y p? = 0.5 × 10−1 , entonces:

1.2. Errores

16

ER =

|0.51×10−1 −0.50×10−1 | |0.51×10−1 |

= 0.19 × 10−1 < 0.5 × 10−1

por tanto p? se aproxima a p con 1 cifra significativa. X El concepto de cifra significativa se ha desarrollado para designar formalmente la confiabilidad de un valor num´erico. Las cifras significativas de un n´ umero son aquellas que pueden utilizarse en forma confiable.

1.2.2

Aritm´ etica de Punto Flotante

Antes de continuar con la parte referente a los errores, es necesario conocer las operaciones aritm´eticas de punto flotante. Estas se resumen enseguida. Si X = Xm ×B Xe y Y = Ym ×B Ye , entonces: X + Y = (Xm × B Xe −Ye + Ym ) × B Ye , con Xe ≤ Ye , X − Y = (Xm × B Xe −Ye − Ym ) × B Ye , con Xe ≤ Ye , X × Y = (Xm × Ym ) × B Xe +Ye y X = XYmm × B Xe −Ye Y Ejemplo 1.2.4 Si X = 0.5 × 102 = 50 y Y = 0.3 × 103 = 300, entonces X + Y = (0.5 × 102−3 + 0.3) × 103 = 0.35 × 103 = 350, X − Y = (0.5 × 102−3 − 0.3) × 103 = −0.25 × 103 = −250, X × Y = (0.5 × 0.3) × 102+3 = 0.15 × 105 = 15000, X 0.5 = 0.3 × 102−3 Y = 1.66 × 10−1 = 0.166 X

1.2.3

Errores de Representaci´ on

Para mostrar algunos de los errores que se pueden presentar en un computador, imaginemos que nuestro computador, adem´as de trabajar con el sistema decimal, supongamos que los n´ umeros se almacenan con una mantisa de 4 d´ıgitos decimales, una caracter´ıstica de 2 d´ıgitos decimales, de los cuales, el primero se usa para almacenar el signo.

1.2. Errores

17 Signo caracter´ıstica ± 1 d´ıgito

mantisa 4 d´ıgitios

Tabla 1.2.1: Almacenamiento del computador imaginario

Por ejemplo, si deseamos almacenar los n´ umeros 8123000 y 0.00456 en nuestro computador imaginario, primero se deben de normalizar los n´ umeros quedando 0.8123×107 y 0.4560×10−2 respectivamente, por tanto los n´ umeros almacenados en nuestro computador quedar´ıan de la siguiente manera: Signo caracter´ıstica + 7 − 2

mantisa 8123 4560

Con ayuda de nuestro computador imaginario, se pueden mostrar ciertos errores que se cometen.

a) Suma de N´ umeros muy Distintos en Magnitud Sean X = 0.003 y Y = 700, que sucede si deseamos almacenar el valor de X + Y para usarlo posteriormente. Primeramente para operar los valores X y Y , estos se deben de almacenar en el computador, para esto, primero se deben de normalizar obteniendo: X = 0.3000 × 10−2 Y = 0.7000 × 103 estos n´ umeros no los puede operar directamente el computador ya que tienen distinta caracter´ıstica, para ello, se deben de desnormalizar antes de efectuar la suma. 0.000003 × 103 +0.700000 × 103 0.700003 × 103 por lo tanto, el valor de la suma almacenado en nuestro computador queda as´ı: Signo caracter´ıstica + 3 y por tanto la suma no se realizo.

mantisa 7000

1.2. Errores

18

b) Resta de N´ umeros Casi Iguales Sean X = 0.3153 y Y = 0.3152, que sucede si ahora deseamos almacenar el valor de X − Y . Se tiene: 0.3153 × 100 −0.3152 × 100 0.0001 × 100 ya que la mantisa del valor de la resta esta desnormalizada, el computador primero la normaliza 0.1000 × 10−3 y la almacena de la siguiente manera: Signo caracter´ıstica − 3

mantisa 1000

hasta aqu´ı no ha habido error alguno, sin embargo el valor de la resta s´olo tiene una cifra significativa y no se puede confiar en su exactitud, ya que si este valor es usado para posteriores operaciones puede generar resultados inesperados. Por ejemplo, supongamos que lo siguiente es parte del c´ogido de un programa: X = (A − B) ∗ C con A = 0.3153 × 100 , B = 0.3152 × 100 y C = 0.1000 × 105 . Al operar, se tiene X = 0.00001 × 105 = 1, siendo correcto este valor, sin embargo, supongamos que el valor de A fue calculado anteriormente y el valor fue A? = 0.3154 × 100 , se tiene que el error absoluto, relativo y porcentual entre el valor de A y el valor de A∗ son: EA = |0.3153 − 0.3154| = 0.0001, EA = 0.00031, ER = |0.3153| ERP = ER × 100 = 0.031% Ahora bien, hagamos uso del valor de A? en lugar de A para obtener el valor de X. Se tiene: X ? = (A? − B) ∗ C = (0.3154 × 100 − 0.3152 × 100 ) ∗ 0.1000 × 105 =2 El error absoluto, relativo y porcentual entre el valor de X y X ? son:

1.2. Errores

19

EA = |1 − 2| = 1, ER = EA = 1, |1| ERP = ER × 100 = 100% En conclusi´on, el error del c´alculo de A? en lugar de A nos genero un error porcentual de 0.031%, sin embargo, esto nos genero un error porcentual del 100% en el resultado final del c´alculo de X.

c) Overflow (desbordamiento o sobreflujo) y Underflow (subflujo) Sean X = 0.7000 × 107 y Y = 0.3000 × 109 , se tiene: 0.7000 × 107 ∗0.3000 × 109 0.2100 × 1016 luego, si se desea almacenar el valor de X ∗ Y en nuestro computador imaginario, este requerir´a de 3 d´ıgitos para la caracter´ıstica siendo el primero para el signo, pero nuestro computador s´olo tiene capacidad para dos d´ıgitos siendo el primero para el signo y por tanto, el valor de X ∗ Y no se puede almacenar en nuestro computador y por ende se interrumpen los c´alculos. A este tipo de error se le conoce como overflow, desbordamiento o sobreflujo. Un sobreflujo surge como consecuencia de una operaci´on artim´etica de dos n´ umero v´alidos cuyo resultado es un n´ umero tan grande que el computador no puede manejar. Ahora, consideremos los valores X = 0.6000 × 10−8 y Y = 0.3000 × 10−5 y hagamos la operaci´on producto. Se tiene: 0.6000 × 10−8 ∗0.3000 × 10−5 0.1800 × 10−13 la caracter´ıstica es -13 y nuevamente como en el caso del overflow, este valor no puede almacenarse en nuestro computador, generandose de esta manera lo que se conoce como un error llamado underflow o subflujo. Un subflujo surge como consecuencia de una operaci´on aritm´etica de dos n´ umeros cuyo resultado es un n´ umero tan peque˜ no que el computador no puede manejar. Tanto el sobreflujo como el subflujo pueden aparecer pueden aparecer a partir de operaciones aritm´eticas que involucran productos o cocientes. El subflujo no es un problema tan serio como lo es el sobreflujo. En el caso del sobreflujo los computadores interrumpen los c´alculos y envian un mensaje de error, mientras que en el caso del subflujo, los computadores trabajan el resultado del valor como cero.

1.2. Errores

20

d) Divisi´ on entre un N´ umero muy Peque˜ no Supongamos que se desea operar: X = A − B/C con valores A = 0.1120 × 109 , B = 0.1000 × 106 y C = 0.9000 × 10−3 . Bajo estas condiciones, X = 0.0009 × 109 = 0.9000 × 106 . Supongamos que tras operaciones anteriores al c´alculo de X, se obtuvo el valor de C ? = 0.9001 × 10−3 en lugar de C. El error absoluto, relativo y porcentual entre el valor de C y C ? son: EA = |0.9000 × 10−3 − 0.9001 × 10−3 | = 0.0001 × 10−3 , EA −4 , ER = |0.9000×10 −3 | = 10 ERP = ER × 100 = 0.01% Ahora bien, hagamos uso del valor de C ? en lugar de C para el c´alculo de X. Se tiene: X ∗ = A − B/C ? = 0.1120 × 10 − 0.1000 × 106 /0.9001 × 10−3 = 0.1000 × 107 9

El error absoluto, relativo y porcentual entre el valor de X y X ? son: EA = |0.9000 × 106 − 0.1000 × 107 | = 100000, EA ER = |0.9000×10 6 | = 0.1111, ERP = ER × 100 = 11.11% El error relativo como consecuencia de haber obtenido el valor de C ? se ha multiplicado 1111 veces. Este error surge como consecuencia de realizar una divisi´on donde el cociente se obtuvo de c´alculos anteriores generandose un pequen˜ no error.

1.2.4

Reducci´ on de Errores

En la secci´on anterior se vieron posibles errores que se generan en un computador. Ahora lo que nos interesa es tratar de reducir lo mejor posible los errores generados. Ejemplo 1.2.5 Consideremos la ecuaci´on cuadr´atica x2 − 101.01x + 0.10101 = 0. Para encontrar las ra´ıces, se recurre a la f´ormula general x=

√ −b± b2 −4ac 2a

1.2. Errores

21

sustituyendo los valores en la f´ormula general para nuestro ejemplo y aplicando una aritm´etica de redondeo a cinco d´ıgitos, se tiene: √ 101.01±

∗

(−101.01)2 −4(0.10101)

x = 2 √ 101.01± 10203−0.40404 = √ 2 = 101.01±2 10203 = 101.01±101.01 2 luego x∗+ = x∗− =

101.01+101.01 2 101.01−101.01 2

= 101.01 y =0

Por tanto las ra´ıces son x∗+ = 101.01 y x∗− = 0. Sin embargo las soluciones verdaderas redondeadas a cinco d´ıgitos decimales son x+ = 101.01 y x− = 0.00100 generandose errores relativos porcentuales de:

ERP =

|x+ −x∗+ | × 100 |x+ | |101.01−101.01| × |101.01|

= 100 = 0% y |x −x∗ | ERP = −|x− |− × 100 = |0.00100−0| × 100 |0.001| = 100% Se observa que el m´etodo fue adecuado para la soluci´on x+ , pero no del todo para x− . En este ejemplo se presenta el problema de restar n´ umeros casi iguales, dando como resultado p´erdida de exactitud. Veamos, que sucede, si en lugar de aplicar la f´ormula general tal y como se efecto para encontrar el valor de x− , se recurre a la racionalizaci´on de t´erminos. Se tiene: x− = = = =

√ −b− b2 −4ac √2a √ −b− b2 −4ac −b+√b2 −4ac 2a −b+ b2 −4ac b2 −(b2 −4ac) √ 2a(−b+ b2 −4ac) 2c √ −b+ b2 −4ac

utilizando la expresi´on anterior con a = 1, b = −101.01 y c = 0.10101 se obtiene:

1.2. Errores

22

x− =

√2(0.10101)

101.01+ 101.012 −4(0.10101) 0.20202 101.01+101.01 0.20202 202.02

= = = 0.001

el cual es el valor verdadero utilizando el m´etodo de redondeo a cinco d´ıgitos. Est´a f´ormula alternativa para calcular una ra´ız peque˜ na de una ecuaci´on cuadr´atica casi siempre produce una respuesta m´as exacta que la f´ormula usual. X Ejemplo 1.2.6 Consideremos el polinomio f (x) = x3 − 3.2x2 + 7.1x + 2.5. Evaluemos x = 5.17 en f (x) utilizando una aritm´etica de redondeo d´ıgitos. La siguiente tabla muestra los resultados intermedios para llegar a la soluci´on final.

Exacto redondeo a tres d´ıgitos

x x2 x3 3.2x2 7.1x 5.17 26.7289 138.188413 85.53248 36.707 5.17 26.7 138 85.4 36.7

luego Exacto f (5.17) = 138.188413 − 85.53248 + 36.707 + 2.5 = 91.862933 Redondeo a tres d´ıgitos f (5.17) = 138 − 85.4 + 36.7 + 2.5 = 91.8 obteni´endose un error relativo porcentual de ERP = |91.862933−91.8| × 100 |91.862933| = 0.06% Como m´etodo alternativo, f (x) se puede escribir de la manera: f (x) = ((x − 3.2)x + 7.1)x + 2.5, esto da como resultado que f (5.17) = 91.9, de donde se tiene ERP = 0.04%. Disminuyendo en 0.02% el error. Este m´etodo siempre disminuye el error, en algunos casos da resultados con un porcentaje mucho mejor. Los polinomios siempre deben de anidarse antes de realizar cualquier evaluaci´on, pues de est´a forma se minimiza el n´ umero de c´alculos aritm´eticos y por tanto los errores.

1.2. Errores

1.2.5

23

Propagaci´ on de Errores

El error de propagaci´on sirve para estimar el error resultante en una funci´on a partir de valores de entrada. El objetivo es estimar como se propagan los errores; por ejemplo, si se suman o multiplican dos n´ umeros que tienen error, como se ve reflejado esto al aplicar estos valores a una funci´on. Cuando se eval´ ua una funci´on en un punto x = a, en general se dispone de un valor aproximado a a: a? , teni´endose como error a = a − a? . Nuestro objetivo es determinar el error de propagaci´on de la funci´on; esto es: f = f (a) − f (a∗ )

Figura 1.2.1: Gr´ afica de una funci´on en la cercan´ıa de a

Nuestro objetivo es determinar la relaci´on existente entre a y f . Si 0 < a  1, entonces f (x) puede aproximarse por la tangente al punto x = a, teni´endose: f 0 (a) ≈

f (a)−f (a∗ ) a−a∗

=

f a

1.2. Errores

24

entonces f ≈ a f 0 (a) ≈ a f 0 (a∗ ) y por lo tanto |f | ≈ |a ||f 0 (a∗ )| de esta manera, se define el error de propagaci´on como: |f | = |x ||f 0 (x∗ )| = |x − x∗ ||f 0 (x∗ )| Ejemplo 1.2.7 Si f (x) = x3 + 1 con x∗ = 3.21 y con a = 0.01, entonces el error de propagaci´on es |f | = |0.01||3(3.21)2 +1| = 0.3191. Puesto que f (3.21) = 34.0761, entonces el valor real se encuentra dentro del intervalo [34.0761−0.3191, 34.0761+0.3191] = [33.7570, 34.3952]. X

1.2.6

N´ umero de Condicion

Definici´ on 1.2.6 Sean f y x∗ una aproximaci´on de x, se define el n´ umero condicionado como: NC =

x∗ f 0 (x∗ ) f (x∗ )

El n´ umero condicionado compara el error relativo de x con x∗ y el de f (x) con f (x∗ ). Esto es, mide la sensibilidad de valores de salida de una funci´on debido a perturbaciones peque˜ nas en los valores de entrada. 1. Si N C = 1, entonces el error relativo de f (x) es igual al error relativo de x. 2. Si N C > 1, entonces el error relativo de f (x) es mayor que el error relativo de x. 3. Si N C < 1, entonces el error relativo de f (x) es menor que el error relativo de x. Si N C  1, se dice que la funci´on esta mal condicionada. Ejemplo 1.2.8 Sean f (x) = tan x y x∗ =

π 2

+ 0.02 π2 . Se tiene:

f 0 (x) = sec2 x =

1 cos2 x

luego ∗

2

∗

x ) N C = x (1/cos tan x∗ = 1.6022(1024.9) −31.998 = −51.319

por tanto la funci´on esta mal condicionada. Esto se debe a que en una vecindad de π/2 la tangente tiende tanto a infinito positivo como a infinito negativo.

25

Cap´ıtulo 2 Algoritmos para dar Soluci´ on a la Ecuaci´ on de Una Variable 2.1 2.1.1

Algoritmos Para B´ usqueda de Ra´ıces Algoritmos

A lo largo de los apuntes, se estar´an trabajando procedimientos como los ya vistos en la secci´on 1.1.1, donde se mostraron procedimientos para convertir un n´ umero de base decimal a base binaria y viceversa. Estos procedimientos son llamados algoritmos. Un algoritmo es una serie de pasos que sirven parta resolver un problema. Sus principales caracter´ısticas son: • Contiene valores de entrada. • contiene una serie finita de pasos. • Contiene valores de salida. • No se permiten ciclos infinitos. • Est´a bien definido. A la descripci´on de los algoritmos por medio de texto parecido al c´odigo de computadora se le llama seudoc´odigo. Un algoritmo se compone de valores de Entrada, Proceso(en el que se describe la serie finita de pasos) y Salida; es de est´a forma en que se describir´an los seudoc´odigos de los algoritmos presentados a lo largo de estos apuntes. Ejemplo 2.1.1 Un algoritmo para calcular

2.1. Algoritmos Para B´ usqueda de Ra´ıces N X

26

xii = x11 + x22 + · · · + xN N

i=1

donde N y los valores x1 , x2 , . . . , xN est´an dados y ordenados se describe enseguida: ENTRADA: N, x1 , x2 , . . . , xN . N X SALIDA: SU M A = xii . i=1

PROCESO: Se utiliza una variable auxiliar aux paso 1: Establecer SU M A = 0 y aux = 1. paso 2: Para i = 1 hasta i = N hacer pasos 3 a 6: paso 3: Para j = 1 hasta j = i hacer paso 4: paso 4: aux = aux ∗ xi . paso 5: SU M A = SU M A + aux. paso 6: Poner aux = 1. paso 7: SALIDA (SUMA). PARAR. X Definici´ on 2.1.1 Un algoritmo se dice que es estable si cambios peque˜ nos en los datos de entrada producen cambios peque˜ nos en los datos de salida. En caso contrario se dice que es inestable. Si un algoritmo es estable s´olo para ciertos datos iniciales se dice que es condicionalmente estable. Los algoritmos que nos interesan trabajar son aquellos que sean estables. En este cap´ıtulos, nos dedicaremos a la b´ usqueda de ra´ıces para una funci´on f (x) = 0. Se describir´an algunos m´etodos que nos ayudara´an a aproximar la ra´ız de una funci´on. El primer m´etodo a tratar es el m´etodo de la bisecci´on.

2.1.2

M´ etodo de la Bisecci´ on

Antes de describir de que se trata el m´etodo de la bisecci´on, recordemos el teorema del valor intermedio. Teorema del Valor Intermedio Si f es una funci´on continua en un intervalo [a, b] y y es un n´ umero que se encuentra entre f (a) y f (b), entonces existe c ∈ [a, b] tal que f (c) = y. El m´etodo de la bisecci´on hace uso del teorema del valor intermedio. Se aplica de la siguiente manera: dada una funci´on f continua en un intervalo [a, b] donde de antemano se conoce que existe una sola ra´ız, entonces, se tiene que [f (a) > 0 y f (b) < 0] o [f (a) < 0 y f (b) > 0]. Sin p´erdida de generalidad, supongamos que se tiene la primera situaci´on. Puesto que f (a) > 0 y f (b) < 0 con f funci´on continua, entonces por el teorema del

2.1. Algoritmos Para B´ usqueda de Ra´ıces

27

; p1 es el punto valor intermedio, existe c ∈ [a, b] tal que f (c) = 0. Sea p1 = a + b−a 2 medio del intervalo [a, b]. Por tricotom´ıa de los n´ umeros reales, se tiene p1 < c o p1 = c o p1 > c. Si p1 = c, entonces p1 es la ra´ız que est´abamos buscando. Sin p´erdida de generalidad, supongamos que p1 < c, entonces se tiene que signo(f (p1 )) = signo(f (a)) > 01 y signo(f (b)) < 0 si y s´olo si signo(f (p1 )) ∗ signo(f (b)) < 0, luego por el teorema del valor 1 intermedio existe c0 ∈ [p1 , b] tal que f (c0 ) = 0. Obs´ervese que c0 = c. Tomamos p2 = p1 + b−p ; 2 esto es, p2 es el punto medio del intervalo [p1 , b] y aplicamos nuevamente el razonamiento anterior. Este proceso se hace repetidamente hasta que encerremos la ra´ız en un intervalo suficientemente peque˜ no (figura 2.1.1).

Figura 2.1.1: M´etodo de la bisecci´on

Algoritmo 2.1.1: M´etodo de la Bisecci´on. Entrada: funci´on f ; intervalo [a, b]; n´ umero m´aximo de iteraciones N0 ; tolerancia T OL. Salida: soluci´on aproximada p o mensaje de error. Proceso: se requieren de dos variables auxiliares f a y f p Paso 1: Tomar i = 1 y f a = f (a) paso 2: Para i = 1 hasta i = N0 hacer pasos 3 A 6: paso 3: Tomar p = a + b−a . 2 paso 4: Hacer f p = f (p) paso 5: Si f (p) = 0 o (b − a)/2 < T OL, entonces SALIDA(p). PARAR. paso 6: Si signo(f (a)) ∗ signo(f (p)) > 0, tomar   −1 si x < 0 1 0 si x = 0 la funci´ on signo se define como signo(x) =  1 si x > 0

2.1. Algoritmos Para B´ usqueda de Ra´ıces

28

a = p y f a = f p, en caso contrario, tomar b = p. paso 7: SALIDA(Procedimiento terminado sin ´exito despu´es de N0 iteraciones). PARAR. Ejemplo 2.1.2 La ecuaci´on f (x) = x3 − 7x2 + 14x − 6 = 0 tiene una ra´ız en el intervalo [0, 1] ya que f (0) = −6 < 0 y f (1) = 2 > 0. Para calcular el valor de p3 se tiene: p1 = 0 + 1−0 2 = 0.5, puesto que f (p1 ) = f (0.5) = −0.625 < 0, entonces p2 = 0.5 + = 0.75,

1−0.5 2

puesto que f (p2 ) = f (0.75) = 0.984 > 0, entonces p3 = 0.5 + 0.75−0.5 2 = 0.625. X

2.1.3

Iteraci´ on de Punto Fijo

Definici´ on 2.1.2 Sea g una funci´on. Un punto p en el dominio de g se dice que es punto fijo si g(p) = p. El m´etodo de punto fijo se basa en la definici´on de punto fijo. Sea f una funci´on, definamos g(x) = x−f (x). Si p es un punto fijo de g, entonces f (p) = g(p)−p = 0, esto es, un punto fijo para g es una ra´ız para f . De esta manera, si se desea encontrar las ra´ıces de una funci´on f , basta con encontrar los puntos fijos de una funci´on g construida a partir de f de tal manera que los puntos fijos de g sean ra´ıces de f . El ejemplo 2.1.3 muestra posibles elecciones de la funci´on g para una funci´on f dada. Ejemplo 2.1.3 Si f (x) = x3 + 2x2 + 5x + 1, entonces las funciones: g1 (x) = x − f (x) g2 (x) = (−1 − x3 − 2x2 )/5 g3 (x) = −1/(x2 + 2x + 5)

2.1. Algoritmos Para B´ usqueda de Ra´ıces

29

son ejemplos de posibles elecciones para la funci´on g. El m´etodo de punto fijo requiere de un punto inicial cercano al valor de la ra´ız y partir de ´este, se comienzan a realizar iteraciones hasta aproximarse cada vez m´as al valor de la ra´ız que se esta buscando. Algoritmo 2.1.2: M´etodo de Punto Fijo. Entrada: funci´on f ; aproximaci´on inicial p0 ; n´ umero m´aximo de iteraciones; tolerancia T OL. Salida: soluci´on aproximada p o mensaje de error. Proceso: Paso 1: Para i = 1 y hasta i = N0 hacer pasos 2 a 4: paso 2: Tomar p = g(p0 ). paso 3: Si |p − p0 | < T OL, entonces SALIDA(p). PARAR. paso 4: hacer p0 = p. paso 7: SALIDA(Procedimiento terminado sin ´exito despu´es de N0 iteraciones). PARAR. Ejemplo 2.1.4 Se desea encontrar una aproximaci´on a una ra´ız de la funci´on f (x) = x4 + 2x2 − x − 3 mediante el m´etodo del punto fijo. Para ello se proponen las siguientes funciones como posibles candidatas para la elecci´on de g. g1 (x) = (3 + x − 2x2 )1/4 4 )1/2 g2 (x) = ( x+3−x 2 x+3 1/2 g3 (x) = ( x2 +2 ) 4 +2x2 +3 g4 (x) = 3x 4x3 +4x−1 la tabla 2.1.1 proporciona los resultados del m´etodo de iteraci´on de punto fijo para las cuatro opciones de g con p0 = 1. De la tabla 2.1.1 se observa que aunque las funciones gi son problemas de punto fijo para la misma funci´on f , difieren ampliamente como m´etodos de aproximaci´on a la soluci´on del problema. X

2.1. Algoritmos Para B´ usqueda de Ra´ıces n 0 1 2 3 4 5 6 7 8 9 10

g1 1 1.189207115 1.080057753 1.149671431 1.107820530 1.133932285 1.118003118 1.127857163 1.121813166 1.125539874 1.123249432

g2 1 1.224744871 0.993666159 1.228568645 0.987506429 1.232183418 0.981585828 1.235563209 0.97596246 1.238689029 0.970684691

30 g3 g4 1 1 1.154700538 1.142857143 1.11642741 1.12448169 1.126052233 1.124123164 1.123638885 1.12412303 1.124244497 1.12412303 1.124092553 1.124130676 1.124121111 1.124123511 1.124122909

Figura 2.1.1: M´etodo de punto fijo para la misma funci´on f y con diferentes funciones g

El siguiente teorema es de gran importancia para la exclusi´on de ciertas funciones gi0 s. Teorema de Punto Fijo Sea g una funci´on continua en un intervalo [a, b] tal que es derivable en (a, b). Si existe una constante k, 0 < k < 1 con |g 0 (x)| ≤ k para cada x ∈ (a, b), entonces para cualquier n´ umero p0 ∈ [a, b], la sucesi´on definida por pn = g(pn−1 ), n ≥ 1, converge al u ´nico punto fijo p en [a, b].

2.1.4

M´ etodo de Newton-Raphson y de la Secante

Consideremos la figura 2.1.2. Deseamos encontrar la aproximaci´on a una ra´ız de f de la siguiente manera. Dado un punto inicial p0 , calculamos la tangente a este punto v´ıa la derivada y v”ia la pendiente de una recta dados dos puntos, siendo estos, los puntos (p0 , f (p0 )) y (p1 , 0) como se muestra en la figura 2.1.2. Se tiene f 0 (p0 ) = de donde

f (p0 ) p0 −p1

2.1. Algoritmos Para B´ usqueda de Ra´ıces

p1 = p0 −

31 f (p0 ) . f 0 (p0 )

Encontrado el valor p1 , se hace el mismo procedimiento para el c´alculo de p2 , obteniendose p2 = p1 −

f (p1 ) . f 0 (p1 )

Se continua haciendo el mismo proceso hasta obtener una aproximaci´on adecuada a la ra´ız buscada.

Figura 2.1.2: M´etodo de Newton-Raphson

Algoritmo 2.1.3: M´etodo de Newton-Raphson. ENTRADA: funci´on f ; aproximaci´on inicial p0 ; n´ umero m´aximo de iteraciones N0 ; tolerancia T OL. SALIDA: soluci´on aproximada p o mensaje de error. PROCESO: paso 1: Para i = 1 hasta i = N0 hacer pasos 2 a 4: paso 2: Hacer p = p0 − f (p0 )/f 0 (p0 ) paso 3: Si |p − p0 | < T OL, entonces SALIDA(p). PARAR. paso 4: Hacer p0 = p. paso 5: SALIDA(Procedimiento terminado sin ´exito despu´es de N0 iteraciones). PARAR. Ejemplo 2.1.5 Sean f (x) = x2 − 6 y p0 = 1. Encontrar el valor de p2 . Se tiene

2.1. Algoritmos Para B´ usqueda de Ra´ıces

p1 = 1 − =1− = 3.5

32 f (1) f 0 (1) 12 −6 2(1)

luego p2 = 3.5 −

f (3.5) f 0 (3.5) 3.52 −6 2(3.5)

= 3.5 − = 2.067. X

El elemento pn en el m´etodo de Newton-Raphson est´a dado por: pn = pn−1 −

f (pn−1 ) f 0 (pn−1 )

para n > 0.

Para n > 1 anteriormente se tienen c´alculados al menos los puntos p0 y p1 . Si se desea evitar el problema de evaluar en la derivada de una funci´on se puede recurrir a calcular a la secante que une a los puntos (pn−1 , f (pn−1 )) y (pn−2 , f (pn−2 )) como una aproximaci´on a la derivada en el punto pn−1 , obteni´endose: pn = pn−1 − = pn−1 − = pn−1 −

f (pn−1 ) f 0 (pn−1 ) f (pn−1 ) f (pn−1 )−f (pn−2 ) pn−1 −pn−2

f (pn−1 )(pn−1 −pn−2 ) f (pn−1 )−f (pn−2 )

Esta ultima f´ormula alternativa al m´etodo de Newton-Raphson recibe el nombre del m´etodo de la secante. Para ello, se requieren dos puntos de aproximaci´on a la ra´ın el lugar de uno como en el m´etodo de Newton-Raphson. Algoritmo 2.1.4: M´etodo de la Secante. ENTRADA: funci´on f ; aproximaci´on inicial p0 y p1 ; n´ umero m´aximo de iteraciones N0 ; tolerancia T OL. SALIDA: soluci´on aproximada p o mensaje de error. PROCESO: se requieren dos variables auxiliares q0 y q1 . paso 1: Tomar q0 = f (p0 ) y q1 = f (p1 ) paso 2: Para i = 2 hasta i = N0 hacer pasos 3 a 5: paso 3: Hacer p = p1 − q1 (p1 − p0 )/(q1 − q0 )

2.2. Ra´ıces de Polinomios

33

paso 4: Si |p − p1 | < T OL, entonces SALIDA(p). PARAR. paso 5: Hacer p0 = p1 , q0 = q1 , p1 = p y q1 = f (p). paso 6: SALIDA(Procedimiento terminado sin ´exito despu´es de N0 iteraciones). PARAR. Ejemplo 2.1.6 Encontrar el valor de p3 con p0 = 3 y p1 = 2 para la funci´on f (x) = x2 − 6. Se tiene: (2)(2−3) p2 = 2 − ff (2)−f (3) = 2.3333

luego p3 = 2.3333 − f (2.3333)(2.3333−2) f (2.3333)−f (2) = 2.4615 X

2.2

Ra´ıces de Polinomios

Las funciones polinomiales son de gran importancia en el estudio del an´alisis num´erico. Su importancia se notar´a en los cap´ıtulos siguientes donde se tratar´a el tema de interpolaci´on. Est´a secci´on esta dedicada a este tipo de funciones especiales: los polinomios. Definici´ on 2.2.1 Un polinomio de grado n es una funci´on P : R → R que tiene la siguiente forma: P (x) = an xn + an−1 xn−1 + · · · + a1 x + a0 con ai valores constantes llamadas coeficientes del polinomio P (x) y an 6= 0 El siguiente teorema es de suma importancia en la teor´ıa algebra´ıca que concierne a los polinomios. Su importancia data en el hecho de que siempre es posible encontrar ra´ıces para un polinomio de grado n ≥ 1 en el campo de los n´ umeros complejos.

2.2. Ra´ıces de Polinomios

34

Teorema Fundamental del Algebra Si P (x) es un polinomio de grado n ≥ 1 con coeficientes reales o complejos, entonces P (x) tiene al menos una ra´ız (posiblemente compleja). Corolario 2.2.1 Si P (x) es un polinomio de grado n ≥ 1 con coeficientes reales o complejos, entonces existen constantes u ´nicas x1 , x2 , . . . , xk (posiblemente complejas), y enteros k positivos m1 , m2 , . . . , mk con Σi=1 mi = n tal que: P (x) = an (x − x1 )m1 (x − x2 )m2 · · · (x − xk )mk Corolario 2.2.2 Si P (x) y Q(x) son polinomios de grado a lo m´as n con P (xi ) = Q(xi ) para x1 , x2 , . . . , xk con k > n, entonces P (x) = Q(x) para todos los valores de x. Definici´ on 2.2.1 Sea P (x) un polinomio. r es ra´ız de multiplicidad m de P (x), si existe un polinomio Q(x) tal que P (x) = (x − r)m Q(x) con Q(r) 6= 0 Algoritmo de la Divisi´ on Dados los polinomios P (x) y Q(x), existen dos u ´nicos polinomios C(x) y R(x) tal que: P (x) = C(x)Q(x) + R(x) con grado(R(x)) < grado(Q(x)) o grado(R(x)) = 0 Teorema del Resto Si Q(x) es el cociente y R(x) es el resto de la divisi´on de un polinomio P (x) por el polinomio x − a aplicando el algoritmo de la divisi´on, entonces R(a) = P (a). Dado un polinomio P (x). Por el algoritmo de la divisi´on para polinomios, se tiene P (x) = Q(x)(x − a) + R(x) con grado(R(x)) = 0. luego P 0 (x) = Q0 (x)(x − a) + Q(x) y por tanto

2.2. Ra´ıces de Polinomios

35 P 0 (a) = Q(a).

El m´etodo de Horner o divisi´on sint´etica es un m´etodo que nos ayuda a encontrar el cociente y el resto al dividir un polinomio por x − a para alg´ un a, y por ende por lo anterior, nos ayuda a encontrar el valor de P 0 (a). M´ etodo de Horner Sea P (x) = an xn + an−1 xn−1 + · · · + a1 x + a0 . Si bn = an y bk = ak + bk+1 x0 para k = n − 1, n − 2, . . . , 1, 0, entonces b0 = P (x0 ). M´as a´ un, si Q(x) = bn xn−1 + bn−1 xn−2 + · · · + b2 x + b1 , entonces P (x) = (x − x0 )Q(x) + b0 . El teorema anterior describe el m´etodo de Horner, sin embargo para llevarlo a la pr´actica se sigue el siguiente algoritmo. Algoritmo 2.2.1: M´etodo de Horner. ENTRADA: grado n del polinomio; coeficientes a0 , a1 , . . . , an ; x0 . SALIDA: y = P (x0 ); z = P 0 (x0 ). PROCESO: paso 1: Tomar y = an y z = an paso 2: Para j = n − 1, n − 2, . . . , 1 tomar y = x 0 y + aj z = x0 z + y paso 3: Hacer y = x0 y + a0 paso 4: SALIDA(y,z). PARAR. Al realizar los c´alculos con el m´etodo de Horner se recurre a una tabla que se construye como en el ejemplo siguiente. Ejemplo 2.2.1 Al aplicar el m´etodo de Horner al polinomio p(x) = 3x4 + 2x3 − x2 + 4x + 2 en x0 = −2, se tiene:

2.2. Ra´ıces de Polinomios coeficiente de x0 = −2

x4 a4 = 3 b4 = 3

36 x3 x2 x1 x0 a3 = 2 a2 = −1 a1 = 4 a0 = 2 b4 x0 = −6 b3 x0 = 8 b2 x0 = −14 b1 x0 = 20 b3 = −4 b2 = 7 b1 = −10 b0 = 22

por lo tanto P (x) = (x + 2)(3x3 − 4x2 + 7x − 10) + 22. X

37

Cap´ıtulo 3 Algoritmos de Interpolaci´ on 3.1

Algoritmos de Interpolaci´ on

Figura 3.1.1: Gr´ afica que muestra puntos aislados tomados de una muestra de alg´ un experimento hipot´etico.

Supongamos que la figura 3.1.1 muestra los datos registrados en un experimento f´ısico que depende del tiempo. El eje de la abscisa muestra los diferentes tiempos en que se llevo a cabo el mismo experimento repetidas veces, y el eje de la ordenada muestra los valores aproximados del resultado del experimento. Este tipo de registros se obtienen a menudo en un laboratorio de f´ısica, puede haber registros para lapsos cortos de tiempo o para lapsos muy largos. Un ejemplo de tales registros, puede ser el de la poblaci´on ya sea a nivel de un pais o la poblaci´on mundial. Por ejemplo, el eje de las abscisas puede representar las mediciones registradas en el a˜ no t0 , t1 , etc., y el eje de la ordenada puede representar la poblaci´on para el respectivo a˜ no de registro. Una vez obtenidos los datos como se muestran en la figura 3.1.1, lo interesante del asunto es determinar el valor que le corresponde a la funci´on para un valor no registrado, esto es, por ejemplo, en la figura se muestran datos registrados para valores t0 , t1 hasta t5 , para

3.1. Algoritmos de Interpolaci´ on

38

estos valores, conocemos el valor exacto que le corresponde a la funci´on f (t), sin embargo que sucede si se desea conocer un valor t0 que se encuentre en el intervalo (t0 , t1 ), este valor no fue registrado. La problem´atica a resolver enseguida ser´a encontrar m´etodos con los cuales sea posible encontrar dicho valor. La idea data en encontrar una funci´on continua que pase por los puntos registrados, de tal manera que sea posible determinar de manera confiable aquellos valores que no fueron registrados. A esto se le conoce como interpolaci´ on. Comenzaremos dicho estudio con el polinomio de Lagrange.

Figura 3.1.2: Ejemplo de Interpolaci´on.

3.1.1

Polinomio de Lagrange

Consideremos el problema de encontrar un polinomio de grado m´aximo n que pase por los n + 1 puntos: (x0 , f (x0 )), (x1 , f (x1 )), . . . , (xn , f (xn )). Sea Ln.k (x) = =

(x−x0 )(x−x1 )...(x−xk−1 )(x−xk+1 )···(x−xn ) (xk −x0 )(xk −x1 )···(xk −xk−1 )(xk −xk+1 )···(xk −xn ) n Y (x − xi ) i=0,i6=k

(xk − xi )

para k = 0, 1, . . . , n

para k = 0, 1, . . . , n

Se define el n-´esimo polinomio interpolante de Lagrange como: Pn (x) = f (x0 )Ln,0 (x) + f (x1 )Ln,1 (x) + · · · + f (xn )Ln,n (x) n X = f (xk )Ln,k (x). k=0

3.1. Algoritmos de Interpolaci´ on

39

Si el grado del polinomio es claro, entonces escribiremos Lk en lugar de Ln,k . Observemos que  Lk (xj ) =

1 si j = k 0 si j 6= k

luego, se tiene: Pn (xk ) = f (xk ) para k = 0, 1, . . . , n. esto es, el polinomio de Lagrange es tal que pasa por cada uno de los puntos (xk , f (xk )). Ejemplo 3.1.1 Para encontrar el polinomio de grado dos para los valores de la siguiente tabla: x f (x) 8.1 16.94410 8.3 17.56492 8.6 18.50515 primeramente se encuentran los valores de L0 (x), L1 (x) y L2 (x), obteniendo: L0 (x) = = = L1 (x) = = = L2 (x) = = =

(x−8.3)(x−8.6) (8.1−8.3)(8.1−8.6) x2 −16.9x+71.38 (−0.2)(−0.5) 1 (x2 − 16.9x + 71.38) 0.1 (x−8.1)(x−8.6) (8.3−8.1)(8.3−8.6) x2 −16.7x+69.66 (0.2)(−0.3) 1 − 0.06 (x2 − 16.7x + 69.66) (x−8.1)(x−8.3) (8.6−8.1)(8.6−8.3) x2 −16.4x+67.23 (0.5)(0.3) 1 (x2 − 16.4x + 67.23) 0.15

por tanto, el polinomio dos de Lagrange est´a dado por: P2 (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + f (x2 )L2 (x) = (16.94410)(10x2 − 169x + 713.8) + (17.56492)(−16.66x2 + 278.222x − 1160.53) + (18.50515)(6.66x2 − 109.224x + 447.751) = 0.054x2 + 2.19x − 4.23 X

3.1. Algoritmos de Interpolaci´ on

40

En el ejemplo 3.1.1 se calculo el polinomio dos de Lagrange, a partir de ´este se pueden encontrar valores aproximados para datos que no se encuentran en la tabla, por ejemplo, se tiene P2 (8.2) = 17.35896, P2 (8.4) = 17.97624 y P2 (8.5) = 18.2865. Enseguida se muestran el m´etodo de Neville seguido del m´etodo de diferencias divididas de Newton que de igual manera sirven para interpolar.

3.1.2

M´ etodo de Neville

Algoritmo 3.1.1: M´etodo de Neville. ENTRADA: valores x, x0 , x1 , . . . , xn ; valores P0,0 = f (x0 ), P1,0 = f (x1 ), . . . Pn,0 = f (xn ). SALIDA: Tabla 3.1.1. con P (x) = P0,n PROCESO: paso 1: Para j = 0, 2, . . . , n hacer para i = 0, . . . , n − j tomar (x−xi )Pi+1,j−1 +(xi+j −x)Pi,j−1 Pi,j = xi+j −xi paso 2: SALIDA(P (x)). PARAR. i 0 1 ··· n

xi x0 x1 ··· xn

Pi,0 P0,0 P1,0 ··· Pn,0

Pi,1 P0,1 P1,1 ···

Pi,2 P0,2 P1,2 ···

··· ··· ··· ···

P i, n P0,n

Tabla 3.1.1: M´etodo de Neville.

√ Ejemplo 3.1.2 Apliquemos √ el m´etodo de Neville para para aproximar el valor de 3 considerando la funci´on f (x) = x y los valores x0 = 0, x1 = 1, x2 = 2, x3 = 4 y x4 = 5. Enseguida se muestra la tabla 3.1.1 para este problema en especifico. i 0 1 2 3 4

xi 0 1 2 4 5

Pi,0 Pi,1 Pi,2 Pi,3 Pi,4 0 3 0.9319 1.4918 1.6180 1 1.6213 1.6785 1.7022 1.4142 1.7071 1.7260 2 1.7640 2.2360

3.1. Algoritmos de Interpolaci´ on

41

Por ejemplo, para calcular los valores P0,1 , P1,1 , P1,2 , se hace uso del paso 1 del algoritmo 3.1.1, teni´endose: +(x1 −x)P0,0 P0,1 = (x−x0 )P1,0 x1 −x0 = (3−0)(1)+(1−3)(0) 1−0 =3 +(x2 −x)P1,0 P1,1 = (x−x1 )P2,0 x2 −x1 = (3−1)(1.4142)+(2−3)(1) 2−1 = 1.6213 +(x3 −x)P1,1 P1,2 = (x−x1 )P2,1 x3 −x1 = (3−1)(1.7071)+(4−3)(1.6213) 4−1 = 1.6785

El paso 2 del algoritmo 3.1.1 muestra el valor del polinomio interpolado y evaluado en el valor de la variable √ x, por tanto, para este ejemplo en especifico, se tiene 1.6180 = P4,4 = P (x = 3) ≈ f (3) = 3. X

3.1.3

Diferencias Divididas de Newton

Sean (x0 , f (x0 )), (x1 , f (x1 )), . . . , (xn , f (xn )) n puntos que se conocen a los que se desea interpolar. Para abreviar denotemos fi = f (xi ). Consideremos el siguiente polinomio de grado n: Pn = a0 + (x − x0 )a1 + (x − x0 )(x − x1 )a2 + · · · + (x − x0 )(x − x1 ) · · · (x − xn−1 )an

(3.1.1)

Si los valores a0i s se escogen adecuadamente de modo que Pn (xi ) = f (xi ) en los n + 1 puntos conocidos, entonces se puede conseguir que Pn sea un polinomio de interpolaci´on. Nuestro objetivo ser´a encontrar los valores adecuados para las a0i s. Para ello, consideremos primeramente la siguiente notaci´on. Se denomina primera diferencia dividida entre los valores xs y xt a el valor: f [xs , xt ] =

ft −fs xt −xs

Obs´ervese que f [xs , xt ] = f [xt , xs ]. Las diferencias de orden superior se definen en t´erminos de diferencias de orden inferior, por ejemplo, la diferencia dividida de segundo orden de los valores x0 y x2 se define como:

3.1. Algoritmos de Interpolaci´ on

42

f [x0 , x1 , x2 ] =

f [x1 ,x2 ]−f [x0 ,x1 ] x2 −x0

y en general, la diferencia dividida de orden n entre los valores x0 y xn se define como: f [x0 , x1 , . . . , xn ] =

f [x1 ,x2 ,...,xn ]−f [x0 ,x1 ,...xn−1 ] . xn −x0

La diferencia dividida de orden cero se define como: f [xs ] = fs . Consideremos nuevamente el polinomio mostrado en la ecuaci´on 3.1.1, Si deseamos que el polinomio interpole los puntos (xi , fi ) para i = 0, . . . , n, entonces se debe de cumplir que Pn (xi ) = fi , obteni´endose entonces: f0 = Pn (x0 ) = a0 de donde a0 = f0 = f [x0 ], continuando, se tiene f1 = Pn (x1 ) = a0 + (x1 − x0 )a1 de donde a1 =

f1 −a0 x1 −x0

=

f1 −f0 x1 −x0

= f [x1 , x0 ].

continuando con este proceso, se tiene en general, que: ai = f [x0 , x1 , . . . , xi ] por tanto Pn (x) = f0 + (x − x0 )f [x0 , x1 ] + (x − x0 )(x − x1 )f [x0 , x1 , x2 ] + · · · + (x − x0 )(x − x1 ) · · · (x − xn−1 )f [x0 , x1 , . . . , xn ].

3.1. Algoritmos de Interpolaci´ on

43

Enseguida se muestra el algoritmo de el m´etodo de diferencias divididas. Algoritmo 3.1.2: M´etodo de Diferencias Divididas. ENTRADA: valores x, x0 , x1 , . . . , xn ; valores f (x0 ), f (x1 ), . . . f (xn ). SALIDA: Tabla 3.1.2 PROCESO: paso 1: Para i = 1, 2, . . . , n hacer para j = 0, . . . , n − i calcular [j,...,i−1] f [j, . . . , n] = f [j+1,...,i]−f y se coloca en la xj+1 −xj columna i de la tabla 3.1.2 paso 2: SALIDA(P (x)). PARAR. xi x0 x1 ··· xn

fi f0 f1 ··· fn

f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [x0 , x1 ] f [x0 , x1 , x2 ] f [x1 , x2 ] f [x1 , x2 , x3 ] ··· ···

··· ··· ··· ···

f [xi , xi+1 , . . . , xi+n ] f [x0 , x1 , . . . , xn ]

Tabla 3.1.2: M´etodo de Diferencias Divididas.

Ejemplo 3.1.3 Apliquemos el m´etodo de √ √diferencias divididas para para aproximar el valor de 3 considerando la funci´on f (x) = x y los valores x0 = 0, x1 = 1, x2 = 2, x3 = 4 y x4 = 5. Enseguida se muestra la tabla 3.1.2 para este problema en especifico. i xi f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ] f [xi , xi+1 , xi+2 , xi+3 , xi+4 ] 0 0 1 -0.2929 0.0631 -0.0115 1 1 0.4142 -0.0404 0.0053 2 1.4142 0.2929 -0.0189 3 2 0.2360 4 2.2360 Por ejemplo, para calcular los valores f [x0 , x1 ] y f [x1 , x2 , x3 ], se hace uso del paso 1 del algoritmo 3.1.2, teni´endose:

3.2. Algoritmos de Ajuste de Curvas

44

−f0 f [x0 , x1 ] = xf11 −x 0 1−0 = 1−0 =1 [x1 ,x2 ] f [x1 , x2 , x3 ] = f [x2 ,xx33]−f −x1 = 0.2929−0.4142 2−1 = 0.0404

Se tiene que el polinomio interpolado es: P4 (x) = f0 +(x−x0 )f [x0 , x1 ]+(x−x0 )(x−x1 )f [x0 , x1 , x2 ]+(x−x0 )(x−x1 )(x−x2 )f [x0 , x1 , x2 , x3 ]+ (x − x0 )(x − x1 )(x − x2 )(x − x3 )f [x0 , x1 , x2 , x3 , x4 ] = (x)(1)+(x)(x−1)(−0.2929)+(x)(x−1)(x−2)(0.0631)+(x)(x−1)(x−2)(x−4)(−0.0115) y por tanto P4 (3) = 1.6902. X

3.2

Algoritmos de Ajuste de Curvas

Supongamos que se se han registrado una cierta cantidad de puntos que pudieron haber sido el resultado de un experimento, y dichos registros de puntos son parecidos a como se muestran en la tabla 3.1.3 a). El objetivo de est´a secci´on es tratar de resolver la problem´atica de encontrar un polinomio de grado n >= 1 que se ajuste m´as a los puntos registrados, esto es, obtener un polinomio que se le parezca m´as a los puntos registrados, de tal manera que sea posible encontrar el valor aproximado de un dato no registrado. La primera parte, ser´a encontrar dicho polinomio de grado n = 1 (figura 3.1.3 b)), enseguida se desarrollar´a de manera general para encontrar un polinomio de grado n > 1.

Figura 3.1.3: a) Registro de puntos. b) Ajuste de puntos a un polinomio de grado 1.

3.2. Algoritmos de Ajuste de Curvas

3.2.1

45

Ajuste por M´ınimos Cuadrados a una Recta

Consideremos la figura 3.1.3 b) en donde se muestra el ajuste de puntos por medio de una recta. Consideremos el punto registrado (xi , yi ), el correspondiente valor de la ecuaci´on de la recta para este valor de la abscisa es a1 xi + a0 . La problem´atica consiste en encontrar los valores de a0 y a1 de tal manera que la distancia entre el valor de yi y a1 xi + a0 sea m´ınima. El m´etodo de m´ınimos cuadrados consiste en determinar la mejor linea de aproximaci´on, cuando el error es la suma de los cuadrados de las diferencias entre los valores de y dados; esto es, consiste en determinar los valores de a0 y a1 de tal manera que el valor de la siguiente diferencia sea m´ınima. m X

[yi − (a1 xi + a0 )]2

(3.1.1)

i=1

con {(xi , yi )}m i=1 los puntos de datos a los que se desea realizar el ajuste. Para que la ecuaci´on 3.1.1 tenga un valor m´ınimo, debemos de tener que:

0=

∂ ∂a0

m m X X [yi − (a1 xi − a0 )]2 = 2 (yi − a1 xi − a0 )(−1) i=1

i=1

y 0=

∂ ∂a1

m X

m X [yi − (a1 xi + a0 )] = 2 (yi − a1 xi − a0 )(−xi )

i=1

i=1

2

Las ecuaciones anteriores se simplifican en las siguientes ecuaciones

a0 m + a1

m X

xi =

i=1

m X yi i=1

y a0

m X i=1

m m X X 2 x i + a1 x i = x i yi i=1

i=1

de donde, la soluci´on al anterior sistema de ecuaciones con a0 y a1 con inc´ognitas es:

3.2. Algoritmos de Ajuste de Curvas

46

m m m m X X X X x2i yi − xi y i xi i=1

a0 =

i=1

i=1

i=1

m m X X m( x2i ) − ( xi )2 i=1

(3.1.2)

i=1

y m m m X X X m xi yi − xi y i

a1 =

i=1 m X

m(

i=1

x2i )

i=1 i=1 m X

−(

xi )

(3.1.3) 2

i=1

Ejemplo 3.1.4 Ajustar los siguientes puntos a la ecuaci´on de una recta. xi yi

1.0 1.1 1.3 1.5 1.9 2.1 1.84 1.96 2.21 2.45 2.94 3.18

Para realizar el ajuste respectivo, se sustituyen los valores directamente en las ecuaciones 3.1.2 y 3.1.3, teni´endose: a0 = (14.17)(14.58)−(22.78)(8.9) 6(14.17)−(8.9)2 = 0.663 a1 = 6(22.78)−(8.9)(14.58) 6(14.17)−(8.9)2 = 1.190 Por lo tanto y = a0 + a1 x = 0.663 + 1.190x es el polinomio que mejor se ajusta a los puntos dados. X En el ejemplo anterior, se observa que y(x = 1.1) = 1.197, y(x = 1.5) = 2.448 y y(x = 2.1) = 3.162, es decir, este polinomio se aproxima a los valores registrados en la tabla correspondientes a los valores de las ordenadas respectivas, por otro lado, si deseamos conocer alguna aproximaci´on de un valor no registrado, ahora podemos recurrir a dicho polinomio, por ejemplo y(x = 1.4) = 2.329. Enseguida, nuestro objetivo ser´a mostrar el ajuste de puntos a un polinomio de grado n ≥ 1.

3.2. Algoritmos de Ajuste de Curvas

3.2.2

47

Ajuste por M´ınimos Cuadrados a un Polinomio de Grado n

La secci´on anterior, se mostr´o el ajuste de puntos a una recta, sin embargo que sucede si dichos puntos no se acoplan a una recta del todo bien, que sucede si en lugar de acoplarse a una recta, se ajustan mejor a un polinomio de segundo grado, o a uno de tercero. La figura 3.1.4 muestra el ajuste de puntos a un polinomio de grado 2. En est´a secci´on se mostrar´an las f´ormulas respectivas para el ajuste de puntos a un polinomio de grado n.

Figura 3.1.4: Ajuste de puntos a un polinomio de grado n.

Supongamos que {(xi , yi )}m i=1 son los puntos que se desean ajustar a un polinomio de la forma Pn (x) = an xn + an−1 xn−1 + · · · + a1 x + a0 con n < m − 1, mediante ajuste de m´ınimos cuadrados. La forma de resolver esta problem´atica es de manera similar a como se hizo en la secci´on anterior al encontrar los valores a0 y a1 mediante la f´ormula 3.1.1. En este caso, se desean encontrar los valores a0 , a1 , . . . , an de tal manera que el valor de la siguiente diferencia sea m´ınima.

E= =

m X i=1

m X

(yi − Pn (xi ))2

i=1

yi2

−2

m X

Pn (xi )yi +

i=1

(3.1.4) m X

(Pn (xi ))2

i=1

m m X n m X n X X X j 2 = yi − 2 ( aj xi )yi + ( aj xji )2 i=1

i=1 j=1

i=1 j=1

De igual manera al caso de la secci´on anterior, para que la ecuaci´on 3.1.4 tenga el m´ınimo valor, se debe de cumplir que

3.2. Algoritmos de Ajuste de Curvas ∂E ∂aj

48

= 0 para j = 0, 1, . . . , n

(3.1.5)

de donde se obtiene un sistema de ecuaciones de n+1 inc´ognitas (siendo las a0i s las inc´ognitas) y n + 1 ecuaciones. Se tiene que el sistema de ecuaciones 3.1.5 es equivalente al siguiente sistema de ecuaciones: a0

m X

x0i + a1

i=1

m X

x1i + a2

i=1

m X

x2i + · · · + an

i=1

m X

xni =

i=1

m X

yi x0i

i=1

m m m m m X X X X X a0 x1i + a1 x2i + a2 x3i + · · · + an xn+1 = yi x1i i

a0

i=1

i=1

m X

m X

i=1

xni

+ a1

i=1

i=1

xn+1 i

+ a2

m X i=1

··· xn+2 i

i=1

+ · · · + an

(3.1.6)

i=1 m X

xin+n

=

i=1

m X

yi xni

i=1

Por tanto, para ajustar los puntos {(xi , yi )}m i=1 a un polinomio de grado n < m − 1 se debe de resolver el sistema de ecuaciones 3.1.6. Ejemplo 3.1.5 Consideremos los puntos dados en el ejemplo 3.1.4 y procedamos a realizar el ajuste de dichos puntos a un polinomio de grado 2. Para ello, se tienen que encontrar constantes a0 , a1 y a2 tales que satisfagan el sistema de ecuaciones 3.1.6 con m = 6 y n = 2, esto es: 6a0 + 8.9a1 + 14.17a2 = 14.58 8.9a0 + 14.17a1 + 24.023a2 = 22.808 14.17a0 + 24.023a1 + 42.8629a2 = 38.0962 Resolviendo el anterior sistema de ecuaciones, se tiene que P2 (x) = a0 + a1 x + a2 x2 = 0.5965 + 1.2553x − 0.0108x2 . X

49

Cap´ıtulo 4 Algoritmos para Derivaci´ on e Integraci´ on Num´ erica El objetivo de est´e cap´ıtulo es mostrar la derivaci´on e integraci´on de funciones por medio ´ de m´etodos num´ericos. Unicamente se mostrar´an las f´ormulas sin el previo desarrollo que conlleva a estas. Si se desea tener un estudio profundo de las f´ormulas, se puede revisar [1]. Comenzaremos con la parte referente a la derivaci´on seguida de integraci´on simple e integraci´on compuesta.

4.1

Algoritmos para Derivaci´ on

Supongamos que se tiene una funci´on f ∈ C n+1 (I) con I un intervalo. Sea {x0 , x1 , . . . , xn } una sucesi´on de n + 1 puntos distintos en I. La f´ormula 4.1.1 recibe el nombre de f´ormula de (n+1) puntos para aproximar f 0 (xj ). 0

f (xj ) =

m X k=0

f (xk )L0k (xj ) +

f (n+1) (ε(xj )) m Πk=0,k6=j (xj − xk ) (n + 1)!

(4.1.1)

para alguna ε(x) ∈ I y con Lk (x) el k-´esimo polinomio de Lagrange. De la ecuaci´on 4.1.1 se concluye que f 0 (x).

m X

f (xk )L0k (xj ) es una aproximaci´on al valor de

k=0

La figura 4.1.1 muestra la interpretaci´on gr´afica que se tiene al usar la ecuai´on 4.1.1 en el caso que se consideren tres puntos. Obs´ervese que el punto a encontrar la aproximaci´on de su derivada se encuentra en el medio de los otros dos puntos.

4.1. Algoritmos para Derivaci´ on

50

Figura 4.1.1: f´ormula de tres puntos.

Existen dos casos especiales del uso de la ecuaci´on 4.1.1, una de ellas se muestra en la figura 4.1.1, donde se hace uso de la ecuaci´on 4.1.1 para el caso de tres puntos; y la otra cuando se consideran cinco puntos. Enseguida se muestran las f´ormulas para el caso de tres y cinco puntos cuyas f´ormulas se les conoce como f´ormula de tres puntos y f´ormula de cinco puntos respectivamente. En el caso de la f´ormula de tres puntos se muestran dos variantes. f´ormulas de tres puntos: f 0 (x0 ) ≈ f 0 (x0 ) ≈

1 [−3f (x0 ) + 4f (x0 2h 1 [f (x0 + h) − f (x0 2h

+ h) − f (x0 + 2h)] − h)]

(4.1.2) (4.1.3)

f´ormula de cinco puntos: f 0 (x0 ) ≈

1 [f (x0 12h

− 2h) − 8f (x0 − h) + 8f (x0 + h) − f (x0 + 2h)]

Ejemplo 4.1.1 Considere la siguiente tabla de valores con f (x) = x 1.5 1.6 1.7 1.8 1.9

(4.1.4)

√ x+1

f (x) 1.5811 1.6124 1.6431 1.6733 1.7029

1 0 0 Se tiene que f 0 (x) = 2(x+1) 1/2 , por lo tanto f (1.7) = 0.3042. Al aproximar f (1.7) mediante las f´ormulas de tres y cinco puntos se obtienen los siguientes resultados.

F´ormula de tres puntos:

4.2. Algoritmos de Integraci´ on

51

f´ormula 4.1.2 con h = 0.1: f 0 (1.7) = f´ormula 4.1.3 con h = 0.1: f 0 (1.7) =

1 [−3f (1.7) + 4f (1.8) − f (1.9)] 2(0.1) 1 [f (1.8) − f (1.6)] = 0.3043 2(0.1)

= 0.3041

F´ormula de cinco puntos: f´ormula 4.1.4 con h = 0.1: f 0 (1.7) =

4.2

1 [f (1.5) 12(0.1)

− 8f (1.6) + 8f (1.8) − f (1.9)] = 0.3042. X

Algoritmos de Integraci´ on

Consideremos la figura 4.1.2. Supongamos que deseamos calcular el a´rea bajo la curva de la funci´on f (x) dentro de un intervalo [a, b]. La soluci´on a esta problem´atica es bien conocida y es mediante la integraci´on de la funci´on f con extremo inferior el valor de a y extremo superior el valor de b. Sin embargo, en est´a secci´on se resolver´a el problema de integrar una funci´on f (x) sobre un intervalo [a, b] pero desde el punto de vista del an´alisis num´erico.

Figura 4.1.2: Integraci´ on mediante f´ormula cerrada de Newton-Cotes.

Consideremos nuevamente la figura 4.1.2. Supongamos que se desea calcular la integral sobre el intervalo [a, b]. Para ello, primeramente hagamos una partici´on del intervalos [a, b] de la siguiente manera. Consideremos la partici´on {x0 , x1 , x2 , . . . , xn } del intervalo [a, b] compuesta por (n + 1) puntos de tal forma que x0 = a, x1 = a + h, x2 = a + 2h, . . . , xi = a + ih, . . . xn = b donde h = (b − a)/n. Considerando est´a partici´on, se pueden obtener f´ormulas para aproximar el valor de la integral. Puesto que se incluyen los extremos del intervalo [a, b] como elementos de la partici´on, las f´ormulas concernientes al c´alculo de la integral bajo est´a situaci´on se les conoce como f´ormulas cerradas de (n+1) puntos de NewtonCotes. Las f´ormulas cerradas de Newton-Cotes para el caso n = 1, n = 2 y n = 3 reciben un nombre en particular; dichas f´ormulas se muestran enseguida. n = 1: Regla del Trapecio: R x1 x0

f (x)dx ≈ h2 [f (x0 ) + f (x1 )]

4.2. Algoritmos de Integraci´ on

52

n = 2: Regla de Simpson: R x2 x0

f (x)dx ≈ h3 [f (x0 ) + 4f (x1 ) + f (x2 )]

n = 3: Regla de tres octavos de Simpson: R x3 f (x)dx ≈ 3h [f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )] 8 x0 La figura 4.1.3 muestra la regla del trapecio aplicada a una funci´on f (x). Para este caso, la partici´on u ´nicamente cuenta con tres elementos de lo cuales dos son los extremos del intervalo a integral, de est´a manera, se observa que si se desea realizar una mejor exactitud de dicha integral, entonces se deben de considerar una mayor cantidad de puntos sobre el intervalo. Para la regla del trapecio y la regla de Simpson, existen sus variantes respectivas, donde el n´ umero de elementos de la partici´on es a consideraci´on de quien hace uso de dichas f´ormulas, esto es, para obtener un mayor grado de precisi´on considerando la regla del trapecio y la regla de Simpson, se tienen la regla compuesta del trapecio y la regla compuesta de Simpson, donde a diferencia de las reglas anteriores se consideran una cantidad arbitraria de (n + 1) puntos de la partici´on considerando que xi = a + ih con h = (b − a)/n (figura 4.1.2).

Figura 4.1.3: F´ormula del Trapecio.

Regla compuesta del Trapecio: Rb a

f (x)dx ≈

h [f (a) 2

+2

n−1 X

f (xj ) + f (b)]

j=1

Regla compuesta de Simpson: (n/2)−1

Rb a

f (x)dx ≈

h [f (a) 3

+2

X j=1

n/2 X f (x2j ) + 4 f (x2j−1 ) + f (b)] j=1

4.2. Algoritmos de Integraci´ on

53

Si en lugar de considerar los nodos x0 = a y xn = b de la partici´on como hasta ahora lo hemos estado haciendo con las f´ormulas cerradas de Newton-Cotes, ahora consideramos la partici´on formada por {x0 , x1 , x2 , . . . , xn } con x0 = a + h y xi = x0 + ih = a + (i + 1)h, donde h = (b − a)/(n + 2), esto implica que xn = b − h. Por tanto, si deseamos incluir los extremos del intervalo en la partici´on, estos se deben de anexar con una numeraci´on acorde, teni´endose x−1 := a y xn+1 := b (figura 4.1.4).

Figura 4.1.3: Integraci´ on mediante f´ormula abierta de Newton-Cotes.

A las f´ormulas de integraci´on para el caso de est´a partici´on, se les conoce como f´ormulas abiertas de Newton-Cotes. F´ormulas abiertas de Newton-Cotes para: n = 0: R x1 x−1

f (x)dx ≈ 2hf (x0 )

n=1 R x2 x−1

f (x)dx ≈

3h [f (x0 ) 2

+ f (x1 )]

n=2 R x3 x−1

f (x)dx ≈

4h [2f (x0 ) 3

− f (x1 ) + 2f (x2 )]

R 0.5 2 Ejemplo 4.1.2 Se tiene que 0 x−4 dx = 2ln|x − 4||00.5 = 2[ln(3.5) − ln(4)] = −0.2670. Resolviendo la integral mediante las f´ormulas de Newton-Cotes cerradas abiertas, se obtiene: F´ormulas cerradas de Newton-Cotes: Regla del trapecio(h = (b − a) = 0.5):

4.2. Algoritmos de Integraci´ on R 0.5 0

2 dx x−4

=

0.5 [f (0) 2

+ f (0.5)] = −0.2678

Regla de Simpson (h = (b − a)/2) = 0.25): R 0.5 2 dx = 0.25 [f (0) + 4f (0.25) + f (0.5)] = −0.2670 x−4 3 0 Regla de tres octavos de Simpson(h = (b − a)/3 = 0.1666): R 0.5 0

2 dx x−4

=

3(0.1666) [f (0) 8

+ 3f (0.1666) + 3f (0.333) + f (0.5)] = −0.2669

F´ormulas abiertas de Newton-Cotes para: n = 0(h = (b − a)/(2) = 0.25) : R 0.5 2 dx = 2(0.25)f (0.25) = −0.2666 x−4 0 n = 1(h = (b − a)/3) = 0.1666 R 0.5 0

2 dx x−4

=

3(0.1666) [f (0.1666) 2

+ f (0.3333)] = −0.2667

n = 2(h = (b − a)/4 = 0.125) R 0.5 0

2 dx x−4

=

4(0.125) [2f (0.125) 3

− f (0.25) + 2f (0.375)] = −0.2670. X

54

55

Cap´ıtulo 5 Algoritmos Num´ ericos Utilizados en Algebra Lineal En este cap´ıtulo se muestran algoritmos para resolver sistemas de ecuaciones lineales de la forma: a11 x1 + a12 x2 + · · · + a1n xn = b1 (5.1) a21 x1 + a22 x2 + · · · + a2n xn = b2 ··· an1 x1 + an2 x2 + · · · + ann xn = bn donde las aij son coeficientes y las xi son las inc´ognitas. Dado un conjunto de ecuaciones como en 5.1, con coeficientes dados, se desea encontrar los valores de los x0i s tales que cumplan el sistema de ecuaciones.

5.1

Eliminaci´ on Gaussiana con Sustituci´ on hacia Atr´ as

Sea Ej la representaci´on de la ecuaci´on j del sistema de ecuaciones 5.1, esto es: Ej : aj1 x1 + aj2 x2 + · · · + ajn xn = bj , entonces las operaciones fundamentales para un sistema de ecuaciones como 5.1, se muestran enseguida: 1. Si c ∈ R − {0} es una constante, entonces cEj multiplica toda la ecuaci´on j por c. Est´a nueva ecuaci´on se sustituye por Ej . Esta operaci´on se denota mediante cEj → Ej . 2. La operaci´on Ei + cEj → Ei significa que la ecuaci´on j se multiplica por la constante c y se suma con la ecuaci´on i sustituyendo esta ultima ecuaci´on por el resultado de lo anterior.

5.1. Eliminaci´ on Gaussiana con Sustituci´ on hacia Atr´ as

56

3. La opreaci´on Ei ↔ Ej significa que las ecuaciones i y j se intercambian de posici´on. Para representar un sistema de ecuaciones, a menuda suele utilizarse su representaci´on matricial sobre el conjunto de coeficientes. Consideremos el sistema de ecuaciones 5.1. Sean   a11 a12 · · · a1n  a21 a22 · · · a2n    A =  .. .. ..   . . .  an1 an2 · · · ann y   b1  b2    b =  ..   .  bn se define la matriz aumentada de A con b   a11   a [A, b] =  . 21  ..  an1

denotada como [A, b] por  . a12 · · · a1n .. b1  .  a22 · · · a2n .. b2  .. .. .. ..  . . . .   .. an2 · · · ann . bn

Si [A, b] es la matriz aumentada del sistema de ecuaciones 5.1 y si Rj representa el j-´esimo rengl´on de la matriz, entonces las operaciones fundamentales para la matriz aumentada son: 1. Si c ∈ R − {0} es una constante, entonces cRj multiplica toda la fila j por c. Est´a nueva ecuaci´on se sustituye por Rj . Esta operaci´on se denota mediante cRj → Rj . 2. La operaci´on Ri + cRj → Ri significa que la fila j se multiplica por la constante c y se suma con la fila i sustituyendo esta ultima fila por el resultado de lo anterior. 3. La opreaci´on Ri ↔ Rj significa que las filas i y j se intercambian de posici´on. Consideremos la matriz aumentada [A, b] del sistema de ecuaciones 5.1. Mediante operaciones fundamentales, esta matriz puede llevarse a otra equivalente de la forma:   .. a a · · · a1n . a1,n+1  11 12  ..  0  a · · · a . a 22 2n 2,n+2  [A, b] ≡   .  .. .. .. ..  ..  . . . . 0 0 · · · 0 ann an,n+n tal que aij = 0 si i < j. Esta ultima matriz es equivalente al sistema de ecuaciones siguiente:

5.1. Eliminaci´ on Gaussiana con Sustituci´ on hacia Atr´ as

57

a11 x1 + a12 x2 + · · · + a1n xn = a1,n+1 a22 x2 + · · · + a2n xn = a2,n+1 ··· ann xn = an,n+1 de donde xn = xn−1 =

an,n+1 , ann an−1,n+1 −an−1,n xn an−1,n−1

y en general xi =

ai,n+1 −ai,n xn −···−ai,i+1 xi+1 aii

para i = n − 1, n − 2, . . . , 2, 1. Al m´etodo anterior para la soluci´on del sistema de ecuaciones, se le conoce como eliminaci´on gaussiana con sustituci´on hacia atr´as. Ejemplo 5.1.1 Resolver el sistema de ecuaciones x1 + 2x2 + 3x3 = 1 4x1 + 5x2 + 6x3 = −2 7x1 + 8x2 + 10x3 = 5 utilizando el m´etodo de eliminaci´on gaussiana con sustituci´on hacia atr´as. La matriz aumentada del sistema es: 

 .. . 1   ..  [A, b] =  4 5 6 . −2   .. 7 8 10 . 5 1 2

3

el cual, mediante operaciones fundamentales, es equivalente a la matriz   .. 1 2 3 . 1    [A, b] ≡  0 1 2 ... 2   .. 0 0 1 . 10 de donde, x3 = 10/1 = 10, x2 = 2(−18) = 7. X

2−2x3 1

= 2 − 2(10) = −18 y x1 =

1−3x3 −2x2 1

= 1 − 3(10) −

5.2. Algoritmos Basados en T´ ecnicas Iterativas

5.2

58

Algoritmos Basados en T´ ecnicas Iterativas

En est´a secci´on se mostrar´an m´etodos iterativos para la resoluci´on de sistemas de ecuaciones como en 5.1. para ello, es necesario, primeramente definir algunos conceptos previos.

5.2.1

Normas de Vectores y Matrices

Definici´ on 5.2.1 Dado un vector x ∈ Rn , una norma en Rn es una funci´on k • k : Rn → R tal que: 1. kxk ≥ 0 para cada x ∈ Rn . 2. kxk = 0 si y s´olo si x = 0. 3. kcxk = ckxk para toda c ∈ R y x ∈ Rn . 4. kx + yk ≤ kxk + kyk para todo x, y ∈ Rn . Proposici´ on 5.2.1 Si x = (x1 , x2 , . . . , xn ) ∈ Rn es un vector, entonces n X kxk2 = [ x2i ]1/2 y kxk∞ = max1≤i≤n |xi | i=1

son normas sobre Rn . A las normas k•k2 y k•k∞ de la proposici´on anterior se les conoce como norma euclideana y norma infinito respectivamente. Ejemplo 5.2.1 Si x = (2, −1, 3) ∈ R3 , entonces kxk2 = (22 + (−1)2 + 32 )1/2 = (14)1/2 = 3.7416 y kxk∞ = max1≤i≤3 |xi | = max1≤i≤3 {|2|, | − 1|, |3|} = 3. Definici´ on 5.2.2 Sean x y y dos vectores en Rn . Una distancia d en Rn es una funci´on n d : R × Rn → R tal que 1. d(x, y) ≥ 0. 2. d(x, y) = d(y, x). 3. d(x, y) = 0 si y s´olo si x = y 4. d(x, y) ≤ d(x, z) + d(z, y) para cualquier z ∈ Rn Dada una norma, se induce una distancia dada de la siguiente manera: Proposici´ on 5.2.2 Si k • k es una norma en Rn , entonces la funci´on d : Rn × Rn → R definida como

5.2. Algoritmos Basados en T´ ecnicas Iterativas

59

d(x, y) = kx − yk es una distancia. Definici´ on 5.2.3 Sea An el conjunto de todas las matrices cuadradas de tama˜ no n × n. Una norma sobre An , es una funci´on k • k : An → R tal que para cada A ∈ An : 1. kAk ≥ 0. 2. kAk = 0 si y s´olo si A es la matriz cero de tama˜ no n × n. 3. kcAk = |c|kAk 4. kA + Bk ≤ kAk + kBk 5. kABk ≤ kAkkBk Enseguida se presentan los algoritmos iterativos de Jacobi y de Gauss-Seidel que manejan la norma de vectores.

5.2.2

Algoritmo de Jacobi

Dado un sistema de ecuaciones como en 5.1, primero se despeja la variable xi en la respectiva ecuaci´on Ei obteniendo el siguiente sistema de ecuaciones equivalente: 12 x2 + · · · + aa1n xn + ab111 x1 = aa11 11 21 x2 = aa22 x1 + · · · + aa2n xn + ab222 22 ··· an1 an2 xn = ann x1 + ann x2 + · · · an,n−1 xn−1 + ann

(5.2)

bn ann

enseguida se propone una posible soluci´on al sistema de ecuaciones como en 5.2. Sea x0 = (x01 , x02 , . . . , x0n ) una propuesta de soluci´on. Ahora se sustituyen los valores de x0 en la ecuaci´on 5.2 obteniendo un nuevo vector de soluciones x1 = (x11 , x12 , . . . , x1n ) a partir de x0 de la siguiente manera: + · · · + aa1n x0 + ab111 11 n + · · · + aa2n x0 + ab222 22 n ··· an1 0 an2 0 1 xn = ann x1 + ann x2 + · · · an,n−1 x0n−1 + ann x11 = x12 =

a12 0 x a11 2 a21 0 x a22 1

bn ann

Enseguida se sustituyen los valores del vector x1 en la anterior ecuaci´on para obtener un nuevo vector de soluciones x2 de la siguiente manera

5.2. Algoritmos Basados en T´ ecnicas Iterativas + · · · + aa1n x1 + ab111 11 n + · · · + aa2n x1 + ab222 22 n ··· an1 1 an2 1 2 xn = ann x1 + ann x2 + · · · an,n−1 x1n−1 + ann x21 = x22 =

y as´ı sucesivamente hasta que asignada previamente.

60

a12 1 x a11 2 a21 1 x a22 1

kxk −xk−1 k kxk k

bn ann

< T OL para alguna k y para alguna tolerancia T OL

Ejemplo 5.2.2 Consideremos el siguiente sistema de ecuaciones 4x1 − x2 + x3 = 7 4x1 − 8x2 + x3 = −21 −2x1 + x2 + 5x3 = 15 de aqu´ı, despejamos cada variable xi de la ecuaci´on Ei para obtener una ecuaci´on de la forma 5.2 x1 = 41 x2 − 14 x3 + 74 x2 = 12 x1 + 18 x3 + 21 8 x3 = 25 x1 − 51 x2 + 3 Proponemos un vector de soluci´on. Por ejemplo, sea x0 = (1, 2, 2), sustituimos en la ecuaci´on anterior para obtener un vector de soluci´on x1 = (x11 , x12 , x13 ) con x11 = 14 x02 − 14 x03 + 74 = 1.75 x12 = 12 x01 + 81 x03 + 21 = 3.375 8 2 0 1 0 1 x3 = 5 x1 − 5 x2 + 3 = 3.00 Nuestro nuevo vector de soluciones es x1 = (1.75, 3.375, 3.00), enseguida se repite el procedimiento anterior tomando ahora en cuenta el nuevo vector soluci´on x1 para obtener un nuevo vector soluci´on x2 . La siguiente tabla muestra la sucesi´on de vectores soluci´on {xk }19 on. k=1 que es donde se estabiliza la soluci´

5.2. Algoritmos Basados en T´ ecnicas Iterativas

61

xk1 xk2 xk3 1.0 2.0 2.0 1.75 3.375 3.0 1.84375 3.875 3.025 1.9625 3.925 2.9625 1.99062500 3.97656250 3.00000000 1.99414063 3.99531250 3.00093750 .. .. .. . . . 15 1.99999993 3.99999985 2.99999993 .. .. .. .. . . . . 19 2.00000000 4.00000000 3.00000000 k 0 1 2 3 4 5 .. .

de donde el vector soluci´on para el sistema de ecuaciones es x19 = (2, 4, 3). X

5.2.3

Algoritmo de Gauss-Seidel

El algoritmo iterativo de Gauss-Seidel es muy parecido al algoritmo de Jacobi, para ello, retomemos la ecuaci´on 5.2 y consideremos nuevamente x1 = aa12 x2 + · · · + aa1n xn + ab111 11 11 21 x1 + · · · + aa2n xn + ab222 x2 = aa22 22 ··· n2 n1 x1 + aann x2 + · · · an,n−1 xn−1 + xn = aann ann

(5.2)

bn ann

Nuevamente se propone un vector soluci´on x0 = (x01 , x02 , . . . , x0n ) y se sustituye en la primera ecuaci´on del anterior sistema de ecuaciones, aquella donde, se encuentra despejada la variable x11 x11 =

a12 0 x a11 2

+ ··· +

a1n 0 x a11 n

+

b1 a11

enseguida se actualiza el vector x0 teniendo como primera entrada x11 de la operaci´on anterior, esto es, x0 = (x11 , x02 , . . . , x0n ), se sustituyen los valores necesarios para encontrar el valor de x12 del sistema de ecuaciones, esto es, se tome el valor de la ecuaci´on anterior, donde se encuentra despejada la variable x12 , esto es: x12 =

a21 1 x a22 1

+ ··· +

a2n 0 x a22 n

+

b2 a22

actualizamos las entradas del vector x0 , teni´endose x0 = (x11 , x12 , x03 , . . . , x0n ), se sigue con el procedimiento anterior hasta que se hayan obtenido todas las entradas del vector de soluci´on correspondiente x1 . Despu´es, se toma a x1 como vector soluci´on y se repite el proceso

5.2. Algoritmos Basados en T´ ecnicas Iterativas

62

anterior para encontrar el valor del vector de soluci´on x2 y as´ı sucesivamente hasta que kxk −xk−1 k < T OL para alguna k y para alguna tolerancia T OL asignada previamente. kxk k Ejemplo 5.2.3 Consideremos nuevamente el sistema de ecuaciones del ejemplo 5.2.2. Consideremos el despeje de la variable xi de la ecuaci´on Ei x1 = 14 x2 − 14 x3 + 74 x2 = 12 x1 + 18 x3 + 21 8 x3 = 25 x1 − 51 x2 + 3 proponemos nuevamente el vector soluci´on x0 = (1, 2, 2). Sustituimos los valores respectivos en la primera ecuaci´on para obtener x11 = 14 x02 − 14 x03 +

7 4

= 1.75,

actualizamos nuestro vector x0 = (1.75, 2, 2) y sustituimos en la ecuaci´on dos de nuestro sistema de ecuaciones obteniendo x12 = 21 x11 + 18 x03 +

21 8

= 3.75

se actualiza nuevamente el vector soluci´on, teniendo x0 = (1.75, 3.75, 2) y finalmente se sustituyen los valores respectivos en la tercera ecuaci´on para obtener x13 = 52 x11 − 15 x12 + 3 = 2.95 en este momento todas las entradas del vector soluci´on forman parte del vector de solcui´on x1 , esto es x0 = (1.75, 3.75, 2.95). Se continua con el mismo procedimiento para obtener las entradas del vector soluci´on x2 y as´ı sucesivamente.La siguiente tabla muestra la sucesi´on de vectores soluci´on {xk }10 on. k=1 que es donde se estabiliza la soluci´ k 0 1 2 3 .. . 8 9 10

xk1 1.0 1.75 1.95 1.995625 .. .

xk2 xk3 2.0 2.0 3.37 2.95 3.96875 2.98625 3.999609375 2.99903125 .. .. . . 1.99999983 3.99999988 2.99999996 1.99999998 3.99999999 3.00000000 2.00000000 4.00000000 3.00000000

de donde el vector soluci´on para el sistema de ecuaciones es x10 = (2, 4, 3). X

63

Cap´ıtulo 6 Algoritmos para dar Soluci´ on Num´ erica a Ecuaciones Diferenciales Ordinarias En lo que sigue, se muestran algoritmos para la resoluc´on de ecuaciones diferenciales ordinarias por medio de ecuaciones diferenciales. Para ello, consideramos el problema de encontrar una soluci´on al problema de valor inicial dy dt

6.1 6.1.1

= f (t, y), con a ≤ t ≤ b y y(a) = α

Algoritmos por M´ etodos de Taylor M´ etodo de Euler

Supongamos que se tiene el problema con valor inicial y 0 = f (t, y), con a ≤ t ≤ b y y(a) = α en n n´ umeros uniformemente espaciados en el intervalo [a, b]. Para resolver el anterior sistema de ecuaciones, tomemos los valores h = (b − a)/n t0 = a w0 = α

6.1. Algoritmos por M´ etodos de Taylor

64

Puesto que el problema se dividi´o en n puntos uniformemente espaciados en el intervalo [a, b], debemos calcular una pareja de puntos (ti , wi ) para cada punto de la partici´on. Para esto, para cada valor 1 ≤ i ≤ n se hace wi = wi−1 + hf (ti−1 , wi−1 ) ti = a + hi En cada paso se muestra el valor de la pareja (ti , wi ). Ejemplo 6.1.1 Aplicar el m´etodo de Euler para aproximar soluciones del siguiente problema de valor inicial y 0 = 1 + (t − y)2 , con 2 ≤ t ≤ 3, y(2) = 1 y h = 0.5. Puesto que h = 0.5 y 2 ≤ t ≤ 3, se tiene n = (b − a)/h = (3 − 2)/0.5 = 2. Se tiene que t0 = a = 2 y w0 = α = 1; y ti = a+hi, wi = wi−1 +hf (ti−1 , wi−1 ) = wi−1 +h(1+(ti−1 −wi−1 )2 ) para i = 1, 2. La siguiente tabla muestra los valores de las parejas (ti , wi ) ti wi 2 1 2.5 2 3 2.625 X

6.1.2

M´ etodo de Taylor de Orden N

Nuevamente se tiene el problema con valor inicial y 0 = f (t, y), con a ≤ t ≤ b y y(a) = α en n n´ umeros uniformemente espaciados en el intervalo [a, b]. Para cada i, ti se calcula de igual manera que en el m´etodo de Euler. Los valores wi se calculan de la siguiente manera w0 = α wi+1 = wi + hT N (ti , wi ) para i = 0, 1, . . . , n − 1 y n−1 N T (ti , wi ) = f (ti , wi ) + h2 f 0 (ti , wi ) + · · · + h n! f n−1 (ti , wi )

6.2. Algoritmo de Runge-Kutta

65

Ejemplo 6.1.2 Apliquemos el m´etodo de Taylor de orden dos para resolver el problema de valor inicial del ejercicio 6.1.1. Se tiene que t0 = 2 y w0 = 1. Luego T N =2 (t0 , w0 ) = f (t0 , w0 ) + h2 f 0 (t0 , w0 ) = 1 + (t0 − w0 )2 + h (2(t − y)(1 − y 0 )) = 1 + (t0 − w0 )2 + h2 (2(t − y)(1 − (1 + (t0 − w0 )2 ))). Sustituyendo se tiene 2 que T 2 (t0 , w0 ) = T 2 (2, 1) = 1.5, sustituyendo en el respectivo valor para w1 , se tiene que w1 = 1.75. La siguiente tabla muestra la pareja de valores (ti , wi ) ti wi 2 1 2.5 1.75 3 2.425 X

6.2 6.2.1

Algoritmo de Runge-Kutta Algoritmo de Runge-Kutta

Para resolver el problema con valor inicial y 0 = f (t, y), con a ≤ t ≤ b y y(a) = α con el m´etodo de Runge-Kutta de orden 2 y orden 4, se divide el intervalo [a, b] en n n´ umeros uniformemente espaciados. Hay dos variantes para el m´etodo de Runge-Kutta de orden 2, uno de los cuales es conocido como m´etodo modificado de Euler y el otro como m´etodo de Heun. Las f´ormulas son muy parecidas a las de la secci´on anterior, de est´a manera solo se presentar´an las f´ormulas de los m´etodos. Para el m´etodo modificado de Euler (Runge-Kutta de orden 2) se tienen las siguientes f´ormulas

wi+1 = wi +

h [f (ti , wi ) 2

w0 = α + f (ti+1 , wi + hf (ti , wi ))] para i = 0, 1, . . . , n − 1

Para el m´etodo de Heun (Runge-Kutta de orden 2) se tienen las siguientes f´ormulas

6.2. Algoritmo de Runge-Kutta

w0 = α wi+1 = wi + h4 [f (ti , wi + 3f (ti + 32 h, wi + 32 hf (ti , wi )))] para i = 0, 1, . . . , n − 1 Finalmente se muestran las f´ormulas para el m´etod de Runge-Kutta de orden cuatro.

wi+1

w0 = α k1 = hf (ti , wi ) k2 = hf (ti + h2 , wi + 12 k1 ) k3 = hf (ti + h2 , wi + 12 k2 ) k4 = hf (ti+1 , wi + k3 ) 1 = wi + 6 (k1 + 2k2 + 2k3 + k4 ) para i = 0, 1, . . . , n − 1.

66

67

Bibliograf´ıa [1] Richard L. Burden y J.Douglas Faires. An´alisis Num´erico. Thomson Learning. 7a edici´on. 2002. [2] Steven C. Chapra y Raymond P. Canale. M´etodods Num´ericos para Ingenieros. Mc Graw Hill. Quinta edici´on. 2007. [3] Curtis F. Gerald y patrick O. Wheatley. An´alisis Num´erico con Aplicaciones. Pearson Education. Sexta edici´on. M´exico 2000. [4] A. Nieves y F. Dominguez. M´etodos Num´ericos Aplicados a la Ingenier´ıa. CECSA. M´exico. [5] William Stallings. Organizaci´on y Arquitectura de Computadores. Dise˜ no para Optimizar Prestaciones. Prentice Hall. Madrid 2000.