• Author / Uploaded
• daniel alain ramon chavez

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 42 43 47 53

. . . . . .

57 58 59 61 63 68 71

3. Integraci´ on Num´ erica 73 3.1. M´etodo del Trapecio . . . . . . . . . . . . . . . . . . . . . . . 74 3.1.1. Simple . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3

´INDICE GENERAL

4

3.2.

3.3. 3.4. 3.5.

3.1.2. Compuesto . . . . . . . . . . . . . . . . . . . . M´etodo de Simpson . . . . . . . . . . . . . . . . . . . . 3.2.1. Simple . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. Compuesto . . . . . . . . . . . . . . . . . . . . 3.2.3. Error en las f´ormulas de Newton–Cotes cerradas 3.2.4. Error en las f´ormulas de Newton–Cotes abiertas 3.2.5. Error en la f´ormulas de tipo interpolatorio . . . Cuadratura Gaussiana . . . . . . . . . . . . . . . . . . 3.3.1. Polinomios Ortogonales . . . . . . . . . . . . . . Cuadratura Gaussiana . . . . . . . . . . . . . . . . . . Ejercicios Propuestos . . . . . . . . . . . . . . . . . . .

4. Sistemas de ecuaciones lineales 4.1. M´etodos Directos . . . . . . . . . . . . . . . . 4.1.1. Eliminaci´on Gaussiana . . . . . . . . . 4.1.2. Pivoteo Parcial . . . . . . . . . . . . . 4.1.3. Pivoteo Total . . . . . . . . . . . . . . 4.1.4. Factorizaci´on LU . . . . . . . . . . . . 4.1.5. Factorizaci´on por matrices elementales 4.1.6. 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 . . . . . . . . . . 5. M´ etodos num´ ericos para problemas 5.1. M´etodos de Taylor . . . . . . . . . 5.2. M´etodos de Runge-Kutta . . . . . . 5.2.1. Runge-Kutta de 2 pasos . . 5.2.2. Runge-Kutta de 3 pasos . . 5.2.3. Runge-Kutta de 4 pasos . . 5.3. M´etodos Multipasos . . . . . . . . 5.4. Ejercicios Propuestos . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

de valor inicial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

75 79 79 83 90 91 91 91 97 97 102

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

105 105 109 114 120 125 128 133 136 137 140 144

. . . . . . .

145 . 149 . 152 . 152 . 153 . 153 . 153 . 153

. . . . . . .

. . . . . . .

6. Pr´ acticas Calificadas 155 6.1. Examen Parcial 2015-2 . . . . . . . . . . . . . . . . . . . . . . 155 6.2. Segunda Pr´actica Calificada . . . . . . . . . . . . . . . . . . . 162 6.3. Cuarta Pr´actica Calificada, 2016-1 . . . . . . . . . . . . . . . 167

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 argumen7

´ 8CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES to, 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 de f 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. Utilice el m´etodo de Punto Fijo para determinar una aproximaci´on con un error absoluto inferior a ×10−5 , del u ´nico cero de la funci´on f (x) = 1 + x + exp(x) = 0

´ 1.5. METODO DEL PUNTO FIJO

37

que se sabe est´a en el intervalo [−2, −1]. Resoluci´ on: Para la funci´on f definimos para todo x ∈ [−2, −1] la siguiente funci´on: g(x) = −1 − ex . Observar que −2 ≤ g(−1) < g(x) < g(−2) ≤ −1, para todo x ∈ [−2, −1], es decir, g([−2, −1]) ⊂ [−2, −1]. Por otro lado: g 0 (x) = −ex ⇒ |g 0 (x)| ≤ e−1 < 1, para todo x ∈ [−2, −1]. Entonces, por el Teorema 1.8 podemos concluir que existe un u ´nico punto fijo xˆ de la funci´on g en el intervalo [−2, −1] y por tanto, xˆ ser´a tambi´en la u ´nica ra´ız de f en el mismo intervalo. Adem´as, el esquema iterativo del punto fijo dado por (1.45) genera una sucesi´on de puntos que convergen a xˆ para cualquier x0 ∈ h−2, −1i. La Tabla 1.10 muestra los resultados de aplicar 1.45 iniciando en x0 = −1.5 ∈ h−2, −1i hasta que |xn − xn−1 | < 10−5 . Tabla 1.10: M´etodo de n xn−1 1 -1.50000000 2 -1.22313016 3 -1.29430749 4 -1.27408761 5 -1.27968604 6 -1.27812461 7 -1.27855922 8 -1.27843818 9 -1.27847188

punto fijo para g(x) = −1 − ex . xn |xn − xn−1 | -1.22313016 0.27686984 -1.29430749 0.07117733 -1.27408761 0.02021989 -1.27968604 0.00559843 -1.27812461 0.00156143 -1.27855922 0.00043461 -1.27843818 0.00012104 -1.27847188 0.00003370 -1.27846250 0.00000939

Ejemplo 1.5.3. ¿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:

´ 38CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES 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: √ √ x2 + 1 2 −1 < 0 ≤ ≤

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

 = 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.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 .

´ 1.5. METODO DEL PUNTO FIJO

39 √

Tabla 1.11: 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.4. ¿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.

´ 40CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES 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.5 ⇒ x1 = g(x0 ) =

√

e−x0 ≈ 0.7788,

por tanto: n > 15.77, es decir, si tomamos n = 16 garantizamos la tolerancia de error pedida. La Tabla 1.12 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.12: 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 . Ejemplo 1.5.5. Considere la ecuaci´on: cos(x) = tan(x), de la cual se sabe que tiene una soluci´on xˆ = 0.6662394321 con un error de 0.9 × 10−9 . (a) Defina un m´etodo de punto fijo (diferente a Newton) y demuestre su convergencia (sin hacer iteraciones) en el intervalo que usted defina.

´ 1.5. METODO DEL PUNTO FIJO

41

(b) ¿Cu´antas iteraciones necesitar´ıa para obtener el error de 0.9 × 10−9 ? comenzando con x0 = 0.6?. (c) Use el m´etodo de Newton con el punto inicial x0 = 0.6 para aproximar la soluci´on. Resoluci´ on: π (a) Consideremos la funci´on g(x) = arctan(cos(x)) para todo 0 ≤ x ≤ . 4 En este invervalo la funci´on g es continua y adem´as cumple √ 2 ≤ cos(x) ≤ 1, 2 por tanto se obtiene 0 < arctan

√ ! 2 π ≤ arctan(cos(x)) ≤ arctan(1) = . 2 4

es decir, g([0, π/4]) ⊂ [0, π/4]. Por otro lado: g 0 (x) = −

sin(x) , 1 + cos2 (x)

π luego,para todo x ∈ h0, i se tiene 4 π  sin(x) 0 |g (x)| = − ≤ | sin(x)| ≤ sin < 1. 1 + cos2 (x) 4 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, π/4] y por tanto, xˆ ser´a tambi´en la u ´nica soluci´on de cos(x) = tan(x) en el mismo intervalo. (b) Iniciando en x0 = 0.6 se obtiene del item anterior x1 = g(x0 ) = arctan(cos(x0 )) = 0.68999971, y que λ = sin(π/4). Luego, a partir de (1.46) se obtiene ” n ” como sigue:   0.9 × 10−9 (1 − λ) ln |x1 − x0 | n> = 56.7. ln(λ) Por tanto, eligiendo n = 57 garantizamos la tolerancia pedida, es decir: |ˆ x − x57 | < 0.9 × 10−9 .

´ 42CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES (c) Usando el m´etodo de Newton se tiene que xk+1 = xk +

f (xk ) , para todo k ≥ 0, f 0 (xk )

donde f (x) = cos(x) − tan(x),

f 0 (x) = − sin(x) − sec2 (x).

Los resultados son mostrados en la Tabla 1.13. Tabla 1.13: M´etodo de k xk 0 0.600000 1 0.669464 2 0.666247 3 0.666239

1.6.

Newton para f (x) = cos(x) − tan(x). f (xk ) f 0 (xk ) 0.14119881 -2.03268565 -0.00722808 -2.24685197 -0.00001733 -2.23609379 -0.00000000 -2.23606798

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)).

(1.48)

Entonces, el problema ahora es encontrar xˆ ∈ Rn tal que F (ˆ x) = 0.

(1.49)

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

´ 1.6. METODOS PARA SISTEMAS DE ECUACIONES NO LINEALES 43

1.6.1.

M´ etodo de Newton

Este m´etodo hace uso del polinomio de Taylor de grado 1 para varias (k) (k) variables. Para ello elegimos un punto x(k) = (x1 , . . . , xn ) ∈ Rn suficientemente cercano a xˆ, es decir, xˆ = x(k) + h(k) para alg´ un h(k) ∈ Rn , luego, para cada i = 1, 2, . . . , n se tiene: 0 = fi (x(k) + h(k) ) = fi (x(k) ) + El Sistema 1.50 escrito  ∂f1 (k) (x ) . . .  ∂x1    ∂f2 (k)  (x ) . . .  ∂x1   .. ..  . .     ∂fn  (x(k) ) . . . ∂x1

∂fi (k) (k) ∂fi (k) (k) (x )h1 + . . . + (x )hn ∂x1 ∂xn

en forma matricial es:  ∂f1 (k) (x )  ∂xn    ∂f2 (k)   h(k)  f1 (x(k) ) 1 (x )  (k) (k)     ∂xn  f2 (x )   h2   .  = − .. ..   ..   . .   h(k) fn (x(k) ) n  ∂fn (k)  (x )  ∂xn

La matriz del Sistema 1.51 dado por  ∂f1 (k) ∂f1 (k) (x ) . . . (x )  ∂x1 ∂xn   ∂f2 (k)  ∂f2 (k)  (x ) . . . (x )  ∂x1 ∂xn  JF (x(k) ) =  .. .. ..  . . .     ∂fn ∂fn (k)  (x(k) ) . . . (x ) ∂x1 ∂xn

(1.50)

    

(1.51)

              

(1.52)

es llamada Jacobiano de la funci´on F dada en (1.48). Luego, usando (1.48), (k) (k) (1.52) y h(k) = (h1 , . . . , hn ) ∈ Rn en el Sistema 1.51 se obtiene JF (x(k) )h(k) = −F (xk )

(1.53)

Resolviendo el Sistema 1.53 se obtiene el valor de h(k) , luego, definimos la siguiente aproximaci´on de xˆ como sigue: x(k+1) = x(k) + h(k) .

(1.54)

´ 44CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES El proceso descrito se resume en el siguiente algoritmo. Algoritmo 4: M´etodo de Newton Entrada: x(0) , T ol. 1 inicio 2 k ← 0; 3 mientras condi¸c˜ao de parada n˜ao ´e satisfeita hacer 4 Obtener h(k) resolviendo el sistema: JF (x(k) )h(k) = −F (xk ) x(k+1) ← x(k) + h(k) ; k ← k + 1; fin mientras devolver x(k) ;

5 6 7 8

fin

