• Author / Uploaded
• Juan Carlos Retuerto Aramburu

Citation preview

C´alculo Num´erico con aplicaciones para las ciencias e ingenier´ıa Angel Ramirez - Irla Mantilla 06 de Abril de 2015

2

´Indice general 1. B´ usqueda de Ra´ıces en Ecuaciones No Lineales 1.1. M´etodo de la bisecci´on . . . . . . . . . . . . . . . 1.1.1. An´alisis de convergencia . . . . . . . . . . 1.1.2. An´alisis de error . . . . . . . . . . . . . . 1.2. M´etodo de la falsa posici´on . . . . . . . . . . . . . 1.2.1. An´alisis de convergencia . . . . . . . . . . 1.2.2. An´alisis de Error . . . . . . . . . . . . . . 1.2.3. M´etodo de la falsa posici´on modificada . . 1.3. M´etodo de Newton . . . . . . . . . . . . . . . . . 1.3.1. An´alisis de convergencia . . . . . . . . . . 1.3.2. An´alisis de error . . . . . . . . . . . . . . 1.4. M´etodo de la secante . . . . . . . . . . . . . . . . 1.5. M´etodo del punto fijo . . . . . . . . . . . . . . . . 1.5.1. An´alisis de convergencia . . . . . . . . . . 1.5.2. An´alisis de error . . . . . . . . . . . . . . 1.6. M´etodos para sistemas de ecuaciones no lineales . 1.6.1. M´etodo de Newton . . . . . . . . . . . . . 1.6.2. M´etodo del punto fijo . . . . . . . . . . . . 1.7. Ejercicios . . . . . . . . . . . . . . . . . . . . . . 2. Interpolaci´ on polinomial 2.1. M´etodo de coeficientes indeterminados 2.2. M´etodo de Lagrange . . . . . . . . . . 2.3. M´etodo de Newton . . . . . . . . . . . 2.4. M´etodo de Diferencias Divididas . . . . 2.5. An´alisis de error . . . . . . . . . . . . . 2.6. Fen´omeno Runge . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . . . . .

7 8 10 13 15 17 19 20 23 26 28 29 31 32 35 40 40 40 40

. . . . . .

43 44 44 46 48 53 56

3. Integraci´ on Num´ erica 59 3.1. M´etodo del Trapecio . . . . . . . . . . . . . . . . . . . . . . . 60 3.1.1. Simple . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3

´INDICE GENERAL

4 3.1.2. Compuesto . . . . . . . 3.2. M´etodo de Simpson . . . . . . . 3.2.1. Simple . . . . . . . . . . 3.2.2. Compuesto . . . . . . . 3.3. Cuadratura Gaussiana . . . . . 3.3.1. Polinomios Ortogonales . 3.4. Ejercicios Propuestos . . . . . . 4. Sistemas de ecuaciones lineales 4.1. M´etodos Directos . . . . . . . . 4.1.1. Eliminaci´on Gaussiana . 4.1.2. Factorizaci´on LU . . . . 4.1.3. Factorizaci´on QR . . . . 4.2. M´etodos Indirectos . . . . . . . 4.2.1. M´etodo de Jacobi . . . . 4.2.2. M´etodo de Gauss-Seidel 4.2.3. M´etodo de relajaci´on . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

. . . . . . . .

. . . . . . .

61 65 65 69 73 73 74

. . . . . . . .

77 77 77 77 79 82 82 82 82

5. M´ etodos num´ ericos para problemas de valor inicial 83 5.1. M´etodos de Taylor . . . . . . . . . . . . . . . . . . . . . . . . 87 5.2. Ejercicios Propuestos . . . . . . . . . . . . . . . . . . . . . . . 89

Agradecimientos No podr´ıa haber escrito este pedazo de libro en espa˜ nol sin la inestimable ayuda de ...

6

´INDICE GENERAL

Cap´ıtulo 1 B´ usqueda de Ra´ıces en Ecuaciones No Lineales Esta cap´ıtulo se centra en el problema de encontrar un cero de una funci´on real de variable real. En primer lugar, algunos de los aspectos te´oricos y num´ericos de encontrar ceros en R pueden ser convenientes para el mismo problema para funciones de Rn en R. Por lo tanto, nuestro problema es hallar x tal que: f (x) = 0

(1.1)

donde f : Rn → R es una funci´on no lineal por lo general, veremos que condiciones debe de satisfacer f para aproximar una ra´ız de ella. Es com´ un encontrarnos con funciones f : D ⊂ R → R no lineales para las cuales deseamos encontrar sus ra´ıces, que normalmente, o son dif´ıciles o incluso hasta imposible calcularlas de forma anal´ıtica y exacta. Con el objetivo de paliar este problema, los m´etodos num´ericos a estudiar en este cap´ıtulo, nos ayudar´an a encontrar la ra´ız de la funci´on de una forma tan aproximada como lo deseemos, para ello, se estudiar´an m´etodos num´ericos adecuados y su respectiva implementaci´on en un computador. Entonces, antes de buscar las ra´ıces de una funci´on dada, primero debemos de garantizar la existencia de la ra´ız as´ı como su unicidad en un intervalo dado, para que as´ı puedan ser aplicados algunos de los m´etodos que se estudiar´an a continuaci´on. Por tanto, recordemos el siguiente teorema b´asico del c´alculo diferencial. Teorema 1.1 Sea f : [a, b] ⊂ R → R una funci´on continua la cual satisface que f (a)f (b) < 0, entonces existe c ∈ ha, bi tal que f (c) = 0. El teorema anterior garantiza la existencia de al menos una ra´ız para la funci´on, para garantizar la unicidad debemos de utilizar alg´ un otro argumento, 7

´ 8CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES por ejemplo, que la funci´on sea mon´otona o que tenga derivada positiva en el intervalo. Un uso adecuado del Teorema 1.1 permitir´a deducir el M´etodo de la Bisecci´on que se detalla a continuaci´on.

1.1.

M´ etodo de la bisecci´ on

Este m´etodo surge de forma natural al aplicar adecuadamente el Teorema 1.1, pues, b´asicamente consiste en reducir el intervalo inicial a la mitad en cada iteraci´on, de modo que se consigue aislar la ra´ız en cada nuevo intervalo. La elecci´on del nuevo intervalo se basa en el Teorema 1.1. Es decir: Dada una funci´on f y dado un intervalo inicial [a0 , b0 ] que satisfacen las con0 , obteniendo diciones del Teorema 1.1, calculamos el punto medio c0 = a0 +b 2 as´ı dos nuevos subintervalos I1 = [a0 , c0 ] e I2 = [c0 , b0 ]. El nuevo intervalo [a1 , b1 ] ser´a aquel que siga cumpliendo la condici´on del Teorema 1.1. Repetimos el proceso hasta la iteraci´on n y obtenemos un intervalo final [an , bn ] m´as peque˜ no que el original que contiene a la ra´ız.

Figura 1.1: Esquema del m´etodo de la bisecci´on

A continuaci´on presentamos el Algoritmo del m´etodo de la bisecci´on:

Algoritmo 1 (M´ etodo de la bisecci´ on.)

´ ´ 1.1. METODO DE LA BISECCION

9

Paso 1: Ingresar una funci´on continua f . Paso 2: Ingresar intervalo inicial [a0 , b0 ]. Paso 3: Verificar que f (a0 )f (b0 ) < 0, caso contrario regresar al Paso 2. a0 + b 0 y f (c0 ). Paso 4: Calcular c0 = 2 Paso 5: Hacer k = 0. Si f (a0 )f (c0 ) < 0 entonces: ak+1 = ak , bk+1 = ck y f (bk+1 ) = f (ck ), sino: ak+1 = ck , bk+1 = bk y f (ak+1 ) = f (ck ). Y el nuevo intervalo ser´a [ak+1 , bk+1 ]. Paso 6: Mientras |ck+1 − ck | > T olerancia hacer ak+1 + bk+1 ck+1 = 2 Calcular f (ck+1 ). k = k + 1; Si f (a0 )f (c0 ) < 0 entonces: ak+1 = ak , bk+1 = ck y f (bk+1 ) = f (ck ), sino: ak+1 = ck , bk+1 = bk y f (ak+1 ) = f (ck ). Y el nuevo intervalo ser´a [ak+1 , bk+1 ]. Paso 7: La aproximaci´on de la ra´ız viene dada por: ak+1 + bk+1 . ck+1 = 2 Ejemplo: Aproximar una ra´ız para la funci´on f (x) = 2 arctan(x − 3) − 0.02x2 . Resoluci´ on Primero observamos que la funci´on f es continua en todo R. Considerando el intervalo [0, 4], se observa que f (0)f (4) < 0, entonces, por el Teorema 1.1 podemos concluir que existe al menos una ra´ız r de f en el intervalo [0, 4]. Para garantizar unicidad, observemos que: f 0 (x) =

2 − 0.04x. 1 + (x − 3)2

Pero, para todo x ∈ [0, 4] se cumplen: 0.2 ≤

2 ≤ 2, 1 + (x − 3)2

0 < 0.2 − 0.16 ≤ 0.2 − 0.04x,

´ 10CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES es decir: f 0 (x) > 0 para todo x ∈ [0, 4], por tanto, la funci´on f es creciente y as´ı, la ra´ız es u ´nica en el intervalo inicial [0, 4]. Garantizada la existencia y unicidad de una ra´ız r en [0, 4], procedemos a aproximarla mediante el m´etodo de la bisecci´on, los cuales se muestran en la Tabla 1.1. Tabla 1.1: M´et. de la bisecci´on para f (x) = 2 arctan(x − 3) − 0.02x2 . i

ai

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0 2 3 3 3 3 3.0625 3.09375 3.09375 3.09375 3.09375 3.09570313 3.09570313 3.09570313 3.09594727 3.09606934 3.09613037 3.09613037

bi

ci

f (c)

4 2 -1.6507963 4 3 -0.1800000 4 3.5 0.6822952 3.5 3.25 0.2787073 3.25 3.125 0.0533975 3.125 3.0625 -0.0627405 3.125 3.09375 -0.0044722 3.125 3.109375 0.0245197 3.109375 3.1015625 0.0100371 3.1015625 3.09765625 0.0027857 3.09765625 3.09570313 -0.0008425 3.09765625 3.09667969 0.0009718 3.09667969 3.09619141 0.0000647 3.09619141 3.09594727 -0.0003889 3.09619141 3.09606934 -0.0001621 3.09619141 3.09613037 -0.0000487 3.09619141 3.09616089 0.0000080 3.09616089 3.09614563 -0.0000203

|ci+1 − ci | < 10−5 F F F F F F F F F F F F F F F F F E

De la Tabla 1.1 podemos observar que: |c17 − c16 | = 0.00000203 < 10−5 , entonces, tenemos que: r ≈ c17 = 3.09614563.

1.1.1.

An´ alisis de convergencia

Hasta ahora, s´olo hemos aplicado el m´etodo de la bisecci´on de forma mec´anica y pr´actica, pero tal vez a´ un no nos hemos preguntado si el m´etodo de la bisecci´on funciona siempre. Es objetivo de esta secci´on confirmar y demostrar el porqu´e funciona este m´etodo.

´ ´ 1.1. METODO DE LA BISECCION

11

Figura 1.2: M´etodo gr´afico para estimar el valor de la ra´ız r. Recordando que el m´etodo de la bisecci´on en cada iteraci´on reduce el intervalo a la mitad, entonces, en la iteraci´on n se cumple que: 1 |bn − an | = |bn−1 − an−1 | 2 lo anterior implica que: |bn − an | =

1 |b0 − a0 | 2n

cuando n → ∞ se tiene que l´ım |bn − an | = 0, as´ı, podemos concluir que: n→∞

l´ım (bn − an ) = 0

n→∞

(1.2)

Por otro lado, observamos que los intervalos generados en cada iteraci´on satisfacen la siguiente relaci´on: a0 ≤ a1 ≤ a2 ≤ · · · an ≤ b n ≤ · · · ≤ b 3 ≤ b 2 ≤ b 1 ≤ b 0 es decir, la sucesi´on {an }∞ on creciente y acotada por b0 , por n=0 es una sucesi´ tanto, es una sucesi´on convergente. Del mismo modo, la sucesi´on {bn }∞ n=0 es

´ 12CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES

Figura 1.3: Este gr´afico es nuevo una sucesi´on decreciente y acotada por a0 , entonces, tambi´en es una sucesi´on convergente. As´ı, a partir de la ecuaci´on (1.2) se obtiene: 0 = l´ım (bn − an ) = l´ım bn − l´ım an n→∞

n→∞

n→∞

es decir: l´ım bn = l´ım an .

n→∞

(1.3)

n→∞

La ecuaci´on (1.4) nos dice que los extremos de los intervalos [an , bn ] van a coincidir en el l´ımite, sea r este l´ımite com´ un. Por otro lado, como f es una funci´on continua, se cumple tambi´en lo siguiente: l´ım f (an ) = f ( l´ım an ) = f (r) = f ( l´ım bn ) = l´ım f (bn ).

n→∞

n→∞

n→∞

n→∞

Recordemos que los extremos del intervalo en el m´etodo de la bisecci´on son elegidos tal que f (an )f (bn ) < 0, entonces, tomando nuevamente l´ımite obtenemos: 0 ≥ l´ım f (an )f (bn ) = l´ım f (an ) l´ım f (bn ) = f (r)f (r) ≥ 0. n→∞

n→∞

n→∞

de donde obtenemos: f (r) = 0, es decir, r es ra´ız de la funci´on f . Esto nos permite concluir que el m´etodo de la bisecci´on siempre converge a una ra´ız de una funci´on continua.

´ ´ 1.1. METODO DE LA BISECCION

1.1.2.

13

An´ alisis de error

El an´alisis de convergencia realizado en la secci´on anterior. nos dice que despu´es de realizar un n´ umero infinito de veces el m´etodo de la bisecci´on, este nos dar´a el valor exacto de la ra´ız r de una funci´on coninua f , siempre que el intervalo inicial [a, b] satisfaga f (a)f (b) < 0. Lamentablemente un proceso infinito no puede ser realizado por ning´ un computador y mucho menos hacerlo manualmente. Lo que garantiza el an´alisis de convergencia y podemos rescatar, es que cuando mayor n´ umero de iteraciones realicemos, la aproximaci´on xn estar´a m´as cerca de la ra´ız, entonces es natural la siguiente pregunta: ¿Ser´a posible determinar apriori el n´ umero de iteraciones de forma que alcancemos una tolerancia de error ingresada al inicio del proceso?. Veamos que la pregunta anterior tiene una respuesta afirmativa, para ello consideremos una funci´on f y un intevalo inicial [a0 , b0 ] que satisfacen las condiciones del Teorema 1.1. Sea r la ra´ız tal que r ∈ [a0 , b0 ], despu´es de n iteraciones se ha reducido el intervalo inicial al intervalo [an , bn ] y la aproxian + b n . Por tanto, el error en la maci´on a la ra´ız r en ese intervalo es cn = 2 iteraci´on n viene dado por: en = |r − cn | (1.4) Por ser cn el punto medio del intervalo [an , bn ], el cual contiene a la ra´ız r, podemos concluir que: 1 1 1 |r − cn | ≤ (bn − an ) = 2 (bn−1 − an−1 ) = ... = n+1 (b0 − a0 ) 2 2 2 Dada una tolerancia inicial tol, tal que deseamos que en < tol, entonces, bastar´a exigir que: 1 (b0 − a0 ) < tol, n+1 2 tomando logaritmos para despejar n, se obtiene:   b 0 − a0 log tol n> − 1. (1.5) log(2) As´ı, la relaci´on (1.5) nos permite determinar apriori un n´ umero m´ınimo de iteraciones a realizar para alcanzar una tolerancia de error deseada. Ejemplo 1.1.1. Sea f (x) = x3 − 2x−1 + 3. ¿Cu´antas iteraciones se necesita para aproximar una ra´ız con una tolerancia de error de 10−5 ?. Resoluci´ on

´ 14CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES La funci´on f tiene como dominio R\{0}. Adem´as: f (−2) = −4 y f (−1) = 4, es decir, f (−1)f (−4) < 0 y como la funci´on f es continua en [−2, −1], entonces, por el Teorema 1.1 concluimos que existe una ra´ız r ∈ h−2, −1i. Para garantizar la unicidad, vemos que: f 0 (x) = 3x2 +

2 > 0 para todo x ∈ [−2, −1], x2

