JavierChavarria

2. Resoluci´on num´erica de ecuaciones no lineales Sea f : R → R. Dada la ecuaci´on f(x) = 0 1. ¿Tiene alguna soluci´on? 2. ¿Cu´antas? 3. ¿C´omo las calculamos? 2.1. Existencia y unicidad de soluciones Teorema 2.1 (Bolzano) Sea f : [a, b] → R continua tal que f(a)f(b) < 0. Entonces existe c ∈ (a, b) tal que f(c) = 0. Teorema 2.2 (Rolle) Sea f : [α, β] → R una funci´on continua en [α, β] y derivable en (α, β). Si f(α) = f(β), entonces existe ξ ∈ (α, β) tal que f′(ξ) = 0. Corolario 2.3 Si f′(x) 6= 0 en (a, b), entonces la c del teorema de Bolzano es ´unica. 2.2. M´etodo del punto fijo. Cambiamos el problema. Vamos a resolver: x = g(x) Las soluciones de este problema se llaman puntos fijos de g. Teorema 2.4 Si g ∈ C[a, b] y a ≤ g(x) ≤ b ∀x ∈ [a, b], entonces existe un punto fijo de g: existe c ∈ [a, b] tal que c = g(c). Teorema 2.5 Si adem´as g es derivable en (a, b) y |g′(x)| < 1 ∀x ∈ (a, b), entonces el punto fijo es ´unico. La iteraci´on de punto fijo se define como x0 ∈ R (semilla) xn+1 = g(xn) para n ≥ 0. Generamos as´ı de manera recursiva una sucesi´on (xn). 2 Resoluci´on num´erica de ecuaciones no lineales Teorema 2.6 Si g es continua y l´ım n→∞ xn = c, entonces c = g(c). Ejemplo 2.1 Vamos a resolver la ecuaci´on x = e−x con dos cifras de precisi´on: x0 = = 0.50 x1 = e−0.50 = 0.61 . . . x2 = e−0.61... = 0.55 . . . x3 = e−0.55... = 0.58 . . . x4 = e−0.58... = 0.56 . . . x5 = e−0.56... = 0.57 . . . x6 = e−0.57... = 0.56 . . . x7 = e−0.56... = 0.57 . . . x8 = e−0.57... = 0.57 . . . 0.45 0.5 0.55 0.6 0.65 0.45 0.5 0.55 0.6 0.65 exp(−x) x La soluci´on es c = 0.57 . . .. Nota: Puede que haya punto(s) fijo(s) pero que la iteraci´on de punto fijo no converja. Por ejemplo, para resolver x = x2 que tiene dos puntos fijos c = 0 y c = 1. Si arrancamos en x0 = 2, se tiene que x1 = 4, x2 = 16, x3 = 32, ... Claramente xn → +∞. Definici´on 2.7 Una funci´on g : [a, b] → R se dice contractiva si existe 0 < k < 1 tal que |f(x1) − f(x2)| ≤ k|x1 − x2| para todo x1, x2 ∈ [a, b]. Ejemplo 2.2 a) La funci´on f(x) = x/2 + 3 es contractiva en cualquier intervalo: |f(x1) − f(x2)| = |x1/2 + 3 − (x2/2 + 3)| = |x1/2 − x2/2| = 1 2|x1 − x2|. Luego la constante de contractividad en este caso es k = 1/2. b) La funci´on f(x) = |x|/2 es contractiva en cualquier intervalo. Por las propiedades del valor absoluto tenemos: |f(x1) − f(x2)| = |x1|/2 − |x2|/2 ≤ 1 2|x1 − x2|. Teorema 2.8 Si g : [a, b] → R es contractiva con constante de contractividad 0 < k < 1 en [a, b] y adem´as a ≤ g(x) ≤ b ∀x ∈ [a, b], entonces la iteraci´on de punto fijo converge hacia la ´unica soluci´on c del problema para cualquier semilla x0 ∈ [a, b]. Adem´as se tienen los siguientes resultados sobre el error: |c − xn| ≤ kn|c − x0|, |c − xn| ≤ kn 1 − k |x1 − x0|, |c − xn| ≤ k 1 − k |xn − xn−1| 2.3 M´etodos de aproximaci´on de ra´ıces. 3 En la pr´actica, para comprobar que una funci´on derivable g(x) cumple la condici´on de mapeo a ≤ g(x) ≤ b ∀x ∈ [a, b] no hay que comprobar todos los puntos del intervalo. Basta con hallar todos los puntos xi ∈ (a, b) tales que g′(xi) = 0 y comprobar que a ≤ g(xi) ≤ b, a ≤ g(a) ≤ b y a ≤ g(b) ≤ b. El siguiente resultado nos proporciona una forma de ver si una funci´on es contractiva. Teorema 2.9 Si g ∈ C1[a, b] y existe 0 < k < 1 tal que |g′(x)| ≤ k < 1 ∀x ∈ [a, b], entonces g es contractiva en [a, b]. En la pr´actica, para comprobar si una funci´on derivable 2 veces g(x) es contractiva, habr´a que hallar los puntos pi ∈ (a, b) tales que g′′(pi) = 0 y ver que |g′(pi)| < 1 para todos ellos, |g′(a)| < 1 y |g′(b)| < 1. En este caso se puede tomar k = m´ax{|g′(a)|, |g′(b)|, |g′(pi)|}. 2.3. M´etodos de aproximaci´on de ra´ıces. M´etodo de bisecci´on: Se basa simplemente en el Teorema de Bolzano. Generamos 3 sucesiones (an), (bn) y (cn) de la siguiente manera: a0 = a, b0 = b Para n ≥ 0 • cn+1 = an + bn 2 • si f(cn+1) · f(bn) < 0, tomamos an+1 = cn+1, bn+1 = bn • si f(cn+1) · f(an) < 0, tomamos an+1 = an, bn+1 = cn+1 Por ejemplo: las sucesiones generadas por el m´etodo de bisecci´on para resolver la ecuaci´on x − cos(x) = 0 en [0, π/3] con XTOL = 5 · 10−3, son n a b c fa fb fc n n n+1 ====================================== 0 0.00 1.05 0.52 -1.000 +0.547 -0.342 1 0.52 1.05 0.79 -0.342 +0.547 +0.078 2 0.52 0.79 0.65 -0.342 +0.078 -0.139 3 0.65 0.79 0.72 -0.139 +0.078 -0.032 4 0.72 0.79 0.75 -0.032 +0.078 +0.023 5 0.72 0.75 0.74 -0.032 +0.023 -0.005 6 0.74 0.75 0.74 -0.005 +0.023 +0.009 7 0.74 0.74 0.74 -0.005 +0.009 +0.002 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 4 Resoluci´on num´erica de ecuaciones no lineales Teorema 2.10 Sea f : [a, b] → R continua tal que f(a)f(b) < 0. Entonces la sucesi´on cn generada por el m´etodo de bisecci´on converge hacia alg´un c ∈ (a, b) tal que f(c) = 0. Adem´as |cn − c| ≤ bn − an = b − a 2n . M´etodo de regula falsi: En el m´etodo de bisecci´on no se usa ninguna informaci´on sobre el valor de la funci´on, s´olo sobre su signo. Podemos mejorar eso cambiando la forma de tomar cn+1. Ahora va a ser el punto de corte con el eje OX de la recta que pasa por los puntos (a, f(a)) y (b, f(b)): cn+1 = bn − f(bn) bn − an f(bn) − f(an) Las sucesiones generadas por el m´etodo de regula falsi para x − cos(x) = 0 en [0, π/3] con XTOL = 5 · 10−3 son n a b c fa fb fc n n n+1 ====================================== 0 0.00 1.05 0.68 -1.000 +0.547 -0.103 1 0.68 1.05 0.74 -0.103 +0.547 -0.006 2 0.74 1.05 0.74 -0.006 +0.547 -0.000 Teorema 2.11 Sea f : [a, b] → R continua tal que f(a)f(b) < 0. Entonces la sucesi´on cn generada por el m´etodo de regula falsi converge hacia alg´un c ∈ (a, b) tal que f(c) = 0. Adem´as, existe una constante 0 < k < 1 y n0 ∈ N, tal que |cn − c| ≤ k|cn−1 − c| ∀n ≥ n0. Velocidad de convergencia: Antes de seguir estudiando m´etodos, vamos a dar una definici´on que nos ayude a distinguir los m´etodos “r´apidos” de los m´etodos lentos. Definici´on 2.12 Sea (xn) una sucesi´on tal que xn → x. Diremos que el orden de convergencia de xn hacia x es p > 0 si existe k > 0 y n0 ∈ N tal que |x − xn+1| ≤ k|x − xn|p ∀n ≥ n0 Si p = 1, y 0 < k < 1 la convergencia se dice lineal y k se llama factor de convergencia. Si p = 2, la convergencia se dice cuadr´atica. Si l´ım n→∞ x − xn+1 |x − xn| = 0 la convergencia se dice superlineal. 2.3 M´etodos de aproximaci´on de ra´ıces. 5 En la pr´actica, una de las maneras de comprobar si el orden de convergencia de una sucesi´on es p es calcular el l´ımite l´ım n→∞ |x − xn+1| |x − xn|p = ℓ Si p > 1 y ℓ < ∞, el orden de convergencia es p. Si p = 1 y ℓ = 0, la convergencia es superlineal. Si p = 1 y 0 < ℓ < 1, la convergencia es lineal. Ejemplo 2.3 a) xn = 1/n → x = 0, pero la convergencia no es ni siquiera lineal porque l´ım n→∞ |0 − 1/(n + 1)| |0 − 1/n|1 = l´ım n→∞ 1/(n + 1) 1/n = 1 b) xn = 1/2n → x = 0 linealmente porque l´ım n→∞ 1/2n+1 |1/2n|1 = 1/2 Converger linealmente es converger muy despacio. En cada iteraci´on aumentamos el n´umero de decimales exactos en una cantidad fija que es aproximadamente log10(1/k). Para 1/n, log10(1) = 0, as´ı que no podemos decir nada sobre los decimales exactos ganados en cada iteraci´on. Para 1/2n, log10(2) ≈ 0.3, as´ı que m´as o menos cada 3 iteraciones ganamos 1 decimal exacto. n | 1/n | 1/2^n ----+------+---------- 1 | 1.00| 0.50000 2 | 0.50| 0.25000 3 | 0.33| 0.12500 4 | 0.25| 0.06250 5 | 0.20| 0.03125 6 | 0.17| 0.01563 7 | 0.14| 0.00781 8 | 0.13| 0.00391 9 | 0.11| 0.00195 10 | 0.10| 0.00098 c) El m´etodo de bisecci´on converge, pero la convergencia cn → c no es ni siquiera lineal en general, ya que puede darse el caso de que el error |cn − c| crezca en algunas iteraciones en vez de disminuir. Sin embargo, la convergencia bn − an → 0 s´ı es lineal. d) xn = 1/22n → x = 0 cuadr´aticamente porque l´ım n→∞ |1/22n+1 | (1/22n)2 = l´ım n→∞ 22n2 22n2 = 1 Converger cuadr´aticamente es converger bastante r´apido. Aproximadamente en cada iteraci ´on duplicamos el n´umero de decimales exactos. n | 1/2^(2^n) ----+-------------------- 1 | 0.250000000000000 6 Resoluci´on num´erica de ecuaciones no lineales 2 | 0.062500000000000 3 | 0.003906250000000 4 | 0.000015258789063 5 | 0.000000000232831 6 | 0.000000000000000 M´etodo de la secante: Observa que al hacer regula falsi para resolver x − cos(x) = 0 en [0, π/3], siempre nos hemos apoyado en el punto b original, que est´a lejos de la soluci´on. Para hacer la iteraci´on 2 lo mejor ser´ıa haberse apoyado en los puntos de abscisas 0.68 y 0.74, que son los que mejor informaci´on tienen sobre la soluci´on por ser los ´ultimos que se han encontrado. En eso consiste el m´etodo de la secante: Tomar x0, x1 dos semillas (preferentemente que cumplan la condici´on de Bolzano de cambio de signo) Para n ≥ 1 definimos xn+1 = xn − f(xn) xn − xn−1 f(xn) − f(xn−1) En este caso la sucesi´on generada es: n x fx ================== 0 0.0000 -1.0000 1 1.0472 +0.5472 2 0.6768 -0.1027 3 0.7354 -0.0062 4 0.7391 +0.0001 Como el ejemplo es muy sencillo sale el mismo trabajo que para regula falsi (antes calcul´abamos c1, c2, c3, y aqu´ı x2, x3, x4). Teorema 2.13 Sea δ > 0 y f ∈ C2(c − δ, c + δ), tal que f(c) = 0 y f′(c) 6= 0. Existe ε > 0 tal que si x0, x1 ∈ (c − ε, c + ε), la sucesi´on xn generada por el m´etodo de la secante converge hacia c, y el orden de convergencia es al menos φ = (1 + √5)/2 ≈ 1.618. Adem´as l´ım n→∞ |c − xn+1| |c − xn| = 1 2 f′′(c) f′(c) −1 M´etodo de Newton Cuando en el m´etodo de la secante los puntos consecutivos est´an muy juntos, la recta secante a la curva en esos puntos tiende hacia la tangente, y la pendiente de esta recta es la derivada en el punto. El m´etodo de Newton se escribe pues Tomar x0 una semilla. 2.3 M´etodos de aproximaci´on de ra´ıces. 7 Para n ≥ 0 definimos xn+1 = xn − f(xn) f′(xn) . De nuevo x − cos(x) = 0 con XTOL = 5 · 10−3. n x fx dfx ========================== 0 0.0000 -1.0000 +1.0000 1 1.0000 +0.4597 +1.8415 2 0.7504 +0.0189 +1.6819 3 0.7391 +0.0000 +1.6736 4 0.7391 +0.0000 +1.6736 Teorema 2.14 Sea δ > 0 y f ∈ C2(c − δ, c + δ), tal que f(c) = 0 y f′(c) 6= 0. Existe ε > 0 tal que si x0 ∈ (c−ε, c+ε), la sucesi´on xn generada por el m´etodo de la secante converge hacia c, y el orden de convergencia es al menos cuadr´atico. Adem´as l´ım n→∞ |c − xn+1| |c − xn|2 = 1 2 f′′(c) f′(c) Algunas consideraciones: 1. Newton es un m´etodo de punto fijo para la funci´on g(x) = x − f(x) f′(x) 2. Newton y la secante son m´etodos de car´acter local: convergen bajo ciertas condiciones si la semilla est´a “cerca” de la soluci´on. 3. Newton falla si f′(xn) = 0. Newton y la secante convergen linealmente si f′(c) = 0. 4. Para parar estos algoritmos, se pueden usar varios tests de paradas y combinaciones AND y OR de ellos: Controlar el error absoluto: |xn − xn−1| < XTOL Controlar el error relativo: 2|xn − xn−1| |xn| + |xn−1| < RTOL Buscamos que f(c) = 0, as´ı que un posible criterio es |f(xn)| < FTOL 8 Resoluci´on num´erica de ecuaciones no lineales Comparaci´on de los m´etodos. Si intentamos resolver x−cos(x) = 0 con XTOL = 5 ·10−16, el n´umero de iteraciones es M´etodo iteraciones Punto fijo 89 Bisecci´on 51 Regula Falsi 13 Secante 7 Newton 6 Claramente, Newton y Secante son los m´etodos m´as r´apidos, pero ¿cu´al es mejor?. El orden de convergencia del m´etodo de Newton es mayor (2 > 1.618), pero en cada iteraci´on del m´etodo de la secante hay que evaluar una vez la funci´on, mientras que en cada iteraci´on del m´etodo de Newton hay que evaluar la funci´on y su derivada. Si evaluar f′(x), lleva el mismo trabajo que evaluar f(x), tendremos que cada iteraci´on del m´etodo de Newton lleva el mismo trabajo que 2 iteraciones del m´etodo de la Secante. Pero al hacer 2 iteraciones el m´etodo de la secante, tendremos que |xn+2 − c| ≤ k|xn+1 − c|1.618 ≤ k􀀀k|xn − c|1.6181.618 ≤ ˜k|xn − c|2.618 ya que ((1+√5)/2)2 ≈ 2.618. En este caso, con el mismo trabajo, hemos multiplicado el n´umero de decimales exactos por 2.6, as´ı que es preferible el m´etodo de la secante. En general se tiene la siguiente “regla”: Si (trabajo de evaluar f′(x) ) > 44%(trabajo de evaluar f(x) ) es preferible la secante. En caso contrario es preferible Newton. Combinaci´on de Newton y bisecci´on: Newton es r´apido, pero poco fiable si en alguna iteraci´on f′(xn) ≈ 0, ya que tendremos que dividir por algo cercano a cero, con lo cual aumenta el error de redondeo y el resultado se hace muy grande. Bisecci´on es lento pero seguro. He aqu´ı un algoritmo que combina lo mejor de cada uno de estos m´etodos. a0 = a, b0 = b, c0 = b Para n ≥ 0 • cn+1 = cn − f(cn) f′(cn) • Si cn+1 6∈ (an, bn) ◦ % Newton ha ido mal. Bisecci´on ◦ cn+1 = an + bn 2 • fin si • si f(cn+1) · f(bn) < 0, tomamos an+1 = cn+1, bn+1 = bn • si f(cn+1) · f(an) < 0, tomamos an+1 = an, bn+1 = cn+1 2.4 Ejercicios 9 2.4. Ejercicios Ejercicio 2.1 M´etodo del punto fijo. 1. Demuestra que x = cos(x) tiene una ´unica soluci´on en [0, π/3], que la iteraci´on de punto fijo converge en ese intervalo y calcula el (peor) ritmo de convergencia k. 2. Con una calculadora calcula k32. Sale 0.01. 32 es el n´umero de iteraciones que hay que hacer para tener 2 cifras de precisi´on en el peor de los casos. 3. Calcula la soluci´on con una calculadora tomando como semilla 0.75 hasta obtener 2 cifras de precisi´on. Cuenta el n´umero de iteraciones necesario. (Escribe 0.75 en la calculadora y dale a la tecla del coseno hasta que las dos primeras cifras se repitan.) Ejercicio 2.2 M´etodo de bisecci´on: a) ¿Cu´antas iteraciones de bisecci´on son necesarias para garantizar que el error absoluto cometido es m´as peque˜no que un ε > 0 dado? b) Si sabemos que la aproximaci´on dada por cn tiene k decimales exactos, ¿cu´antas iteraciones hay que hacer para ganar otro decimal? Ejercicio 2.3 Probar que 1/2n2 converge superlinealmente, pero que su orden de convergencia no puede ser ning´un p > 1. Escribe una tabla para n = 1, . . . , 7. Ejercicio 2.4 a) Escribe los algoritmos de bisecci´on, regula falsi, secante y Newton para resolver la ecuaci´on x2 − 2 = 0. b) Ejec´utalos a mano (puedes usar una calculadora) hasta lograr una precisi´on absoluta de 8 decimales. c) Demuestra que el orden de convergencia de la sucesi´on generada en este caso por el m´etodo de Newton es 2. 2.5. Ejercicios complementarios Ejercicio 2.5 Alg´un ejercicio de punto fijo con funci´on m´as complicada Ejercicio 2.6 La orden fzero de Matlab usa una combinaci´on de bisecci´on y secante (en algunos casos tambi´en usa el m´etodo de interpolaci´on cuadr´atica inversa, que no hemos visto este curso). Escribe un algoritmo que combine bisecci´on y secante. 10 Resoluci´on num´erica de ecuaciones no lineales Ejercicio 2.7 (Ning´un m´etodo es el mejor de todos) Sea M > (1 + √5)/2 y sea f(x) = M −Mx si 0 ≤ x ≤ 1 + 1/M −1 si 1 + 1/M ≤ x ≤ M. a) Dibuja la funci´on y calcula la soluci´on de f(x) = 0 en [0,M]. b) ¿Para que valores de x0 es aplicable el m´etodo de Newton?. ¿Cu´antas iteraciones tarda en converger para una precisi´on absoluta εa? c) ¿Se puede aplicar el m´etodo de la secante con x0 = 0, x1 = M? d) ¿Cu´antas iteraciones nb del m´etodo de bisecci´on se necesitan para lograr una precisi´on absoluta dada εa tomando a0 = 0 y b0 = M? e) Aplicamos el m´etodo de Regula Falsi tomando a0 = 0 y b0 = M. Dibuja esquem´aticamente unas cuantas iteraciones suponiendo M >> 1. 1. Probar que si an−1 = 0 y bn−1 > 1 + 1/M, entonces cn = Mn+1/(M + 1)n. 2. Calcular la primera iteraci´on nr en que tiene que cnr < 1 + 1/M. 3. Calcular el error absoluto y el relativo entre dos iteraciones consecutivas con 1 < n ≤ nr: ea = |cn − cn−1|, er = ea/|cn| 4. Explicar qu´e ocurre si para M = 108 tomamos una tolerancia relativa de εr = 10−6. 5. Probar que cnr > 1. Deducir que por tanto anr = 0 y que para este problema el m´etodo de Regula Falsi converge exactamente en nr + 1 iteraciones. f ) Sea ε = 2−52 y M ≥ 20. ¿Es mejor el m´etodo de bisecci´on o el de Regula Falsi?