Ejemplo 1.6.1. 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 e iniciando en x(0) = (1, 1.5). Resoluci´ on: Definimos: 1 f1 (x1 , x2 ) = x1 + ln(x1 x2 ) − 2 = 0, 2 f2 (x1 , x2 ) = x1 x22 + 5x2 − 3. y siendo F (x) = (f1 (x), f2 (x) donde x = (x1 , x2 ), se calcula el Jacobiano de la funci´on F :   1 1  1 + 2x1  2x2 JF (x) =   2 x2 2x1 x2 + 5 Iniciando con x(0) = (1, 1.5) en el Algoritmo 4 obtenemos:     −0.7972674 1.5 1/3 (0) (0) F (x ) = , JF (x ) = . 6.75 2.25 8. Resolviendo el sistema     1.5 1/3 −0.7972674 (0) h =− 2.25 8. 6.75

´ 1.6. METODOS PARA SISTEMAS DE ECUACIONES NO LINEALES 45 obtenemos que (0)

h

 =

0.7669457 −1.0594535



lo que implica x

(1)

=x

(0)

(0)

+h

 =

1.7669457 0.4405465

 .

Evaluamos el criterio de parada kh(0) k = 1.3079172 > 10−5 . Por tanto, repetimos el proceso hasta conseguir que kh(k) k < 10−5 o alcanzar un n´ umero m´aximo de iteraciones sin conseguir el criterio de parada. Las iteraciones para este ejemplo son resumidas en la Tabla 1.14 y se consigui´o la convergencia para k = 5.

k 0 1 2 3 4 5

Tabla 1.14: M´etodo de Newton para el Ejemplo 1.6.1. x(k) F (x(k) ) kh(k) k (k) (k) x2 x1 f1 (x(k) ) f2 (x(k) ) 1 1.5 -0.79726745 6.75 0 1.766946 0.440547 -0.35829764 -0.45433645 1.30791722 1.99078 0.503213 -0.00832713 0.020177 0.23244147 1.999988 0.499999 -0.00001557 -0.00000913 0.00975263 2 0.5 0 0 0.00001178 2 0.5 0 0 0

An´ alisis de convergencia Teorema 1.9. Sea F : Rn → Rn una funci´on de clase C 2 en un conjunto abierto y convexo D ⊂ Rn tal que x? ∈ D. Supongamos que JF (x? ) exista y que adem´as existen constantes positivas R, C y L tal que : kJF (x? )k ≤ C,

(1.55)

y kJF (x) − JF (y)k ≤ L kx − yk

para todo

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

(1.56)

Entonces existe r > 0 tal que para cualquier punto inicial x0 ∈ B(x? , r) la sucesi´on generada por el Algoritmo 4 converge para la soluci´on x? con: kxk+1 − x? k ≤ kxk − x? k2 .

(1.57)

´ 46CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES Demostraci´ on: Afirmamos que existe r > 0 tal que JF (x0 ) es invertible. En efecto, para cualquier x0 ∈ B(x? , R), de la hip´otesis se tiene que: kI − JF−1 (x? JF (x0 ))k = kJF−1 (x? )(J(x? ) − JF (x0 ))k, ≤ kJF−1 (x? )k kJF (x? ) − JF (x0 )k, ≤ CLkx? − x0 k.

(1.58)

Siguiendo [3], de demostrar´a que el radio espectral de la matriz kI−JF−1 (x? )JF (x0 )k es menor que 1, donde: ρ(P ) = ´ınf (kP k). k.k

  1 , obteniendo para cualquier x0 ∈ En (1.58) elegimos r = min R, 2CL B(x? , r): ρ(I − JF−1 (x? )JF (x0 )) ≤ kI − JF−1 (x? )JF (x0 )k ≤

1 . 2

(1.59)

Entonces la matriz I − (I − JF−1 (x? )JF (x0 )) = JF−1 (x? )JF (x0 ) es invertible, adem´as: 1 . (1.60) k(JF−1 (x? )JF (x0 ))−1 k ≤ −1 ? 1 − kI − JF (x )JF (x0 )k 1 > 0 tal que para cualquier x0 ∈ B(x? , r) se tiene Por tanto, existe r = 2CL que JF (x0 ) es invertible. Probaremos ahora que para cualquier x0 ∈ B(x? , r) se tiene que la sucesi´on generada por el Algoritmo 4 converge a x? . En efecto, a partir de (1.59) se cumple que:

kJF−1 (x? )k ≤ 2 kJF−1 (x? )k ≤ 2C. −1 ? 1 − kI − JF (x )JF (x0 )k

(1.61)

Como JF−1 (x0 ) = JF−1 (x0 )JF (x? )JF−1 (x? ), entonces de (1.60) y (1.61) obtenemos: kJF−1 (x? )k ≤ 2C. (1.62) kJF−1 (x0 )k ≤ 1 − kI − JF−1 (x? )JF (x0 )k Usando ahora el polinomio de Taylor alrededor de x0 , [?], se tiene que: F (x? ) = F (x0 ) + JF (x0 )(x? − x0 ) + r(x? − x0 ),

(1.63)

donde r(x? − x0 ) es el resto, el cual satisface: l´ım ?

x0 →x

r(x? − x0 ) = 0. kx? − x0 k2

(1.64)

´ 1.6. METODOS PARA SISTEMAS DE ECUACIONES NO LINEALES 47 De (1.64) y como F (x? ) = 0, obtenemos: kx1 − x? k ≤ kx0 − JF−1 (x0 ) − x? + JF−1 (x0 )F (x? )k, ≤ kJF−1 (x0 )k kF (x0 ) − F (x? ) + JF (x0 )(x0 − x? )k, ≤ CL kx0 − x? k2 .

(1.65)

Por tanto, (1.65) prueba (1.57) para k = 0. Adem´as, como x0 ∈ B(x? , r) podemos concluir de (1.65) que x1 ∈ B(x? , r). Asumiendo ahora que el teorema es v´alido para k y probando de um modo similar para k + 1, usando inducci´on se concluye la demostraci´on. Del Teorema 1.9 podemos deducir que uma desventaja del M´etodo de Newton, es que la existencia de la vecindad donde el M´etodo converge es peque˜ na. Algunas variaciones del M´etodo de Newton para resolver un sistema de ecuaciones no lineales pueden ser vistas en [3, 4, 5, 6, 7, 8].

1.6.2.

M´ etodo del punto fijo

Teorema 1.10. Sea D ⊂ Rn un conjunto cerrado y acotado. Sea G : D → Rn una funci´on continua. Supongamos que G(D) ⊂ D, entonces existe un punto fijo xˆ para G. M´as a´ un, si kJG(x)k ≤ λ < 1 para todo x ∈ D, entonces el punto fijo es u ´nico y el proceso iterativo dado por xn+1 = G(xn ),

para todo

n ≥ 0,

(1.66)

converge al u ´nico punto fijo xˆ para cualquier punto inicial x0 ∈ D. Es m´as, en la iteraci´on n se cumple tambi´en kˆ x − xn k ≤

λn kx1 − x0 k. 1−λ

(1.67)

Como consecuencia de (1.67) se tiene un estimado del n´ umero de iteraciones ” n ” para alcanzar una cierta tolerancia TOL. El c´alculo de ” n ” viene estimado por   TOL (1 − λ) ln kx1 − x0 k n≥ . (1.68) ln(λ) Ejemplo 1.6.2. Considere el sistema de ecuaciones no lineales: x21 − 10x1 + x22 + 8 = 0 x1 x22 + x1 − 10x2 + 8 = 0

´ 48CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES (a) Determinar una regi´on D ⊂ R2 que contiene una u ´nica soluci´on α = (α1 , α2 ) del sistema. (b) Proponga un m´etodo de punto fijo convergente para determinar la soluci´on α = (α1 , α2 ) ∈ D del sistema. Demuestre la convergencia sin iterar. ¿Se puede iniciar en (0, 0)? En caso afirmativo, estimar el n´ umero de −5 iteraciones para alcanzar una tolerancia menor a 10 . (c) Realice 2 iteraciones con el m´etodo propuesto y compare con el m´etodo de Newton. En ambos casos, inicie en (0, 0) (especifique cual de los m´etodos est´a usando). Resoluci´ on: (a) El sistema dado puede escribirse de la forma siguiente: x1 =

x21 + x22 + 8 10

x2 =

x1 x22 + x1 + 8 10

Por tanto, podemos considerar la siguiente funci´on de iteraci´on:   x21 + x22 + 8   10  G(x1 , x2 ) =    x1 x22 + x1 + 8 10 Considere la siguiente regi´on: D = {(x1 , x2 ) ∈ R2 / 0 ≤ x1 ≤ 1, 0 ≤ x2 ≤ 1}. Para todo (x1 , x2 ) ∈ D se cumple x21 + x22 + 8 ≤ 9, esto implica: 0 ≤ G1 (x1 , x2 ) ≤ 0.9. 0 ≤ x1 + 8 ≤ x1 x22 + x1 + 8 ≤ 10, esto implica: 0 ≤ G2 (x1 , x2 ) ≤ 1. por tanto se concluye que (G1 (x1 , x2 ), G2 (x1 , x2 )) ∈ D. Nos falta examinar la matriz jacobiana de G dada por  x  x  JG(x1 , x2 ) =  

1

2

5

5

x22 + 1 10

  2x1 x2  10

´ 1.6. METODOS PARA SISTEMAS DE ECUACIONES NO LINEALES 49 de donde  kJG(x1 , x2 )k∞ = m´ax

x1 + x2 x22 + 1 + 2x1 x2 , 5 10

 .

Se presentan dos casos: x1 + x2 2 , entonces, kJG(x1 , x2 )k∞ < < 1. 5 5 x22 + 1 + 2x1 x2 = , entonces 10

Si kJG(x1 , x2 )k∞ = Si kJG(x1 , x2 )k∞

kJG(x1 , x2 )k∞ =

(x1 + x2 )2 + 1 5 (x1 + x2 )2 + 1 − x21 ≤ < . 10 10 10

De los casos anteriores podemos concluir que   2 5 1 kJG(x1 , x2 k∞ < m´ax , = < 1, 5 10 2 1 as´ı, tenemos que λ = < 1. 2 Por tanto, del Teorema 1.10 podemos concluir que la funci´on G admite un u ´nico punto fijo xˆ ∈ D. (b) Por el item anterior, la funci´on G satisface las condiciones del Teorema 1.10, por tanto, la sucesi´on   (k+1) x  = G(x1(k) , x(k)  1 para todo k ≥ 0, 2 ), (k+1) x2 es decir  

 

(k+1)

x1

(k+1)

x2

(k) x1

2



(k) x2

2

+8 +   10  =   (k)  (k) 2 (k)  x1 x2 + x1 + 8 

       

10 converge al u ´nico punto fijo xˆ ∈ D para cualquier punto inicial x(0) ∈ D. Por otro lado, como x(0) = (0, 0) ∈ D, entonces la sucesi´on converge al u ´nico punto fijo xˆ. Luego, para estimar el n´ umero de iteraciones ” n ” tenemos:     0 0.8 (0) (1) x = ⇒x = ⇒ kx(1) − x(0) k∞ = 0.8, 0 0.8

´ 50CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES y del item anterior se tiene que λ = 0.5. Por tanto, usando (1.68) obtenido a partir del Teorema 1.10 tenemos que   −5 10 × 0.5 ln 0.8 ⇒ n ≥ 17.287712. n≥ ln(0.5) Esto u ´ltimo implica que debemos elegir n = 18 como el n´ umero de iteraciones para alcanzar una tolerancia menor a 10−5 . (c) Empezamos con el m´etodo del Punto Fijo usando el proceso iterativo descrito en el item anterior. Por tanto, si x(0) = (0, 0)T , entonces   0.8 (1) (0) x = G(x ) = 0.8   116    125  0.928 (2) (1)   = . x = G(x ) =  0.9312 582  625 Ahora usamos el m´etodo de Newton siguiendo el Algoritmo 4, para ello definimos:   2 x1 − 10x1 + x22 + 8 F (x1 , x2 ) = x1 x22 + x1 − 10x2 + 8 luego JF (x1 , x2 ) =

2x1 − 10

2x2

x22 + 1

2x1 x2 − 10

!

Iniciando en x(0) = (0, 0)T se tiene:     8 −10 0 (0) (0) F (x ) = , JF (x ) = , 8 1 −10 usando ahora (1.53)–(1.54) resulta     4 4    5   5  0.8 (0) (1) (0) (0)     h = ⇒x =x +h = = . 0.88 22  22  25 25 Para la segunda iteraci´on se tiene:     884 42 44 −  625   5 25   , JF (x(1) ) =  , F (x(1) ) =   1936   1109 1074  − 3125 625 125

´ 1.6. METODOS PARA SISTEMAS DE ECUACIONES NO LINEALES 51 usando ahora (1.53)–(1.54) resulta     0.1918 0.9918 (1) (2) (1) (1) h =− ⇒x =x +h = . 0.1117 0.9917

Tabla 1.15: Criterio de parada: kh(k) k∞ = kx(k) − x(k−1) k∞ < 10−5 . k Punto Fijo Newton (k) (k) (k) (k) (k) x1 x2 kh k∞ x1 x2 kh(k) k∞ 0 0.000000 0.000000 0.000000 0.000000 1 0.800000 0.800000 0.800000 0.800000 0.880000 0.880000 2 0.928000 0.931200 0.131200 0.991787 0.991712 0.191787 3 0.972832 0.973270 0.044832 0.999975 0.999969 0.008257 4 0.989366 0.989435 0.016534 1.000000 1.000000 0.000031 5 0.995783 0.995794 0.006417 1.000000 1.000000 0.000000 6 0.998319 0.998321 0.002536 7 0.999328 0.999329 0.001010 8 0.999732 0.999732 0.000403 9 0.999893 0.999893 0.000161 10 0.999957 0.999957 0.000064 11 0.999983 0.999983 0.000026 12 0.999993 0.999993 0.000010 13 0.999997 0.999997 0.000004

Ejemplo 1.6.3. Considere el sistema de ecuaciones no lineales: x21 + 16x22 − 9 = 0 50 x1 + x22 − 25 = 0 3 (a) Determine una regi´on D ⊂ R2 que contiene una u ´nica soluci´on del sistema. (b) Proponga un esquema iterativo de punto fijo y demuestre la convergencia sin iterar. ¿Se puede iniciar en (0, 0)? En caso afirmativo, estimar el n´ umero de iteraciones para alcanzar una tolerancia menor a 10−5 . (c) Realice las iteraciones con el esquema propuesto y compare con el m´etodo de Newton.

´ 52CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES Resoluci´ on: (a) El sistema dado puede escribirse de las siguiente manera. p 9 − x21 x2 = 4 (25 − x22 ) 50 Por lo tanto podemos considerar la siguiente funci´on de iteraci´on   (25 − x22 )  3  50   G(x1 , x2 ) =  p   9 − x21  4 x1 = 3

considere la siguiente region: D={(x1 ; x2 )R2 /0 ≤ x1 ≤ 2, 0 ≤ x2 ≤ 4 } 0 ≤ x2 ≤ 2 9 ≤ 25 − x22 ≤ 25 3(25−x2 ) 27 ≤ 50 2 ≤ 32 50 27 donde: 0 ≤ 50 ≤ x1 ≤ 32 ≤ 4 la segunda componente: 0 ≤ x1√≤ 2 √

9−x21 ≤ 34 16 √ donde:0 ≤ 45 ≤ 5 4

≤

√

9−x21 16

≤

3 4

≤2

entonces: G(D) ⊂ D (b) Empezamo con el metodo de punto fijo usando el proceso iterativo por tanto: x(0) = (α1 ; α2 ) 3(25 − x22 ) x1 = 50 p 9 − x21 x2 = 4 Ahora usamos el metodo de Newton para ello definimos:  2  x1 + 16x22 − 9 F (x1 ; x2 ) = 50 x + x22 − 25 3 1 ! 3(25−x22 ) G(x1 ; x2 ) = √ 50 2 9−x1 4

1.7. EJERCICIOS

53

JG(x1 ; x2 ) =

0 − √x1 4

( ||JG(x1 ; x2 )||∞ = max

9−x21

2 − 3x 25 0

!

||x1 || 3||x2 || ; p 25 4 9 − x21

)

(JG(xn1 ; xn2 )).hn = F (xn1 ; xn2 )

Definici´ on 1.2. Sea G una funci´on que define un m´etodo iterativo y que G(p) = p, G0 (p) = 0, G00 (p) = 0, . . . , G(m−1) (p) = 0 y G(m) 6= 0. Entonces se dice que el m´etodo de punto fijo asociado tiene orden de convergencia m. Ejemplo 1.6.4. Se proponen los siguientes m´etodos para calcular 211/3 . Clasifique por orden bas´andose para ello en la rapidez de convergencia y suponiendo que p0 = 1 en cada caso. 20pn−1 + (a) pn =

21 p2n−1

21

.

p3n−1 − 21 . (b) pn = pn−1 − 3p2n−1 (c) pn = pn−1 −  (d) pn =

1.7.

21 pn−1

p4n−1 − 21pn−1 . p2n−1 − 21

1/2 .

Ejercicios

1. 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. 2. 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.

´ 54CAP´ITULO 1. BUSQUEDA DE RA´ICES EN ECUACIONES NO LINEALES 3. 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 . 4. 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, . . ..

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

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

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

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

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

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

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

f ) f (x) = x2 − e−x .

6. Determine con un error absoluto menor de 0.001 la soluci´on de la ecuaci´on x − cos(x) = 0. 7. 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. 8. 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. 9. 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. 10. Use el m´etodo de Newton para obtener la menor ra´ız positiva de las siguientes ecuaciones: a) x/2 − tan(x) = 0. b) 2cos(x) − exp(x)/2 = 0. c) x5 − 6 = 0. 11. 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.

1.7. EJERCICIOS

55

12. 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. 13. 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. 14. Aplique el m´etodo de Newton-Raphson a la ecuaci´on x3 −2x2 −3x+10 = 0 con x0 = 1.9. Justifique los resultados obtenidos. 15. 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 π?. 16. Use el m´etodo del punto fijo para encontrar una ra´ız de f (x) = x3 + 4x2 − 10 en el intervalo [1, 2]. 17. Use el m´etodo del punto fijo para encontrar una ra´ız de f (x) = x − x4/5 − 2.

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

Cap´ıtulo 2 Interpolaci´ on polinomial Considere un conjunto de n + 1 puntos de la forma (xi , yi ) para i = 0, 1, . . . , n. El objetivo de este cap´ıtulo es encontrar un polinomio Pn de grado menor o igual a n tal que cumpla la siguiente condici´ on de interpolaci´ on: Pn (xi ) = yi ,

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

(2.1)

El 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.2)

Como Pn debe de satisfacer (2.1), entonces: a0 + a1 xi + a2 x2i + ... + an xni = yi , El conjunto de ecuaciones dado por (2.3) al queda como sigue:    1 x0 x20 x30 · · · xn0   1 x1 x2 x3 · · · xn   1 1 1     ..   .  2 3 n 1 xn xn xn · · · xn

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

(2.3)

ser escrito de forma matricial a0 a1 a2 .. . an





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

y0 y1 y2 .. .

      

(2.4)

yn

donde observamos que la matriz asociada al sistema lineal (2.4) es conocida como matriz de Vandermonde de orden n+1, cuyo determinante se calcula como: 1 x0 x2 x3 · · · xn 0 0 0 1 x1 x2 x3 · · · x n Y 1 1 1 (xj − xi ), (2.5) = .. . 1≤i≤j≤n+1 1 xn x2n x3n · · · xnn 57

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

58

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

Del an´alisis realizado en la secci´on anterior, ya es posible calcular los coeficientes del polinomio de interpolaci´on Pn resolviendo el sistema lineal (2.4). Ejemplo 2.1.1. Dados los puntos (−1, 2), (0, 1), (−2, 3), (1, 0). Calcule el polinomio de interpolaci´on. Resoluci´ on: Dados los 4 puntos (esto es, n = 3), procedemos a calcular los coeficientes ai (i = 0, 1, 2, 3) del polinomio de interpolaci´on P3 resolviendo el sistema lineal (2.4) dado por: 

 1 (−1) (−1)2 (−1)3 a0 2 3  1   0 0 0    a1  1 −2 (−2)2 (−2)3   a2 1 (1) (1)2 (1)3 a3





 2   1   =  ,   3  0

(2.6)

resolviendo el sistema anterior se obtiene: a0 = 1, a1 = −1, a2 = 0, a3 = 0, entonces, el polinomio de interpolaci´on P3 es: P3 (x) = 1 − x. La Figura 2.1 muestra los puntos dados con el polinomio de interpolaci´on P3 (x).

´ 2.2. METODO DE LAGRANGE

59

Figura 2.1: P3 (x) vs (xi , yi )

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: Li (x) =

n Y (x − xj ) (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

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

60

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.1). Por tanto, Pn es el polinomio de interpolaci´on dado en la forma de Lagrange. 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 (3.80):

L0 (x) =

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

L1 (x) =

(x + 1)x (x + 1)(x − 0) = (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.2 muestra los puntos dados con el polinomio de interpolaci´on P2 (x).

´ 2.3. METODO DE NEWTON

61

Figura 2.2: 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)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

62

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

´ 2.4. METODO DE DIFERENCIAS DIVIDIDAS

63

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.3.

Figura 2.3: 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)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

64

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)

´ 2.4. METODO DE DIFERENCIAS DIVIDIDAS

65

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.4: 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−xi )(x−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.

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

66

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

∆(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.5: Tabla de diferencias divididas para el ejemplo dado.

De la Figura 2.5 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.6 muestra el polinomio de interpolaci´on P3 con los puntos dados.

´ 2.4. METODO DE DIFERENCIAS DIVIDIDAS

67

Figura 2.6: 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.5 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

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

68

Figura 2.7: 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)

´ 2.5. ANALISIS DE ERROR

69

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!

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

70

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 P4 (x) = − +0.1836x+0.6450x x − −0.3191x x − x− −0.0505x x − x− 8 4 4 2 4 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.8 compara los valores de f y P4 para los puntos de la partici´on considerada.

Figura 2.8: An´alisis de error.

´ 2.6. FENOMENO RUNGE

71

En la Figura 2.8 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.9 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)

´ POLINOMIAL CAP´ITULO 2. INTERPOLACION

72

1

0.8

0.6

0.4

0.2

0

−0.2

−0.4 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

Figura 2.9: Fen´omeno Runge.

0.4

0.6

0.8

1

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). 73

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

74

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

75

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

76