es decir, f es una funci´on creciente en ese intervalo. As´ı, podemos garantizar que existe una u ´nica ra´ız r ∈ h−2, −1i. Para calcular el n´ umero de iteraciones necesarias para obtener una tolerancia −5 de error de 10 , hacemos uso de la relaci´on (1.5), donde a0 = −2, b0 = −1, es decir: 5 − 1 ≈ 15.6 n > log(2) es decir, debemos tomar n = 16 iteraciones. Los resultados pueden observarse en la Tabla 1.2. Tabla 1.2: Met. de la bisecci´on para f (x) = x3 − 2x−1 + 3. i

ai

0 -2 1 -2 2 -1.75 3 -1.625 4 -1.625 5 -1.625 6 -1.625 7 -1.625 8 -1.621094 9 -1.619141 10 -1.618164 11 -1.618164 12 -1.618164 13 -1.618042 14 -1.618042 15 -1.618042 16 -1.618042

bi

ci

f (c)

-1 -1.5 -1.5 -1.5 -1.5625 -1.59375 -1.609375 -1.617188 -1.617188 -1.617188 -1.617188 -1.617676 -1.61792 -1.61792 -1.617981 -1.618011 -1.618027

-1.5 -1.75 -1.625 -1.5625 -1.59375 -1.609375 -1.617188 -1.621094 -1.619141 -1.618164 -1.617676 -1.61792 -1.618042 -1.617981 -1.618011 -1.618027 -1.618034

0.958333 -1.216518 -0.060246 0.465303 0.206715 0.074296 0.007292 -0.02641 -0.009542 -0.001121 0.003086 0.000983 -0.000069 0.000457 0.000194 0.000063 -0.000003

´ ´ 1.2. METODO DE LA FALSA POSICION

1.2.

15

M´ etodo de la falsa posici´ on

Dada una funci´on f : [a, b] → R continua en [a, b] y tal que f (a)f (b) < 0. Iniciando en el intervalo [x0 , x1 ] = [a, b], la idea es generar una sucesi´on de intervalos [xk−1 , xk ] (k ≥ 1) tal que la ra´ız de f se encuentre en ´el. Con este objetivo, nada nos impide considerar la recta secante que pasa por los puntos (xk−1 , f (xk−1 )) y (xk , f (xk )) dada por y − f (xk ) =

f (xk−1 ) − f (xk ) (x − xk ) xk−1 − xk

(1.6)

Figura 1.4: M´etodo gr´afico para estimar el valor de la ra´ız r. De la Figura 1.4 observamos que la intersecci´on de la recta secante con el eje x es el punto (xk+1 , 0), que reemplazando en (1.6) nos da el valor de xk+1 : xk−1 − xk xk+1 = xk − f (xk ) (1.7) f (xk−1 ) − f (xk ) Para elegir el nuevo intervalo [xk−1 , xk ], seguimos la idea del m´etodo de la bisecci´on, es decir, elegimos aquel que satisface f (xk−1 )f (xk ) < 0

´ 16CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES y la nueva aproximaci´on de la ra´ız xˆ es calculada de forma an´aloga al c´alculo de xk+1 . El proceso termina cuando la aproximaci´on xk+1 satisfaga la condici´on f (xk+1 ) = 0 o |xk+1 − xk | < Tolerancia. Algoritmo 1: Algoritmo del M´etodo de la Secante Entrada: Intervalo [a, b] tal que f (a)f (b) < 0. Tolerancia: tol. 1 inicio 2 k ← 1; 3 ak ← a; 4 bk ← b; b k − ak ; 5 xk ← ak − f (ak ) f (bk ) − f (ak ) 6 mientras |f (xk )| > tol hacer 7 si f (ak )f (xk ) < 0 entonces 8 ak+1 ← ak ; 9 bk+1 ← xk ; 10 en otro caso 11 ak+1 ← xk ; 12 bk+1 ← bk . 13 fin si 14 k ← k + 1; b k − ak . 15 xk ← ak − f (ak ) f (bk ) − f (ak ) 16 fin mientras 17 devolver xk . 18 fin Ejemplo 1.2.1. Sea f (x) = x2 ln(x) − x. D´e un intervalo donde exista una ra´ız de la funci´on f y aprox´ımela con una tolerancia de error de 10−5 . Resoluci´ on La funci´on f tiene como dominio x > 0. Adem´as: f (0.5) = −0.6931472 y f (2) = 0.6931472, es decir, f (0.5)f (2) < 0 y como la funci´on f es continua en [0.5, 2], entonces, por el Teorema 1.1 concluimos que existe una ra´ız r ∈ h0.5, 2i. Para garantizar la unicidad, para todo x ∈ [0.5, 2] se cumple que: f 0 (x) = 2x ln(x) + x − 1, f 00 (x) = 2 ln(x) + 3. El punto cr´ıtico de f es x = 1 donde f 0 (1) = 0. La Tabla 1.3 muestra el an´alisis de signo para esbozar la gr´afica de la funci´on f dada en la Figura 1.5.

´ ´ 1.2. METODO DE LA FALSA POSICION

17

Tabla 1.3: An´alisis de signo para f = x2 ln(x) − x. x1 f 0 (x) + 00 f (x) + + 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Figura 1.5: f (x) = x2 ln(x) − x. La Figura 1.5 garantiza que f tiene una u ´nica ra´ız en el intervalo [0.5, 2], por tanto, inicializando en [a0 , b0 ] = [0.5, 2] la Tabla 1.4 muestra las iteraciones hasta que se cumple que |f (x)| ≤ 10−5 .

1.2.1.

An´ alisis de convergencia

En la iteraci´on k + 1 el error para el m´etodo de la regla falsa se obtiene xk−1 − xk − xˆ (1.8) ek+1 = xk+1 − xˆ = xk − f (xk ) f (xk−1 ) − f (xk ) reordenando (1.8) 



 f (xk ) − f (ˆ x)   ek+1 = xk − xˆ −   f (xk−1 ) − f (xk )  . xk−1 − xk

(1.9)

´ 18CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES Tabla 1.4: Met. de la Falsa Posici´on para f (x) = x2 log(x) − x. i

ai

bi

ci

f (ci )

0 1 2 3 4 5 6 7 8

0.500000 0.500000 1.198490 1.638086 1.742621 1.760031 1.762733 1.763148 1.763211

2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000 2.000000

1.198490 1.198490 1.638086 1.742621 1.760031 1.762733 1.763148 1.763211 1.763221

-0.938416 -0.938416 -0.313789 -0.056052 -0.008799 -0.001353 -0.000207 -0.000032 -0.000005

Hacemos uso de la siguiente notaci´on f [a, b] :=

f (b) − f (a) b−a

para reordenar la expresi´on dada en (1.9)     f [xk , xˆ] f [xk−1 , xk ] − f [xk , xˆ] ek+1 = ek 1 − = ek f [xk−1 , xk ] f [xk−1 , xk ] ahora acomodamos convenientemente (1.11) para obtener   f [xk−1 , xk , xˆ] ek+1 = ek ek−1 f [xk−1 , xk ]

(1.10)

(1.11)

(1.12)

donde se ha usado la notaci´on siguiente f [a, b, c] :=

f [a, b] − f [b, c] . a−c

(1.13)

Usando (1.10) y (1.13) para xk−1 , xk , xˆ (en ese orden), podemos definir el siguiente polinomio P2 (x) = f (xk−1 ) + f [xk−1 , xk ](x − xk−1 ) + f [xk−1 , xk , xˆ](x − xk−1 )(x − xk ) el cual satisface P2 (xk−1 ) = f (xk−1 ), P2 (xk ) = f (xk ) y P2 (ˆ x) = f (ˆ x). Definimos la funci´on g : [xk−1 , xk ] → R como sigue g(x) = f (x) − P2 (x).

(1.14)

´ ´ 1.2. METODO DE LA FALSA POSICION

19

De (1.14) se observa g(xk−1 ) = g(xk ) = g(ˆ x) = 0, por tanto, se puede aplicar el Teorema de Rolle para concluir que existe ηk ∈ hxk−1 , xk i tal que g 00 (ηk ) = 0. Esto implica f 00 (ηk ) = f [xk−1 , xk , xˆ]. 2

(1.15)

Por otro lado, del Teorema del Valor Medio, existe ξk ∈ hxk−1 , xk i tal que f 0 (ξk ) = f [xk−1 , xk ].

(1.16)

Usando (1.15) y (1.16) en (1.12) se obtiene ek+1 =

f 00 (ηk ) ek ek−1 . 2f 0 (ξk )

(1.17)

Con el objetivo de acotar f 00 y f 00 , vamos a imponer que f ∈ C 2 ([x0 , x1 ]) y f 0 (x) 6= 0 para todo x ∈ [x0 , x1 ], luego, para k suficientemente grande podemos considerar lo siguiente, [1]: ek+1 = M ek ek−1 .

(1.18)

f 00 (ˆ x) . De la definici´on de ek dada en (1.8) podemos concluir a 0 2f (ˆ x) partir de (1.18) que el m´etodo de la Regla Falsa converge para puntos x0 , x1 suficientemente cerca de xˆ. donde M =

1.2.2.

An´ alisis de Error

Para el an´alisis de error, consideramos la igualdad (1.18) y definimos la sucesi´on yk = ln(M ek ), que reemplazando en (1.18) resulta yk+1 = yk + yk−1 ,

(1.19)

considerando yk = rk para resolver (1.19), se obtiene que √ √ 1+ 5 1− 5 r1 = , y r2 = = 1 − r1 . 2 2 por tanto, la soluci´on de (1.19) viene dada por yk = Ar1k + Br2k .

(1.20)

A partir de (1.20) se obtiene la siguiente relaci´on yk − r1 yk−1 = Br2k−1

(1.21)

´ 20CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES tomando l´ımite l´ım (yk − r1 yk−1 ) = l´ım Br2k−1 = 0.

k→∞

k→∞

(1.22)

usando la definici´on de yk en (1.22) y de las propiedades de logaritmos se obtiene   1−r1 ek l´ım (yk − r1 yk−1 ) = l´ım M = 0, (1.23) 1 k→∞ k→∞ erk−1 y dado que ek es una sucesi´on convergente, resulta que  l´ım

k→∞

ek 1 erk−1



= M r1 −1 .

(1.24)

El resultado obtenido en (1.24) nos permite afirmar que el m´etodo de Regla Falsa es superlineal.

1.2.3.

M´ etodo de la falsa posici´ on modificada

Se observa en la Tabla 1.4 que el extremo derecho del intervalo en cada iteraci´on permanece constante. En general, cuando la funci´on es convexa (c´oncava) creciente o decreciente, entonces al aplicar el m´etodo de falsa posici´on, uno de los extremos del intervalo inicial permanece constante. Esta caracter´ısca impide que la convergencia sea r´apida y que tomemos como criterio de parada que la longitud del intervalo [xk−1 , xk ] sea suficientemente peque˜ no para alg´ un k > 0. Con el objetivo de superar estas dificultades, se hace una modificaci´on al M´etodo de Falsa Posici´on. Cuando se observe que alg´ un extremo del intervalo permanece constante en dos iteraciones (al menos), entonces, se aplica el Algoritmo de Illinois, el cual considera un factor α de la imagen del extremo que permance constante para el c´alculo del siguiente punto. Este algoritmo es descrito en detalle en [2] y presentado a continuaci´on. El objetivo de la introducci´on del factor α es el prevenir que alg´ un extremo permanezca constante.

´ ´ 1.2. METODO DE LA FALSA POSICION

21

Algoritmo 2: M´etodo de la Falsa Posici´on Modificada Par´ ametros: α = 0.5. Entrada: Puntos x0 y x1 tal que f (x0 )f (x1 ) < 0. Tolerancias: Tol, δ. 1 inicio 2 k ← 1; x1 − x 0 ; 3 x ← x1 − f (x1 ) f (x1 ) − f (x0 ) 4 mientras |f (x)| > Tol y |x − x1 | > δ hacer 5 si f (x)f (x1 ) < 0 entonces 6 x0 ← x1 ; 7 f0 ← f (x0 ) 8 en otro caso 9 f0 ← αf (x0 ); 10 fin si 11 k ← k + 1; 12 x1 ← x; x1 − x0 13 x ← x1 − f (x1 ) f (x1 ) − f (x0 ) 14 fin mientras 15 devolver x. 16 fin Si aplicamos el m´etodo de la Regla Falsa Modificada al Ejemplo 1.2.1 se obtienen los resultados dados en la Tabla 1.5.

Tabla 1.5: M´etodo de la Regla Falsa Modificada i ai bi ci f (ci ) 0 1 2 3 4 5 6

0.500000 2.000000 2.000000 1.638086 1.800302 1.800302 1.763123

2.000000 1.198490 1.638086 1.800302 1.759541 1.763123 1.763317

1.198490 1.638086 1.800302 1.759541 1.763123 1.763317 1.763223

-0.938416 -0.313789 0.105309 -0.010146 -0.000276 0.000260 -0.000000

La interpretaci´on gr´afica del M´etodo de la Regla Falsa Modificada para el Ejemplo 1.2.1 se observa en la Figura 1.6.

´ 22CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-1 0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-1 0.4

0.6

0.8

(a) Paso 1 0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

0.6

0.8

1

1.2

1.2

1.4

1.6

1.8

2

1.4

1.6

1.8

2

(b) Paso 2

0.8

-1 0.4

1

1.4

1.6

1.8

2

-1 0.4

0.6

0.8

(c) Paso 3

1

1.2

(d) Paso 4

Figura 1.6: Interpretaci´on gr´afica del M´etodo de la Falsa Posici´on Modificada. Consideremos ahora la funci´on siguiente dada en [2] f (x) = 2xe−n + 1 − 2e−nx ,

n = 1, 2, . . .

(1.25)

En la Tabla 1.6 se resumen los n´ umeros de iteraciones al aplicar el M´etodo de la Bisecci´on (MB), M´etodo de la Regla Falsa (MRF) y el M´etodo de la Regla Falsa Modificada (MRFM) para la funci´on f dada en (1.25) para cuatro valores diferentes de n. n 1 5 15 20

MB 50 50 50 50

MRF 19 19 19 19

MRFM 8 9 11 10

Tabla 1.6: N´ umero de iteraciones. Criterio de parada: |f (x)| ≤ 10−19 .

´ 1.3. METODO DE NEWTON

1.3.

23

M´ etodo de Newton

El m´etodo de Newton es un m´etodo cl´asico y uno de los m´as r´apidos para aproximar una ra´ız de una funci´on no lineal f : [a, b] ⊂ R → R, Para deducir el m´etodo de Newton, supongamos que xn es una buena aproximaci´on de la ra´ız xˆ de f , es decir, | xˆ − xn | es bien peque˜ no, f ∈ C 2 ([a, b]) y f 0 (xn ) 6= 0; el polinomio de Taylor para f alrededor de xn viene dado por: f (x) = f (xn ) + (x − xn )f 0 (xn ) +

(x − xn )2 00 f (x ) 2

(1.26)

donde x est´a entre x y xn . Como hemos supuesto que xˆ est´a bien cerca de xn se tiene en (1.26) con x = xˆ lo siguiente: f (ˆ x) = f (xn ) + (ˆ x − xn )f 0 (xn ) +

(ˆ x − xn )2 00 f (xˆ ) 2

(1.27)

pero f (ˆ x) = 0, y desde que | xˆ − xn | es bien peque˜ no, entonces el t´ermino 2 (ˆ x − xn ) es mucho menor, de aqu´ı en (1.27): 0 ≈ f (xn ) + (ˆ x − xn )f 0 (xn )

(1.28)

de (1.28) obtenemos: xˆ ≈ xn −

f (xn ) f 0 (xn )

(1.29)

la nueva aproximaci´on dada por (1.29) ser´a el nuevo t´ermino de la sucesi´on, es decir: f (xn ) , ∀ n ≥ 0. (1.30) xn+1 = xn − 0 f (xn ) En resumen, el m´etodo de Newton comienza con una buena aproximaci´on inicial x0 y a partir de ´el se genera una sucesi´on de puntos {xn }∞ n=0 definida por (1.30). Geom´etricamente significa que empezando en el punto inicial x0 que no est´e muy lejos de la ra´ız, se efect´ ua un desplazamiento a lo largo de la tangente hacia su intersecci´on con el eje ”x”, y se toma ´este como la siguiente aproximaci´on. Esto se realiza hasta que valores de xn sucesivos est´en suficientemente pr´oximos o el valor de la funci´on est´e suficientemente cerca de cero. Esta interpretaci´on geom´etrica puede observarse en la Figura 1.7.

´ 24CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES

(a) Paso 1

(b) Paso 2

(c) Paso 3

(d) Paso 4

Figura 1.7: Interpretaci´on geom´etrica del M´etodo de Newton. La rapidez del m´etodo de Newton se ve ensombrecida debido a que es necesario conocer f 0 (x), porque en (1.30) es necesario calcular en cada iteraci´on los valores de f (xn ) y f 0 (xn ). La relaci´on (1.30) sugiere que debemos verificar que f 0 (xn ) no est´e muy pr´oxima de cero, de lo contrario el m´etodo no va a funcionar. Ejemplo 1.3.1. Aplique el m´etodo de Newton para encontrar un cero de: π f (x) = x − 0.8 − 0.2 sen(x), x ∈ [0, ] 2 Resoluci´ on: Tenemos: f (x) = x − 0.8 − 0.2 sen(x) f 0 (x) = 1 − 0.2 cos(x) π x0 = 4 y aplicando la relaci´on (1.30) codificado en el programa Newton.sce(Ver ap´endice)se obtiene:

´ 1.3. METODO DE NEWTON i 0 1 2 3

xi 0.785398 0.967121 0.964335 0.964334

25 f (xi ) -0.156023 0.002470 0.000001 0.000000

f 0 (xi ) 0.858579 0.886466 0.886007 

Ejemplo 1.3.2. Use el m´etodo de Newton en la ecuaci´on x2 = N para deducir el algoritmo para obtener la ra´ız cuadrada de un n´ umero positivo N . Resoluci´ on: Resolver x2 = N es hallar un cero de f (x) = x2 − N ⇒ f 0 (x) = 2x y por el m´etodo de Newton se obtiene:   x2n − N 1 N f (xn ) = xn − = xn + xn+1 = xn − 0 f (xn ) 2xn 2 xn 

Ejemplo 1.3.3. Hallar aproximadamente una soluci´on de la siguiente ecuaci´on: x3 + 4x2 − 10 = 0. Resoluci´ on: Para hallar una soluci´on de la ecuaci´on dada, es equivalente a encontrar una ra´ız de f (x) = x3 + 4x2 − 10, la cual es continua y diferenciable en todo R. Adem´as, observemos que: f (1) = −5,

f (2) = 14,

entonces, f (1)f (2) < 0

as´ı, el Teorema 1.1 garantiza la existencia de al menos una soluci´on en [1, 2]. Por otro lado: f 0 (x) = 3x2 + 8x > 0 para todo x ∈ [1, 2], es decir, la funci´on f : [1, 2] → R es creciente en [1, 2] y as´ı la ra´ız xˆ en ese intervalo es u ´nica. Aplicamos ahora la relaci´on (1.30) codificado en el programa Newton.sce(Ver ap´endice) para obtener los resultados mostrados en la Tabla 1.7 con un criterio de parada dado por |xn+1 − xn | < 10−5 :

´ 26CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES Tabla 1.7: M´et. de Newton para f (x) = x3 + 4x2 − 10. n

xn

f (xn )

f 0 (xn )

xn

f (xn )

f 0 (xn )

0 1 2 3 4 5

2.0000000 1.5000000 1.3733333 1.3652620 1.3652300 1.3652300

14.0000000 2.3750000 0.1343455 0.0005285 0.0000000 0.0000000

28.0000000 18.7500000 16.6448000 16.5139172 16.5133991 16.5133991

1.0000000 1.4545455 1.3689004 1.3652366 1.3652366

-5.0000000 1.5401953 0.0607197 0.0001088 0.0000000

11.0000000 17.9834711 16.5728681 16.5135057 16.5133991

En la Tabla 1.7 observamos que |x5 − x4 | < 10−5 cuando empezamos en x0 = 2 y |x4 − x3 | < 10−5 cuando empezamos en x0 = 1. Adem´as el criterio de parada permite garantizar que la ra´ız es exacta hasta el quinto decimal. 

1.3.1.

An´ alisis de convergencia

El error absoluto en la iteraci´on n + 1 viene dado por: en+1 = xˆ − xn+1 , y considerando xn+1 la aproximaci´on dado por el m´etodo de Newton en la relaci´on (1.30), se obtiene:   f (xn ) en f 0 (xn ) + f (xn ) . (1.31) = en+1 = xˆ − xn − 0 f (xn ) f 0 (xn ) Por otro lado, suponiendo que f ∈ C 2 ([a, b]), podemos desarrollar el polinomio de Taylor de segundo orden alrededor de xn , es decir: f 00 (ξ) f (x) = f (xn ) + f (x)(x − xn ) + (x − xn )2 2 0

donde ξ es un valor entre x y xn . Desde que se supone que xˆ est´a muy pr´oximo de xn , podemos reemplazar en la ecuaci´on anterior: 0 = f (ˆ x) = f (xn ) + f 0 (x)(ˆ x − xn ) +

f 00 (ξ) (ˆ x − xn ) 2 , 2

donde ξ es un valor entre xˆ y xn , por tanto: f (xn ) + f 0 (x)(ˆ x − xn ) = −

f 00 (ξ) f 00 (ξ) 2 (ˆ x − xn )2 = − e 2 2 n

(1.32)

´ 1.3. METODO DE NEWTON

27

Reemplazando en (1.31), concluimos:   f 00 (ξ) en+1 = − 0 e2 , 2f (xn ) n si adem´as se tiene que f 0 (x) 6= 0 para todo x ∈ [a, b], entonces: f 00 (x) 2 |en+1 | ≤ M en donde M = m´ax − 0 x∈[a,b] 2f (x) Teorema 1.2 Sea f : [a, b] ⊂ R → R tal que f y f 0 son continuas. Sea xˆ ∈ ha, bi tal que f (ˆ x) = 0 y f 0 (ˆ x) 6= 0, si adem´as existen R, L > 0 tales que: |f 0 (x) − f 0 (y)| ≤ L|x − y| para todo

x, y ∈ B(ˆ x, R).

(1.33)

entonces ∃ δ > 0 tal que si |x0 − xˆ| < δ, entonces la sucesi´on {xn }∞ n=0 dada por (1.30) converge a xˆ para cualquier punto inicial x0 ∈ B(ˆ x, δ). Demostraci´ on: 0 Como f es continua y f 0 (ˆ x) 6= 0, entonces existe R1 > 0 tal que f 0 (x) 6= 0 para todo x ∈ B(ˆ x, R1 ). Definamos:   1 c0 = sup < ∞. (1.34) f 0 (x) x∈B(ˆ x,R1 ) Para R2 = m´ın{R, R1 } definamos: g(x) = x −

f (x) f 0 (x)

para todo x ∈ B(ˆ x, R2 ).

observando que g(ˆ x) = xˆ. Adem´as: g(x) − g(ˆ x)

= = = =

|g(x) − g(ˆ x)|

≤ (1.33)

≤

f (x) − xˆ f 0 (x) f (x) (x − xˆ) − 0 f (x) 1 [f (ˆ x) − f (x) − f 0 (x)(ˆ x − x)] f 0 (x) Z  1 1 0 0 (f (x + t(ˆ x − x)) − f (x)) dt (ˆ x − x) f 0 (x) Z0 1 1 |f 0 (x + t(ˆ x − x)) − f 0 (x)| dt |(ˆ x − x)| |f 0 (x)| 0 Z 1 L |t(ˆ x − x)| dt|ˆ x − x| |f 0 (x)| 0

x−

´ 28CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES utilizando (1.34) se tiene: |g(x) − g(ˆ x)| ≤

Lc0 |ˆ x − x|2 . 2

(1.35)

  2 L c0 δ Definamos δ < m´ın R2 , , entonces α = < 1. Por tanto, dado L c0 2 x ∈ B(ˆ x, δ) se tiene: |g(x) − g(ˆ x)| ≤

L c0 δ |ˆ x − x| = α δ < δ =⇒ g(x) ∈ B(ˆ x, δ). 2

(1.36)

As´ı, dado x0 ∈ B(ˆ x, δ), por (1.30) y (1.36) implican que: |x1 − xˆ| = |g(x0 ) − g(ˆ x)| < δ =⇒ x1 ∈ B(ˆ x, δ) y procediendo de forma inductiva: xn ∈ B(ˆ x, δ). Por otro lado, de (1.36) tambi´en se tiene lo siguiente: |xn+1 − xˆ| ≤ α |xn − xˆ| ≤ αn |x0 − xˆ| para todo n ≥ 0 es decir, cuando n → ∞ la relaci´on anterior implica que xn → xˆ.

1.3.2.

An´ alisis de error

El teorema siguiente nos informa que es posible encontrar una cota para el error en el m´etodo de Newton. Teorema 1.3 Sean f, f 0 , xˆ, δ como en el Teorema 1.2, entonces existe M > 0 tal que si M δ < 1, se cumple que: |xn+1 − xn | ≤ M |xn − xˆ|2 ,

y,

(1.37)

n

(M δ)2 |xn − x| ≤ . M

(1.38)

Demostraci´ on: De la demostraci´on del Teorema 1.2 se obtiene la desigualdad (1.35). Entonc0 L ces, considerando x0 ∈ B(ˆ x, δ) y M = se tiene que: 2 |xn+1 − xˆ| ≤ M |xn − xˆ|2 c0 Lδ = α < 1, demostrando as´ı la desigualdad (1.37). 2 En la desigualdad (1.37), multiplicando por M a ambos lados se obtiene: observe que M δ =

M |xn − xˆ| ≤ (M |xn−1 − xˆ|)2

´ 1.4. METODO DE LA SECANTE

29

procediendo de forma recursiva: n

(M δ)2 |xn − xˆ| ≤ M probando as´ı la cota de error dada por (1.38).

Los teoremas 1.2 y 1.3 nos dicen que el m´etodo de Newton tiene como principal desventaja la elecci´on del punto inicial x0 , el cual debe de estar cerca de la ra´ız xˆ de la funci´on f . Sin embargo, cuando esto ocurre, la relaci´on (1.37) nos dice que el m´etodo de Newton converge a la soluci´on de forma cuadr´atica, es decir, el n´ umero de decimales exactos se duplica en cada iteraci´on. Otra desventaja del m´etodo de Newton es el c´alculo de la derivada y posterior evaluaci´on en cada iteraci´on.

1.4.

M´ etodo de la secante

Una de las desventajas del m´etodo de Newton es la necesidad de que la funci´on f sea derivable y que no se anule en cada iteraci´on. El M´etodo de la secante descrito a continuaci´on no hace uso de la derivada de la funci´on. Para ello, dado dos puntos iniciales x0 y x1 , se considera la recta secante que pasa por los puntos (x0 , f0 ), (x1 , f1 ) dada por Algoritmo 3: Algoritmo del M´etodo de la Secante Entrada: x0 , x1 , tol. 1 inicio 2 k ← 1; 3 mientras condi¸c˜ao de parada n˜ao ´e satisfeita hacer fk (xk − xk−1 ) 4 xk+1 ← xk − ; fk − fk−1 5 xk−1 ← xk ; 6 xk ← xk+1 ; 7 k ← k + 1; 8 fin mientras 9 devolver xk ; 10 fin

´ 30CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES

(a) Paso 1

(b) Paso 2

(c) Paso 3

(d) Paso 4

Figura 1.8: Interpretaci´on geom´etrica del M´etodo de la Secante. Ejemplo 1.4.1. Sea el polinomio P (x) = x3 − 2x2 + 4x − 4. Se pide: 1. Demostrar anal´ıticamente que es la u ´nica ra´ız del polinomio. 2. Encontrar una ra´ız de dicho polinomio con una precisi´on de 10−6 . Resoluci´ on: 1. Para demostrar que existe una u ´nica ra´ız, en este caso, basta analizar la derivada del polinomio  2 24 4 + > 0 para todo x ∈ R. P (x) = 3x − 4x + 4 = 3 x − 6 9 0

2

Por tanto, el polinomio es una funci´on creciente en todo su dominio y as´ı, intersecta al eje x s´olo una vez, garantizando as´ı la unicidad de la ra´ız.

´ 1.5. METODO DEL PUNTO FIJO

31

2. Para aproximar la ra´ız, consideramos como aproximaciones iniciales los puntos x0 = 3 y x1 = 6. Las iteraciones realizadas se muestran en la Tabla 1.8 con el criterio de parada |xi+1 − xi | < 10−5 . Tabla 1.8: M´etodo de la Secante para P (x) = x3 − 2x2 + 4x − 4.

1.5.

i

xi−1

xi

xi+1

|xi+1 − xi |

f (xi+1 )

1 2 3 4 5 6 7 8 9

3.000000 6.000000 2.653061 2.407525 1.797483 1.514003 1.346433 1.301010 1.295733

6.000000 2.653061 2.407525 1.797483 1.514003 1.346433 1.301010 1.295733 1.295598

2.65306122 2.40752500 1.79748326 1.51400332 1.34643306 1.30101021 1.29573267 1.29559810 1.29559774

3.34693878 0.24553623 0.61004174 0.28347994 0.16757027 0.04542284 0.00527755 0.00013457 0.00000036

11.20896905 7.99218686 2.53561237 0.94200877 0.20089249 0.02091146 0.00051994 0.00000138 0.00000000

M´ etodo del punto fijo

Primero debemos de saber que significa que un punto xˆ sea un punto fijo para una funci´on g. Definici´ on 1.1 Dada una funci´on g : [a, b] ⊂ R → R, decimos que un punto xˆ ∈ [a, b] es un punto fijo de g cuando satisface lo siguiente: g(ˆ x) = xˆ. De forma natural surge la pregunta: ¿Qu´e relaci´on tiene el punto fijo de una funci´on g con la b´ usqueda de una ra´ız para una funci´on f ?. Pues veamos, si xˆ es punto fijo de g entonces g(ˆ x) = xˆ, es decir, xˆ es una ra´ız para la funci´on f (x) = g(x) − x. Rec´ıprocamente, si r es una ra´ız de una funci´on f , entonces f (r) = 0, es decir, r es punto fijo de la funci´on g(x) = x − f (x). Por tanto, buscar una ra´ız para una funci´on f es equivalente a buscar el punto fijo para alguna funci´on g adecuada. Los teoremas siguientes dan condiciones necesarias para que una funci´on g tenga punto fijo. Teorema 1.4 Sea g : [a, b] ⊂ R → [a, b] una funci´on continua. Entonces existe xˆ ∈ ha, bi tal que g(ˆ x) = xˆ, es decir, existe un punto fijo para la funci´on g.

´ 32CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES Demostraci´ on: Definamos: f : [a, b] ⊂ R → R del modo que sigue: f (x) = x − g(x). De donde observamos: f (a) = a − g(a) < 0,

porque g(a) ∈ [a, b],

f (b) = b − g(b) > 0,

porque g(b) ∈ [a, b],

es decir, f (a) < 0 < f (b) y f es una funci´on continua, entonces, por el Teorema 1.1, concluimos que existe xˆ ∈ ha, bi tal que f (ˆ x) = 0, y de la definici´on de f se tiene que xˆ satisface g(ˆ x) = xˆ, que es lo quer´ıamos probar. Teorema 1.5 Sea g : R → R una funci´on continua y acotada. Entonces existe xˆ tal que g(ˆ x) = xˆ, es decir, existe un punto fijo para la funci´on g. Demostraci´ on: Como g es una funci´on acotada, entonces, existe k > 0 tal que |g(x)| < k para todo x ∈ R, es decir: −k < g(x) < k. Definamos: f : [−k, k] ⊂ R → R del modo que sigue: f (x) = x − g(x), y observamos que: f (k) = k − g(k) < 0, f (−k) = −k − g(−k) > 0,

porque g(k) ∈ [−k, k], porque g(−k) ∈ [−k, k],

es decir, f (k) < 0 < f (−k), as´ı, por el Teorema 1.1 existe xˆ ∈ h−k, ki tal que f (ˆ x) = 0. Por tanto, xˆ es un punto fijo para g.

1.5.1.

An´ alisis de convergencia

Los teoremas anteriores garantizan la existencia de un punto fijo punto fijo. El siguiente resultado nos permite asegurar la unicidad del punto fijo, y va m´as all´a, porque nos dice como calcular el punto fijo. Teorema 1.6 Sea g : R → R acotada y diferenciable, tal que existe λ ∈ h0, 1i que satisface: |g 0 (x)| ≤ λ < 1 para todo x ∈ R. Entonces, existe un u ´nico ∞ punto fijo xˆ de la funci´on g. Adem´as, la sucesi´on de puntos {xn }n=0 generada por: xn = g(xn−1 ) para todo n ≥ 1, (1.39) converge al u ´nico punto fijo xˆ para cualquier punto inicial x0 ∈ R.