[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

77

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

78 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

79

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

80 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

81

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

82

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

83

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

84

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

85

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

ESC = −

2 X f (iv) (ξi )

i=1

90

h5 (iv) h = −n f (ξ) para alg´ un ξ ∈ ha, bi. (3.39) 180 5

o de forma equivalente ESC = −n

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

−x Ejemplo R 1 3.2.1. Sea f (x) = 1 + e sin(4x). Aproximar el valor de la integral 0 f (x) dx usando la regla trapezoidal y Simpson simple. Compare los resultados obtenidos en cada aproximaci´on con el resultado exacto de la integral. Z 1 21e − 4 cos (4) − sin (4) f (x) dx = = 1.3082506046426 . . . 17e 0

Resoluci´ on: Dada nuestra funci´on f (x) = 1 + e−x sin(4x). Usando el M´etodo de Trapecio Simple: Z 1 1−0 {f (0) + f (1)}. f (x)dx = 2 0 Reemplazando: f (0) = 1 f (1) = 1 + e−1 sin (4) Z 1 1 f (x) dx = 1 + 1 + e−1 sin (4) ≈ 0.860793960475 2 0 Usando el M´etodo de Simpson Simple:     Z 1 1 1−0 f (x)dx = f (0) + 4f + f (1) . 3 2 0 Reemplazando: f (0)  = 1  1 f = 1 + e−0.5 sin(2), 2 f (1) = 1 + e−1 sin (4), Z

1

f (x)dx = 0

1−0 3

C´alculo de errores:

{1 + 4 [1 + e−0.5 sin (2)] + 1 + e−1 sin (4)} ≈ 2.64255166454.

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

86

1.308250604642 − 0.860793960475 ≈ 0.342026705418 Error (T S) = 1.308250604642 1.30825060464 − 2.64255166454 ≈ 1.0199124351 Error (SS) = 1.30825060464 Ejemplo 3.2.2. 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 de 10−2 al usar el MSC. Resoluci´ on: Aplicando el M SS a f (x) = x2 ln(x) se tiene M SS(f ) =

(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: 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 n≥

4

M (e − 1)5 ≈ 2.01980, 180 tol

as´ı, debemos tomar n = 4 porque debe ser un n´ umero par. La Tabla 3.3 muestra los c´alculos respectivos.

´ 3.2. METODO DE SIMPSON

87

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

Ejemplo 3.2.3. Determinar el n´ umero de sumandos necesarios, en las f´ormulas compuesta del Trapecio y de Simpson, para calcular con 6 cifras decimales exactas la siguiente integral. Z 2 ln (x)dx 1

Resoluci´ on: Para el c´alculo de errores o tolerancia tenemos: T ol = 10−6 (2) 2 f (θ) |ET (f )| = h (b−a) 12 |ES (f )| = Adem´as h =

b−a n

= n1 .

Pasamos a obtener las derivadas:

h4 (b−a) 180

(4) f (θ)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

88

f (1) (x) =

f (x) = ln (x);

f (3) (x) =

2 x3

;

1 x

;

f (2) (x) = − x12

f (4) (x) = − x64

Acotamos las derivadas respectivas: Usando el m´etodo del Trapecio. 1≤

x

≤2

⇒

−1 ≤ − x12 ≤ − 14

1 4

⇒

≤ f (2) (x) ≤ 1

Reemplazando: 1 n2 (12)

(1) f (2) (x) ≤ 10−6

1 n2 (12)

;

≤ 10−6

n = 288.67 ≈ 289

;

Usando el m´etodo de Simpson. 1≤

x

≤2

⇒

−6 ≤ − x64 ≤ − 61

Reemplazando: (4) 1 f (x) ≤ 10−6 (1) 4 n (180)

;

1(6) n4 (180)

⇒

≤ 10−6

1 6

≤ f (4) (x) ≤ 6.

;

n = 13.51 ≈ 14 Z

Ejemplo 3.2.4. 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.2. METODO DE SIMPSON

89

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 donde M = e4 decir:

p

(119)2 + (120)2 = 9227.0874. Ahora, acotamos ES(f ), es 24 (iv) ≤ 8M ≤ 10−4 f (ξ) |ES(f )| = 2 45n4 4 180n

de donde obtenemos: r n≥

4

8M 104 = 63.640785 45

por tanto, para el m´etodo de Simpson tenemos que n = 64. Ejemplo 3.2.5. Un cuerpo est´a sometido a la fuerza conservativa F (x) = x − x1 y se desplaza desde la posici´on x = 1 hasta x = 1.8. Determinar el trabajo realizado utilizando la regla de Simpson compuesta con paso h = 0.2. Calcular el error relativo y estimar un n´ umero de intervalos para obtener un −5 error de 10 . Aproximar al quinto decimal. Resoluci´ on: El trabajo pedido es:  2  r=1.8 Z 1.8 Z 1.8 1 x W (x) = F (x)·dx = (x− )·dx = − ln(x) = 0.5322133 ≈ 0.53221 x 2 1 1 r=1 Aproximamos el trabajo por el m´etodo de Simpson compuesto:

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

90 I=

0.2 (F (1) + 4F (2) + 2F (1.4) + 4F (1.6) + F (1.8)) = 0.5321693 ≈ 0.53217 3

El error relativo es: ER = −0.0000019 Para estimar el n´ umero de subintervalos, se tiene que:

|f

(iv)

120 (x)| = 5 ≤ 24, x

Y el error por el m´etodo de Simpson compuesto viene dado por: 4 h (b − a)5 (iv) 24(b − a)5 (iv) ≤ 10−5 |E| = (b − a)f (ζ) = f (ζ) ≤ 180 180n4 180n4 De donde se obtiene n ≥ 8.13, es decir, podemos tomar n = 9.

3.2.3.

Error en las f´ ormulas de Newton–Cotes cerradas

Tomando los nodos x0 = a < x1 < . . . < xn−1 < xn = b equidistantes: h=

b−a n

se tiene lo siquiente: Si n es par: b

Z

f (x)dx = a

n X i=0

h(n+3 (n+2) Ai f (xi ) + f (ξ) (n + 2)!

n

Z

t2 (t − 1) . . . (t − n)dt.

0

(3.41) Si n es par: Z

b

f (x)dx = a

n X i=0

h(n+2 (n+1) Ai f (xi ) + f (ξ) (n + 1)!

Z

n

t(t − 1) . . . (t − n)dt. 0

(3.42)

3.3. CUADRATURA GAUSSIANA

3.2.4.

91

Error en las f´ ormulas de Newton–Cotes abiertas

Tomando los nodos x−1 = a < x0 < x1 < . . . < xn−1 < xn < xn+1 = b equidistantes: b−a h= n+2 se tiene lo siquiente: Si n es par: Z

b

f (x)dx = a

n X i=0

h(n+3 (n+2) Ai f (xi )+ f (ξ) (n + 2)!

n+1

Z

t2 (t−1) . . . (t−n)dt.

−1

(3.43) Si n es par: Z

b

f (x)dx = a

n X i=0

h(n+2 (n+1) Ai f (xi ) + f (ξ) (n + 1)!

Z

n+1

t(t − 1) . . . (t − n)dt. −1

(3.44)

3.2.5.

Error en la f´ ormulas de tipo interpolatorio

El error en las f´ormulas de tipo interpolatorio tiene una expresi´on de la forma: Rf ((a, b)) = k(b − a)m+2 f (m+1) (ξ),

para todo f ∈ C m+1 ((a, b)), (3.45)

donde m es el orden de exactitud de la f´ormula.

3.3.

Cuadratura Gaussiana

Consideremos la siguiente integral Z w(x)f (x) dx

(3.46)

I

donde I = [a, b] es un intervalo y w(x) ≥ 0 para todo x ∈ I. La funci´on w : R → R es llamada funci´on de peso. Los m´etodos de integraci´on num´erica de la secci´on anterior motivan una f´ormula general para aproximar (3.46) de la siguiente forma Q(f ) = w1 f (x1 ) + . . . + xn f (xn )

(3.47)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

92

denominada f´ ormula de cuadratura de ” n ” puntos. Los valores wi (i = 1, . . . , n) son llamados pesos y los puntos xi i = 1, . . . , n son llamados nodos. Para aproximar el valor de (3.46) usando (3.47) es necesario calcular los pesos wi y los nodos xi para todo i = 1, . . . , n. Es decir, tenemos 2n variables. Una forma de calcular estas ” 2n ” variables es exigir que (3.47) aproxime (3.46) de forma exacta para un cierto tipo de funciones. En particular, por simplicidad, podemos exigir que (3.47) sea exacta para todo polinomio P2n−1 de grado menor o igual a ” 2n − 1 ”, es decir: Z (3.48) Q(P ) = w(x)P (x) dx para todo P ∈ P2n−1 I

Con este objetivo, basta que (3.48) sea exacta para la base {1, x, x2 , . . . , x2n−1 }, obteni´endose el siguiente sistema:

              

1 x1 .. . .. . xi1 .. . .. .

1 x2

xi2

x2n−1 x2n−1 1 2

··· 1 · · · xn .. .. . . .. .. . . · · · xin .. .. . . .. .. . . 2n−1 · · · xn

Z



              

w1 w2 .. . wi .. . wn

1 dx

  Z I   w(x)x dx     I   ..   . = Z     w(x)xi dx   I  ..  .  Z  w(x)x2n−1 dx 

                

(3.49)

I

del Sistema (3.49) se observan las siguientes complicaciones: El sistema (3.49) es no lineal. La matriz asociada al Sistema 3.49 no es cuadrada. Los nodos xi (i = 1, . . . , n) deben ser distintos. Una forma de encontrar la soluci´on ser´ıa aplicar el m´etodo de Newton para encontrar una soluci´on del Sistema 3.49, sin embargo, encontrar el punto inicial no siempre es f´acil. Como casos particulares, veamos los siguientes ejemplos: Ejemplo 3.3.1. Determine una f´ormula cuadratura de 2 puntos Q(f ) = A1 f (x1 ) + A2 f (x2 ) para aproximar Z 1 f (x)dx −1

3.3. CUADRATURA GAUSSIANA

93

Resoluci´ on: En Q(f ) tenemos 4 variables {x1 , x2 , A1 , A2 } a determinar. Para esto basta imponer que Q(f ) sea exacta para polinomios P3 de grado menor igual a 3, es decir: Z 1 f (x)dx = Q(f ) para todo P3 ∈ P3 . −1

Con este fin, basta que sea exacta para la base {1, x, x2 , x3 }, luego: Z 1 1 dx = 2 = A1 + A2 (3.50) f (x) = 1, ⇒ −1 Z 1 x dx = 0 = A1 x1 + A2 x2 (3.51) f (x) = x, ⇒ −1 Z 1 2 2 f (x) = x ⇒ xdx == = A1 x21 + A2 x22 (3.52) 3 −1 Z 1 3 f (x) = x ⇒ xdx = 0 = A1 x31 + A2 x32 (3.53) −1

De (3.51) en (3.53) resulta: A2 x2 (x2 +x1 )(x2 − x1 ) = 0.

(3.54)

Para que Q(f ) sea una f´ormula de cuadratura de dos puntos, entonces A1 6= 0, A2 6= 0 y x1 6= x2 . Por tanto, de (3.54) se tiene: x2 (x2 + x1 ) = 0 lo cual implica los siguientes dos casos: x2 = 0

∨

x2 = −x1 .

I Si x2 = 0, entonces de (3.51) se obtiene x1 = 0 y esto contradice el hecho de obtener una f´ormula de cuadratura de 2 puntos. I Si x2 = −x1 , esto implica en (3.52): 2 x21 (A1 + A2 ) = , 3 y usando (3.50) se concluye: √ √ 3 3 x1 = ∨ x1 = − 3 3 En cualquier caso, se obtiene de (3.50) y (3.51): A1 = A2 = 1.

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

94

√

Por lo tanto , en cualquiera de los casos Q(f ) = f (− √ √ Z 1 3 3 ) + f( ) f (x)dx ≈ f (− 3 3 −1

√ 3 3 ) + f ( ), 3 3

y as´ı:

Una aplicaci´on sencilla de Q(f ) ser´a aproximar el valor de la integral cuando: f (x) = cos (x). √ √ Z 1 3 3 cos (x)dx = 1.68294 y Q(f ) = cos (− ) + cos ( ) = 1.67582. 3 3 −1

Ejemplo 3.3.2. Sea f una funci´on continua. Use el m´etodo de cuadratura sobre I = [−1, 1] y w(x) = 1 para todo x ∈ I, para obtener una f´ormula de cuadratura de 3 puntos. Resoluci´ on: La f´ormula de cuadratura de 3 puntos tiene la forma: Z 1 f (x) dx = w1 f (x1 ) + w2 f (x2 ) + w3 f (x3 ) −1

Hacemos exacta la f´ormula para la base de los polinomios de grado menor o igual a 5, es decir:  β = 1, x1 , x2 , x3 , x4 , x5 Obtenemos el siguiente sistema de ecuaciones: f (x) = 1, ⇒ w1 + w2 + w3 = 2 f (x) = x, ⇒ w1 x1 + w2 x2 + w3 x3 = 0 2 f (x) = x2 ⇒ w1 x21 + w2 x22 + w3 x23 = 3 f (x) = x3 ⇒ w1 x31 + w2 x32 + w3 x33 = 0 2 f (x) = x4 , ⇒ w1 x41 + w2 x42 + w3 x43 = 5 5 5 5 5 f (x) = x , ⇒ w1 x1 + w2 x2 + w3 x3 = 0

(3.55) (3.56) (3.57) (3.58) (3.59) (3.60)

Multiplicamos la ecuaci´on (3.56) por x21 , para luego restarla de (3.58), resultando: w2 x32 + w3 x33 − w2 x2 x21 − w3 x3 x21 = 0. (3.61) Multiplicamos la ecuaci´on (3.58) por x21 , para luego restarla de (3.60), resultando: w2 x52 + w3 x53 − w2 x32 x21 − w3 x33 x21 = 0. (3.62)

3.3. CUADRATURA GAUSSIANA

95

Multiplicamos la ecuaci´on (3.61) por x22 , para luego restarla de (3.62) obtenemos: w3 x3 (x3 − x1 )(x3 + x1 )(x2 − x3 )(x2 + x3 ) = 0. (3.63) Desde que queremos una f´ormula de tres puntos, entonces w3 6= 0, x3 6= x1 y x3 6= x2 , luego (3.63) se reduce a: x3 (x3 + x1 )(x3 + x2 ) = 0.

(3.64)

Por tanto, de (3.64) resta analizar los siguientes casos: x3 = 0

∨

x3 = −x1

∨

x3 = −x2 .

Veamos cada caso: Si x3 = 0, entonces x1 6= 0 y x2 6= 0. Adem´as, el sistema inicial (3.55)(3.60) se reduce a: w1 + w2 + w3 = 2 w 1 x1 + w 2 x2 = 0 2 w1 x21 + w2 x22 = 3 3 3 w 1 x1 + w 2 x2 = 0 2 w1 x41 + w2 x42 = 5 5 5 w 1 x1 + w 2 x2 = 0

(3.65) (3.66) (3.67) (3.68) (3.69) (3.70)

Multiplicando (3.66) por x22 y restando (3.68) se obtiene: w1 x1 (x2 − x1 )(x2 + x1 ) = 0 Desde que buscamos f´ormula de 3 puntos, entonces: w1 6= 0 y x1 6= x2 , por tanto se obtiene x1 = −x2 (3.71) Reemplazando en (3.66): w1 = w2 . Considerando (3.67) y (3.69) se obtiene: x2 = −x1 = Que reemplazando en (3.67): w2 = w1 =

q

3 5

√

=

15 5

5 9

Reemplazando en (3.65) se obtiene: w3 = 89 . De este modo, la f´ormula de cuadratura de 3 puntos queda como:  √  √  Z 1 5 15 5 15 8 f (x) dx = f − + f + f (0). (3.72) 9 5 9 5 9 −1

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

96 En los otros casos:

x3 = −x1

∨

x3 = −x2 ,

se obtiene la misma regla de cuadratura dada por (3.72).

R1 Ejemplo 3.3.3. Para aproximarla integral I(f ) = 0 xf (x)dx, tal que f ∈ C 2 ([0.1]), consideremos la f´ormula de cuadratura Q(f ) = Af (0) + Bf 0 (1) donde f 0 representa la derivada de f . Determine los pesos A y B de modo que la f´ormula sea exacta cuando f sea un polinomio de grado menor o igual a 1. Resoluci´ on: Debemos calcular A y B tal que I(f ) = Q(f ) donde f es un polinomio de grado menor o igual a 1. Para f (x) = 1 se tiene f 0 (x) = 0, entonces de I(f ) = Q(f ) resulta: 1

Z

xdx = A(1) + B(0), 0

1 esto da A = . 2 Para f (x) = x se tiene f 0 (x) = 1, entonces de I(f ) = Q(f ) resulta: Z

1

x2 dx = A(0) + B(1),

0

1 esto da B = . 3 Por tanto: 1 1 Q(f ) = f (0) + f 0 (1). 2 3

Una de las principales desventajas del m´etodo de cuadratura es tener que resolver el sistema no lineal dado por (3.49). Entonces, una forma de reducir el c´alculo es poder elegir los nodos de una forma adecuado de modo que el grado de precisi´on de la f´ormula de cuadratura no sea afectada.

3.4. CUADRATURA GAUSSIANA

3.3.1.

97

Polinomios Ortogonales

1. Polinomios de Legendre: I = [−1, 1], w(x) = 1, L0 (x) = 1, L1 (x) = x, n 2n + 1 xLn (x) − Ln−1 (x) ∀n ≥ 1. Ln+1 (x) = n+1 n+1 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 = h−1, 1i, 1 w(x) = √ , 1 − x2 T0 (x) = 1, T1 (x) = x, Tn+1 (x) = 2xTn (x) − Tn−1 (x) ∀n ≥ 1.

3.4.

Cuadratura Gaussiana

Queremos aproximar la siguiente integral: Z w(x)f (x) dx I

(3.73)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

98

mediante una f´ormula de cuadratura de ” n ” puntos de la forma: Q(f ) = w1 f (x1 ) + . . . + wn f (xn )

(3.74)

Estos m´etodos consideran conocidos los nodos {xi }ni=1 , donde ´estos son las ra´ıces de un polinomio ortogonal Pn de grado n. Por tanto, resta encontrar los n pesos w1 , w2 , . . . , wn . Para esto, exigimos que (3.74) calcule el valor exacto de (3.73) para la base {1, x, x2 , . . . , xn−1 }, obteni´endose el siguiente sistema: Z

 

1 1  x1 x2  .  .  .  i xi2  x1  .  .. xn−1 xn−1 1 2

··· ··· .. .

1 xn .. .

··· .. .

xin .. .

        

· · · xnn−1

w1 w2 .. . wi .. . wn

w(x) dx

  ZI   w(x)x dx     I   ..   . = Z     w(x)xi dx   I  ..  .  Z  w(x)xn−1 dx 

                

(3.75)

I

del Sistema (3.75) se observan las siguientes simplificaciones en relaci´on al Sistema (3.49): El sistema (3.75) es lineal. La matriz asociada al Sistema 3.75 es cuadrada. Desde que los nodos xi (i = 1, . . . , n) son distintos por ser ra´ıces de polinomios ortogonales, la soluci´on del Sistema (3.75) es u ´nica. Ejemplo 3.4.1. Aproximar la siguiente integral: Z 3 x e Sen(x) dx 1 + x2 0 empleando el m´etodo de Gauss-Legendre de 3 nodos, es decir, los nodos son las ra´ıces del polinomio ortogonal de Legendere donde ´estos son dados por la f´ormula de Rodr´ıgues mostrado a continuaci´on: Pn (x) =

1 dn ((x2 − 1)n ) 2n n! dxn

3.4. CUADRATURA GAUSSIANA

99

Resoluci´ on: A partir de la f´ormula de Rodr´ıgues se calcula el polinomio de Legendre de grado 3 para luego obtener los 3 nodos de la regla de cuadratura: 1 1 P3 (x) = (5x3 − 3x) = x(5x2 − 3) 2 2 De donde las ra´ıces de P3 son: r r 3 3 = −0.7746, x2 = 0, x3 = = 0.7746. x1 = − 5 5 Para obtener los pesos w1 , w2 , w3 usamos el Sistema (3.75) para n = 3, dado por:      1 1 1 w1 2  x1 x2 x3   w 2  =  0  2 x21 x22 x23 w3 3 Resolviendo el sistema se obtiene: w1 = 0.5556,

w2 = 0.8889,

w3 = 0.5556.

Por tanto, la regla de cuadratura de Gauss-Legendre de 3 puntos tiene la forma: Z 1 f (x)dx = w1 f (x1 ) + w2 f (x2 ) + w3 f (x3 ) (3.76) −1

Luego, para un intervalo arbitrario [a, b] se tiene:   Z b Z b−a 1 b−a a+b f (x)dx = f t+ dt. 2 2 2 a −1

(3.77)

x

Sen(x) Por tanto, usando f (x) = e 1+x en (3.77) se obtiene: 2   Z 3 x Z e Sen(x) 3 1 3 3 dx = f t+ dt 1 + x2 2 −1 2 2 0        3 3 3 3 3 3 3 f x1 + +f x2 + +f x3 + = 2 2 2 2 2 2 2

= 2.8633.  Por tanto, por construcci´on podemos concluir que la f´ormula de cuadratura gaussiana (3.74) es exacta para cualquier polinomio Sn−1 de grado ” n − 1 ”, es decir: Z w(x)Sn−1 (x) dx = Q(Sn−1 ). I

(3.78)

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

100

Con la elecci´on hecha para los nodos, es natural preguntarse: ¿ Q(f ) seguir´a siendo exacta para los polinomios P˜2n−1 de grado 2n − 1?. Al dividir P˜2n−1 entre Pn , del algoritmo de la divisi´on resulta: P˜2n−1 (x) = qn−1 (x)Pn (x) + rn−1 (x)

(3.79)

Recordando que los xi son ra´ıces de Pn para todo i = 1, ..., n, entonces de (3.79) se obtiene P˜2n−1 (xi ) = qn−1 (xi )Pn (xi ) + rn−1 (xi ) = rn−1 (xi ).

(3.80)

Usando (3.79) y (3.80) se obtiene lo siguiente: Z b Z b ˜ w(x)P2n−1 (x)dx = w(x)[qn−1 (x)Pn (x) + rn−1 (x)]dx a

a b

Z

w(x)P˜2n−1 (x)dx =

a

b

Z w(x)qn−1 (x)Pn (x) +

a b

Z

Z

w(x)P˜2n−1 (x)dx =

Z

b

w(x)rn−1 (x)dx a

b

w(x)rn−1 (x)dx = Q(rn−1 ) a

a b

Z

w(x)P˜2n−1 (x)dx = w1 rn−1 (x1 ) + w2 rn−1 (x2 ) + ... + wn rn−1 (xn )

a

Z

b

w(x)P˜ 2n−1 (x)dx = w1 P˜2n−1 (x1 ) + w2 P˜2n−1 (x2 ) + ... + wn P˜2n−1 (xn )

a

Z

b

w(x)P˜ 2n−1 (x)dx = Q(P˜2n−1 ).

a

Esto u ´ltimo indica que la f´ormula de cuadratura de n puntos Q(f ) dada por (3.74) sigue siendo exacta para cualquier polinomio de grado 2n − 1 cuando los nodos {xi }ni=1 son las ra´ıces un polinomio ortogonal de grado n. Definici´ on 3.1. Una regla de cuadratura tiene grado de precisi´on ”p”si es exacta para polinomios de grado menor o igual a p. Ejemplo 3.4.2. Considere la f´ormula de cuadratura de dos puntos dada por: √ ! √ ! 3 3 +f . Q(f ) = f − 3 3 Determine su grado de precisi´on. Resoluci´ on:

3.4. CUADRATURA GAUSSIANA

Z

1

−1

1dx = x|1−1 = 2;

1 x2 xdx = = 0; 2 −1 −1 1 Z 1 x3 2 2 x dx = = ; 3 −1 3 −1 1 Z 1 x4 3 x dx = = 0; 4 −1 −1 1 Z 1 2 x5 4 = ; x dx = 5 −1 5 −1 Z

101

1

Q(1) = 1 + 1 = 2 √ √ 3 3 + =0 Q(x) = − 3 3 1 1 2 + = 3 3 3 ! √ 3 3 3 + Q(x ) = − 3 √ !4 3 Q(x4 ) = − + 3 Q(x2 ) =

√ !3 3 =0 3 √ !4 3 2 = 3 9

Por tanto, Q(f ) tiene grado de precisi´on igual a 3.



Ejemplo 3.4.3. Determine el grado de precisi´on de la siguiente f´ormula de cuadratura: Q(f ) =

7 16 7 1 1 f (−1) + f (0) + f (1) + f 0 (−1) − f 0 (1), 15 15 15 15 15

sobre el intervalo [−1, 1]. Resoluci´ on: f (x) = 1,

0

Z

1

f (x) = 0,

1dx = 2;

Q(1) = 2

xdx = 0;

Q(x) = 0

−1

f (x) = x,

0

Z

1

f (x) = 1, −1

Z

1

2 2 x2 dx = ; Q(x2 ) = 3 3 −1 Z 1 3 0 2 f (x) = x , f (x) = 3x , x3 dx = 0; Q(x3 ) = 0 2

0

f (x) = x , f (x) = 2x,

−1

Z

1

2 2 x4 dx = ; Q(x4 ) = 5 5 −1 Z 1 f (x) = x5 , f 0 (x) = 5x4 , x5 dx = 0; Q(x5 ) = 0 f (x) = x4 , f 0 (x) = 4x3 ,

−1

f (x) = x6 , f 0 (x) = 6x5 ,

Z

1

2 2 x5 dx = ; Q(x6 ) = 7 15 −1

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

102

Por tanto, Q(f ) tiene grado de precisi´on igual a 5.



Ejemplo 3.4.4. Si Q(f ) = w1 f (x1 ) + ... + wn f (xn ) es una f´ormula de cuadratura donde los xi son ra´ıces de alg´ un polinomio ortogonal Pn , entonces Q(f ) tiene grado de precisi´on 2n − 1. 

3.5.

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

4

f (x)dx.

determine un valor aproximado de la integral 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

1

exp(−x) dx x

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).

3.5. EJERCICIOS PROPUESTOS 5. Encuentre la f´ormula Z

103

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(

πx ) 2

104

´ NUMERICA ´ CAP´ITULO 3. INTEGRACION

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

b1 b2 .. . bn−1 bn

      

que escrito en forma matricial: Ax = b

(4.1)

donde A ∈ Rn×n , x ∈ Rn y b ∈ Rn . Definici´ on 4.1. Dado el sistema lineal Ax = b, definimos la matriz aumentada M asociada al sistema lineal de la forma siguiente: M = (A | |b)

4.1.

(4.2)

M´ etodos Directos

Definici´ on 4.2. Dada una matriz A ∈ Rn×n , definimos como operaciones elementales fila para la matriz A a cualquiera de las siguientes operaciones: 1. Intercambiar la fila i con la fila j, denotado por Fi ↔ Fj . 2. Asignar a la fila i la misma fila i pero multiplicada por un n´ umero no nulo λ. Esto es denotado por Fi ← λFi . 105

106

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

3. Asignar a la fila i la misma fila i y sum´andole λ veces la fila j donde λ 6= 0. Esto es denotado por Fi ← Fi + λFj . 

−1 2 −3 1  2 −1 2 5 Ejemplo 4.1.1. Considere la matriz A =   1 −3 1 −1 2 1 −1 2 la matriz original A realice las siguientes operaciones elementales: F1 ← −4F1 y F4 ← F4 − 2F2 . Resoluci´ on

  . Sobre  F2 ↔ F3 ,

Aplicando la operaci´on elemental F2 ↔ F3 a la matriz A se obtiene:   −1 2 −3 1  1 −3 1 −1     2 −1 2 5  2 1 −1 2 Aplicando la operaci´on elemental F1 ← −4F1 a la matriz A se obtiene:   4 −8 12 −4  2 −1 2 5     1 −3 1 −1  2 1 −1 2 Aplicando la operaci´on elemental F4 ← F4 − 2F2 a la matriz A se obtiene:   −1 2 −3 1  2 −1 2 5     1 −3 1 −1  −2 3 −5 −8 

Definici´ on 4.3. Una matriz E ∈ Rn×n es llamada matriz elemental cuando ella es obtenida al aplicar una operaci´on elemental fila a la matriz identidad In . Ejemplo 4.1.2. Considere la matriz identidad I4 . Son matrices elementales:

´ 4.1. METODOS DIRECTOS

107

Operaci´on elemental F1 ↔ F3 . Entonces: 

1  0 E=  1 0

0 1 0 0

1 0 0 0

 0 0   0  1

Operaci´on elemental F3 ← 2F3 . Entonces: 

1  0 E=  0 0

0 2 0 0

0 0 1 0

 0 0   0  1

Operaci´on elemental F2 ← F2 − 3F4 . Entonces: 

1  0 E=  0 0

0 1 0 0

 0 0 0 −3   1 0  0 1 

Definici´ on 4.4. Una matriz L es llamada Triangular Inferior cuando tiene la forma: 

l11 l21

0 l22

0 0

···

0 0

 0 0      0 

   .. L= .   ln−1,1 ln−1,2 ln−1,3 · · · ln−1,n−1 ln,1 ln,2 ln,3 · · · ln,n−1 ln,n Cuando se tiene un sistema lineal de la forma: Lx = b

(4.3)

donde L es una matriz triangular inferior cuyos elementos lii 6= 0 (i = 1, . . . , n). Para calcular la soluci´on x se usa el Algoritmo 5 descrito a continuaci´on.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

108

Algoritmo 5: Sustituci´on Progresiva Entrada: Ingresar una matriz triangular inferior L ∈ Rn×n . 1 inicio 2 para i ← 1 a n hacer i−1 X bi − lij xj lii fin para devolver Soluci´on del sistema lineal x = (x1 , . . . , xn ).

4 5 6

j=1

xi ←

3

fin

Definici´ on 4.5. Una matriz U es llamada Triangular Superior cuando tiene la forma: 

u11 u12 · · · u1,n−2 u1,n−1 u1,n  0 u22 u2,n−2 u2,n−1 u2,n   . .. U =   0 0 0 un−1,n−1 un−1,n 0 0 ··· 0 0 un,n

      

Cuando se tiene un sistema lineal de la forma: Ux = b

(4.4)

donde U es una matriz triangular superior cuyos elementos uii 6= 0 (i = 1, . . . , n). Para calcular la soluci´on x se usa el Algoritmo 6 descrito a continuaci´on. Algoritmo 6: Sustituci´on Regresiva Entrada: Ingresar una matriz triangular superior U ∈ Rn×n . 1 inicio 2 para i ← n a 1 hacer n X bi − uij xj xi ←

3

uii fin para devolver Soluci´on del sistema lineal x = (x1 , . . . , xn ).

4 5 6

j=i+1

fin

´ 4.1. METODOS DIRECTOS

4.1.1.

109

Eliminaci´ on Gaussiana

Dado el sistema lineal Ax = b, el m´etodo consiste en aplicar operaciones elementales fila a la matriz aumentada M asociada al sistema lineal de forma tal que la matriz A sea transformada a una matriz triangular superior. Ejemplo 4.1.3. Resuelva gaussiana.  −1 3  1 −1   3 1 3 −1

el sistema lineal siguiente mediante eliminaci´on  −1 1 x1  x2 2 −1   −2 2   x3 1 2 x4

Resoluci´ on: Consideremos la matriz  −1  1 M = (A | b) =   3 3





 −8   4  =    4 . 16

aumentada  3 −1 1 −8 −1 2 −1 4   1 −2 2 4  −1 1 2 16

Sobre M realizamos las siguientes operaciones elementales: F2 ← F1 + (1)F1 ,

F3 ← F3 + (3)F1 ,

F4 ← F4 + (3)F1 ;

obteni´endose la siguiente matriz aumentada:  −1 3 −1 1 −8  0 2 1 0 −4 M =  0 10 −5 5 −20 0 8 −2 5 −8

   

Sobre el M obtenido realizamos las siguientes operaciones elementales:     10 8 F3 ← F3 + − F2 , F4 ← F4 + − F2 2 2 obteni´endose la siguiente matriz aumentada:   −1 3 −1 1 −8  0 2 1 0 −4   M =  0 0 −10 5 0  0 0 −6 5 8 Sobre el M obtenido realizamos la siguiente operaci´on elemental:   6 F4 ← F4 + − F3 10

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

110 obteni´endose:



−1  0 M =  0 0

 3 −1 1 −8 2 1 0 −4   0 −10 5 0  0 0 2 8

Por tanto, el sistema resultante es:



−1  0   0 0

 3 −1 1 x1   2 1 0   x2 0 −10 5   x3 0 0 2 x4

 −8   −4  =    0  8 



Este sistema tiene la forma del Sistema 4.4 y as´ı podemos aplicar el Algoritmo 6. Esto es:

i = 4 : x4 =

b4 u44 b3 −

u33 b2 −

4 X

u22 b1 −

4 X j=2

u11

=

3 − (5(4)) =2 −10

=

−4 − (1(2) + 0(4)) = −3 2

=

−8 − (3(−3) − 1(2) + 1(4)) =1 −1

u2j xj

j=3

i = 2 : x2 =

8 =4 2

u3j xj

j=4

i = 3 : x3 =

i = 1 : x1 =

4 X

=

u1j xj

As´ı se obtiene la soluci´on x = (1, −3, 2, 4). El proceso para triangulizar una matriz es descrito a continuaci´on:



´ 4.1. METODOS DIRECTOS

111

Algoritmo 7: Proceso de Triangulizaci´on de una matriz Entrada: Ingresar una matriz A ∈ Rn×n . 1 inicio 2 f lag ← 1; 3 para j ← 1 a n hacer 4 si Ajj 6= 0 entonces 5 para i ← j + 1 a n hacer 6 para k ← 1 a n hacer   Ajk 7 Aik ← Aik − Aij Ajj 8 fin para 9 fin para 10 en otro caso 11 Guardar j 12 f lag ← 0 13 Parar. 14 fin si 15 fin para 16 si f lag == 1 entonces 17 devolver Matriz A en forma triangular superior. 18 en otro caso 19 Falla el proceso de triangulizaci´on en la fila j 20 fin si 21 fin Ejemplo 4.1.4. ¿Es posible el siguiente sistema lineal?  2 2 −1  2 2 3   1 −1 1   −1 8 −2 3 −2 1

usar el m´etodo de eliminaci´on Gaussiana para    3 −1 x1    −4 1    x2      −2 1    x3  =  3 −1   x4   −3 2 x5

Resoluci´ on Consideremos la matriz aumentada  2  2  M = (A | b) =   1  −1 3

−9 11 5 0 7

     

 2 −1 3 −1 −9 2 3 −4 1 11   −1 1 −2 1 5   8 −2 3 −1 0  −2 1 −3 2 7

112

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

Usamos el Algoritmo 7 realizando las siguientes operaciones elementales sobre la primera columna:       1 1 3 F2 ← F1 +(−1)F1 , F3 ← F3 + − F1 , F4 ← F4 + F1 , F5 ← F5 + − F1 ; 2 2 2 obteni´endose la siguiente matriz aumentada:   2 2 −1 3 −1 −9    0 0  4 −7 2 20      3 7 3 19   0 −2  −  2 2 2 2  M =    9 3 9  5   0 9 − − −   2 2 2 2     5 15 7 41  0 −5 − 2 2 2 2 Podemos observar que el m´etodo de eliminaci´on gaussiana no puede continuar dado que m22 = 0.  El ejemplo anterior muestra una clara desventaja del m´etodo de eliminaci´on gausiana para calcular de forma directa la soluci´on del sistema lineal. Por tanto, una condici´on necesaria de este m´etodo es que se verifique en cada paso que los elementos de la diagonal sean distintos de cero. El ejemplo siguiente muestra que esta condici´on no es suficiente para garantizar el c´alculo de la soluci´on exacta. Ejemplo 4.1.5. Considere el sistema lineal siguiente: 10−4 x1 + x2 = 1 x1 + x2 = 3 cuya soluci´on exacta es x = 2.00020002... e y = 0.99979997.... Halle la soluci´on del sistema en un computador donde la aritm´etica en punto flotante usa 3 d´ıgitos en la mantisa y redondeo. Resoluci´ on: La matriz aumentada del sistema viene dado por:   0.100 × 10−3 0.100 × 101 0.100 × 101 M= 0.100 × 101 0.100 × 101 0.300 × 101 Realizamos la operaci´on elemental:  F2 ← F2 −

1 0.1 × 10−3

 F1

´ 4.1. METODOS DIRECTOS

113

es decir: m21 m22

 0.1 × 10−3 = 0.1 × 10 + − 0.1 × 10−3   0.1 × 101 1 = 0.1 × 10 + − 0.1 × 10−3 1 = 0.1 × 10 − 0.1 × 105 1



= 0.1 × 101 − 1 = 0 = 0.1 × 101 − 104 (expresando en punto flotante)

= 0.00001 × 105 − 0.1 × 105

(igualando exponentes)

= (0.00001 − 0.1) × 105

(restando mantisas)

= = −0.09999 × 105 = −0.9999 × 104 (expresando en punto flotante) 4 = −1.000 × 10

m23

(redondeo al tercer d´ıgito)

0.1 × 101 = 0.3 × 101 + − 0.1 × 10−3 = 0.3 × 101 − 0.1 × 105



= 0.3 × 101 − 104 (expresando en punto flotante)

= 0.00003 × 105 − 0.1 × 105

(igualando exponentes)

= (0.00003 − 0.1) × 105

(restando mantisas)

= = −0.09997 × 105 = −0.9997 × 104 (expresando en punto flotante) = −1.000 × 104

(redondeo al tercer d´ıgito)

Por lo que la matriz aumentada M queda de la forma siguiente:  −4  10 1 1 M= 0 −104 −104 Por tanto, el sistema resultante es:  −4     10 1 x1 1 = 0 −104 x2 −104 aplicamos ahora el Algoritmo 6 al sistema anterior, resulta: i = 2 : x2 =

b2 u22 b1 −

i = 1 : x1 = obteniendo la soluci´on x = (0, 1).

2 X j=2

u11

=

−104 =1 −104

=

1−1 =0 10−4

u2j xj



114

4.1.2.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

Pivoteo Parcial

Los inconvenientes observados en el Ejemplo 4.1.4 y Ejemplo 4.1.5 pueden ser superados si consideramos la siguiente variaci´on al m´etodo de eliminaci´on gaussiana descrito a continuaci´on. Algoritmo 8: Proceso de pivoteo parcial Entrada: Ingresar una matriz A ∈ Rn×n . 1 inicio 2 para j ← 1 a n hacer 3 maxc ← |Ajj |; 4 para i ← j + 1 a n hacer 5 si |Aij | > maxc entonces 6 maxc ← |Aij |; 7 p←i 8 fin para 9 Intercambiar las filas j y p; 10 para i ← j + 1 a n hacer 11 Haciendo ceros los elementos de cada fila i en la columna j para k ← 1 a n hacer   Ajk 12 Aik ← Aik − Aij Ajj 13 fin para 14 fin para 15 fin para 16 devolver Matriz triangular superior U y vector b 17 fin Ejemplo 4.1.6. Consideremos el sistema del Ejemploo 4.1.4 cuya matriz aumentada es dada por:    M = (A | b) =   

 2 2 −1 3 −1 −9 2 2 3 −4 1 11   1 −1 1 −2 1 5   −1 8 −2 3 −1 0  3 −2 1 −3 2 7

El m´aximo en valor absoluto de la primera columna de A1 es m51 = 3. Por tanto, realizamos la operaci´on elemental: F1 ↔ F5

´ 4.1. METODOS DIRECTOS

115

obteniendo:    M =  

7 3 −2 1 −3 2 2 2 3 −4 1 11 1 −1 1 −2 1 5 −1 8 −2 3 −1 0 2 2 −1 3 −1 −9

     

Ahora hacemos cero los elementos mi1 (i = 2, 3, 4, 5) mediante las operaciones elementales:         1 1 2 2 F1 , F3 ← F3 + − F1 , F4 ← F4 + F1 , F5 ← F5 + − F1 F2 ← F2 + − 3 3 3 3 resultando: 

3 −2

1

−3

2

7

  1 19 7  0 10 −2 −  3 3 3 3    1 2 1 8  M =  0 − 3 3 −1 3 3   5 1 7 22   0 − 2 −  3 3 3 3   10 41 5 7 0 − − 5 − 3 3 3 3

                

22 El m´aximo en valor absoluto de la primera columna de A2 es m42 = . Por 3 tanto, realizamos la operaci´on elemental: F2 ↔ F4 obteniendo: 

3 −2

1

−3

2

7

  7  0 22 − 5 2 − 1  3 3 3 3    1 2 1 8  M =  0 − 3 3 −1 3 3   10 7 1 19   0 −2 −  3 3 3 3   10 5 7 41 0 − 5 − − 3 3 3 3

                

116

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

Ahora hacemos cero los elementos mi2 (i = 3, 4, 5) mediante las operaciones elementales:       1 10 10 F2 , F4 ← F4 + − F2 , F5 ← F5 + − F2 , F3 ← F3 + 22 22 22 resultando: 

3 −2

1

−3

2

7

  1 7  0 22 − 5 2 −  3 3 3 3    10 7 61 13  − M = 0 0 22 11 22 22   34 58 32 2   0 0 − −  11 11 11 11   162 10 45 24 − 0 0 − − 11 11 11 11

                

34 El m´aximo en valor absoluto de la primera columna de A3 es m43 = . Por 11 tanto, realizamos la operaci´on elemental: F3 ↔ F4 obteniendo: 

3 −2

1

−3

2

7

  1 7  0 22 − 5 2 −  3 3 3 3    34 32 2 58  − − M = 0 0 11 11 11 11   13 10 7 61   0 0 −  22 11 22 22   10 45 24 162 0 0 − − − 11 11 11 11

                

Ahora hacemos cero los elementos mi3 (i = 4, 5) mediante las operaciones elementales:     13 10 F4 ← F4 + − F3 , F5 ← F5 + F3 , 68 34

´ 4.1. METODOS DIRECTOS

117

resultando: 

3 −2

1

−3

2

7

  1 7  0 22 − 5 2 −  3 3 3 3    34 58 32 2  M =  0 0 11 − 11 − 11 11   6 6 30   0 0 0 −  17 17 17   55 38 224 0 0 0 − − 17 17 17

                

55 El m´aximo en valor absoluto de la primera columna de A4 es m54 = . Por 17 tanto, realizamos la operaci´on elemental:

F4 ↔ F5

obteniendo: 

3 −2

1

−3

2

7



   1 7   0 22 − 5  2 −   3 3 3 3      34 32 2 58   0 0  − − M =  11 11 11 11     55 38 224    0 0  − − 0  17 17 17     6 6 30  0 0 0 − 17 17 17 Ahora hacemos cero los elementos mi4 (i = 5) mediante las operaciones elementales:   6 F5 ← F5 + F4 , 55

118 resultando:

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES



3 −2

1

−3

2

7



   1 7   0 22 − 5  2 −   3 3 3 3      32 2 58  34  0 0  − − M =  11 11 11 11     55 38 224    0 0  0 − −  17 17 17     6 18  0 0 0 0 55 55 aplicamos ahora el Algoritmo 6 al sistema anterior, resulta: i = 5 : x5 =

b5 u55 b4 −

5 X

u4j xj

j=5

i = 4 : x4 =

= −2

u44 b3 −

5 X

u3j xj

j=4

i = 3 : x3 =

=0

u33 b2 −

5 X

u2j xj

j=3

i = 2 : x2 =

=1

u22 b1 −

i = 1 : x1 =

=3

5 X

u1j xj

j=2

u11

= −1

obteniendo la soluci´on x = (−1, 1, 0, −2, 3).



Ejemplo 4.1.7. Consideremos el sistema del Ejemploo 4.1.5 cuya matriz aumentada es dada por:   0.100 × 10−3 0.100 × 101 0.100 × 101 M= 0.100 × 101 0.100 × 101 0.200 × 101 Realizamos pivoteo parcial sobre A1 mediante la operaci´on elemental: F1 ↔ F2

´ 4.1. METODOS DIRECTOS

119



0.100 × 101 0.100 × 101 0.300 × 101 M= 0.100 × 10−3 0.100 × 101 0.100 × 101 Realizamos la operaci´on elemental:
F2 ← F2 + −0.1 × 10−3 F1



es decir: m21 = 0.1 × 10−3 + (−0.1 × 10−3 ) (0.1 × 101 ) 0.1 × 10−3 − 0.01 × 10−2

= 0.1 × 10−3 − 0.1 × 10−3 = 0

m22 = 0.1 × 101 + (−0.1 × 10−3 ) (0.1 × 101 )

= 0.1 × 101 − 0.01 × 10−2

= 0.1 × 101 − 0.1 × 10−3

(expresando en punto flotante)

= 0.1 × 101 − 0.00001 × 101

(igualando exponentes)

= (0.10000 − 0.00001) × 101 = 0.09999 × 101 (restando mantisas) = 0.9999 × 100

(expresando en punto flotante)

= 1

(redondeo al tercer d´ıgito)

m23 = 0.1 × 101 + (−0.1 × 10−3 ) (0.3 × 101 )

= 0.1 × 101 − 0.03 × 10−2

= 0.1 × 101 − 0.3 × 10−3

(expresando en punto flotante)

= 0.1 × 101 − 0.00003 × 101

(igualando exponentes)

= (0.10000 − 0.00003) × 101 = 0.09997 × 101 (restando mantisas) = 0.9997 × 100

(expresando en punto flotante)

= 1

(redondeo al tercer d´ıgito)

Por lo que la matriz aumentada M queda de la forma siguiente:   1 1 3 M= 0 1 1 Por tanto, el sistema resultante es:      x1 3 1 1 = 0 1 x2 1 aplicamos ahora el Algoritmo 6 al sistema anterior, resulta: b2 1 i = 2 : x2 = = =1 u22 1 b1 − i = 1 : x1 =

2 X j=2

u11

u2j xj =

3−1 =2 1

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

120

obteniendo la soluci´on x = (2, 1).

4.1.3.



Pivoteo Total

Definici´ on 4.6. Dada una matriz A ∈ Rn×n , definimos como operaciones elementales columna para la matriz A a cualquiera de las siguientes operaciones: 1. Intercambiar la columna i con la columna j, denotado por Ci ↔ Cj . 2. Asignar a la columna i la misma columna i pero multiplicada por un n´ umero no nulo λ. Esto es denotado por Ci ← λCi . 3. Asignar a la columna i la misma columna i y sum´andole λ veces la columna j donde λ 6= 0. Esto es denotado por Ci ← Ci + λCj .   −1 2 −3 1  2 −1 2 5   Ejemplo 4.1.8. Considere la matriz A =   1 −3 1 −1 . Sobre 2 1 −1 2 la matriz original A realice las siguientes operaciones elementales: C1 ↔ C3 , C2 ← −5C2 y C4 ← C4 + 2C1 . Resoluci´ on Aplicando la operaci´on elemental C1 ↔ C3 a la matriz A se obtiene:   −3 2 −1 1  2 −1 2 5     1 −3 1 −1  −1 1 2 2 Aplicando la operaci´on elemental C2 ← −5C2 a la matriz A se obtiene:   −3 −10 −1 1  2 5 2 5     1 15 1 −1  −1 −5 2 2 Aplicando la operaci´on elemental C4 ← C4 + 2C1 a la matriz A se obtiene:   −1 2 −3 −1  2 −1 2 9     1 −3 1 1  −2 1 −1 6

´ 4.1. METODOS DIRECTOS

121 

Ejemplo 4.1.9. Consideremos el sistema del Ejemploo 4.1.4 cuya matriz aumentada es dada por:   2 2 −1 3 −1 −9  2 2 3 −4 1 11     5  M = (A | b) =  1 −1 1 −2 1   −1 8 −2 3 −1 0  7 3 −2 1 −3 2 Primero, denotemos por Ind el vector de ´ındices de las variables xi (i = 1, 2, 3, 4, 5), es decir: Ind = ( 1 2 3 4 5 ). Tener en cuenta, cuando realizamos una operaci´on elemental columna, entonces cambia el orden de los elementos del vector Ind. El m´aximo elemento de A1 en valor absoluto es dado por m42 = 8. Por tanto, realizamos las operaciones elementales: F1 ↔ F4 ,

C1 ↔ C2 ,

luego: Ind = ( 2  8 −1  2 2  −1 1 M =   2 2 −2 3

1 3 4 5 ), −2 3 −1 0 3 −4 1 11 5 1 −2 1 −1 3 −1 −9 1 −3 2 7

     

Ahora hacemos cero los elementos mi1 (i = 2, 3, 4, 5) mediante las operaciones elementales:   2 F2 ← F2 + − F1 , 8   1 F3 ← F3 + F1 , 8   2 F4 ← F4 + − F1 , 8   2 F5 ← F5 + F1 , 8

122

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

resultando: 

8 −1 −2

3

−1

0



    7 19 5   0 9 − 11   4 2 4 4       7 3 13 7   0 5 − M =  8 4 8 8     9 1 9 3    0 − − −9    4 2 4 4     11 1 9 7 0 − 7 4 2 4 4 19 El m´aximo elemento de A2 en valor absoluto es dado por m24 = − . Por 4 tanto, realizamos las operaciones elementales: C2 ↔ C4 , luego: Ind = ( 2 4 3 1 5 ), 

8

3

−2 −1 −1

   0 − 19 7  4 2    13 3  M = 0 −8 4   9 1   0 −  4 2   9 1 0 − 4 2

0



  11      7 7 5   8 8   3 9  − −9   4 4   11 7 7 4 4 9 4

5 4

Ahora hacemos cero los elementos mi2 (i = 3, 4, 5) mediante las operaciones elementales:   13 F3 ← F3 + − F2 , 38   9 F4 ← F4 + F2 , 19   9 F5 ← F5 + − F2 , 19

´ 4.1. METODOS DIRECTOS

123

resultando: 

8

3

−2

  7  0 − 19  4 2    17  0 0 − M = 38   22   0 0  19   22 0 0 − 19

−1

−1

9 4

5 4

11

2 19

17 38

47 38

0

63 3 72 − − 19 19 19 32 19

22 19

34 19

                

63 El m´aximo elemento de A3 en valor absoluto es dado por m44 = . Por 19 tanto, realizamos las operaciones elementales: F3 ↔ F4 ,

C3 ↔ C4 ,

luego: Ind = ( 2 4 1 3 5 ), 

8

3

   0 − 19  4     0 M = 0     0 0    0 0

−1

−2

−1

9 4

7 2

5 4

63 19

22 19

−

0 11

3 72 − 19 19

2 17 − 19 38

17 38

47 38

22 32 − 19 19

22 19

34 19

                

Ahora hacemos cero los elementos mi3 (i = 4, 5) mediante las operaciones elementales:   2 F4 ← F4 + − F3 , 63   32 F5 ← F5 + − F3 , 63

124

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

resultando: 

8

3

   0 − 19  4     0 M = 0     0 0    0 0

−1

−2

−1

9 4

7 2

5 4

63 19

22 19

−

0 11

72 3 − 19 19

0

−

61 126

19 42

19 14

0

−

110 63

26 21

26 7

                

110 . Por El m´aximo elemento de A4 en valor absoluto es dado por m54 = − 63 tanto, realizamos las operaciones elementales: F4 ↔ F5 , luego: 

8

3

   0 − 19  4     0 M = 0     0 0    0 0

−1

−2

−1

9 4

7 2

5 4

63 19

22 19

−

0 11

3 72 − 19 19

0

−

110 63

26 21

26 7

0

−

61 126

19 42

19 14

                

Ahora hacemos cero los elementos mi4 (i = 5) mediante las operaciones elementales:   61 F4 , F5 ← F5 + − 220 resultando el vector de ´ındices de las variables: Ind = ( 2 4 1 3 5 ),

´ 4.1. METODOS DIRECTOS

125

y la matriz aumentada queda:   0 8 3 −1 −2 −1     19 9 7 5  0 −  11   4 4 2 4      63 22 3 72    0 − − M = 0 19 19 19 19      110 26 26    0  0 0 −  63 21 7     6 18  0 0 0 0 55 55 aplicamos ahora el Algoritmo 6 al sistema anterior, resulta: i = 5 : xInd(5) =

b5 u55 b4 −

u44 b3 −

5 X

u33 b2 −

5 X

u22 b1 −

5 X

⇒ x3 = 0

= −1 ⇒ x1 = −1

u2j xInd(j)

j=3

i = 2 : xInd(2) =

=0

u3j xInd(j)

j=4

i = 3 : xInd(3) =

⇒ x5 = 3

u4j xInd(j)

j=5

i = 4 : xInd(4) =

i = 1 : xInd(1) =

5 X

=3

= −2 ⇒ x4 = −2

u1j xInd(j)

j=2

u11

=1

obteniendo la soluci´on x = (−1, 1, 0, −2, 3).

4.1.4.

⇒ x2 = 1 

Factorizaci´ on LU

Ejemplo 4.1.10. Consideremos la siguiente matriz   0 1 A= 2 3

(4.5)

126

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

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

= = = =

l11 u11 l11 u12 l21 u11 l21 u12 + l22 u22

(4.7) (4.8) (4.9) (4.10)

De (4.7) se obtiene: l11 = 0 o u11 = 0. Si l11 = 0 entonces contradecimos (4.8). Si u11 = 0 entonces contradecimos (4.9). Por tanto, podemos concluir que no existen matrices L y U tal que A = LU , a pesar que la matriz A es no singular.  Ejemplo 4.1.11. Consideremos la siguiente matriz   0 10 4 A= 0 2 1  2 3 1 observamos que es decir  0  0 A= 2

(4.11)

det(A) = 4 6= 0. Deseamos calcular L y U tal que A = LU ,     u11 u12 u13 10 4 l11 0 0 2 1  =  l21 l22 0   0 u22 u23  3 1 l31 l32 l33 0 0 u33

(4.12)

De (4.12) obtenemos el siguiente sistema de ecuaciones: 0 = l11 u11 , 10 = l11 u12 , 4 = l11 u13

(4.13)

De (4.13) podemos observar que l11 6= 0 y as´ı u11 = 0. Otra vez, a partir de (4.13) se observa que: 0 = l21 u11 2 = l31 u11

(4.14)

Como u11 = 0 se puede observar que se contradice la segunda ecuaci´on de (4.14). Por tanto, la matriz A no admite descomposici´on LU a pesar que A es una matriz no singular. 

´ 4.1. METODOS DIRECTOS

127

El ejemplo anterior nos motiva a establecer condiciones sobre una matriz A tal que pueda existir la factorizaci´on LU . Con este fin, damos primero las siguientes definiciones. Definici´ on 4.7 (Menor Principal). Dada una matriz A = (aij ) ∈ Rn×n , el menor complementario de un elemento aij es el determinante de la submatriz de A que se obtiene al eliminar la fila i y la columna j. 1. Los menores complementarios de una matriz A = (aij ) ∈ Rn×n son denotados por Aij . 2. Los menores complementarios Aii (i = 1, . . . , n) son llamados menores principales.

Ejemplo 4.1.12. Halle los menores complementarios de la siguiente matriz: 

 −1 3 1 A =  2 −1 −1  . −2 0 −1 Resoluci´ on Tenemos: A11 A21 A31

−1 −1 2 −1 2 −1 = −2; = = 1, A12 = = −4, A13 = 0 −1 −2 −1 −2 0 3 1 −1 1 −1 3 = 6; = = −3, A22 = = 3, A23 = 0 −1 −2 −1 −2 0 3 −1 1 −1 3 1 = −5. = −2, A32 = = −1, A33 = = −1 −1 2 −1 2 −1

Teorema 4.1. Existencia de descomposici´on LU . FALTA ESCRIBIR Teorema 4.2. Si una matriz A ∈ Rn×n es no singular y admite factorizaci´on LU , entonces esta factorizaci´on es u ´nica.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

128

4.1.5.

Factorizaci´ on por matrices elementales

Factorizaci´ on de Doolitle

Suponiendo que la matriz A admite factorizaci´on LU , para hallar las matrices L y U partimos de:

      

  a11 a12 · · · a1n  a21 a22 a2n     a31 a32 a3n  =  .. . . . ..    . . an1 an2 · · · ann

l11 0 · · · l21 l22 l31 l32 .. ... .

0 0 0 .. .

ln1 ln2 · · · lnn

      

u11 u12 · · · u1n 0 u22 u2n 0 0 u3n .. . . . .. . . 0 0 · · · unn

      

El razonamiento es como sigue, formamos un sistema de ecuaciones manteniendo la fila constante y hacemos variar todas las columnas. Veamos el caso de mantener la fila 1 constante y variamos todas las columnas, as´ı obtenemos:

a11 = l11 u11 a12 = l11 u12 a13 = l11 u13 .. .

(4.15)

a1n = l11 u1n

Del sistema (4.15) podemos observar que es posible calcular los valores u1j para todo j = 1, . . . , n si asignamos cualquier valor lii 6= 0.

´ 4.1. METODOS DIRECTOS

129

Algoritmo 9: Factorizaci´on LU por el m´etodo de Doolitle Entrada: Ingresar una matriz A ∈ Rn×n . Asignar valores no nulos para lii (i = 1, . . . , n). 1 inicio 2 para i ← 1 a n hacer 3 para j ← 1 a i − 1 hacer j−1 X aij − lik ukj lij ←

4

ujj

;

fin para para j ← i a n hacer i−1 X aij − lik ukj

5 6

uij ←

7

k=1

lii

;

fin para fin para devolver Matrices L y U .

8 9 10 11

k=1

fin

Ejemplo 4.1.13. Encuentre la factorizaci´on LU por el m´etodo de Doolittle de la siguiente matriz   2 2 1 A = 4 7 2 2 11 5 Resoluci´ on: Primero analizamos los menores principales de la matriz A para garantizar la existencia de la descomposici´on LU . Se debe de cumplir Aii 6= 0 para todo i = 1, 2, 3. Veamos: 7 2 2 1 2 2 = 6 6= 0 A11 = = 20 6= 0, A22 = = 8 6= 0, A33 = 11 6 2 5 4 7 Por tanto, por el Teorema 4.1, existe la factorizaci´on LU de la matriz A. Vamos a calcular las matrices L y U usando el m´etodo descrito de Doolitle, para esto, partimos de 

    2 2 1 l11 0 0 u11 u12 u13 4 7 2 = l21 l22 0   0 u22 u23  2 11 5 l31 l32 l33 0 0 u33

130

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES Columna Fila 1 Columna Columna Columna Fila 2 Columna Columna Columna Fila 3 Columna Columna

1 2 3 1 2 3 1 2 3

2 = l11 u11 u11 = 2 l11 =1 2 = l11 u12 =⇒ u12 = 2 1 = l11 u13 u13 = 1 4 = l21 u11 l21 = 2 l22 =1 7 = l21 u12 + l22 u22 =⇒ u22 = 3 2 = l21 u13 + l22 u23 u23 = 0 2 = l31 u11 l31 = 1 l33 =1 11 = l31 u12 + l32 u22 =⇒ l32 = 3 5 = l31 u13 + l32 u23 + l33 u33 u33 = 1

Por tanto, las matrices L y U son: 

 1 0 0 L = 2 1 0 , 1 3 1

  2 2 1 U = 0 3 0 0 0 4

Factorizaci´ on de Crout Algoritmo 10: Factorizaci´on LU por el m´etodo de Crout Entrada: Ingresar una matriz A ∈ Rn×n . Asignar valores no nulos para uii (i = 1, . . . , n). 1 inicio 2 para j ← 1 a n hacer 3 para i ← 1 a j − 1 hacer i−1 X aij − lik ukj uij ←

4

lii

;

fin para para i ← j a n hacer j−1 X aij − lik ukj

5 6

lij ←

7

k=1

ujj

;

fin para fin para devolver Matrices L y U .

8 9 10 11

k=1

fin

Ejemplo 4.1.14. Encuentre la factorizaci´on LU por el m´etodo de Crout de

´ 4.1. METODOS DIRECTOS

131

la siguiente matriz   −1 1 −1 A= 0 2 1  1 0 −1 Resoluci´ on: Primero se garantiza la existencia de la descomposici´on LU verificando que Aii 6= 0 para todo i = 1, 2, 3. Veamos: 2 1 −1 −1 −1 1 = −2. A11 = = −2, A22 = = 2, A33 = 0 −1 1 −1 0 2 Por tanto, por el Teorema 4.1, existe la factorizaci´on LU de la matriz A. Vamos a calcular las matrices L y U usando el m´etodo de Crout, para esto, partimos de      −1 1 −1 l11 0 0 u11 u12 u13  0 2 1  = l21 l22 0   0 u22 u23  1 0 −1 l31 l32 l33 0 0 u33 Fila Fila Fila Fila Fila Fila Fila Fila Fila

1 −1 = l11 u11 2 Columna 1 0 = l21 u11 3 1 = l31 u11 1 1 = l11 u12 2 Columna 2 2 = l21 u12 + l22 u22 3 0 = l31 u12 + l32 u22 1 −1 = l11 u13 2 Columna 3 4 = l21 u13 + l22 u23 3 −1 = l31 u13 + l32 u23 + l33 u33

Por tanto, las matrices L y  −1  0 L= 1

u11 =1

=⇒

u22 =1

=⇒

u33 =1

=⇒

l11 = −1 l21 = 0 l31 = 1 u12 = −1 l22 = 2 l32 = 1 u13 = 1 u23 = 2 l33 = −4

U son:  0 0 2 0 , 1 −4



 1 −1 1 U = 0 1 2  0 0 1

Factorizaci´ on de Cholesky Definici´ on 4.8. Una matriz A ∈ Rn×n es llamada Definida Positiva si y s´olo si para cualquier x ∈ Rn \{0} se cumple que xT Ax > 0.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

132

Ejemplo 4.1.15. Considerar la matriz   4 −1 −1 1 A = −1 4 −1 1 4 (a) Sin calcular, demostrar que existe la descomposicion de Cholesky de A. (b) Determinar la descomposicion de Cholesky de A. Resoluci´ on: (a) 

     4 −1 −1 λ 0 0 4 − λ −1 −1 1  −  0 λ 0  =  −1 4 − λ 1 . (A−λI) = −1 4 −1 1 4 0 0 λ −1 1 4−λ Calculando el determinante: 4 − λ −1 −1 −1 4 − λ 1 = −λ3 + 12λ2 − 45λ + 54 = −(λ − 3)2 (λ − 6) = 0 −1 1 4 − λ de donde se obtiene λ = 3 (multiplicidad 2) y

λ = 6.

Como todos los valores propios son positivos entonces A es una matriz definida positiva. (b) 

    4 −1 −1 l11 0 0 l11 l21 l31 1  = l21 l22 0   0 l22 l32  A = LLT = −1 4 −1 1 4 l31 l32 l33 0 0 l33 de donde obtenemos: 4 −1 −1 4 1 4

= l11 2 = 2l21 = 2l31 = 41 + l√22 2 = 14 + 215 l32 15 = 14 + 100 + l33 2

=⇒ =⇒ =⇒ =⇒ =⇒ =⇒

l11 l21 l31 l22 l32 l33

=2 = 1/2 = −1/2 √ = √215 15 = 10 √ = 3 510

´ 4.1. METODOS DIRECTOS

133

Por tanto: 

 2 0 √0  L = −1/2 √ 15/2 √0 −1/2 15/10 3 10/5

4.1.6.

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): " # 2 n−1 X X A = ||u1 ||q1 ||u2 ||q2 + hq1 , v2 iq1 ||u3 ||q3 + hqi , v3 iqi · · · ||un ||qn + hqi , vn iqi i=1

i=1

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

134

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.16. 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

135

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.16) −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

136

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.17)

donde M es conocida como Matriz de iteraci´ on y d es llamado Vector de iteraci´ on. El sistema (4.17) genera una secuencia de puntos xk dados por xk = M xk−1 + d,

k ≥ 1,

(4.18)

que convergen a la soluci´on del sistema lineal equivalente (4.1). Teorema 4.3. Si para alguna norma matricial k · k subordinada a una norma vectorial se cumple que kM k < 1, entonces el esquema iterativo (4.18) converge a la soluci´on del sistema (4.17). Demostraci´ on: Sea xˆ soluci´on exacta del sistema (4.17) y xk un punto generado por el esquema iterativo dado por (4.18). kˆ x − xk k = kM xˆ + d − (M xk−1 + d)k = kM (ˆ x − xk−1 )k, desde que k · k es una norma subordinada a una norma vectorial, entonces: kˆ x − xk k ≤ kM kkˆ x − xk−1 k, realizando este mismo procedimiento ” k ” veces, se obtiene: kˆ x − xk k ≤ kM kk kˆ x − x0 k,

(4.19)

y dado que kM k < 1, entonces cuando k → ∞ se tiene que xk → xˆ. Definici´ on 4.9. Definimos el radio espectral de una matriz M ∈ Mn como sigue: ρ(M ) = m´ax{|λ| / λ es un autovalor de la matriz M }. Teorema 4.4. Dado el esquema iterativo (4.18), se cumple que converge a la soluci´on exacta del sistema (4.17) si y s´olo si ρ(M ) < 1. Demostraci´ on:

´ 4.2. METODOS INDIRECTOS

137

(⇒) Los puntos xk de la secuencia dada por (4.18) convergen a la soluci´on exacta del sistema (4.17). Supongamos ahora que ρ(M ) ≥ 1, entonces existe λ autovalor de la matriz M tal que |λ| ≥ 1, y existe un vector no nulo v tal que M v = λv. Consideremos el punto inicial x0 = xˆ + v, entonces: kˆ x − xk k = kM k (ˆ x − x0 )k = |λ|k k(ˆ x − x0 )k → ∞ cuando k → ∞, lo cual es una contradicci´on. Por tanto, ρ(M ) < 1. (⇐) Para aproximar la soluci´on exacta del sistema lineal Ax = b consideremos la siguiente descomposici´on para la matriz A: A=D+L+U

(4.20)

donde 

a11 0 · · · 0 0  0 a22 0 0   . .. D=   0 0 an−1,n−1 0 0 0 ··· 0 an,n 

0 a21

0 0

0 0

···

    ,  

0 0

   .. L=  .   an−1,1 an−1,2 an−1,3 · · · 0 an,1 an,2 an,3 · · · an,n−1 

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

 0 0      0  0

(4.21)

(4.22)

    .  

(4.23)

Los m´etodos iterativos a estudiar se obtienen eligiendo un sistema lineal equivalente de una forma conveniente a partir de las matrices D, L, U . Veremos esto en detalle en cada una de las siguientes secciones.

138

4.2.1.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

M´ etodo de Jacobi

Reemplazando la descomposici´on de A = D +L+U en el sistema Ax = b: (D + L + U )x = b siempre que D sea invertible, podemos obtener: x = −D−1 (L + U )x + D−1 b.

(4.24)

el esquema de punto fijo asociado a (4.24) es xk+1 = −D−1 (L + U )xk + D−1 b

(4.25)

por tanto, la matriz de iteraci´on asociada al m´etodo de Jacobi es MJ = −D−1 (L + U )

(4.26)

y el vector de iteraci´on es dJ = D−1 b. De la definici´on de las matrices D, L y U se tiene que

D−1

 1  a11   0    =    0   0 

0 ··· 1 a22 .. .

0 a21 .. .

0

0

0

1

0 0

0

0

an−1,n−1 ··· a12 0

1

0

       ,     

y

an,n

· · · a1,n−1 a2,n−1 ...

a1,n a2,n .. .

   L+U =   an−1,1 an−1,2 0 an−1,n an,1 an,2 · · · an,n−1 0

      

reemplazando las expresiones de D−1 y L+U en la expresi´on (4.25) se obtiene para cada i = 1, 2, . . . , n la siguiente relaci´on: bi − xk+1 = i

n X j=1,j6=i

aii

aij xki .

(4.27)

´ 4.2. METODOS INDIRECTOS

139

Ejemplo 4.2.1. Considere el sistema de ecuaciones lineales: −8x + 5y − z = −1.3 2x − 5y + z = 0.7 3x − 3y + 7z = 3.8 (a) Calcular la matriz asociada al m´etodo de Jacobi. (b) Verifique que el sistema es convergente. (c) Calcule aplicando el m´etodo de Jacobi, una aproximaci´on de la soluci´on del sistema con un error inferior a 0.02. Resoluci´ on: (a) Las matrices D,  −8 0 D =  0 −5 0 0

L y U son:    0 0 0 0 0 , L =  2 0 0 , 7 3 −3 0



 0 5 −1 U =  0 0 1 . 0 0 0

Usando la relaci´on (4.26) se calcula la matriz asociada al m´etodo iterativo de Jacobi:  −1   5 1 1  0 −8 8   −8 0 0        0 5 −1     1 1   2     2 0 1 = − − MJ = −  0 − 0  0 − .  5   5 5  3 −3 0         3 1 3 − 0 0 0 7 7 7 (b) Para garantizar la convergencia del m´etodo, basta observar que la matriz A es diagonal dominante. (c) (toda esta parte falta dado por  laexpresi´on 0 0  0  X = 0  −1  1  x1  0 10     1   1  x2  =  0    200      −47 x13 0 4000

corrigir)El esquema iterativo de Jacobi viene (4.25), por tanto, usando como punto inicial

   −3 19  0  x   10     1   10   1.9    −87   0   −249   −1.245  +   x2  =   200  200    3.47575     13903  489  x03 4000 4000

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

140

  0     2    x2  =  0         x23 0 



x21

x31







0

    3    x2  =  0       3 0 x3   0     4    x2  =  0         x43 0 

x41



−1 10 1 200 −47 4000

   −3 19  1  x   10     1   10   0.981775     −87   1   −249   −2.763176    x2 + =   200  200   3.915289    13903  489  x13 4000 4000

−1 10

−3 10

1 200

−87 200

−47 4000

489 4000

−1 10 1 200 −47 4000



x21

 

    2     x2 +       2 x3

19 10 −249 200 13903 4000





1.001730



       =  −2.961966        3.986861

   19 −3  3  x   10     1   10   1.000138     −87   3   −249   −2.994094    x2 + =   200  200   3.997946    13903  489  x33 4000 4000

Mientras que la soluci´on exacta del sistema es   1 x∗ =  −3  4 Notar que con 4 iteraciones el error es menor a 0.02

4.2.2.

M´ etodo de Gauss-Seidel

Reemplazando (4.20) en (4.1) se obtiene: x = (D − L)−1 U x + (D − L)−1 b.

(4.28)

el esquema de punto fijo asociado a (4.24) es xk+1 = (D − L)−1 U xk + (D − L)−1 b

(4.29)

´ 4.2. METODOS INDIRECTOS

141

De la definici´on de las matrices D, L y U se tiene que Ejemplo 4.2.2. Considere el sistema de ecuaciones lineales: 20x + 2y + 6z = 38 x + 20y + 9z = −23 2x − 7y − 20z = −57 (a) Calcular la matriz asociada al m´etodo de Gauss–Seidel para la resoluci´on aproximada del sistema. (b) Verifique que el sistema es convergente. (c) Calcule aplicando el m´etodo de Gauss–Seidel, una aproximaci´on de la soluci´on del sistema con un error inferior a 0.02. Resoluci´ on: (a) Las matrices D, L y U son:     20 0 0 0 0 0 D =  0 20 0  , L =  −1 0 0  , 0 0 −20 −2 7 0



 0 −2 −6 U =  0 0 −9  0 0 0

La matriz asociada al m´etodo iterativo de Gauss–Seidel viene dado por: MGS = (D − L)−1 U reemplazando se obtiene −1 10

 

MGS

 0  0 2 6   0 0 9 =  0  0 0 0   0

−1 

20 0 0  1 20 0  = 2 −7 20



1 200 −47 4000

 −3 10    −87   200   489  4000

(b) Para garantizar la convergencia del m´etodo, basta observar que la matriz A es diagonal dominante. (c) El esquema iterativo de Gauss–Seidel viene dado por: x(k+1) = MGS x(k) + d

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

142



 0 Usamos como punto inicial X 0 =  0  0     −3 −1 19   1   x01 x1  0 10  10   10            1.9     1     −249  1 −87    x2  =  0 x02    =  −1.245      +  200  200 200       3.47575      13903   0 1 −47 489 x3 x3 0 4000 4000 4000   0     2    x2  =  0         x23 0 



x21

x31







0

    3    x2  =  0       x33 0   0     4    x2  =  0         x43 0 

x41



−1 10 1 200 −47 4000

   19 −3  1  x   10     1   10   0.981775     −87   1   −249   −2.763176    x2 + =   200  200   3.915289      1 489 13903 x3 4000 4000

−1 10

−3 10

1 200

−87 200

−47 4000

489 4000

−1 10 1 200 −47 4000



x21

 

    2     x2 +       x23

19 10 −249 200 13903 4000





1.001730



       =  −2.961966        3.986861

   19 −3  3  x   10     1   10   1.000138    −87   3   −249   −2.994094  +   x2  =   200  200    3.997946     13903  489  x33 4000 4000

Mientras que la soluci´on exacta del sistema es   1 x∗ =  −3  4 Notar que con 4 iteraciones el error es menor a 0.02

´ 4.2. METODOS INDIRECTOS

143

Ejemplo 4.2.3. Considere el siguiente polinomio caracter´ıstico P (λ) = −λ3 + 1.08λ − 0.432. Suponga que P (λ) es el polinomio caracter´ıstico de la matriz de iteraci´on M de alguno de los m´etodos iterativos estudiados (GaussJacobi o Gauss-Seidel), cuando es aplicado al sistema lineal Ax = b, responda las siguientes preguntas: (a) Indigue, justificando adecuadamente, ¿cu´al es el m´etodo estudiado? (b) Averigue si el m´etodo en  a  b (c) Sabiendo que M = b

estudio es convergente.  b b c b , determine los valores de a, b, c y d. b d

(d) Suponiendo que los elementos de la diagonal de la matriz A del sistema son todos iguales a 5, determine los elementos restantes de A. Resoluci´ on: (a) Veamos primero el caso del m´etodo iterativo de Gauss-Seidel.  −1  0 a12 a13 a11 0 0 = −(D + L)−1 U = −  a21 a22 0   0 0 a23  . 0 0 0 a31 a32 a33 

MGS

Como la matriz inversa de una matriz triangular inferior es otra matriz triangular inferior, entonces MGS puede expresarse como: −1   (1)    a11 0 0 0 a12 a13 0 m12 m13   MGS = −  a(1) a122 0   0 0 a23  =  0 m22 m23  . 21 (1) (1) (1) 0 0 0 0 m31 m33 a31 a32 a33 Ahora calculamos su polinomio caracter´ıstico: PGS (λ) = det(MGS − λI) = −λ(λ − m22 )(λ − m33 ) + m32 m23 λ. Se puede observar que PGS no tiene t´ermino independiente, as´ı, P (λ) no es el polinomio caracter´ıstico asociado al m´etodo de Gauss-Seidel. (b) bla (c) bla

144

4.2.3.

CAP´ITULO 4. SISTEMAS DE ECUACIONES LINEALES

M´ etodo de relajaci´ on

De la descomposici´on A = D + L + U , considere un par´ametro no nulo w ∈ R y la siguiente descomposici´on:     1 w−1 A= D+L + D+U , (4.30) w w reemplazando (4.30) en el sistema lineal Ax = b se obtiene:     w−1 1 D+L + D+U x = b, w w de donde tenemos el siguiente sistema equivalente:     1 w−1 D+L x=− D + U x + b, w w por tanto, se tiene el siguiente esquema de punto fijo: −1    −1  w−1 1 1 D+L D+U x+ D+L b, x=− w w w de donde resulta el siguiente esquema iterativo de punto fijo:  −1    −1 1 w−1 1 k k+1 x =− b, D+L D+U x + D+L w w w

(4.31)

(4.32)

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 | 145

(5.3)

´ ´ 146CAP´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.

147

´ 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)

´ ´ 148CAP´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

149

´ 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) para alg´ un entero positivo n. Definimos el polinomio de Taylor de la funci´on f alrededor del punto x0 ∈ I del modo siguiente:

Pn,x0 (x) =

n X f (i) (x0 ) i=0

i!

(x − x0 )i (5.9)

f (n) (x0 ) Pn,x0 (x) = f (x0 ) + f 0 (x0 )(x − x0 ) + ... + (x − x0 )n n! Dados x, x0 ∈ I, existe s entre x y x0 tal que f (t) = Pn,x0 (t) + Rn+1,x0 (t)

(5.10)

´ ´ 150CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL Teorema 5.1. Sea I un intervalo, x0 ∈ I, f : I → R y dado n ∈ N existen f, f 0 , f 00 , ..., f (n) , f (n+1) en I ◦ . Entonces, para cualquier x ∈ I con x 6= x0 existe un punto c en el intervalo abierto de extremos x0 y x tal que: f (x) − Pn,x0 (x) = f (x) −

n X f (i) (x0 ) i=0

i!

(x − x0 )i =

f (n+1) (c) (x − x0 )n+1 . (n + 1)!

La expresi´on: f (n+1) (ξ) (x − x0 )n+1 (n + 1)! es conocido como Resto en la forma de Lagrange. Rn+1,x0 (x) =

(5.11)

Demostraci´ on: Sean x0 , x ∈ I valores fijos y tal que x0 6= x. Sin p´erdidad de generalidad, podemos suponer que x0 < x. Consideremos el intervalo J = [x0 , x], por tanto J ⊂ I. Sobre J definamos las siguientes funciones: ϕ:J → R t 7→ ϕ(t) :=

n X f ( k)(t) k=0

k!

(x − t)k

ψ:J → R t 7→ ψ(t) := −(x − t)n+1

(5.12)

(5.13)

Desde que f (k) (0 ≤ k ≤ n) es continua y (x − t)k es una funci´on polin´omica (y por tanto continua e infinitamente derivables), podemos garantizar que ϕ es continua en J. Adem´as, como f (k) (k = 0, . . . , n) son derivables en I ◦ entonces tambi´en son derivables en J ◦ . Por otro lado, podemos observar que ψ es una funci´on infinitamente derivable en J. De lo anterior, tenemos que ϕ y ψ son funciones continuas en J y derivables en J ◦ , as´ı, podemos aplicar el Teorema del Valor Medio Generalizado, esto es, existe c ∈ hx0 , xi tal que: ψ 0 (c)(ϕ(x) − ϕ(x0 )) = ϕ0 (c)(ψ(x) − ψ(x0 ))

(5.14)

A partir de 5.12 y 5.13 calculemos cada expresi´on: ϕ(x) = f (x) ϕ(x0 ) = Pn,x0 (x) n n X X f (k+1) (t) f (k) (t) f (n+1) (t) ϕ0 (t) = (x − t)k − k(x − t)k−1 = (x − t)n k! k! n! k=0 k=1 ψ(x) = 0 ψ(x0 ) = −(x − x0 )n+1 ψ 0 (t) = (n + 1)(x − t)n