´ 1.5. METODO DEL PUNTO FIJO

33

Demostraci´ on: Desde que g es diferenciable en R, entonces es continua en R, y por hip´otesis es acotada, entonces existe k ∈ R, k > 0 tal que |g(x)| ≤ k. Luego, por el Teorema 1.5 concluimos que existe un punto fijo xˆ ∈ h−k, ki de g. Demostremos que el punto fijo es u ´nico. En efecto, supongamos que existe otro punto fijo x∗ de la funci´on g, es decir: g(x∗ ) = x∗ . Sin p´erdida de generalidad podemos suponer que xˆ < x∗ . Podemos aplicar el Teorema del Valor Medio en el intervalo [ˆ x, x∗ ] para obtener un punto c ∈ hˆ x, x∗ i tal que: |g(x∗ ) − g(ˆ x)| = |g 0 (c)(x∗ − xˆ)| ≤ λ|(x∗ − xˆ)| Desde que xˆ y x∗ son puntos fijos de g, la relaci´on anterior implica: |x∗ − xˆ| ≤ λ|x∗ − xˆ|, es decir: 0 ≤ (1 − λ)|x∗ − xˆ| ≤ 0

⇒

(1 − λ)|x∗ − xˆ| = 0

pero λ < 1, entonces x∗ = xˆ. Por tanto, el punto fijo es u ´nico. ∞ Resta probar la convergencia de la sucesi´on {xn }n=0 hacia el u ´nico punto fijo xˆ. Veamos: |xn − xˆ| = |g(xn−1 ) − g(ˆ x)| = |g 0 (cn )(xn−1 − xˆ)|

(1.40)

para alg´ un cn entre xn y xˆ. Por hip´otesis, existe λ ∈ h0, 1i tal que |g 0 (x)| ≤ λ < 1, entonces, la ecuaci´on (1.40) implica: |xn − x∗ | ≤ λ |xn−1 − xˆ|.

(1.41)

Repitiendo n veces el proceso de la ecuaci´on (1.40) , la relaci´on (1.41) permite observar lo siguiente: |xn − xˆ| ≤ λn |x0 − xˆ| Como λ ∈ h0, 1i, entonces l´ım λn = 0, por tanto, se satisface: n→∞

l´ım xn = xˆ,

n→∞

as´ı hemos probado la convergencia.

(1.42)

´ 34CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES

(a) Paso 1

(b) Paso 2

(c) Paso 3

(d) Paso 4

Figura 1.9: Interpretaci´on geom´etrica del M´etodo de Punto Fijo. Ejemplo 1.5.1. Aproxime la ra´ız para la siguiente funci´on f (x) =

1 arctan(x − 2) − x. 3

Resoluci´ on: Observar que el dominio de f es R. Recordando que aproximar una ra´ız def es equivalente a aproximar el punto fijo para una funci´on g adecuada. Podemos considerar por ejemplo: g(x) =

1 arctan(x − 2). 3

Primero verifiquemos que se cumplen las condiciones del Teorema 1.6. En efecto, el dominio de g es R, y es acotada en su dominio porque: 1 π |g(x)| = arctan(x − 2) < para todo x ∈ R, 3 6

´ 1.5. METODO DEL PUNTO FIJO

35

adem´as: 1 ≤ 1 = λ < 1 para todo x ∈ R. |g (x)| = 2 3(1 + (x − 2) ) 3 0

Por tanto, se satisfacen las condiciones del Teorema 1.6 y as´ı concluimos que existe un u ´nico punto fijo xˆ de la funci´on g, el cual tambi´en es ra´ız de la funci´on f . En la Tabla 1.9 se muestran los resultados despu´es de aplicar el m´etodo del punto fijo a la funci´on g y se ha detenido cuando |xi+1 − xi | < 10−5 . Se muestran los resultados para dos puntos iniciales: x0 = 0 y para x0 = 50, para ambos casos se da la convergencia hacia la soluci´on. Tabla 1.9: Met. del punto fijo para f (x) = arctan(x − 2)/3 − x i

xi

0 0 1 -0.36905 2 -0.390459 3 -0.39153 4 -0.391583 5 -0.391585 6

1.5.2.

g(xi )

f (xi )

xi

-0.36905 -0.390459 -0.39153 -0.391583 -0.391585 -0.391586

-0.36905 -0.021409 -0.001071 -0.000053 -0.000003 0.000000

g(xi )

f (xi )

50 0.516655 -49.483345 0.516655 -0.325876 -0.842532 -0.325876 -0.388248 -0.062372 -0.388248 -0.39142 -0.003172 -0.39142 -0.391577 -0.000158 -0.391577 -0.391585 -0.000008 -0.391585 -0.391586 0.000000

An´ alisis de error

Hasta ahora hemos analizado la convergencia del m´etodo del punto fijo. Similar al m´etodo de la bisecci´on, para el m´etodo del punto fijo se puede establecer tambi´en una cota para el error absoluto, entre la soluci´on exacta y la aproximada. Teorema 1.7 Sean g, λ, {xn }∞ n=0 y x0 establecidos como en el Teorema 1.6 y sea xˆ el punto fijo de g. Entonces, en la iteraci´on n se cumple: |ˆ x − xn | ≤

λn |x1 − x0 |. 1−λ

Demostraci´ on: De la definici´on de xn en (1.39), xˆ punto fijo de g y de la condici´on que cumple λ para g 0 , se cumple: |xn − xn−1 | = |g(xn−1 ) − g(xn−2 )| = |g 0 (c∗ )||xn−1 − xn−2 | ≤ λ|xn−1 − xn−2 |

´ 36CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES para alg´ un c∗ entre xn−1 y xn−2 . Repitiendo el proceso anterior n−1 veces, se obtiene la siguiente desigualdad: |xn − xn−1 | ≤ λn−1 |x1 − x0 |

(1.43)

Consideremos ahora m, n ∈ N, m > n, entonces: |xm − xn | = ≤ ≤ ≤ ≤ ≤

|xm − xm−1 + xm−1 − · · · − xn+1 + xn+1 − xn | |xm − xm−1 | + |xm−1 − xm−2 | + · · · + |xn+1 − xn | λm−1 |x1 − x0 | + λm−2 |x1 − x0 | + · · · + λn |x1 − x0 | (λm−1 + λm−2 + · · · + λn )|x1 − x0 | λn (λm−n−1 + λm−n−2 + · · · + λ2 + λ + 1)|x1 − x0 | λn (1 + λ + · · · + λm−n−1 + λm−n + · · · + λm )|x1 − x0 |

ahora, fijando n y haciendo m → ∞, entonces xm → xˆ y obtenemos: |ˆ x − xn | ≤

λn |x1 − x0 | 1−λ

(1.44)

Por lo general, una funci´on s´olo est´a definida en un intervalo [a, b] ⊂ R, cuando suceda ello, el siguiente resultado da condiciones suficientes similares a los presentados en el Teorema 1.6. Teorema 1.8 Sea g : [a, b] ⊂ R → [a, b] ⊂ R continua en [a, b] y diferenciable en ha, bi, tal que existe λ ∈ h0, 1i que satisface: |g 0 (x)| ≤ λ < 1 para todo x ∈ ha, bi. Entonces, existe un u ´nico punto fijo xˆ ∈ ha, bi de la funci´on g. Adem´as, la sucesi´on de puntos {xn }∞ n=0 generada por: xn = g(xn−1 )

para todo

n ≥ 1,

(1.45)

converge al u ´nico punto fijo xˆ para cualquier punto inicial x0 ∈ ha, bi. Adem´as, en la iteraci´on n se sigue cumpliendo: |ˆ x − xn | ≤

λn |x1 − x0 |. 1−λ

(1.46)

La demostraci´on del Teorema 1.8 es muy similar a la presentada en el Teorema 1.6. El siguiente ejemplo muestra una aplicaci´on del Teorema 1.8. Ejemplo 1.5.2. ¿Cu´antas√iteraciones ser´an necesarias para aproximar una x2 + 1 ra´ız de la funci´on f (x) = − x con una tolerancia de error menor a 5 10−5 ?. Resoluci´ on:

´ 1.5. METODO DEL PUNTO FIJO

37

Veamos primero √ si la funci´on f tiene al menos una ra´ız. Para ello, considex2 + 1 remos g(x) = definida en el intervalo |x| ≤ 1. En este intervalo se 5 tiene: √ −1 < 0 ≤

√ x2 + 1 2 ≤

10−5 (1 − λ) |x1 − x0 | log(λ)

 = 12.8,

es decir, si tomamos n = 13 garantizamos la tolerancia de error pedida. Si cambiamos el punto inicial a x0 = 0.5, se tiene que: x1 = g(x0 ) ≈ 0.224, por tanto: n > 11.7, es decir, si tomamos n = 12 iteraciones garantizamos la tolerancia de error pedida. La Tabla 1.10 muestra los resultados despu´es de aplicar el m´etodo de Punto Fijo para g, los c´alculos se han hecho hasta el n´ umero de iteraciones estimado para cada punto inicial x0 .

´ 38CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES √ Tabla 1.10: M´et. punto fijo para f (x) = i

xi

g(xi )

0 -0.500000 0.223607 1 0.223607 0.204939 2 0.204939 0.204157 3 0.204157 0.204125 4 0.204125 0.204124 5 0.204124 0.204124 6 0.204124 0.204124 7 0.204124 0.204124 8 0.204124 0.204124 9 0.204124 0.204124 10 0.204124 0.204124 11 0.204124 0.204124 12 0.204124 0.204124 13 0.204124 0.204124

f (xi )

xi

0.723607 -0.018668 -0.000782 -0.000031 -0.000001 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

0.500000 0.223607 0.204939 0.204157 0.204125 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124

x2 + 1 − x. 5 g(xi ) f (xi )

0.223607 0.204939 0.204157 0.204125 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124 0.204124

-0.276393 -0.018668 -0.000782 -0.000031 -0.000001 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

Podemos observar que en ambos casos se cumple que |x6 − x5 | < 10−5 , es decir, se alcanz´o la tolerancia pedida mucho antes del valor estimado para el n´ umero de iteraciones n. Ejemplo 1.5.3. ¿Cu´antas iteraciones ser´an necesarias para aproximar una ra´ız de la funci´on f (x) = x2 − e−x con una tolerancia de error menor a 10−5 ?. Resoluci´ on: Veamos primero √ si la funci´on f tiene al menos una ra´ız. Para ello, consideremos g(x) = e−x definida en el intervalo [0, 1]. En este intervalo se tiene: 0≤

√

e−1 ≤ g(x) ≤ 1

es decir, se cumple: g(x) ⊂ [0, 1] para todo x ∈ [0, 1]. Adem´as: √ e−x 1 |g 0 (x)| = ≤ < 1 para todo x ∈ h0, 1i. 2 2 Entonces, por el Teorema 1.8 podemos concluir que existe un u ´nico punto fijo xˆ de la funci´on g en el intervalo [0, 1] y por tanto, xˆ ser´a tambi´en la u ´nica ra´ız de f en el mismo intervalo.

´ 1.5. METODO DEL PUNTO FIJO

39

El n´ umero de iteraciones necesarias para alcanzar una tolerancia de error menor a 10−5 , se calcula usando (1.46), donde: 1 λ= , 2

y si

x0 = 0 ⇒ x1 = g(x0 ) =

√

ex0 ≈ 0.7788,

por tanto: n > 15.77, es decir, si tomamos n = 16 garantizamos la tolerancia de error pedida. La Tabla 1.11 muestra los resultados despu´es de aplicar el m´etodo de Punto Fijo para g, los c´alculos se han hecho hasta el n´ umero de iteraciones estimado para cada punto inicial x0 .

Tabla 1.11: M´et. punto fijo para f (x) = x2 − e−x . i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

xi

g(xi )

f (xi )

0.50000000 0.77880078 0.67746297 0.71267379 0.70023667 0.70460470 0.70306752 0.70360810 0.70341794 0.70348483 0.70346130 0.70346958 0.70346667 0.70346769 0.70346733 0.70346746

0.77880078 0.67746297 0.71267379 0.70023667 0.70460470 0.70306752 0.70360810 0.70341794 0.70348483 0.70346130 0.70346958 0.70346667 0.70346769 0.70346733 0.70346746 0.70346741

0.27880078 0.10133782 0.03521082 0.01243711 0.00436803 0.00153719 0.00054058 0.00019015 0.00006688 0.00002352 0.00000827 0.00000291 0.00000102 0.00000036 0.00000013 0.00000004

Observar que |x12 − x11 | < 10−5 .

´ 40CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES

1.6.

M´ etodos para sistemas de ecuaciones no lineales

En esta secci´on vamos a considerar el siguiente sistema de ecuaciones no lineales f1 (x1 , . . . , xn ) = 0, f2 (x1 , . . . , xn ) = 0, (1.47) .. . fn (x1 , . . . , xn ) = 0. Considerando ahora la notaci´on vectorial, sea F : Rn → Rn definido por F (x) = (f1 (x), f2 (x), . . . fn (x)). Entonces, el problema ahora es encontrar xˆ ∈ Rn tal que F (ˆ x) = 0.

(1.48)

Para aproximar el valor exacto de xˆ haremos uso de dos m´etodos num´ericos cl´asicos que describiremos a continuaci´on.

1.6.1.

M´ etodo de Newton

Este m´etodo hace uso del polinomio de Taylor para varias variables

1.6.2.

1.7.

M´ etodo del punto fijo

Ejercicios

1. Utilice el m´etodo de Punto Fijo para determinar una aproximaci´on con un error absoluto inferior a 5×10−5 , del u ´nico cero de la funci´on f (x) = 1 + x + exp(x) = 0 que se sabe est´a en el intervalo [−2, −1]. 2. Localice gr´aficamente las ra´ıces de f (x) = 0, siendo f (x) =| x | − exp(x). Use el m´etodo de Bisecci´on para determinar un valor aproximado para los ceros de f con un error que no exceda a 0.15. 3. Aplique el m´etodo de Newton a f (x) = x2 −q (donde q > 0). Demuestre que si xn tiene despu´es del punto decimal k d´ıgitos correctos, entonces xn+1 tendr´a despu´es del punto decimal al menos 2k−1 d´ıgitos correctos, con tal que q > 0.006 y k ≥ 1.

1.7. EJERCICIOS

41

4. En el m´etodo de la Secante demuestre que si xn → q conforme n → ∞ y si f 0 (q) 6= 0, entonces q es un cero de f . 5. Aproxime la soluci´on para el sistema de ecuaciones no lineales siguientes:   x + 1 ln(x x ) − 2 = 0 1 2 1 2  x1 x22 + 5x2 − 3 = 0 con un error menor a 10−5 . 6. Compare el n´ umero de iteraciones que realizan el M´etodo de la Bisecci´on, el M´etodo de la Regla Falsa y el M´etodo de la Regla Falsa Modificada al ser aplicadas a las siguientes funciones, [2]: a) f (x) = x2 − (1 − n)n ,

n = 1, 2, 3, . . ..

2

b) f (x) = (1 + (1 − n) )x − (1 − nx)2 ,

n = 1, 2, 3, . . ..

Use n = 2, 5, 15, 20. 7. Localice gr´aficamente los ceros de las siguientes funciones a) f (x) = 4cos(x) − exp(2x).

c) f (x) = x ln(x).

1 −

b) f (x) = x/2 − tan(x).

d ) f (x) = 2x − 3x.

e) f (x) = x3 + x − 1000. f ) f (x) = x2 − e−x .

8. Determine con un error absoluto menor de 0.001 la soluci´on de la ecuaci´on x − cos(x) = 0. 9. Considere la funci´on f (x) = x2 /2 + x(ln(x) − 1). Obtener los puntos cr´ıticos con el auxilio del m´etodo de las secantes. 10. Considere la ecuaci´on f (x) = 4x3 + 3x − 1 = 0. Intente obtener una ra´ız real de la ecuaci´on usando el m´etodo de las secantes. Use x0 = 0 y x1 = 1 como estimativas iniciales de la ra´ız buscada. Haga 5 iteraciones de forma ordenada. 11. Considere la ecuaci´on f (x) = 2x4 + 3x2 − 10 = 0. Realice 3 iteraciones del m´etodo de las secantes iniciando en x0 = 0.5 y x1 = 1.2. 12. Use el m´etodo de Newton para obtener la menor ra´ız positiva de las siguientes ecuaciones:

´ 42CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES a) x/2 − tan(x) = 0. b) 2cos(x) − exp(x)/2 = 0. c) x5 − 6 = 0. 13. Usando el m´etodo de Regula Falsi resuelva la ecuaci´on x + tan(x) = 0 partiendo del intervalo [1.9, 2.1] iterando hasta que la soluci´on tenga una precisi´on de 0.001. Repetir el c´alculo pero usando el m´etodo de Newton partiendo de x = 1.7 hasta alcanzar una precisi´on de 0.0001. 14. Considere la ecuaci´on 2x − cos(x) = 3. Se pide a) Demostrar que esta ecuaci´on tiene una s´ola ra´ız. b) Determinar el valor de la ra´ız mediante el m´etodo de NewtonRaphson, partiendo del punto x = 0 y haciendo las iteraciones necesarias para que el error sea menor que 0.0001. 15. Considere la funci´on f (x) = exp(x) − 4x2 . Se pide: a) Localizar graficamente las ra´ıces de la funci´on f . b) Considere el intervalo I = [−1, 5]. Realizar dos iteraciones del m´etodo de la bisecci´on y escoja el punto medio del u ´ltimo intervalo obtenido como aproximaci´on inicial para el m´etodo de Newton. Aplique el m´etodo de Newton hasta alcanzar una precisi´on de 10−2 . Comparando con la localizaci´on de los zeros realizada en el item anterior, identifique cual ra´ız fue obtenida en este proceso y justifique por que la convergencia fue para esta ra´ız. 16. Aplique el m´etodo de Newton-Raphson a la ecuaci´on x3 −2x2 −3x+10 = 0 con x0 = 1.9. Justifique los resultados obtenidos. 17. El valor de π puede ser obtenido como soluci´on de las siguientes ecuaciones sin(x) = 0 y cos(x) + 1 = 0. a) Aplique el m´etodo de Newton con x0 = 3 y precisi´on 10−2 en cada caso y compare los resultados. Justifique los resultados obtenidos. b) ¿Puede aplicarse el m´etodo de la bisecci´on en la resoluci´on de la ecuaci´on cos(x) + 1 = 0 para obtener el valor de π?. 18. Use el m´etodo del punto fijo para encontrar una ra´ız de f (x) = x3 + 4x2 − 10 en el intervalo [1, 2]. 19. Use el m´etodo del punto fijo para encontrar una ra´ız de f (x) = x − x4/5 − 2.

Cap´ıtulo 2 Interpolaci´ on polinomial Dada una funci´on continua f para la cual s´olo conocemos su valor yi para n + 1 puntos discretos xi , es decir: f (xi ) = yi ,

i = 0, 1, 2, ..., n.

(2.1)

El objetivo de este cap´ıtulo es encontrar un polinomio Pn de grado n y establecer condiciones para la funci´on f tal que se cumpla la siguiente condici´ on de interpolaci´ on, es decir: Pn (xi ) = f (xi ),

i = 0, 1, 2, ..., n.

(2.2)

Como primer paso es garantizar la existencia del polinomio Pn , el cual, en forma general tiene la forma: Pn (x) = a0 + a1 x + a2 x2 + ... + an xn

(2.3)

Como Pn debe de satisfacer (2.2), entonces: a0 + a1 xi + a2 x2i + ... + an xni = f (xi ),

i = 0, 1, 2, ..., n.

(2.4)

El conjunto de ecuaciones dado por (2.4) al ser escrito de forma matricial queda como sigue: 

x20 x21

x30 x31

xn0 xn1

1 x0 ···  1 x1 ···   ..  . 2 3 1 xn xn xn · · · xnn





     

43

a0 a1 a2 .. . an





      =    

f (x0 ) f (x1 ) f (x2 ) .. . f (xn )

      

(2.5)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

44

donde observamos que la matriz asociada al sistema lineal (2.5) es conocida como matriz de Vandermonde, cuyo determinante se calcula como: 1 x0 x2 x 3 · · · xn 0 0 0 1 x1 x2 x3 · · · xn 1 1 1 (2.6) = .. . 1 xn x2n x3n · · · xnn podemos observar que para garantizar la existencia y unicidad del polinomio Pn es necesario y suficiente que los n + 1 puntos xi (i = 0, 1, 2, ..., n) sean todos distintos, lo cual ser´a considerado de ahora en adelante salvo se mencione expl´ıcitamente lo contrario. Una vez garantizada la existencia y unicidad del polinomio Pn , ahora llamado Polinomio de interpolaci´ on, procedemos a estudiar algunos m´etodos que nos permiten calcularlo.

2.1.

M´ etodo de coeficientes indeterminados

2.2.

M´ etodo de Lagrange

Este m´etodo fue desarrollado por Lagrange () el cual consiste primero en definir los llamados polinomios de Lagrange dados por: n Y (x − xj ) Li (x) = (xi − xj ) j=0

(2.7)

j6=i

es f´acil observar que todas las Li (i = 0, 1, ..., n) cumplen las siguientes propiedades: Cada Li es un polinomio de grado n.  1, i = j Li (xj ) = 0, i 6= j Calculadas todas las Li , definimos: Pn (x) = f (x0 )L0 (x) + · · · + f (xn )Ln (x) =

n X

f (xj )Lj (x).

(2.8)

j=0

De la forma como es construido el polinomio Pn dado en (2.8), puede observarse que satisface la condici´on de interpolaci´on dado en (2.2). Por tanto, Pn es el polinomio de interpolaci´on dado en la forma de Lagrange.

´ 2.2. METODO DE LAGRANGE

45

Ejemplo 2.2.1. Dados los puntos (−1, 1), (2, −1), (0, −1). Calcule el polinomio de interpolaci´on en la forma de Lagrange. Resoluci´ on: Dados los 3 puntos, procedemos a calcular L0 , L1 y L2 seg´ un (2.7):

L0 (x) =

(x − 2)x (x − 2)(x − 0) = (−1 − 2)(−1 − 0) 3

L1 (x) =

(x + 1)(x − 0) (x + 1)x = (2 + 1)(2 − 0) 6

L2 (x) =

(x + 1)(x − 2) (x + 1)(x − 2) =− (0 + 1)(0 − 2) 2

luego, el polinomio de interpolaci´on en la forma de Lagrange (2.8) viene dado por:

P2 (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + f (x2 )L2 (x) P2 (x) = L0 (x) − L1 (x) − L2 (x) P2 (x) =

(x − 2)x (x + 1)x (x + 1)(x − 2) − + 3 6 2

La Figura 2.1 muestra los puntos dados con el polinomio de interpolaci´on P2 (x).

46

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

Figura 2.1: P2 (x) vs f (xi )

Dados los puntos (x0 , f (x0 )), ..., (xn , f (xn )) es relativamente f´acil calcular Pn en la forma de Lagrange, pero, ¿qu´e sucede si agregamos a los puntos iniciales el punto (xn+1 , f (xn + 1))?. Ante esta situaci´on, el m´etodo de Lagrange no es ventajoso, porque tenemos que rehacer todos los c´alculos para obtener el polinomio Pn+1 . El m´etodo de Newton, descrito en la siguiente secci´on, resuelve este problema.

2.3.

M´ etodo de Newton

Dados los puntos {(xi , yi )}ni=0 y construimos el polinomio de interpolaci´on Pn . A este conjunto de puntos vamos a agregar un nuevo punto (xn+1 , yn+1 ) y deseamos calcular el polinomio de interpolaci´on Pn+1 aprovechando que ya conocemos Pn . Dado que Pn+1 es de grado n + 1 y Pn es de grado n, luego se cumple: Pn+1 (x) = Pn (x) + Qn+1 (x)

(2.9)

´ 2.3. METODO DE NEWTON

47

para alg´ un polinomio Qn+1 de grado n+1. Reemplazando x = xi (i = 0, ..., n) en (2.9), se obtiene: Pn+1 (xi ) = Pn (xi ) + Qn+1 (xi )

⇒

Qn+1 (xi ) = 0.

(2.10)

Es decir, (2.10) nos dice que x0 , x1 , ..., xn son las n + 1 ra´ıces del polinomio Qn+1 . Por tanto: Qn+1 (x) = an+1 θn+1 (x) (2.11) para alguna constante an+1 ∈ R y donde θn+1 (x) = (x − x0 )(x − x1 )...(x − xn ). Reemplazamos en (2.10) para x = xn+1 y obtenemos: yn+1 = Pn (xn+1 ) + an+1 θn+1 (xn+1 ) de ah´ı: an+1 =

yn+1 − Pn (xn+1 ) θn+1 (xn+1 )

(2.12)

Ejemplo 2.3.1. Para los puntos (−1, 1), (2, −1), (0, −1) del ejemplo anterior, calcule el nuevo polinomio de interpolaci´on despu´es de agregar el punto (−2, −2). Resoluci´ on: Procedemos de la forma siguiente: Para (x0 , y0 ), definimos P0 (x) = y0 = 1. Para (x0 , y0 ), (x1 , y1 ) calculamos a1 =

y1 − P0 (x1 ) 2 =− , (x1 − x0 ) 3

luego: 2 P1 (x) = P0 (x) + a1 (x − x0 ) = 1 − (x + 1). 3 Para (x0 , y0 ), (x1 , y1 ), (x2 , y2 ) calculamos a2 =

y2 − P1 (x2 ) 2 = , (x2 − x1 )(x2 − x0 ) 3

luego: 2 2 P2 (x) = P1 (x) + a2 (x − x0 )(x − x1 ) = 1 − (x + 1) + (x + 1)(x − 2). 3 3

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

48

Para (x0 , y0 ), (x1 , y1 ), (x2 , y2 ), (x3 , y3 ) calculamos a3 =

y3 − P2 (x3 ) 19 = , (x3 − x2 )(x3 − x1 )(x3 − x0 ) 24

luego: P3 (x) = P2 (x) + a3 (x − x0 )(x − x1 )(x − x2 ) 2 2 19 P3 (x) = 1 − (x + 1) + (x + 1)(x − 2) + (x + 1)(x − 2)x. 3 3 24 Los resultados pueden observarse en la Figura 2.2.

Figura 2.2: yi vs P3 (x)

2.4.

M´ etodo de Diferencias Divididas

El valor de an+1 dado en (2.12) es llamada F´ ormula de diferencia dividida de orden n + 1 y es denotada por f [x0 , x1 , ..., xn , xn+1 ] := an+1

(2.13)

´ 2.4. METODO DE DIFERENCIAS DIVIDIDAS

49

y por convenci´on se acepta que: f [x0 , x0 ] := a0 = y0

(2.14)

As´ı, el polinomio de interpolaci´on en la forma de Newton puede ser escrito como sigue: n+1 i−1 X Y Pn+1 (x) = f [x0 , x0 ] + ai (x − xj ) (2.15) i=1

j=0

y usando la notaci´on (2.13) se tiene: Pn+1 (x) = f [x0 , x0 ] +

n+1 X

f [x0 , ..., xi ]

i=1

i−1 Y (x − xj )

(2.16)

j=0

Una ventaja de la notaci´on usada en (2.13) es que permite obtener cada coeficiente dado en (2.16) de una forma pr´actica, que mostraremos a continuaci´on: Consideremos los puntos: (xi , yi ), (xi+1 , yi+1 ), · · · (xi+n , yi+n ), (xi+n+1 , yi+n+1 ),

(2.17)

y consideremos los polinomios siguientes: Pn interpola a (xi , yi ), (xi+1 , yi+1 ), · · · (xi+n , yi+n ), es decir Pn (xi+j ) = yi+j

j = 0, 1, 2, ..., n

Qn interpola a (xi+1 , yi+1 ), · · · (xi+n , yi+n ), (xi+n+1 , yi+n+1 ), es decir Qn (xi+j ) = yi+j

j = 1, 2, 3, ..., n + 1.

Por tanto, de (2.16) los polinomios Pn y Qn tienen la siguiente forma: Pn (x) = f [xi , xi ] +

n X

f [xi , xi+1 , ..., xi+j ](x − xi )...(x − xi+j−1 )

(2.18)

j=1

Qn (x) = f [xi+1 , xi+1 ] +

n+1 X

f [xi+1 , xi+2 , ..., xi+j ](x − xi+1 )...(x − xi+j−1 )

j=2

(2.19)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

50

Ahora, a partir de (2.18) y (2.19) construimos un polinomio Hn+1 de grado n + 1 que interpola a todos los puntos dados en (2.17), el cual es dado por Hn+1 (x) =

x − xi xi+n+1 − x Pk (x) + Qk (x) xi+n+1 − xi xi+n+1 − xi

(2.20)

De la unicidad del polinomio de interpolaci´on y usando nuevamente (2.16), el polinomio Hn+1 tambi´en puede ser escrito como: Hn+1 (x) = f [xi , xi ] +

n+1 X

f [xi , xi+1 , ..., xi+j ](x − xi )(x − xi+1 )...(x − xi+j−1 )

j=1

(2.21) procedemos a igualar los coeficientes de grado n + 1 en las expresiones de Hn+1 dadas en (2.20) y (2.21), obteniendo la diferencia dividida de orden n + 1 del modo siguiente: f [xi , xi+1 , ..., xi+n+1 ] =

f [xi+1 , xi+2 , ..., xi+n+1 ] − f [xi , xi+1 , ..., xi+n ] (2.22) xi+n+1 − xi

Por simplicidad, consideremos los puntos (xi , yi ), (xi+1 , yi+1 ), (xi+2 , yi+2 ) y (xi+3 , yi+3 ), podemos representar la expresi´on dada en (2.22) en una forma m´as pr´actica en la siguiente tabla: ∆(0)

∆(1)

∆(2)

xi

yi

xi+1

yi+1

f [xi , xi+1 ]

xi+2

yi+2

f [xi+1 , xi+2 ]

f [xi , xi+1 , xi+2 ]

xi+3

yi+3

f [xi+2 , xi+3 ]

f [xi+1 , xi+2 , xi+3 ]

∆(3)

f [xi , xi+1 , xi+2 , xi+3 ]

Figura 2.3: Tabla de diferencias divididas. donde cada ∆(j) , (j = 1, 2, 3) representa la diferencia dividida de orden j. Luego, el polinomio de interpolaci´on en la forma de Newton es:

P3 (x) = f [xi , xi ]+f [xi , xi+1 ](x−xi )+f [xi , xi+1 , xi+2 ](x−xi )(x−xi+1 )+f [xi , xi+1 , xi+2 , xi+3 ](x− Ejemplo 2.4.1. Considere los puntos (1, 1/2), (0, −1/2), (3/2, 1) y (2, 2). Use el m´etodo de diferencias divididas para calcular el polinomio de interpolaci´on en la forma de Newton.

´ 2.4. METODO DE DIFERENCIAS DIVIDIDAS

51

Resoluci´ on: Construimos la tabla seg´ un la Figura 2.3:

∆(0)

∆(1)

∆(2)

1

1 2

0

−

2 3

1

9 4

−

2

2

3 4

−

1 2

∆(3)

1 15 4 3 4

3

Figura 2.4: Tabla de diferencias divididas para el ejemplo dado.

De la Figura 2.4 podemos construir el polinomio de interpolaci´on en la forma de Newton, dado por:

  1 15 2 P3 (x) = + (x − 1) − (x − 1)x + 3(x − 1)x x − . 2 4 3

La Figura 2.5 muestra el polinomio de interpolaci´on P3 con los puntos dados.

52

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

Figura 2.5: Polinomio de interpolaci´on por diferencias divididas.

Una ventaja del m´etodo de diferencias divididas es que nos permite encontrar tambi´en el polinomio de interpolaci´on para un subconjunto de los puntos dados, por ejemplo, a partir de la Figura 2.4 podemos obtener el polinomio de interpolaci´on P2 para los puntos (0, −1/2), (2/3, 1) y (2, 2), el cual es dado por:

  1 9 3 2 P2 (x) = − + x − x x − . 2 4 4 3

´ 2.5. ANALISIS DE ERROR

53

Figura 2.6: Diferencias divididas para un subconjunto de puntos.

2.5.

An´ alisis de error

Dados los puntos (x0 , y0 ), (x1 , y1 ), ..., (xn , yn ), vamos a suponer que los valores yi (i = 0, ..., n) son im´agenes de cada xi para alguna funci´on continua f , es decir: yi = f (xi ) para todo i = 0, 1, ..., n. Por los puntos dados, podemos construir un polinomio de interpolaci´on Pn que satisface: Pn (xi ) = yi . Ahora, dado un punto x tal que x 6= xi para todo i = 0, 1, ..., n y: a := m´ın {xi } ≤ x ≤ m´ax {xi } =: b i=1,...,n

i=1,...,n

y al aproximar el valor f (x) usando Pn es l´ogico que se comete un error. La pregunta es, se puede saber que tan buena es la aproximaci´on? es decir, qu´e se puede decir del error: E(x) = f (x) − Pn (x)

(2.23)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

54

Para responder esta pregunta, fijando x definamos g : [a, b] → R como sigue: g(t) = f (t) − Pn (t) −

θn+1 (t) [f (x) − Pn (x)] θn+1 (x)

(2.24)

donde θn+1 (t) = (t − x0 )(t − x1 )...(t − xn−1 )(t − xn ). Observemos que si t = xi (i = 0, 1, 2, ..., n) en (2.24), entonces g(xi ) = 0. As´ı mismo, si t = x, tambi´en implica g(x) = 0. Es decir, {x0 , ..., xn , x} son n + 2 ra´ıces de la funci´on g. Luego, exigiendo que f sea n+1 veces derivable, podemos aplicar el Teorema de Rolle continuamente y obtenemos que existe ξ ∈ (a, b) tal que g (n+1) (ξ) = 0. As´ı, a partir de (2.24) se tiene: g (n+1) (t) = f (n+1) (t) −

(n + 1)! E(x) θn+1 (x)

Por tanto: g (n+1) (ξ) = f (n+1) (ξ) −

(n + 1)! E(x) = 0 θn+1 (x)

lo que implica E(x) =

f (n+1) (ξ) θn+1 (x) (n + 1)!

(2.25)

La relaci´on dada en (2.25) nos da una expresi´on para el error que se comete en alg´ un punto x ∈ (a, b). Observar que E(xi ) = 0 para todo i = 0, 1, ..., n. x8 3cos(2x) − . Calcule una cota del error 84 8 cuando se aproxima f (0.65) cuando se considera 5 puntos uniformemente espaciados en [0, 1]. Resoluci´ on: Consideremos la partici´on uniforme: Ejemplo 2.5.1. Sea f (x) =

x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75, x4 = 1. Es decir, n = 4, luego, el error (2.25) viene dado por: E(x) =

f (v) (ξ) (x − x0 )(x − x1 )(x − x2 )(x − x3 )(x − x4 ) 5!

´ 2.5. ANALISIS DE ERROR

55

para alg´ un ξ ∈ (0, 1). Calculando f (v) se obtiene: f (v) (x) = 12sen(2x) + 80x3 ,

para todo x ∈ [0, 1].

Luego: |f (v) (x)| ≤ 92, por tanto:

(v) f (ξ) (x − x0 )(x − x1 )(x − x2 )(x − x3 )(x − x4 ) |E(x)| = 5! 92 |E(0.65)| ≤ |(0.65)(0.40)(0.15)(−0.15)(−0.35)| = 0.0016. 120 El polinomio de interpolaci´on P4 es:         3 1 1 1 1 1 P4 (x) = − +0.1836x+0.6450x x − −0.3191x x − x− −0.0505x x − x− x 8 4 4 2 4 2 Luego, los valores son: f (0.65) = −0.0999327,

P4 (0.65) = −0.1002010

por tanto, el error exacto es: |E(0.65)| = |f (0.65) − P4 (0.65)| = 0.0002683 < 0.0016. La Figura 2.7 compara los valores de f y P4 para los puntos de la partici´on considerada.

Figura 2.7: An´alisis de error.

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

56

En la Figura 2.7 podemos observar la buena aproximaci´on de P4 en relaci´on a f dentro del intervalo de interpolaci´on [0, 1] y tal vez nos hace suponer que la aproximaci´on mejorar´a conforme aumentemos el n´ umero de puntos en la partici´on. En la siguiente secci´on mostramos un ejemplo en el cual observaremos que nuestra suposici´on est´a errada.

2.6.

Fen´ omeno Runge

Consideremos la funci´on siguiente:

f (x) =

1 , 1 + 25x2

para todo x ∈ [−1, 1].

(2.26)

Consideremos tambi´en una partici´on regular en el intervalo [−1, 1]

−1 = x0 < x1 < · · · < xn−1 < xn = 1 (n es impar).

(2.27)

La Figura 2.8 muestra la gr´afica de la funci´on f compar´andola con los polinomios de interpolaci´on P5 , P9 y P11 . Podemos observar que la aproximaci´on en los extremos del intervalo cada vez empeora conforme aumenta el valor del grado del polinomio de interpolaci´on. En general, se prueba que:

 l´ım

n→∞

 m´ax |f (x) − Pn (x)| = ∞.

−1≤x≤1

(2.28)

´ 2.6. FENOMENO RUNGE

57

1

0.8

0.6

0.4

0.2

0

−0.2

−0.4 −1

−0.8

−0.6

−0.4

−0.2

0

Figura 2.8: Fen´omeno Runge.

0.2

0.4

0.6

0.8

1

58

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

Cap´ıtulo 3 Integraci´ on Num´ erica Es bien conocido el Teorema Fundamental del C´alculo, el cual establece que: Teorema 3.1 Sea f : [a, b] → R una funci´on continua. Luego: Si definimos F : [a, b] → R como sigue: Z x F (x) = f (t) dt, a

entonces se cumple que F 0 (x) = f (x)

para todo

x ∈ [a, b].

Tal funci´on F es llamada primitiva o antiderivada de la funci´on f . Si una funci´on G es una primitiva de f en [a, b], entonces: Z b f (x) dx = G(b) − G(a). a

Por tanto, dada f : [a, b] → R una funci´on continua, para calcular Z b f (x) dx

(3.1)

a

el Teorema 3.1 establece que s´olo es necesario hallar una primitiva G de la funci´on f . En la pr´actica, hallar la funci´on G no siempre es posible, es conocido que para ciertas funciones f es imposible hallar una primitiva y as´ı no es posible usar el Teorema 3.1. Nuestro objetivo en el presente cap´ıtulo ser´a aproximar el valor de la integral dada en (3.1). 59

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

60

3.1.

M´ etodo del Trapecio

Para una funci´on f que satisface las condiciones del Teorema 3.1 y que adem´as f (x) ≥ 0 para todo x ∈ [a, b], sabemos que la integral dada en (3.1) viene dada por el a´rea bajo la curva de f . Esta ´area puede ser aproximada intuitivamente por el a´rea del trapecio tal y como se muestra en la Figura 3.1, es decir   Z b f (a) + f (b) . (3.2) f (x) dx ≈ (b − a) 2 a

´ Figura 3.1: Area del trapecio bajo la gr´afica de f . En la siguiente secci´on obtendremos de forma anal´ıtica la expresi´on dada en (3.2) con su respectivo t´ermino de error.

3.1.1.

Simple

Siendo f : [a, b] → R una funci´on continua y consideramos x0 = a, x1 = b. Con la finalidad de utilizar interpolaci´on polinomial, exigimos que f sea dos veces derivable. Luego, se cumple f (x) = P1 (x) + E1 (x) donde P1 (x) = f0 +

(x − x0 ) (f1 − f0 ) h

(3.3)

´ 3.1. METODO DEL TRAPECIO

61

y la expresi´on para el error viene dado por f 00 (zx ) (x − x0 )(x − x1 )dx 2

E1 (x) =

para alg´ un zx ∈ hx0 , x1 i. Luego: Z Z x1 P1 (x)dx =

x1

x0

x0



 (x − x0 ) f0 + (f1 − f0 ) h

(3.4)

h (f0 + f1 ), = 2 Z x1 E1 (x) dx ET S = x0 x1

f 00 (zx ) = (x − x0 )(x − x1 ) dx 2 x0 Z f 00 (ξ) x1 = (x − x0 )(x − x1 )dx 2 x0 Z

(3.5)

f 00 (ξ) 3 h, 12 donde se ha observado que g(x) = (x − x0 )(x − x1 ) no cambia de signo para todo x ∈ [x0 , x1 ], por tanto, para usar el teorema del valor medio para integrales exigimos que f 00 sea continua en [x0 , x1 ] y as´ı obtenemos la expresi´on para el error ET S. Luego, integrando (3.3) y reemplazando en ´el los resultados obtenidos en (3.4) y (3.5) se obtiene Z x1 f 00 (ξ) 3 h h (3.6) f (x)dx = (f0 + f1 ) − 2 12 x0 = −

luego, a partir de (3.6), se denomina m´ etodo del trapecio simple para la funci´on f a la siguiente igualdad: Z x1 h f (x)dx ≈ M T S(f ) = (f0 + f1 ) (3.7) 2 x0

3.1.2.

Compuesto

Es natural pensar que usar el m´etodo del trapecio simple para una funci´on f da en general aproximaciones muy alejadas del valor exacto de la integral, como puede observarse en la Figura 3.1. Con el objetivo de mejorar la aproximaci´on, es natural pensar en considerar una partici´on regular del intervalo

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

62

[a, b], es decir, dado un n´ umero natural n se construyen n subintervalos del modo siguiente: a = x0 < x1 < ... < xn−1 < xn = b b−a xn − x0 = y xi = x0 + ih (i = 0, ..., n). n n Luego aplicamos el MTS en cada subintervalo Ii = [xi , xi+1 ], (i = 0, ..., n−1). Veamos el procedimiento:

donde h =

Z

xn

f (x)dx =

n−1 Z X

x0

xi+1

f (x)dx =

xi

i=0

n−1 X

MTSi (f ) +

i=0

n−1 X

ETSi

i=0

donde Z

xi+1

MTSi (f ) =

f (x) dx = xi

h (fi + fi+1 ), 2

i = 0, ..., n − 1

y el error total es ETC =

n−1 X i=0

ETSi = −

n−1 00 X f (ξi ) i=0

12

h3

(3.8)

donde cada ξi ∈ Ii (i = 0, ..., n − 1). El m´ etodo del trapecio compuesto es definido como: M T C(f ) =

n−1 X i=0

M T Si (f ) =

n−1 X h i=0

2

(fi + fi+1 )

(3.9)

realizando la suma en (3.9), el m´etodo del trapecio compuesto queda ! n−1 X h (3.10) M T C(f ) = f0 + 2 fi + fn 2 i=1 A continuaci´on vamos a analizar el error para el m´etodo del trapecio compuesto, ETC, con el objetivo de encontrar una forma m´as simple de acotar. Sabemos que f 00 es continua en el intervalo cerrado y acotado [a, b], por tanto, existen x∗ , xˆ ∈ [a, b] donde f 00 alcanza su m´ınimo y m´aximo respectivamente, es decir f 00 (x∗ ) ≤ f 00 (x) ≤ f 00 (ˆ x) para todo x ∈ [a, b] (3.11) en particular, para ξi ∈ Ii ⊂ I, (i = 0, ..., n − 1) tambi´en se cumple f 00 (x∗ ) ≤ f 00 (ξi ) ≤ f 00 (ˆ x)

(3.12)

´ 3.1. METODO DEL TRAPECIO

63

sumando para i = 0, ..., n − 1 en (3.12) se obtiene 00

∗

n f (x ) ≤

n−1 X

f 00 (ξi ) ≤ n f 00 (ˆ x)

(3.13)

i=0

es decir

n−1 X

f 00 (x∗ ) ≤

f 00 (ξi )

i=0

≤ f 00 (ˆ x) n dado que f 00 es una funci´on continua en [a, b], entonces podemos aplicar el Teorema del Valor Intermedio, es decir, existe ξ ∈ ha, bi tal que n−1 X

f 00 (ξ) =

f 00 (ξi )

i=0

n

(3.14)

Reemplazando adecuadamente (3.14) en la expresi´on para el error del m´etodo del trapecio compuesto dado por (3.8), obtenemos ETC = −

n−1 00 X f (ξi )

h3 00 f (ξ) para alg´ un ξ ∈ ha, bi. 12

(3.15)

h3 00 (b − a)3 00 f (ξ) = − f (ξ) para alg´ un ξ ∈ ha, bi. 12 12n2

(3.16)

i=0

12

h3 = −n

o de forma equivalente ETC = −n

Ejemplo 3.1.1. Aplice el MTS a la siguiente integral Z e x2 ln(x) dx. 1

Luego, estime el n´ umero de subintervalos n para alcanzar una aproximaci´on de 10−2 al usar el MTC. Resoluci´ on: Aplicando el MTS a f (x) = x2 ln(x) se tiene     f (1) + f (e) 0 + e2 M T S(f ) = (e − 1) = (e − 1) = 6.34824. 2 2 Ahora, para estimar el n´ umero de subintervalos n, debemos acotar f 00 seg´ un lo sugiere la f´ormula dada en (3.16). A partir de f obtenemos: f 00 (x) = 2 ln(x) + 3

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

64 y como 1 ≤ x ≤ e, entonces

|f 00 (x)| ≤ 2 ln(e) + 3 = 5. Por tanto, de (3.16) se obtiene (e − 1)3 00 (e − 1)3 |ET C| = f (ξ) ≤ 5 12n2 12n2 y queremos una tolerancia de error de 10−2 , por tanto, exigimos que: (e − 1)3 ≤ 10−2 5 12n2 de donde se obtiene r n≥

5 × 102 (e − 1)3 ≈ 14.53. 12

Por tanto, debemos tomar n = 15 como n´ umero de subintervalos para alcanzar la tolerancia pedida. La Tabla 3.1 muestra los c´alculos respectivos.

Tabla 3.1: Resultados del MTC para f (x) = x2 ln(x). i xi f (xi ) Peso Peso f (xi ) 0 1.0000 0.00000000 1 0.00000000 1 1.1146 0.13472274 2 0.26944547 2 1.2291 0.31163516 2 0.62327032 3 1.3437 0.53330895 2 1.06661791 4 1.4582 0.80208584 2 1.60417169 5 1.5728 1.12011539 2 2.24023078 6 1.6873 1.48938419 2 2.97876839 7 1.8019 1.91173888 2 3.82347776 8 1.9164 2.38890456 2 4.77780912 9 2.0310 2.92249983 2 5.84499966 10 2.1455 3.51404919 2 7.02809838 11 2.2601 4.16499338 2 8.32998675 12 2.3746 4.87669810 2 9.75339620 13 2.4892 5.65046148 2 11.30092296 14 2.6037 6.48752041 2 12.97504082 15 2.7183 7.38905610 1 7.38905610

´ 3.2. METODO DE SIMPSON

65

De la Tabla 3.1 se obtiene M T C(f ) = 4.58238800, y realizando una simple integraci´on por partes se puede obtener en este caso el valor exacto de la integral, es decir Z e 2e3 1 + = 4.57456376, f (x) dx = 9 9 1 por tanto, se puede calcular el error exacto Z e f (x) dx − M T C(f ) = 0.00782 < 0.01. |ET C| = 1

3.2.

M´ etodo de Simpson

En la secci´on anterior, hemos usado interpolaci´on polinomial con dos puntos, a saber, los extremos del intervalo [a, b]. Con esos puntos hemos deducido el m´etodo del trapecio simple para luego generalizar realizando una partici´on regular de [a, b] y obtener as´ı el m´etodo del trapecio compuesto. Siguiendo con esta idea, en el intervalo [a, b] consideramos los puntos a = x0 < x1 =

a+b < x2 = b, 2

y definiendo h = x2 −x1 = x1 −x0 . Con los puntos x0 , x1 y x2 podemos aplicar interpolaci´on polinomial para obtener otro m´etodo de integraci´on num´erica, que ser´a descrito en la siguiente secci´on.

3.2.1.

Simple

Dado que tenemos los puntos (x0 , f0 ), (x1 , f1 ) y (x2 , f2 ), usamos interpolaci´on polinomial mediante el m´etodo de coeficientes indeterminados para obtener P2 , polinomio de segundo grado, que se puede expresar de la siguiente forma P2 (x) = p + qx + rx2 (3.17) que satisface P2 (x0 ) = f0 ,

P2 (x1 ) = f1 ,

y P2 (x2 ) = f2 .

(3.18)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

66 Se cumple

f (x) = P2 (x) + E(x)

(3.19)

donde E(x) representa el error cometido en la interpolaci´on. Integrando Z

x2

Z

x2

x2

E(x) dx

P2 (x) dx +

f (x) dx = x0

Z

(3.20)

x0

x0

Definimos el M´ etodo de Simpson simple para la funci´on f a la siguiente expresi´on Z

x2

M SS(f ) =

P2 (x) dx.

(3.21)

x0

lo cual es representado geom´etricamente en la Figura 3.2.

´ Figura 3.2: Area del trapecio bajo la gr´afica de f .

Vamos a usar las relaciones dadas en (3.17) y (3.18) para encontrar una

´ 3.2. METODO DE SIMPSON

67

expresi´on de M SS(f ) dado en (3.21). Z

x2

P2 (x) dx

M SS(f ) = x0

  x x2 x3 2 = px + q + r 2 3 x0 r q = p(x2 − x0 ) + (x22 − x20 ) + (x32 − x30 ) 2 3 =

(x2 − x0 ) [6p + 3q(x2 + x0 ) + 2r(x22 + x2 x0 + x20 )] 6

=

(x2 − x0 ) [(p + qx2 + rx22 ) + (p + qx0 + rx20 ) + 2q(x0 + x2 ) + 6

+ r(x22 + 2x2 x0 + x20 ) + 4p] " #    2 (x2 − x0 ) x0 + x2 x0 + x2 = f0 + f2 + 4q + 4r + 4p 6 2 2 =

(x2 − x0 ) [f0 + f2 + 4(qx1 + rx21 + p)] 6

=

(x2 − x0 ) 2h [f0 + 4f1 + f2 ] = [f0 + 4f1 + f2 ] 6 6

As´ı el M´etodo de Simpson simple dado en (3.21) queda expresado de la siguiente forma M SS(f ) =

h [f0 + 4f1 + f2 ] 3

(3.22)

La expresi´on del error en la interpolaci´on viene dado por ESS =

f 000 (ξ) (x − x0 )(x − x1 )(x − x2 ) 6

(3.23)

A diferencia del m´etodo del trapecio simple, para el m´etodo de Simpson simple no es simple determinar el error para la aproximaci´on de la integral a partir de (3.23). Para este objetivo vamos a usar otro m´etodo que pasamos a describir a continuaci´on. Primero consideramos el intervalo [−h, h] y la

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

68

siguiente funci´on que describe el error para cada valor de h: Z h E(h) = f (x) dx − M SS(f ) −h

Z

(3.24)

h

h f (x) dx − [f (−h) + 4f (0) + f (h)] 3 −h

=

A partir de (3.24) se calculan las derivadas de E que se muestran en la Tabla 3.2. Tabla 3.2: Derivadas de la funci´on de error E dada por (3.24) Z

h

h f (x) dx − [f (−h) + 4f (0) + f (h)] 3 −h

E(h)

=

E 0 (h)

=

2 2 4 h f (h) + f (−h) − f (0) − [−f 0 (−h) + f 0 (h)] 3 3 3 3

E 00 (h)

=

1 0 1 h f (h) − f 0 (−h) − [f 00 (−h) + f 00 (h)] 3 3 3

E 000 (h)

h = − [−f 000 (−h) + f 000 (h)] 3

1 h E (iv) (h) = − [−f 000 (−h) + f 000 (h)] − [f (iv) (−h) + f (iv) (h)] 3 3 E (v) (h)

2 2 h = − f (iv) (−h) − f (iv) (h) − [−f (v) (−h) + f (v) (h)] 3 3 3

De la Tabla 3.2 se observa que E(0) = 0,

E 0 (0) = 0,

E 00 (0) = 0,

4 E 000 (0) = 0, E (iv) (0) = 0, E (v) (0) = − f (iv) (0). 3 Usando E (iv) de la Tabla 3.2 y acomodando como sigue  000   (iv)  2h f (h) − f 000 (−h) f (h) + f (iv) (−h) (iv) E (h) = − + 3 2h 2

(3.25)

(3.26)

Exigiendo que f (iv) exista y sea continua en [−h, h], entonces podemos aplicar el Teorema del Valor Medio y Teorema del Valor Intermedio en (3.26) para garantizar la existencia de ξ1 , ξ2 ∈ h−h, hi tal que E (iv) (h) = −

2h (iv) [f (ξ1 ) + f (iv) (ξ2 )] 3

(3.27)

´ 3.2. METODO DE SIMPSON

69

acomodando de forma adecuada (3.27) con la finalidad de aplicar nuevamente el Teorema del Valor Intermedio, se tiene   4h f (iv) (ξ1 ) + f (iv) (ξ2 ) (iv) (3.28) E (h) = − 3 2 nuevamente, como f (iv) es continua, entonces en (3.28) existe ξ ∈ h−h, hi tal que 4h (3.29) E (iv) (h) = − f (iv) (ξ) 3 Ahora, integrando E (iv) de (3.29) y usando (3.25) se obtiene E(h) = −

f (iv) (ξ)h5 90

(3.30)

Por tanto, a partir de (3.19) junto a (3.21) y (3.30) se tiene la siguiente igualdad Z x2 h f (4) (ξ)h5 f (x) dx = [f (x0 ) + 4f (x1 ) + f (x2 )] − (3.31) 3 90 x0

3.2.2.

Compuesto

De una forma an´aloga al m´etodo del trapecio compuesto, vamos a considerar una partici´on regular del intervalo [a, b], es decir: a = x0 < x1 < ... < xn−1 < xn = b, b−a donde h = . Dado que para aplicar el m´etodo de Simpson simple necen sitamos 3 puntos: xi−1 , xi , xi+1 , el n´ umero de intervalos n est´a restringido a que sea un n´ umero natural par. Veamos el procedimiento: n

Z

xn

f (x)dx = x0

n

2 Z X

i=1

x2i

f (x)dx =

x2i−2

2 X

MSSi (f ) + ESC

i=1

donde Z

x2i

M SSi (f ) =

f (x) dx = x2i−2

h (f2i−2 + 4f2i−1 + f2i ), 3

i = 1, ...,

n 2

n

ESC = −

2 X h5

i=1

90

f (4) (ξi )

(3.32)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

70

El M´ etodo de Simpson Compuesto es definido como: n

M SC(f ) =

2 X

n

M SSi (f ) =

i=1

2 X h

i=1

3

(f2i−2 + 4f2i−1 + f2i )

(3.33)

realizando la suma en (3.33), el m´etodo de Simpson compuesto queda   n n −1 2 2 X X h M SC(f ) = f0 + 4 f2i−1 + 2 f2i + fn  (3.34) 3 i=1 i=1 De un modo an´alogo para el an´alisis de error para el m´etodo del trapecio compuesto, vamos a encontrar una forma m´as adecuada para el error para el m´etodo de Simpson compuesto dado por (3.32). Usando (3.30) y sabiendo que f (iv) es continua en el intervalo cerrado y acotado [a, b], por tanto, existen x∗ , xˆ ∈ [a, b] donde f (iv) alcanza su m´ınimo y m´aximo respectivamente, es decir f (iv) (x∗ ) ≤ f (iv) (x) ≤ f (iv) (ˆ x) para todo x ∈ [a, b]

(3.35)

n en particular, para ξi ∈ [x2i−2 , x2i ] ⊂ [a, b], (i = 1, ..., ) tambi´en se cumple 2 f (iv) (x∗ ) ≤ f 00 (ξi ) ≤ f (iv) (ˆ x) sumando para i = 1, ...,

(3.36)

n en (3.36) se obtiene 2 n

2 X n (iv) ∗ n f (x ) ≤ f (iv) (ξi ) ≤ f (iv) (ˆ x) 2 2 i=1

es decir

(3.37)

n

2 f (iv) (x∗ ) ≤

2 X

f (iv) (ξi )

i=1

≤ f (iv) (ˆ x)

n

dado que f (iv) es una funci´on continua en [a, b], entonces aplicamos de nuevo el Teorema del Valor Intermedio, es decir, existe ξ ∈ ha, bi tal que n

2 f (iv) (ξ) =

2 X

f (iv) (ξi )

i=1

n

(3.38)

´ 3.2. METODO DE SIMPSON

71

Reemplazando adecuadamente (3.38) en la expresi´on dado por (3.32), obtenemos n

ESC = −

2 X f (iv) (ξi )

i=1

90

h5 = −n

h5 (iv) f (ξ) para alg´ un ξ ∈ ha, bi. (3.39) 180

o de forma equivalente ESC = −n

h5 (iv) (b − a)5 (iv) f (ξ) = − f (ξ) para alg´ un ξ ∈ ha, bi. (3.40) 180 180n4

Ejemplo 3.2.1. Aplice el MSS a la siguiente integral Z e x2 ln(x) dx. 1

Luego, estime el n´ umero de subintervalos n para alcanzar una aproximaci´on −2 de 10 al usar el MSC. Resoluci´ on: Aplicando el M SS a f (x) = x2 ln(x) se tiene (e − 1) (f (1) + 4f ((1 + e)/2) + f (e)) = 4.57135. 6 Ahora, para estimar el n´ umero de subintervalos n, debemos acotar f (iv) seg´ un lo sugiere la f´ormula dada en (3.40). A partir de f obtenemos: M SS(f ) =

f (iv) (x) = −

2 , x2

y como 1 ≤ x ≤ e, entonces |f (iv) (x)| ≤ M = 2. Por tanto, de (3.40) se obtiene (e − 1)5 (iv) |ESC| = f (ξ) ≤ M 180n4

(e − 1)5 180n4

y queremos una tolerancia de error de tol = 10−2 , por tanto, exigimos que: (e − 1)5 ≤ tol M 180n4 de donde se obtiene r

M (e − 1)5 ≈ 2.01980, 180 tol as´ı, debemos tomar n = 4 porque debe ser un n´ umero par. n≥

4

La Tabla 3.3 muestra los c´alculos respectivos.

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

72

Tabla 3.3: Resultados del MSC para f (x) = x2 ln(x). i xi f (xi ) Peso Peso f (xi ) 0 1.00000 0.00000 1 0.00000 1 1.42957 0.73036 4 2.92142 2 1.85914 2.14337 2 4.28673 3 2.28871 4.33717 4 17.34869 4 2.71828 7.38906 1 7.38906

De la Tabla 3.3 se obtiene M SC(f ) = 4.57434, y sabemos que el valor exacto de la integral es Z e 2e3 1 f (x) dx = + = 4.57456, 9 9 1 por tanto, se puede calcular el error exacto Z e |ESC| = f (x) dx − M SC(f ) = 0.00023 < 0.01. 1

Z Ejemplo 3.2.2. Determine el n´ umero de intervalos n para aproximar

2

e2x sen(3x) dx

0

con una exactitud de 10−4 usando el m´etodo del trapecio compuesto y m´etodo de Simpson compuesto. Resoluci´ on: Sea f (x) = e2x sen(3x), luego: f 0 (x) f 00 (x) f 000 (x) f (iv) (x)

= = = =

e2x (2sen(3x) + 3cos(3x)) e2x (−5sen(3x) + 12cos(3x)) e2x (−46sen(3x) + 9cos(3x)) e2x (−119sen(3x) − 120cos(3x))

Para el m´etodo del trapecio se tiene que: ET (f ) = −

b − a 2 00 h f (ξ) para alg´ un ξ ∈ ha, bi. 12

Por tanto, debemos acotar f (00) (x) para todo 0 ≤ x ≤ 2. Por dato, obtenemos lo siguiente: p |e2x | ≤ e4 , | − 5sen(3x) + 12cos(3x)| ≤ (5)2 + (12)2 ,

3.3. CUADRATURA GAUSSIANA

73

entonces: |f (00) (x)| ≤ M p donde M = e4 (5)2 + (12)2 = 709.77595. Ahora, acotamos ET (f ), es decir: 22 (00) 2M |ET (f )| = 2 f (ξ) ≤ 2 ≤ 10−4 12n2 3n de donde obtenemos: r n≥

2M 104 = 2175.2792, 3

por tanto, para el m´etodo del trapecio tenemos que tomar n = 2176. Para el m´etodo de Simpson compuesto, debemos acotar f (iv) (x) para todo 0 ≤ x ≤ 2. De esto, tenemos lo siguiente: p |e2x | ≤ e4 , | − 119sen(3x) − 120cos(3x)| ≤ (119)2 + (120)2 , entonces: |f (iv) (x)| ≤ M p donde M = e4 (119)2 + (120)2 = 9227.0874. Ahora, acotamos ES(f ), es decir: 24 8M (iv) f (ξ) ≤ ≤ 10−4 |ES(f )| = 2 4 180n 45n4 de donde obtenemos: r n≥

4

8M 104 = 63.640785 45

por tanto, para el m´etodo de Simpson tenemos que n = 64.

3.3. 3.3.1.

Cuadratura Gaussiana Polinomios Ortogonales

1. Polinomios de Legendre: I = [−1, 1], w(x) = 1, L0 (x) = 1, L1 (x) = x, 2n + 1 n Ln+1 (x) = xLn (x) − Ln−1 (x) ∀n ≥ 1. n+1 n+1

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

74

2. Polinomios de Laguerre: I = [0, ∞i, w(x) = e−x , L0 (x) = 1, L1 (x) = 1 − x, n (2n + 1 − x) Ln (x) − Ln−1 (x) ∀n ≥ 1. Ln+1 (x) = n+1 n+1 3. Polinomios de Hermite: I = h−∞, +∞i, 2 w(x) = e−x , H0 (x) = 1, H1 (x) = 2x, Hn+1 (x) = 2xHn (x) − 2nHn−1 (x) ∀n ≥ 1. 4. Polinomios de Tchebyshev de primera especie: I = [−1, 1], 1 , 1 − x2 T0 (x) = 1, T1 (x) = x, Tn+1 (x) = 2xTn (x) − Tn−1 (x) ∀n ≥ 1. w(x) = √

3.4.

Ejercicios Propuestos

1. Dados los siguientes puntos de una funci´on y = f (x): x y

0 0.78

1 2.04

2 3.71

3 4.11

4 3.89 Z

determine un valor aproximado de la integral

4

f (x)dx. 0

a) Por la regla del Trapecio b) Por la regla de Simpson 2. Cu´antos intervalos ser´an precisos para aproximar la integral con un error menor que 5×10−4 , donde: Z 2 exp(−x) dx x 1

3.4. EJERCICIOS PROPUESTOS

75

a) Por la regla del Trapecio b) Por la regla de Simpson 3. Aproximar la integral: Z 0