´ 5.1. METODOS DE TAYLOR

151

Reemplazando estos resultados en 5.14 se obtiene: (n + 1)(x − c)n (f (x) − Pn,x0 (x)) =

f (n+1) (c) (x − c)n (x − x0 )n+1 n!

de lo anterior obtenemos: f (n+1) (c) Rn+1,x0 (x) = f (x) − Pn,x0 (x) = (x − x0 )n+1 , (n + 1)!

c ∈ hx0 , xi.

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: P4 (x) = 2 +

x − 1 (x − 1)2 (x − 1)3 5 (x − 1)4 − + − 4 64 512 16384

luego P4 (5) =

179 = 2.79688, 64

´ ´ 152CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL encuanto el valor exacto es √ 5 = 2.23607.

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. 5.2.1.

M´ etodos de Runge-Kutta Runge-Kutta de 2 pasos

Estos presentan la siguiente forma: K1 = f (ti , yi ) K2 = f (ti + αh, yi + βK1 h) yi+1 = yi + h(aK1 + bK2 )

(5.15)

donde se cumple las siguientes relaciones: a+b = 1 1 αb = 2 1 βb = 2 Entre los casos particulares m´as conocidos tenemos:

(5.16)

´ 5.3. METODOS MULTIPASOS

153

1 1. M´ etodo de Heun: Se presenta cuando b = , resultando: 2 K1 K2