1

4 dx 1 + x2

utilizando el m´etodo de Romberg, con una tolerancia de 5 × 10−9 y 20 iteraciones. R1 4. Deduzca la f´ormula de Newton - Cotes para 0 f (x) dx, usando polinomios de interpolaci´on de Lagrange con nodos 0, 31 , 23 , 1. Use este resultado para evaluar la integral cuando f (x) = sen(πx). 5. Encuentre la f´ormula Z

1

f (x) dx ≈ A0 f (0) + A1 f (1) 0

que resulta exacta para todas las funciones de la forma f (x) = a exp(x) + bcos( Z

πx ) 2

b

f (x) dx apoy´andose en la

6. Determine la regla de integraci´on para a

regla de cuadratura gaussiana: Z 1 1 1 f (x) dx ≈ f (− √ ) + f ( √ ) 3 3 −1

76

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

Cap´ıtulo 4 Sistemas de ecuaciones lineales A continuaci´on consideremos el siguiente sistema lineal:     x1 a11 a12 · · · a1,n−1 a1,n     a21 a22 a2,n−1 a2,n    x2      ..    ...  .  =        an−1,1 an−1,2 an−1,n−1 an−1,n   xn−1   xn an,1 an,2 · · · an,n−1 an,n

b1 b2 .. . bn−1 bn

      

que escrito en forma matricial: Ax = b

(4.1)

donde A ∈ Rn×n , x ∈ Rn y b ∈ Rn .

4.1.

M´ etodos Directos

4.1.1.

Eliminaci´ on Gaussiana

4.1.2.

Factorizaci´ on LU

Recordemos lo siguiente: Una matriz L es llamada Triangular  a11 0 0  a21 a22 0   L=   an−1,1 an−1,2 an−1,3 an,1 an,2 an,3 77

Inferior cuando tiene la forma:  ··· 0 0 0 0    ..  .  · · · an−1,n−1 0  · · · an,n−1 an,n

78

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES Una matriz U es llamada Triangular Superior cuando tiene la forma: 

a11 a12 · · · a1,n−2 a1,n−1 a1,n  0 a22 a2,n−2 a2,n−1 a2,n   . .. U =   0 0 0 an−1,n−1 an−1,n 0 0 ··· 0 0 an,n

      

Consideremos la siguiente matriz  A=

0 1 1 0

 (4.2)

observamos que det(A) = −1 6= 0. Deseamos calcular L y U tal que A = LU , es decir      0 1 l11 0 u11 u12 A= = (4.3) 1 0 l21 l22 0 u22 De (4.3) obtenemos el siguiente sistema de ecuaciones: 0 1 1 0

= = = =

l11 u11 l11 u12 l21 u11 l21 u12 + l22 u22