= f (ti , yi ) = f (ti + h, yi + K1 h) h yi+1 = yi + (K1 + K2 ) 2 Tambi´en es llamado M´ etodo de Euler Mejorado.

(5.17)

2. M´ etodo de Euler Modificado: Se presenta cuando b = 1, resultando: K1 = f (ti , yi ) K2 = f (ti + h/2, yi + K1 h/2) yi+1 = yi + hK2

5.2.2.

(5.18)

Runge-Kutta de 3 pasos

Estos presentan la siguiente forma: K1 K2 K3

= f (ti , yi ) = f (ti + h/2, yi + hK1 /2) = f (ti + h, yi + 2hK2 − hK1 ) h = yi + (K1 + 4K2 + K3 ) 6

yi+1

5.2.3.

(5.19)

Runge-Kutta de 4 pasos

Estos presentan la siguiente forma: K1 K2 K3 K4 yi+1

= = = =

f (ti , yi ) f (ti + h/2, yi + K1 h) f (ti + h/2, yi + K2 h) f (ti + h, yi + K3 h) h = yi + (K1 + 2K2 + 2K3 + K4 ) 6

5.3.

M´ etodos Multipasos

5.4.

Ejercicios Propuestos

(5.20)

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 .

´ ´ 154CAP´ITULO 5. METODOS NUMERICOS PARA PROBLEMAS DE VALOR INICIAL 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. 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: 1 y(x) = 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

Cap´ıtulo 6 Pr´ acticas Calificadas 6.1.

Examen Parcial 2015-2

Problema 1: Conteste verdadero o falso en cada caso con la justificaci´on: a. El coeficiente del t´ermino de mayor grado del polinomio de interpolaci´on que cumple f(0)=1; f(1)=2; f’(1)=3 es 2. b. Al calcular −10π + 6e − 3/62 usando aritm´etica de tres d´ıgitos con redondeo el error relativo es 3.6 × 10−3 . c. La iteraci´o√ n de punto fijo xi+1 = xi + 0.25(xi 2 − 3) converge siempre a 3 si se escoge un punto inicial x0 adecuado. d. Si f (x) = xn−1 para n ≥ 1 entonces si {x1 , x2 , ..., xn+1 } son distintos puntos de interpolaci´on entonces f [x1 , ..., xn ] = 1. e. El m´etodo de Newton puede considerarse como un caso l´ımite del m´etodo de la falsa posici´on. Resoluci´ on: a. Verdadero. Teniendo tres datos podemos considerar que el polinomio es de grado 2 y por el m´etodo de coeficientes indeterminados se tiene: P (X) = ax2 +bx+c y de las condiciones iniciales establecemos: P (0) = c = 1 y de las otras a + b + c = 2 y tambi´en 2a + b = 3 que nos da las ecuaciones: a + b + 1 = 2 y 2a + b = 3 155

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

156

de donde se obtiene a = 2; b = −1. Por tanto, el coeficiente buscado es a = 2. b. Verdadero. −10π + 6e − 3/62 = −31.4 + 16.3 − 0.0484 = −15.1 y | − 15.1 − (−15.154622)| n = −15.154622 tenemos: er = = 3.6 × 10−3 . | − 15.154622| c. Falso. El MPF con g(x) = x + 0.25(x2 − 3) entonces g 0 (x) = 1 +√0.5x y como se tendr´ıa que tener convergencia los xi ser´ıan cercanos a 3 y por tanto g 0 (x) > 0 lo cual no asegura la convergencia. d. Verdadero. Podemos usar la f´ormula f [x1 , ..., xn ] =

f (n−1) (a) y dife(n − 1)!

renciando f (x) = xn−1 se produce lo siguiente: f 0 (x) = (n − 1)xn−2 f 00 (x) = (n − 1)(n − 2)xn−3 .. . f (n−1) (x) = (n − 1)! f (n) (x) = 0 y por tanto: f [x1 , ..., xn ] =

f (n−1) (a) = 1. (n − 1)!

e. Falso. Es caso l´ımite del m´etodo de la secante pero no del m´etodo de la falsa posici´on ya que este requiere que se cumpla f (a)f (b) < 0. Problema 2: Se pide: −3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) es una apro2h ximaci´on para f 0 (x0 ).

a. Demuestre que

b. Dar la expresi´on para el t´ermino principal del error. c. Si f (x) = xsen(x) + x2 cos(x) y h = 0.1, establecer el valor de f 0 (1.2). Resoluci´ on:

6.1. EXAMEN PARCIAL 2015-2

157

a. Desarrollamos en serie de Taylor cada una de las expresiones presentes en −3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) E= 2h como sigue: −3f (x0 ) = −3f (x0 )   h2 f 00 (x0 ) h3 f 000 (x0 ) 0 4f (x0 + h) = 4 f (x0 ) + hf (x0 ) + + + ... 2 3! (2h)2 f 00 (x0 ) (2h)3 f 000 (x0 ) + + ... −f (x0 + 2h) = − f (x0 ) + (2h)f 0 (x0 ) + 2! 3! Por tanto al sumar para obtener E: E = f 0 (x0 ) −

2h2 000 f (x0 ) + . . . 6

(6.1)

por tanto: f 0 (x0 ) ≈

−3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) . 2h

b. A partir de (6.1) podemos concluir que el t´ermino principal del error es 2 2h 000 2h2 000 T = f (x0 ) = |f (x0 )|. 6 6 c. De acuerdo a la expresi´on para f 0 (x0 ), considerando h = 0.1, el valor aproximado de f 0 (1.2) es: f 0 (1.2) ≈

−3f (1.2) + 4f (1.3) − f (1.4) 2 · (0.1)

f 0 (1.2) ≈

−3(1.64024) + 4(1.70470) − 1.71277 0.2

f 0 (1.2) ≈ 0.92655. Problema 3: Usando el MN y el MB determine la ra´ız m´axima de f (x) = x3 − 6x2 + 11x − 6.1 , con un error de 10−3 . Para MN usar x0 = 3 y para MB usar [a, b] = [3, 3.5].

!

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

158

Resoluci´ on: Usando el m´etodo de Newton se debe de usar para todo k ≥ 0 la siguiente relaci´on: xk+1 = xk −

f (xk ) f 0 (xk )

xk+1 = xk −

x3k − 6x2k + 11xk − 6.1 . 3x2k − 12xk + 11

Los resultados usando el m´etodo de Newton son mostrados en la Tabla 6.1. i 0 1 2 3

xi 3.00000000 3.05000000 3.04669556 3.04668053

error 3.00000000 0.05000000 0.00330444 0.00001503

Tabla 6.1: Resultados del m´etodo de Newton

Para el m´etodo de la bisecci´on, iniciando en [a, b] = [3, 3.5], se verifica que f (a)f (b) < 0, por tanto, los nuevos intervalos en cada iteraci´on i se muestran en la Tabla 6.2. i 0 1 2 3 4 5 6 7 8

a 3.000000 3.000000 3.000000 3.000000 3.031250 3.031250 3.039063 3.042969 3.044922

c 3.250000 3.125000 3.062500 3.031250 3.046875 3.039063 3.042969 3.044922 3.045898

b 3.500000 3.250000 3.125000 3.062500 3.062500 3.046875 3.046875 3.046875 3.046875

f (a) -0.100000 -0.100000 -0.100000 -0.100000 -0.034540 -0.034540 -0.017238 -0.008444 -0.004012

f (c) 0.603125 0.198828 0.036963 -0.034540 0.000445 -0.017238 -0.008444 -0.004012 -0.001786

Tabla 6.2: Resultados del m´etodo de la bisecci´on

f (b) 1.775000 0.603125 0.198828 0.036963 0.036963 0.000445 0.000445 0.000445 0.000445

error 0.500000 0.250000 0.125000 0.062500 0.031250 0.015630 0.007810 0.00391 0.00195

6.1. EXAMEN PARCIAL 2015-2

159

Problema 4: El dinero necesario para pagar la cuota correspondiente a un cr´edito hipotecario a inter´es fijo se suele estimar mediante la denominada A[1 − (1 + i)−n ] ’ecuaci´on de la anualidad ordinaria’ Q = en donde i Q es la cantidad (en d´olares) pedida del pr´estamo, A es la cuota que debe pagar el beneficiario del pr´estamo, i es la tasa de inter´es (en tanto por 1) fijado por la entidad bancaria que concede el pr´estamo y n es el n´ umero de periodos durante los cuales se realizan pagos de la cuota. Una pareja desea adquirir una vivienda mediante un pr´estamo de 150000 d´olares a pagar semestralmente durante un plazo de 15 a˜ nos. La cantidad m´axima que ellos desean pagar es de 900 d´olares mensuales, calcular, usando MN cu´al ser´ıa el tipo m´aximo de inter´es al que pueden negociar su pr´estamo con las entidades bancarias, Use el MN con i0 = 0.03. Resoluci´ on: Puesto que el pago es semestral, en 15 a˜ nos realizaran un total de 30 cuotas dem´as dado que pueden pagar 900 euros al mes, cada semestre podr´a afrontar el pago de 5400 euros. Ello hace que la ecuaci´on de la anualidad ordinaria quede: 5400 [1 − (1 + i)−30 ] 150000 = i 5400 [1 − (1 + i)−30 ] = 0 f (i) = 150000 − i " # −30 5400 1 − (1 + i) f 0 (i) = − 30(1 + i)−31 . i i Por lo que el m´etodo de Newton con i0 = 0.03 nos conducir´a a la siguiente sucesi´on de valores: k 0 1 2 3 4

ik 0.03000000 0.00227545 0.00449705 0.00503585 0.00503892

error 0.03000000 0.03227545 0.00677249 0.00053881 0.00000307

Por tanto, el inter´es al que pueden acceder es de i = 0.5 % aproximadamente.

160

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

Problema 5: Sea f : R → R una funci´on continua y acotada (es decir, ∃M/|f (x)|M, ∀x ∈ R). Demostrar que la ecuaci´on f (x) − x = 0, tiene al menos una ra´ız real.

Resoluci´ on: Como f (x) es acotada se tiene que existe M tal que −M < f (x) < M para todo x ∈ R, entonces definimos la funci´on h : [−M, M ] → R del modo siguiente: h(x) = f (x) − x. Observe que h es una funci´on continua y que adem´as h(−M ) = f (x) + M > 0 y h(M ) = f (M ) − M < 0. Es decir, h(M ) < 0 < h(−M ) y como h es continua, podemos aplicar el Teorema del Valor Intermedio para garantizar que existe c ∈< −M, M > tal que h(c) = 0, y esto completa la demostraci´on. Problema 6: Deduzca la expresi´on para determianr el n´ umero m´ınimo de iteraciones que se debe realizar en el M´etodo de la Bisecci´on, para asegurar que se logre un error menor a ε, si el intervalo donde se empieza la iteraci´on es [a, b].

Resoluci´ on: De la teor´ıa de la f´ormula para el error en el m´etodo de la Bisecci´on se b−a establece como sigue: en = |xn − n| ≤ n+1 de donde se debe cumplir 2 b−a en = |xn − n| ≤ n+1 < ε y para hallar el n´ umero de pasos necesarios se 2 b−a b−a tendr´a n+1 < ε que equivale a 2n+1 > , de donde tomando logaritmos 2 ε b−a naturales: (n + 1)ln(2) > ln( ) que equivale finalmente a: ε

b−a ) ε −1 ln(2)

ln( n>

6.1. EXAMEN PARCIAL 2015-2

161

Problema 7: Asuma que un polinomio de interpolaci´on pasa por los puntos: x f (x)

0 22

15 24

18 37

22 25

24 123

Construir la tabla de diferencias dividas y el polinomio de interpolaci´on de Newton.

Resoluci´ on: Para Newton construimos la tabla de diferencias divididas:

k 0 1 2 3

xk 0 15 18 22

f (xk ) 22 24 37 25

f [xk , xk+1 ]

f [xk , xk+1 , xk+2 ]

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

2/15 13/3 -3

7/30 -22/21

-269/4620

Por tanto el polinomio de Newton estar´a dado por:

N (x) = 22 +

2 7 269 x + x(x − 15) − x(x − 15)(x − 18). 15 30 4620

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

162

6.2.

Segunda Pr´ actica Calificada

Problema 1: Conteste verdadero o falso en cada caso con la justicaci´on: a. Si para un metodo de punto fijo se tiene x = g(x) = 1 + 0.5 sen(x) entonces la iteraci´on de PF converge a una soluci´on, sin importar el punto inicial x(0) . b. Si g : [0, 1] → [0, 1] y est´a definida por g(x) = |g(x)| < 8(3)9 0.5 y la serie de punto fijo converge.

1 (1+x2 )

entonces

1

c. Si se toma f (x) = 3x 3 ) con x(0) > 0 entonces la sucesi´on formada usando el M´etodo de Newton converge a la ra´ız. d. Si se toma los puntos de interpolaci´on {(15, 24), (18, 37), (22, 25)} entonces una aproximaci´on a f 00 (18) es -2.096. e. Dado n + 1 puntos no colineales el polinomio de interpolaci´on tiene siempre n ra´ıces reales. Resoluci´ on : a. Verdadero. x = g(x) = 1 + 21 sin(x) derivando g(x) → g 0(x) = 12 sin(x) 1 1 =⇒ −1 ≤ sin(x) ≤ 1 =⇒ 0 ≤ | sin(x)| ≤ 2 2 =⇒ |g 0(x) | ≤

1 ≤ 1 ←− 2

converge. b. Verdadero. g : [0, 1] → [0, 1] g(x) =

1 1+x2

−→ g 0(x) = −1(+x2 )−2 2x =

−2x 1+x2

9 =⇒ |g 0(x) | ≤ √ < 1 ←− 3 converge. 1

1

c. Falso. f(x) = 3x 3 g(x) = x −

3x 3 x

−2 3

−→ g(x) = x − 3x = −2x

−→ g 0(x) = −2 > 1

´ 6.2. SEGUNDA PRACTICA CALIFICADA

163

−→ |g 0(x) | = 2 > 1 No converge Si escogemos un x0 > 0 el |g(x) 0 | > 1 por lo tanto no se puede asegurar que converge. x0 = 15 x1 = 18 x2 = 82 d. V o F?. Como f (x0 ) = 24 f (x1 ) = 37 f (x2 ) = 25 f (x) =

24(24 − 18)(x − 22) 37(x − 15)(x − 22) 24(x − 18)(x − 15) + + (−3)(−7) (3)(−4) (4)(7) −22 2 817 x + x + cte 21 4 −22 ∗ 2 −→ f (18)00 ≈ −2.096 −→ f (x)00 = 21 f (x) =

e. Falso. Por (NM) puntos pasa un unico polinomio de grado ”n”por lo tanto el numero de raices tendra que ser n-raices ≤ n; no necesariamente tendra ” n ” ra´ıces. Problema 2: Una banda met´alica flexible de 10 m de longitud se dobla en forma de arco circular al asegurar los extremos entre si por medio de un cable de 8 pies de longitud. Use el m´etodo de Newton para aproximar el radio r del arco circular. Use como punto inicial r(0) = 4 con al menos 4 cde. Resoluci´ on: Del gr´afico 10 = 2θr,

sin(θ) =

4 r

4 → θ = arcsin( ) r 4 ⇒ arcsin( = 5 r 4 5 5 4 ⇒ arcsin( ) = ⇒ sin( ) = r r r r 5 r 5 r = csc( ) ⇒ f(r) = − csc( ) 4 r 4 r Por el m´etodo de Newton: rnH = rn −

rn 4

− csc( r5n )

( 14 + csc( r5n ) cot( r5n ) r52 ) n

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

164

csc( r5n ). cot( r5n ) r5n + csc( r5n )

rnH =

( 14 + csc( r5n ) cot( r5n ) r52 ) n

⇒ r0 = 4 r0 = 4 r5 = 4.363385 r1 = 4.14.. r6 = 4.38113085 r2 = 4.240266 r7 = 4.393291 r3 = 4.298663 r8 = 4.401660281 r4 = 4.337317 r9 = 4.407436912 ⇒ r = 4.4204 Problema 3: Usando la tabla mostrada a continuaci´on, construir el polinomio de interpolaci´on de Newton y calcular el valor aproximado de f 000 (0.7) (tercera derivada). x 0.6 0.8 0.9 1.1 f (x) 1.16454 1.51736 1.68333 1.99121 Resoluci´ on: x f (x)

0.6 1.16484

0.8 1.51736

0.9 1.68333

1.1 1.99121

000

f (α) f [x0 , x1 , x2 , x3 ] = 8! 000

f (α) = 6xf [x0 , x1 , x2 , x3 ] n 0 1 2 3

x1 0.6 0.8 0.9 1.1

f(x1 ) 1.16464 1.51736 1.68333 1.99121

f[x0 , x1 ] 1.7636 1.6597 1.5394

f[x0 , x1 ; x2 ] -0.346333 -0.411

→ f 000 (0.7) = 6(−0.129334) = −0.776004.

´ 6.2. SEGUNDA PRACTICA CALIFICADA

165

Problema 4: Asuma que se tienen los puntos de interporlaci´on x(0) , x(1) , x(2) igualmente espaciados por una cantidad h y tambien si |f (α)| ≤ M , determine una expresi´on (en t´erminos de h y M ) para el m´aximo error que se produce en el polinomio de interpolaci´on. Resoluci´ on: |f 000 (α)| ≤ M sabemos que: n

f (n+1) (α) Y Error = (x − xi ) (n + 1)! i=0 pero: x − x0 ≤ 2h x − x1 ≤ 2h x − x2 ≤ 2h

=⇒ |(x − x0 )(x − x1 )(x − x2 )| ≤ 8h3 f 000 (α) =⇒ Error = (x − x0 )(x − x1 )(x − x2 ) 3! 000 f (α) 8h3 M =⇒ (x − x0 )(x − x1 )(x − x2 ) ≤ . 3! 3! Problema 5: Halle el valor de x con una aproximacion de dos decimales en el cual la curva y = e−x ln(x) tiene un punto de inflexi´on. Use un m´etodo de punto fijo con x0 = 1.5. Resoluci´ on: Los puntos de inflexi´on de la curva y = e−x ln(x) se hallan con la segunda derivada igualando a cero, por tanto: ⇒ y 00 (x) = e−x ln(x) − de lo anterior se obtiene:

e−x e−x e−x − − 2 =0 x x x 2

1

x = e( x + x2 ) ,

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

166

por tanto, el esquema de punto fijo resulta ( x2 +

xn+1 = e

n

1 ) x2 n

Iniciando en x0 = 1.5 obtenemos: x1 = 5.9166935 x2 = 1.4428 ; x3 = 6.46602 x4 = 1.395464; x5 = 7.005 x7 = 1.357780 Para varias iteraciones se tiene: x = 2.0590296 ⇒ x ≈ 2.05. Problema 6: 2 3 Sea el polinomio P3 (x) R 2 = a0 + a1 x + a2 x + a3 x tal que P3 (0) = 1 , P3 (2) = 3 y tambien 0 P3 (x)dx = 4. Hallar el valor de P3 (1)

Resoluci´ on De los datos sobre P3 se obtiene: P3 (0) = 1 = a0 + a1 0 + a2 02 + a3 03 → a0 = 1 P3 (x) = 1 + a1 x + a2 x2 + a3 x3 P3 (2) = 3 = 1 + a1 2 + a2 22 + a3 23 → a1 + 2a2 + 4a3 = 1 Z 2 (1 + a1 x + a2 x2 + a3 x3 ) = 4 0 2

x 2 x3 x4 |0 + a2 |20 + a3 |20 = 4 → 3a1 + 4a2 + 6a3 = 3, 2 3 4 despu´es de resolver el sistema se obtiene: 2 + a1

a1 = 2a3 + 1 a2 = −3a3 Piden: P3 (1) = a0 + a1 + a2 + a3 = 1 + (2a3 + 1) + (−3a3 ) + a3 = 2, → P3 (1) = 2.

´ 6.3. CUARTA PRACTICA CALIFICADA, 2016-1

6.3.

167

Cuarta Pr´ actica Calificada, 2016-1

Problema 1: Conteste Verdadero ´o Falso en cada caso con la justificaci´on. En el caso de la ecuaci´on diferencial se refiere a y 0 = f (x, y) y(x0 ) = y0 dy a. Dado 3 + 5y 2 = sen(x); y(0.3) = 5, el valor de y(0.6) usando dx ME con h = 0.3 es 7.4704. b. Se puede afirmar que el ME aproxima la soluci´on mediante una l´ınea recta. c. El siguiente esquema es un m´ etodo de RK valido:yi+1 = yi + k2 h 1 1 Donde k1 = f (xi , yi ) k2 = f xi + h, yi + k1 h 2 2 d. El m´etodo de integraci´on de Simpson genera una f´ormula de Adams-Bashforth. e. En los m´etodos PC la formula correctora puede ser un M´etodo de Taylor de orden dos. Resoluci´ on : senx − 5y 2 ME: y(x + h) = y(x) + 3   sen(0.3) − 5(5)2 hf (x, y) y(0.3 + 0.3) = y(0.6) = y(0.3) + h y(0.6) = 3 −7.47045

a. Falso. 3y 0 + 5y 2 = senx → y 0 =

b. Verdadero. El m´etodo de Euler: y(x+h)=y(x)+hf(x,y) y(x + h) − y(x) = f (x, y) = y 0 h Es una aproximaci´on lineal. c. Verdadero. RK-2 dice yi+1 = gi + ak1 + bk2 a + b = 1 bα = 1/2 bβ = 1/2 k1 = hf (xi , yi ) k2 = hf (xi + αh, yi + βk1 ) → a = 0 → b = 1rightarrowα = 1/2, β = 1/2 Queda:

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

168

yi+1 = yi + k2 h hk1 h Donde: k1 = f (xi , yi ) k1 = f (xi + , yi + ) 2 2 d. Falso. El m´etodo de Simpson es: b

Z

f (x)dx = I = a

h [f (a) + 2f (a + h) + f (b)] 3

h2 e. Falso. y(x + h) = y(x) + hf (x, y) + y 00 donde: y 00 = fx + fy f No 2 usamos un punto adicional, entonces no hay forma que el m´etodo de Taylor corrija la soluci´on. Problema 2: Si y 00 − 2y 0 + 2y = e2t sen(t) 0 ≤ t ≤ 1 y(0) = −0.4, y 0 (0) = −0.6. Establecerla como un sistema de ecuaciones diferenciales y realice un paso con h = 0.1 usando M RK − 4. Dar todos los detalles de la soluci´on. Resoluci´ on y 00 − 2y 0 + 2y = e2t sen(t) z1 = y z2 = z10 = y 0 z20 = y 00 z10 = z2 → z20 − 2z2 + 2z1 = e2t sent → z20 = 2z2 − 2z1 + e2t sent Hacemos:   z Z= 1 z2     z1 (0) −0.4 = z2 (0) −0.6  0   z1 z2 Z = 0 = 2t = F (t, Z) z2 e sent + 2z2 − 2z1 0

Usando RK4: 1 Zn+1 = Zn + [k1 + 2k2 + 2k3 + k4 ] 6

´ 6.3. CUARTA PRACTICA CALIFICADA, 2016-1

169

Donde: k1 = hF (tn , Zn ) k2 = hF (tn + h/2, Zn +) k3 = hF (tn + h/2, Zn + k2 /2) k4 = hF (tn + h, Zn + k3 ) 1era iteraci´ on: 1 Z(0.1) = Z(0) + (k1 + 2k2 + 2k3 + k4 ) 6      −0.4 −0.06 k1 = 0.1F 0, = −0.6 0.4      −0.35 −0.055 k2 = 0.1F 0.05, = −0.55 −0.034476      −0.4275 −0.61724 k3 = 0.1F 0.05, = −0.617238 −0.032424      −0.41724 −0.063242 k4 = 0.1F 0.05, = −0.632424 0.0880002   −0.64462 Z(0.1) = −0.663957 z1 (0.1) = y1 (0.1) = −0.64462 z2 (0.1) = y10 (0.1) = −0.663957 Problema 3: Una pelota a 1200K se enfr´ıa a una temperatura ambiental de 300K. Asumiendo que el calor es perdido solo debido a la radiaci´on, la ecuaci´on diferencial para la temperatura de la pelota est´a dθ dada por: = −2.2067x10−12 (θ4 − 81x108 ) donde θ esta en K y dt t en segundos. Hallar la temperatura en t=480 segundos usando RK-2 con a=1/2 y asumiendo un paso de h=240 segundos. Resoluci´ on RK-2 con a = 1/2 → b = 1/2, α = 1, β = 1

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

170 1 θn+1 = θn + (k1 + k2 ) 2 k1 = hf (tn , θn )

k2 = hf (tn + h, θn + k1 ) 1era iteraci´ on: 1 θ(240) = θ(0) + (k1 + k2 ) 2 k1 = −1093.91 k2 = 4.22272 θ(240) = 1200 +

−1093.91 + 4.22272 = 655.156 2

2da iteraci´ on: 1 1 θ(480) = θ(240) + (k1 + k2 ) = 655.156 + (k1 + k2 ) 2 2 k1 = −43.2842 k2 = −48.4444 θ(480) = 584.267 Problema 4: Z Hallar el valor de

8

6x3 dx considerando como una ecuaci´on di-

5

ferencial y usando el ME con h = 1.5. Resoluci´ on Z 8 6x3 dx = I = y(8) − y(5) y(x + 8) = y(x) + I donde: y 0 = 6x3 = 5

f (x, y) → ME:

´ 6.3. CUARTA PRACTICA CALIFICADA, 2016-1

171

y(x + h) = y(x) + hf (x, y)

y(5 + 1.5) = y(5) + 1125 y(5 + 1.5 + 1.5) = y(5 + 1.5) + 1.5x6x(5 + 1.5)3 y(8) = y(5) + 1125 + 1.5x6x(6.5)3 y(8) − y(5) = I = 3596.63 Problema 5: Establecer el m´etodo de Taylor de orden 2 para resolver: 2 3y = e−5x , y(0) = 7 asumiendo que h = xi+1 − xi .

dy + dx

Resoluci´ on dy e−5x − 3y 2 = e−5x − 3y → y 0 = dx 2 h2 MT-2: y(x + h) = y(x) + hf (x, y) + y 00 yi+1 = yi + h 2   h2 9ye−5x − 13e−5x Donde: 2 4

e−5x − 3y 2



−1 − 3(7) 2





+

h = xi+1 − xi  1era iteraci´ on: x = 0 → y = 7 y(0 + h) = 7 + h   h2 9.7.1 − 13 2 4

y(h) = 7 − 11h +

25h2 4

+

´ CAP´ITULO 6. PRACTICAS CALIFICADAS

172 Problema 6:

En el siguiente sistema resolver aplicando el M´etodo de Gauss Seidel. Use como punto inicial a (0, 0, 0). Inferir la soluci´on exacta.      1 −3 −1 x1 −4 5 −1 1  x2  =  2  1 1 4 x3 −1 Resoluci´ on Ordenando la matriz para que sea una matriz diagonal dominante: 5x1 − x2 + x3 = 2 x1 − 3x2 − x3 = −4 x1 + x2 + 4x3 = −1 Luego, el esquema iterativo de Gauss-Seidel es: (k)

(k+1) x1

(k+1)

(k+1) x2

=

−x1

(k+1)

=

−x1

(k+1)

x3

(k)

2 + x2 − x3 = 5

(k)

− 4 + x3 −3 (k+1)

− x2 4

1era iteraci´ on: 2 4 1 x1 = x2 = x3 = − 5 3 4 2da iteraci´ on: x1 = 0.716667 x2 = 1.38333 x3 = −0.683333 3era iteraci´ on: dx1 = 0.813333 x2 = 1.8 x3 = −0.77499 Podemos inferir: x1 = 1 x2 = 2 x3 = −1

−1

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. [3] A. Quarteroni, R. Sacco, and F. Saleri, Numerical mathematics, vol. 37. Springer Science & Business Media, 2010. [4] J. Herskovits and S. R. Mazorche, “A feasible directions algorithm for nonlinear complementarity problems and applications in mechanics,” Structural and Multidisciplinary Optimization, vol. 37, no. 5, pp. 435– 446, 2009. [5] A. Ram´ırez, S. Mazorche, J. Herskovits, and G. Chapiro, “An interior point algorithm for mixed complementarity nonlinear problems,” June 2017. Submitted. [6] A. Ram´ırez, S. Mazorche, and G. Chapiro, “Numerical simulation of an in situ combustion model formulated as mixed complementarity problem,” Revista Interdisciplinar de Pesquisa em Engenharia, vol. 2, no. 17, pp. 172–181, 2016. [7] G. Chapiro, A. E. Gutierrez, J. Herskovits, S. R. Mazorche, and W. S. Pereira, “Numerical solution of a class of moving boundary problems with a nonlinear complementarity approach,” Journal of Optimization Theory and Applications, vol. 168, no. 2, pp. 534–550, 2016. [8] A. E. R. Guti´errez, “Aplica¸ca˜o do m´etodo de complementaridade n˜ao linear para o estudo de combust˜ao de oxigˆenio in situ,” 2013.

173