(4.4) (4.5) (4.6) (4.7)

De (4.4) se obtiene: l11 = 0 o u11 = 0. Si l11 = 0 entonces contradecimos (4.5). Si u11 = 0 entonces contradecimos (4.6). Por tanto, podemos concluir que no existen matrices L y U tal que A = LU , a pesar que la matriz A es no singular. El ejemplo anterior nos motiva a establecer condiciones sobre una matriz A tal que pueda existir la factorizaci´on LU . Antes de estudiar las t´ecnicas que nos permitan calcular la factorizaci´on LU de una matriz A, debemos saber primero qu´e condiciones debe cumplir la matriz A para que exista la factorizaci´on LU .

´ 4.1. METODOS DIRECTOS

79

Factorizaci´ on de Doolitle Factorizaci´ on de Crout Factorizaci´ on de Cholesky

4.1.3.

Factorizaci´ on QR

Processo de ortogonalizaci´ on de Gram-Schmidt Dados los vectores {v1 , v2 , ..., vn } linealmente independientes, el proceso de ortogonalizaci´on de Gram-Schmidt permite construir los vectores {q1 , q2 , ..., qn } que son ortonormales, del modo siguiente: u1 = v1

u1 ||u1 || u2 q2 = ||u2 || u3 q3 = ||u3 || .. . un qn = ||un || q1 =

u2 = v2 − hq1 , v2 iq1 u3 = v3 − hq1 , v3 iq1 − hq2 , v3 iq2 .. . un = vn − hq1 , vn iq1 − · · · − hqn−1 , vn iqn−1 del sistema anterior podemos obtener: v1 v2 v3 v4 .. .

= ||u1 ||q1 = ||u2 ||q2 + hq1 , v2 iq1 = ||u3 ||q3 + hq1 , v3 iq1 + hq2 , v3 iq2 = ||u4 ||q4 + hq1 , v4 iq1 + hq2 , v4 iq2 + hq3 , v4 iq3

vn = ||un ||qn +

n−1 X

hqi , vn iqi

i=1

de donde reordenando adecuadamente se obtiene: A = [v1 | v2 | v3 | v4 | · · · | vn ] reemplazando las expresiones para cada columna vi (i = 1, 2, 3, ..., n): A = ||u1 ||q1 ||u2 ||q2 + hq1 , v2 iq1 "

2 X ||u ||q + hq , v iq 3 3 i 3 i ··· i=1

# n−1 X hqi , vn iqi ||un ||qn + i=1

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

80

que puede escribirse de la forma siguiente:  

||u1 || 0 ··· 0  0 ||u2 || 0  A = [q1 |q2 | · · · |qn ]  . . {z } | . Q 0 0 · · · ||un ||

hq1 , v2 i ||u1 ||

1

   0     0  

1 0

0

0

|

hq1 , v3 i hq1 , vn−1 i ··· ||u1 || ||u1 || hq2 , v3 i hq2 , vn−1 i ||u2 || ||u2 || 1 .. . 0 ··· 1

{z

de esta forma hemos conseguido la factorizaci´on: A = QR donde Q es una matriz ortonormal (QT Q = I) e R es una matriz triangular superior. En la pr´actica se procede calcular Q calcul´andose {q1 , q2 , ..., qn } mediante el proceso de Gram-Schmidt y luego calcular R = QT A. Ejemplo 4.1.1. Considere el siguiente circuito: 15 Ω

ig

6Ω

12 Ω I2

+ −

125 V

I1 20

Ω

5Ω

13 Ω I3

         }

R

2Ω



i0

Figura 4.1: Circuito para las preguntas 3, 4 y 5. Escribir un sistema lineal de 3 ecuaciones en funci´on de I1 , I2 , I3 con el objetivo de calcular ig e i0 . Resoluci´ on:

´ 4.1. METODOS DIRECTOS

81

Usando la ley de Kirchof, se obtiene el siguiente sistema para el circuito mostrado:      13 −6 −5 I1 125  −6 66 −20   I2  =  0  (4.8) −5 −20 25 I3 0 Luego, usando el proceso de ortonormalizaci´on de Gram-Schmidt se tiene el siguiente cuadro: i ui

qi 



 0.8571946 u1 =  −0.3956283  q1 = ||u1 || −0.3296902     15.13913 0.2340525 u2 2 u2 = v2 − hq1 , v1 i =  56.243478  q2 = =  0.8695298  ||u2 || −28.130435 −0.4348993     5.8454677 0.4587339 u3 =  0.2956285  3 u3 = v3 − hq1 , v3 i − hq2 , v3 i =  3.7670792  q3 = ||u3 || 10.677721 0.8379540  13 1 u1 = v1 =  −6  −5

luego, las matrices Q y R son:   0.8571946 0.2340525 0.4587339 Q = [q1 | q2 | q3 ] =  −0.3956283 0.8695298 0.2956285  −0.3296902 −0.4348993 0.8379540   15.165751 −24.66083 −4.6156633 0 64.682637 −29.433341  R = QT A =  0 0 12.74261 Usando factorizaci´on QR el sistema Ax = b se convierte en Rx = ˆb = QT b, es decir:      15.165751 −24.66083 −4.6156633 x1 107.14933  0 64.682637 −29.433341   x2  =  29.256558  0 0 12.74261 x3 57.341744 y luego procedemos por sustituci´on regresiva, es decir: ˆb3 x3 = = 4.5 R33 x2 = x1 =

ˆb2 − R23 x3 = 2.5 R22

ˆb1 − R12 x2 − R13 x3 = 12.5 R11

82

4.2.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

M´ etodos Indirectos

Tambi´en llamados iterativos porque iniciando en un punto x0 ∈ Rn se generan una secuencia de puntos xk que esperamos converjan a la soluci´on del sistema lineal dado en (4.1). Los m´etodos a estudiar est´an basados en el m´etodo del punto fijo porque tienen la forma general siguiente: x = Mx + d

(4.9)

donde M es conocida como Matriz de iteraci´ on y d es llamado Vector de iteraci´ on. Bajo condiciones adecuadas, el sistema (4.9) genera una secuencia de puntos xk dados por xk+1 = M xk + d (4.10) que convergen a la soluci´on del sistema lineal equivalente (4.1).

4.2.1.

M´ etodo de Jacobi

4.2.2.

M´ etodo de Gauss-Seidel

4.2.3.

M´ etodo de relajaci´ on

Cap´ıtulo 5 M´ etodos num´ ericos para problemas de valor inicial Un Problema de Valor Inicial tiene la forma siguiente: 0 y (t) = f (t, y), t0 ≤ t ≤ T y(t0 ) = y0

(5.1)

En este cap´ıtulo estudiaremos t´ecnicas num´ericas que permitan aproximar la soluci´on exacta de (5.1). Para este fin es necesario saber si (5.1) tiene soluci´on y si ´esta es u ´nica. La teor´ıa de ecuaciones diferenciales ordinarias establece el siguiente Teorema de existencia y unicidad de soluci´on para (5.1), conocido son como teorema de Picard, el cual establece que si las funciones f y ∂f ∂y continuas en un rect´angulo cerrado D del plano XY con lados paralelos a los ejes, y el punto t0 es un punto interior de D, entonces existe un n´ umero real positivo H con la propiedad de que el problema (5.1) tiene soluci´on u ´nica en el intervalo | t − t0 |≤ H. Si se debilita la condici´on de que la derivada parcial de f respecto de y sea continua en el rect´angulo y se exige que satisfaga la desigualdad conocida como condici´ on de Lipschitz en la variable y : | f (t, y1 ) − f (t, y2 ) ≤ K | y1 − y2 |

(5.2)

en los puntos del recinto D, se consigue una versi´on m´as fuerte del teorema de Picard, que se conoce como el teorema de Picard - Lindel¨ of cuyo enunciado es como sigue: Si f es continua en t e y en el rect´angulo cerrado D definido por: | t − t0 |≤ a, | y − y0 |≤ b y existe una constante K (constante de Lipschitz) independiente de t, y1 e y2 tal que si (t, y1 ), (t, y2 ) son puntos interiores de D se cumple que: | f (t, y1 ) − f (t, y2 ) |≤ K | y1 − y2 | 83

(5.3)

´ ´ 84CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL (condici´on de Lipschitz), entonces la ecuaci´on diferencial (5.1) tiene soluci´on u ´nica y(t) definida en un entorno del punto t0 (| t − t0 |≤ H), tal que y(t0 ) = y0 , siendo H = m´ın( Mb , a) donde M es la cota superior en D de la funci´on f en D, es decir | f (t, y) |≤ M, ∀ (t, y) ∈ D. Este valor de H define el entorno del punto x0 para el cual existe soluci´on u ´nica para el problema (5.1). H es una cantidad que debemos tener siempre presente, pues de no hacerlo, s´olo obtendremos resultados que no tienen ning´ un sentido. Con este objetivo, primero consideramos una partici´on regular del intervalo [0, T ], es decir: 0 = t0 < t1 < ... < ti < ... < tn−1 < tn = T

(5.4)

tn − t0 para todo i = 0, 1, ..., n − 1. donde h = ti+1 − ti = n Los m´etodos a estudiar est´an basados en el segundo teorema fundamental del c´alculo, del cual sabemos que se cumple: Z ti+1 y 0 (x) dx = y(ti+1 ) − y(ti ) ti

de donde se obtiene: Z

ti+1

y(ti+1 ) = y(ti ) +

y 0 (x) dx

(5.5)

ti

De (5.5) concluimos que para conocer el valor de y(ti+1 ) es necesario conocer el valor de y(ti ) y el valor de la integral. Del cap´ıtulo anterior, conocemos t´ecnicas num´ericas que permiten aproximar el valor de la integral, as´ı, obtenemos valores aproximados para y(ti+1 ) que ser´an denotados por yi+1 (i = 0, ..., n − 1). Adem´as, usaremos la notaci´on fi para representar el valor f (ti , yi ) Una forma sencilla de aproximar la integral bajo la curva de f (t, y(t)) en el intervalo [ti , ti+1 ] es considerar el a´rea del rect´angulo mostrado en la Figura 5.1 con altura fi , es decir: Z ti+1 y 0 (x) dx ≈ (ti+1 − ti )fi ti

Por tanto, la identidad (5.5) da la siguiente igualdad: yi+1 = yi + hfi = yi + hf (ti , yi )

(5.6)

El m´etodo mostrado en (5.6) es llamado M´ etodo de Euler expl´ıcito.

85

´ Figura 5.1: Area del rect´angulo con altura fi .

Nada nos impide considerar ahora el rect´angulo considerado en la Figura 5.2 de altura fi+1 . As´ı, la identidad (5.5) da la siguiente igualdad:

yi+1 = yi + hfi+1 = yi + hf (ti+1 , yi+1 )

donde (5.7) es llamado M´ etodo de Euler impl´ıcito.

(5.7)

´ ´ 86CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL

´ Figura 5.2: Area del rect´angulo con altura fi+1 .

Es m´as, definiendo:

ti+ 1 = 2

ti + ti+1 , 2

Considerando la Figura 5.3:

yi+ 1 = y(ti+ 1 ), 2

2

fi+ 1 = f (ti+ 1 , yi+ 1 ) 2

2

2

´ 5.1. METODOS DE TAYLOR

87

´ Figura 5.3: Area del rect´angulo con altura fi+ 1 . 2

obtenemos el siguiente esquema num´erico: yi+1 = yi + hfi+ 1 = yi + hf (ti+ 1 , yi+ 1 ) 2

2

2

(5.8)

Hasta ahora, los m´etodos descritos por (5.6), (5.7) y (5.8) han sido obtenidos de forma intuitiva y geom´etrica, bas´andonos en el hecho de que la integral de una funci´on positiva coincide con el a´rea bajo la curva. Una deducci´on m´as formal ser´a presentada en la secci´on siguiente.

5.1.

M´ etodos de Taylor

Sea I ⊂ R un intervalo abierto y f : I → R una funci´on tal que existen f 0 , f 00 , ..., f (n) y f (n+1) para alg´ un entero positivo n. Dados x, x0 ∈ I, existe s entre x y x0 tal que f (t) = Pn,x0 (t) + Rn+1,x0 (t)

(5.9)

donde Pn,x0 (t) es llamado polinomio de Taylor de orden n y viene dado por Pn,x0 (t) = f (x0 ) + f 0 (x0 )(t − x0 ) +

f (00) (x0 ) f (n) (x0 ) (t − x0 )2 + ... + (t − x0 )n 2! n! (5.10)

´ ´ 88CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL adem´as Rn+1,x0 (t) =

f (n+1) (s) (t − x0 )n+1 (n + 1)!

(5.11)

conocido como Resto en la forma de Lagrange. Ejemplo 5.1.1. √ Construir la f´ormula de Taylor para el polinomio de grado 4 de f (x) = 3 + x alrededor√del punto 1 y obtener una cota del error cometido al aproximar el valor 5. Resoluci´ on: La f´ormula de Taylor para el polinomipo de grado 4 alrededor del punto 1 viene dado por: 4 X f (i) (1) P4 (x) = (x − 1)i i! i=0 Calculamos las derivadas: √ f (x) = 3+x f 0 (x)

=

f 00 (x)

=

f 000 (x)

=

f (iv) (x) = f (v) (x) =

⇒ f (1) = 2

1 √ 2 3+x 1 − p 4 (3 + x)3 3 p 8 (3 + x)5 15 − p 16 (3 + x)7 105 p 32 (3 + x)9

⇒ f 0 (1) =

1 4

⇒ f 00 (1) = − ⇒ f 000 (1) =

1 32

3 256

⇒ f (iv) (1) = − ⇒ f (v) (1) =

15 2048

105 . 16384

Por tanto: x − 1 (x − 1)2 (x − 1)3 5 (x − 1)4 P4 (x) = 2 + − + − 4 64 512 16384 luego P4 (5) =

179 = 2.79688, 64

encuanto el valor exacto es √ 5 = 2.23607.

5.2. EJERCICIOS PROPUESTOS

89

3

2.8

2.6

2.4

2.2

2

1.8

1.6

1.4

1.2

1 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Figura 5.4: Polinomio de Taylor de orden 1

5.2.

Ejercicios Propuestos

1. Use el teorema de Picard-Lindel¨of para demostrar que el problema de valor inicial: y 0 = tan(y), y(0) = 0 tiene una soluci´on en el intervalo | x |< π4 . 2. Los m´etodos para resolver problemas de valor inicial tambi´en se pueden usar para calcular integrales, tanto definidas como indefinidas. Por ejemplo, podemos calcular: Z 2 exp(−s2 ) ds 0

resolviendo el problema de valor inicial: y 0 = exp(−x2 )

y(0) = 0

con t en el intervalo [0,2]. D´e la soluci´on a este problema usando el m´etodo de la serie de Taylor de orden 4. Al consultar una tabla de la funci´ on error Z t 2 exp(−s2 ) ds erf (t) = √ 2 0 obtenemos y(2) ≈ 0.8820813907.

´ ´ 90CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL 3. Deduzca las f´ormulas de Runge -Kutta de tercer orden: 1 y(x + h) = y(x) + (2k1 + 3k2 + 4k3 ) 9 donde:

   k1 = hf (x, y) k2 = hf (x + 21 , y + 12 k1 )   k = hf (x + 3 , y + 3 k ) 3 4 4 2

Demuestre que para la ecuaci´on diferencial y 0 = x + y concuerda con el m´etodo de la serie de Taylor de orden 3. 4. Calcule la soluci´on de: y 0 = −2xy 2

y(0) = 1

en x = 1 usando h = 0.25 y el m´etodo de Adams - Bashforth - Moulton en conjunci´on con el m´etodo de Runge - Kutta de cuarto orden. Muestre la soluci´on calculada usando cinco d´ıgitos significativos en 0.25, 0.5, 0.75, 1.0. Compare sus resultados con la soluci´on exacta: y(x) =

1 1 + x2

5. Utilice el m´etodo de Runge - Kutta de cuarto orden para aproximar la soluci´on exacta del problema de valor inicial siguiente: x2 y 00 − 2xy 0 + 2y = x3 ln(x), y(1) = 0, y 0 (1) = 0 tomando h = 0.05.

1 ≤ x ≤ 1.5

Bibliograf´ıa [1] D. G. Luenberger, Y. Ye, et al., Linear and nonlinear programming. Springer, third ed., 2008. [2] M. Dowell and P. Jarratt, “A modified regula falsi method for computing the root of an equation,” BIT Numerical Mathematics, vol. 11, no. 2, pp. 168–174, 1971.

91