Source: https://www.slideshare.net/lopeztobal/clculo-races-con-octave
Timestamp: 2018-01-20 17:25:55+00:00

Document:
Método de la Secante by Fabricio Flores 63841 views
Matlab programacion by julio sanchez 46215 views
, En el Paro at Profesor Secundaria
Libro dedicado al cálculo aproximado de raíces de ecuaciones no lineales utilizando Octave: Bisección, Regula, Secante, Pto. Fijo, Newton-Raphson, Wegstein, Müller, Sturm, etc.
Rohatan22
Jose Benedetty , President at trabajo en mi casa haciendo nada xD hahhaha ojala pagarn por eso ¬¬ jaja
1. Resolución De Ecuaciones Métodos Iterativos Cristóbal López Silla - Licenciado en Matemáticas.
2. Índice general Listado de Teoremas, Proposiciones, Corolarios Y Deﬁniciones . . . . . . . . . . . . . . . Listado Ejercicios y Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 8 1. PREFACIO 1.1. Software A Utilizar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Software Utilizado Para Escribir Éste Documento . . . . . . . . . . . . . . . . . . . 1.3. Reﬂexiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 12 13 14 2. MÉTODOS ITERATIVOS PARA RAÍCES REALES 2.1. Método De la Bisección . . . . . . . . . . . . . . . 2.1.1. Introducción . . . . . . . . . . . . . . . . . 2.1.2. Algoritmo Bisección . . . . . . . . . . . . 2.2. Método Regula Falsi . . . . . . . . . . . . . . . . 2.2.1. Introducción . . . . . . . . . . . . . . . . . 2.2.2. Algoritmo Regula-Falsi . . . . . . . . . . 2.3. Método De Punto Fijo . . . . . . . . . . . . . . . 2.3.1. Introducción . . . . . . . . . . . . . . . . . 2.3.2. Convergencia Del Método De Punto Fijo 2.3.3. Algoritmo Del Punto Fijo . . . . . . . . . 2.4. El Método De Wegstein . . . . . . . . . . . . . . . 2.4.1. Introducción . . . . . . . . . . . . . . . . . 2.4.2. Algoritmo De Wegstein . . . . . . . . . . 2.5. Método De Newton . . . . . . . . . . . . . . . . . 2.5.1. Introducción . . . . . . . . . . . . . . . . . 2.5.2. Ejemplo Gráﬁco Del Método De Newton 2.5.3. Algoritmo De Newton . . . . . . . . . . . 2.6. Método De La Secante . . . . . . . . . . . . . . . 2.6.1. Introducción . . . . . . . . . . . . . . . . . 2.6.2. Ejemplo Gráﬁco Método Secante . . . . . 2.6.3. Algoritmo Secante. . . . . . . . . . . . . . 2.7. Métodos Illinois Y Pegasus . . . . . . . . . . . . . 2.7.1. Algoritmo De Illinois . . . . . . . . . . . . 2.7.2. Algoritmo De Pegasus . . . . . . . . . . . 2.8. Método De Steffensen . . . . . . . . . . . . . . . 2.8.1. Introducción . . . . . . . . . . . . . . . . . 15 15 15 15 18 18 19 22 22 22 23 26 26 26 27 27 28 28 31 32 32 34 36 36 38 39 39 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3. CAPÍTULO 0 4 2.8.2. Combinando el método de 2 −Aitken y Punto Fijo Para Steffensen . . . . 2.8.3. Algoritmo De Steffensen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. EJERCICIOS RESUELTOS MÉTODOS ITERATIVOS 4. MÉTODOS PARA POLINOMIOS 4.1. Propiedades De Los Polinomios . . . . . . . 4.2. Método McLaurin De Acotación De Raíces 4.3. Separación De Raíces Reales . . . . . . . . . 4.3.1. Regla De Los Signos De Descartes . 4.3.2. Método De Sturm . . . . . . . . . . . 4.3.3. Método de Sturm . . . . . . . . . . . 4.4. Polinomios Con Coeﬁcientes Racionales . . . . . . . . . . . . . . . . . . . . . . . 40 41 45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 123 125 127 127 129 131 132 5. EJERCICIOS RESUELTOS POLINOMIOS 143 6. CÁLCULO DE RAÍCES COMPLEJAS 6.1. Método De Müller . . . . . . . . . . . . . . . . 6.1.1. Introducción . . . . . . . . . . . . . . . . 6.1.2. Ejemplo Gráﬁco Del Método De Müller 6.1.3. Algoritmo De Müller . . . . . . . . . . . 167 167 167 169 171 Bibliografía . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 RESOLUCIÓN ECUACIONES
4. Índice de ﬁguras 1. 2. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. √ Método Newton f ( x ) = x2 − 2 para aproximar 2 . . . . . . . . . . . . . . . . . . . √ Método Secante f ( x ) = x2 − 2 para aproximar 2 . . . . . . . . . . . . . . . . . . . Gráﬁca polinomio p( x ) = 2x3 − 3x − 1 . . . . . . . . . . . . . . . . Método Bisección p( x ) = 2x3 − 3x − 1 en [1,1.5] . . . . . . . . . . Método Bisección p( x ) = 2x3 − 3x − 1 en [−0.75, − 0.25] . . . . . 2 Función f ( x ) = e− x . . . . . . . . . . . . . . . . . . . . . . . . . . Función f ( x ) = 2 cos( x ) recta y = x . . . . . . . . . . . . . . . . . Función f ( x ) = x − 2 cos( x ) . . . . . . . . . . . . . . . . . . . . . . Tasas De Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . Billar Circular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Función Billar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tasas Convergencia Billar Newton . . . . . . . . . . . . . . . . . . Catenaria x cosh( 1 ) − 0.1 . . . . . . . . . . . . . . . . . . . . . . . x 1 Catenaria . . . . . . . . . . . . . . . . . . . . . . . . . acosh( x+0.1 ) x TCL Catenaria x cosh( 1 ) − 0.1 . . . . . . . . . . . . . . . . . . . . . x 1 TCL Catenaria . . . . . . . . . . . . . . . . . . . . . . acosh( x+0.1 ) x Regula Falsi p( x ) = x10 − 1 en [0.5,1.5] . . . . . . . . . . . . . . . . TC Newton p( x ) = x10 − 1 en [0.5,1.5] . . . . . . . . . . . . . . . . Función Logística f ( x ) = ax (1 − x ) a ∈ [1,4] anatural . . . . . . . TCL Función Logística f ( x ) = 0.5x (1 − x ) . . . . . . . . . . . . . . TCL Función Logística f ( x ) = 3 x (1 − x ) Estimación inicial 0.1 . . 2 TCL Función Logística f ( x ) = 3 x (1 − x ) Estimación inicial 0.375 2 Pto. Fijo Función Logística f ( x ) = 7 x (1 − x ) . . . . . . . . . . . . 2 Errores Pto. Fijo Función Logística f ( x ) = 15 x (1 − x ) . . . . . . . 4 Fractal De Feigenbaum . . . . . . . . . . . . . . . . . . . . . . . . . f ( x ) = 8x − cos x − 2x2 . . . . . . . . . . . . . . . . . . . . . . . . Derivada de g1 ( x ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . Derivada de g3 ( x ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . √ Función De Raabe f ( a) = a ln a − 2a + ln( 2π ) . . . . . . . . . . Problema de la caja . . . . . . . . . . . . . . . . . . . . . . . . . . . Función f ( x ) = x3 − 13x2 + 40x − 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 49 50 51 54 56 58 59 60 63 65 . . . . . . . . . . 66 . . . . . . . . . . 67 . . . . . . . . . . 67 . . . . . . . . . . . . . . . 69 70 73 74 75 75 76 77 78 86 90 91 93 94 95 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 33 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
5. CAPÍTULO 0 6 30. 31. 32. FReynolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 FReynolds2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Savitsky1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 1. 2. Corona Circular Acotación McLaurin . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Polinomio f ( x ) = x4 − 4x3 − x2 + 12x − 6 . . . . . . . . . . . . . . . . . . . . . . . . 139 1. Raíces Unidad De x11 + 1 = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 1. Müller para f ( x ) = 16x4 − 40x3 + 5x2 + 20x + 6 . . . . . . . . . . . . . . . . . . . . 170 RESOLUCIÓN ECUACIONES
6. ECUACIONES CAPÍTULO 0 7 Listado de Teoremas, Proposiciones, Corolarios Y Deﬁniciones 2.1.1. Bolzano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1. Punto Fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5.1. Convergencia De Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.8.1. Convergencia Cuadrática Steffensen . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.8.1. Diferencia Progresiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.0.1. Número Iteraciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.1.1. Ecuación Algebraica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.1.1. Teorema Fundamental Del Álgebra . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.1.1. Galois . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.1.2. División De Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.1.2. Polinomio Divisor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.1.1. El Resto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.1.3. Multiplicidad De Una Raíz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.1.4. MCD Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.2.1. Teorema Acotación McLaurin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 4.2.1. Corolario Acotación McLaurin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 4.2.1. Corona Circular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.3.1. Separación Raíces En Intervalos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.3.1. Cambios De Signo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.3.1. Regla de los signos de Descartes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.3.2. Secuencia De Sturm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.3.1. Multiplicidad Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.3.2. Teorema De Sturm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.3.1. Número Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.4.1. Divisores Término Independiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 4.4.2. Estudio Raíces Racionales En Polinomios . . . . . . . . . . . . . . . . . . . . . . . 133 RESOLUCIÓN ECUACIONES
7. CAPÍTULO 0 8 Listado Ejercicios y Ejemplos 3.0.1. Polinomio Bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.0.2. Coseno Regula-Falsi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.0.3. Coseno Punto Fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.0.4. Raíces Quintas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.0.5. Función Billar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.0.6. Problema Catenaria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.0.7. Raíces Décimas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.0.8. Ecuación Logística . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.0.9. Diodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.0.10. Economía . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.0.11. Coseno Diferentes Métodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.0.12. Identidad De Raabe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.0.13. Problema De La Caja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 3.0.14. Método De Halley . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.0.15. Polinomio Grado Dos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.0.16. Sistema No Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.0.17. Esfera Sumergida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 3.0.18. Problema Partícula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.0.19. Problema Reynolds-Colebrook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.0.20. Método Savitsky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.2.1. Acotación De Raíces De Un Polinomio . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.3.1. Determinamos Cambios De Signo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.3.2. Hallamos Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.3.3. Determinamos Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.3.4. Separamos En Intervalos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 4.4.1. Resolvemos Polinomio De Grado Cuatro . . . . . . . . . . . . . . . . . . . . . . . 133 4.4.2. Otro Más De Hallar Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 4.4.3. Ejemplo Deﬂación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 RESOLUCIÓN ECUACIONES
8. ECUACIONES CAPÍTULO 0 9 5.0.1. Polinomio Grado Tres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.0.2. Polinomio Grado Trece . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 5.0.3. Secuencia Sturm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 5.0.4. Polinomio Grado Nueve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 5.0.5. Polinomio Coeﬁcientes PI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.0.6. Polinomio Grado Cinco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.0.7. Raíces Complejas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.0.8. Acotar Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 5.0.9. Separación De Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 5.0.10. Aproximamos Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 5.0.11. Estudiamos Raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 RESOLUCIÓN ECUACIONES
9. 1 PREFACIO Aquí empieza un proyecto personal de hacer llegar al estudiante universitario la capacidad de aprender Métodos Numéricos utilizando software libre, en nuestro caso he optado por el uso del software Octave. Octave es una alternativa libre y bastante compatible con el software que la inmensa mayoría conocemos, y que me abstengo a mentar. No sólo contamos con Octave, si quieres puedes utilizar como alternativa el software Scilab, de hecho; no creo que los códigos que presento en este trabajo sea muy difícil de poder portarlos a Scilab. El presente texto no es una guía o tutorial de Octave, pienso que los tutoriales llegan a ser demasiado pesados, particularmente aprendo más realizando directamente mis tareas matemáticas. Hoy en día el San Google y algún foro que vaya sobre tus dudas te pueden solucionar la vida, evitando tener que leer eternos manuales con inﬁnidad de comandos. Por descontado, lo principal es tener interés en lo que estes haciendo, con esa premisa se derriban multitud de barreras. Este libro en formato PDF trata sobre diferentes formas y algoritmos para poder calcular las raíces (o ceros) de cualquier ecuación no lineal de una variable. Se divide principalmente en dos partes. Una primera parte contiene los principales métodos iterativos: bisección, regula, punto ﬁjo, Newton-Raphson, etc. Una segunda parte esta dedicada a cómo calcular las raíces de las ecuaciones polinómicas: desde Rufﬁni hasta el método de Sturm. Terminamos con una dedicación especial a cómo hallar raíces complejas aproximadas con el método de Müller. He incluido una gran cantidad de problemas resueltos que he sacado de diferentes libros que tengo en mi posesión y que están debidamente nombrados al ﬁnal del documento en el apartado de bibliografía. A Éste documento ha sido realizado en su integridad bajo software libre, utilizando L TEXcomo código para generar el documento, y distribuciones Gnu/Linux. Todos los códigos realizados aquí los puedes descargar de forma gratuita en el siguiente enlace de GitHub. ENLACE PARA DESCARGAR CÓDIGOS OCTAVE Los códigos son principalmente de guía, un modelo a seguir, el cual tú puedes cambiar a tu gusto o conveniencia. Lo que sí te recomiendo es que no te dediques a copiar y pegar los códigos, porque lo más seguro es que no obtengas experiencia ni con Octave ni realizando el estudio de las 11
10. CAPÍTULO 1 12 raíces. Otra cosa que te recomiendo es que seas crítico con las soluciones aquí expuestas y las que tú mismo/a obtengas, en la justa medida. Realiza las tareas que te dejo pendientes, si un ejercicio está realizado con determinados métodos, puedes probar con otros métodos. 1.1 Software A Utilizar Deﬁnitivamente hay multitud de software que puedes utilizar para calcular las raíces de una ecuación con los métodos propuestos en éste documento. Los que te aconsejo: 1. GNU Octave-> Alternativa libre a Matlab, es el programa principal que se utiliza en este documento. Esta disponible para Mc-OS, Windows, Linux y Android. 2. Octave Source Packages-> Conjunto de paquetes que extienden las capacidades de Octave, altamente recomendable. 3. EPSTK-> Conjunto de comandos extra Octave para una mejor resolución en los gráﬁcos. Viene muy bien para dibujar en coordenadas polares. 4. QtOctave-> Es un GUI basado en QT4 que te permite un uso más amigable de Octave. En Linux lo puedes encontrar en los repositorios de las principales distribuciones. Fue creado hace unos años por un español, que abandonó el proyecto, pero la comunidad Linux sigue dándole soporte porque es muy útil. 5. WxMaxima-> Es un GUI basado en wxWidgets el cual facilita el cálculo simbólico y no simbólico bajo el programa Maxima. Disponible en las 3 principales plataformas. Viene muy bien como apoyo para cálculos simbólicos de raíces y demás en éste documento. Muy recomendable su uso. 6. Sublime Text + SublimeREPL-> Aunque no es código libre sí es muy seguido este editor de textos no enriquecido. Viene muy bien para editar nuestros propios códigos de Octave. El plugin SublimeREPL es un muy buen plugin de Sublime Text que nos permite ejecutar nuestros códigos Octave desde el mismo editor. Viene a ser un muy buen sustituto de QtOctave. 7. Cantor-> Si eres linuxero y usas KDE te puede venir muy bien para trabajar con Maxima el programa Cantor. RESOLUCIÓN ECUACIONES
11. ECUACIONES 1.2 CAPÍTULO 1 13 Software Utilizado Para Escribir Éste Documento A Éste manual ha sido escrito en su totalidad con el lenguaje L TEXutilizando los siguientes editores: LyX TeXMaker TeXStudio Gummi Como documentos de apoyo he utilizado intensamente Edición De Textos Cientíﬁcos Con Latex, y me ha servido de apoyo el foro TEX de Stackexchange. Para dibujar muchas de las gráﬁcas de Octave en éste documento he utilizado Matlab2TikZ. RESOLUCIÓN ECUACIONES
12. CAPÍTULO 1 1.3 14 Reﬂexiones Antes de terminar el prefacio me gustaría que te quedase muy claro lo siguiente, es cierto que con cualquier programa matemático puedes calcular las raíces con tan sólo un comando, sin tener que escribir los algoritmos, pero el objetivo de esta materia en las universidades y en éste documento es que tú conozcas cómo funcionan los diferentes métodos y seas capaz de elegir cuál es el mejor para determinadas ecuaciones no lineales. Además, este tipo de estudio viene muy bien para aﬁanzar conceptos matemáticos, como también a conocer nuevos conceptos. Otra cosa a tener muy en cuenta es que te ayudan a aprender a ser crítico con los resultados que puedas obtener con una calculadora o con un determinado software cientíﬁco. Sin más, te dejo con el contenido de esta especie de manualillo que espero sinceramente que te sea de verdadera ayuda para aprender este campo especíﬁco de las matemáticas aplicadas. RESOLUCIÓN ECUACIONES
13. 2 2.1 MÉTODOS ITERATIVOS PARA RAÍCES REALES Método De la Bisección 2.1.1. Introducción El método de la bisección consiste básicamente en encontrar la raíz de una función en un determinado intervalo. Dicho intervalo lo vamos dividiendo cada paso por la mitad, quedándonos con la mitad que tenga la raíz. Al hacer cada vez más pequeños los intervalos lo que conseguimos es una mejor aproximación a la raíz buscada. El proceso de ir partiendo los intervalos se apoya en el método de la sucesión de intervalos encajados de Cantor. Para cerciorarnos que la raíz se encuentra en un determinado intervalo lo que hacemos en realidad es aplicar el Teorema de Bolzano. Teorema 2.1.1 Bolzano Sea una función cualquiera f : [ a, b] −→ R, tal que f ∈ C ([ a, b]). Entonces: Si f ( a) · f (b) < 0 ⇒ ∃c ∈] a, b[ / f (c) = 0 El Teorema de Bolzano lo que nos dice es que si en un intervalo dado una función continua posee un cambio de signo entonces, podemos asegurar que la función posee al menos una raíz en dicho intervalo, independientemente de la multiplicidad. 2.1.2. Algoritmo Bisección El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 % Metodo Biseccion 15
14. CAPÍTULO 2 2 3 4 5 16 function [] = biseccion(f,a,b,tol,maxiter) filename = 'biseccion.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 6 7 t = a:0.01:b; 8 9 10 fx = inline(f); func = fx(t); 11 12 13 fa = feval (fx,a); fb = feval (fx,b); 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 if fa*fb > 0 disp ('La funcion debe tener signo distinto en los extremos del intervalo') return end iter = 0; x = [ ]; errores = [ ]; imagenes = [ ]; incr = b-a; % warning('off', 'Octave:possible-matlab-short-circuit-operator'); while (incr > tol) && ( iter < maxiter) c = (a+b)/2; x = [x,c]; fc = feval (fx,c); if fc == 0 a = c; b = c; elseif fa*fc < 0 b = c; fb = fc; else a = c; fa = fc; end errores =[errores,incr]; imagenes = [imagenes,fc]; incr = b-a; iter++; end iter2 = 1:iter; sol = c; if incr > tol RESOLUCIÓN ECUACIONES
15. ECUACIONES 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 CAPÍTULO 2 17 disp('Insuficientes iteraciones.') else for k=1:iter fprintf(fid,' %.0ftt %.15ftt %.15fn',iter2(k),x(k),errores(k)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %e', sol,iter,errores(end)); fclose(fid); printf ('La solucion es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de % en', sol,iter,errores(end)); end clf(); subplot(1,3,1); fplot(fx,[a-5,a+5]) set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({f},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,x,imagenes,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('Metodo Biseccion'); legend({f,'Aproximaciones'},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter2,errores,iter2,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en Biseccion'); RESOLUCIÓN ECUACIONES
16. CAPÍTULO 2 90 91 92 93 94 95 18 legend({'Errores','Puntos Errores'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); set (get (gca, 'ylabel'), 'rotation', 360); end biseccion.m En los ejercicios se ve con todo detalle los pasos que se realizan en el método de la bisección, al igual que en los demás métodos. 2.2 Método Regula Falsi 2.2.1. Introducción Éste método trata de mejorar el de la bisección en su velocidad de convergencia, para mejorarlo no tomamos cada vez el punto medio del intervalo [ a,b] sino que se calcula en cada iteración la intersección con el eje de abscisas la recta que pasa por los puntos ( a, f ( a)) y (b, f (b)), siendo f la curva de partida del problema. Consideremos el punto P(c,0) y calculemos la expresión para la recta. Tenemos que el vector director de la recta es: v (b − a, f (b) − f ( a)) De esta forma podemos expresar la ecuación de la recta en su forma continua como: x−a y − f ( a) = b−a f (b) − f ( a) Sustituimos el punto P y despejamos c: c−a − f ( a) − f ( a) · (b − a) − f ( a) · (b − a) = ⇒ c−a = ⇒c= +a⇒ b−a f (b) − f ( a) f (b) − f ( a) f (b) − f ( a) ⇒c= −b · f ( a) + a · f ( a) + a · f (b) − a · f ( a) f (b) − f ( a) c= a · f (b) − b · f ( a) f (b) − f ( a) RESOLUCIÓN ECUACIONES
17. ECUACIONES CAPÍTULO 2 19 Éste c es el que se irá iterando cada vez en el método de la Regula Falsi, es decir, son nuestras aproximaciones para obtener la raíz que se busca de una función o curva f en un intervalo [ a,b]. Para poder iterarlo hemos de darnos cuenta que el valor c es el que hace el papel del punto medio de la bisección, por lo que obtenemos a partir de [ a,b] 2 nuevos subintervalos: [ a,c] y [c,b]. A cada subintervalo le aplicamos el Teorema de Bolzano, y en el subintervalo que se cumpla el teorema le volvemos a aplicar el método de la Regula Falsi, es decir; volvemos a iterar. De esa forma lo que obtenemos es una sucesión ﬁnita de aproximaciones a la raíz que buscamos en un intervalo [ a,b] y que podemos deﬁnirla como {ci }in=0 . Esto es en esencia los pasos que realizamos en éste método, al que hay que añadir los diferentes criterios de parada. 2.2.2. Algoritmo Regula-Falsi El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 function [] = regulafalsi (f,a,b,tol,maxiter) filename = 'regulafalsi.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 5 6 t = a:0.01:b; 7 8 9 fx = inline(f); func = fx(t); 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 fa = feval(fx,a); fb = feval(fx,b); if fa*fb > 0 disp('La funcion debe tener signo distinto en los extremos del intervalo.') return end iter = 0; x = a; errores = [ ]; imagenes = [ ]; incr = b-a; while incr > tol && iter < maxiter c = (a*fb-b*fa)/(fb-fa); x = [x,c]; fc = feval(fx,c); if fc == 0 a = c; b = c; elseif fa*fc < 0 b = c; fb = fc; else RESOLUCIÓN ECUACIONES
18. CAPÍTULO 2 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 20 a=c; fa=fc; end incr = abs(x(end)-x(end-1)); errores =[errores,incr]; iter++; end iter2 = 1:1:iter; fc1 = feval(fx,x); imagenes = [imagenes,fc1]; sol = c; if incr > tol disp('Insuficientes iteraciones.') else for k=1:iter fprintf(fid,' %.0ftt %.15ftt %.15fn',iter2(k),x(k),errores(k)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %e',sol,iter,errores(end)); fclose(fid); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.nCon un error de %en', sol, iter,errores(end)); end 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 clf(); subplot(1,3,1); fplot(fx,[a-5,b+5]) set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({f},'location', 'southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,x,imagenes,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('Metodo Regula Falsi'); legend({f,'Aproximaciones'},'location','southoutside'); RESOLUCIÓN ECUACIONES
19. ECUACIONES 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 CAPÍTULO 2 legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter2,errores,iter2,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en Regula Falsi'); legend({'Errores','Puntos Errores'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); set (get (gca, 'ylabel'), 'rotation', 360); end regulafalsi.m RESOLUCIÓN ECUACIONES 21
20. CAPÍTULO 2 2.3 22 Método De Punto Fijo Deﬁnicion 2.3.1 Punto Fijo Dada una función f : D ⊆ R −→ R y un punto x ∗ ∈ D, decimos que x ∗ es un punto ﬁjo si f ( x ∗ ) = x ∗ , es decir; x ∗ es solución de la ecuación f ( x ) = x ; ∀ x ∈ D. 2.3.1. Introducción En éste método consideramos una función f : D ⊆ R −→ R cualquiera a la que queremos calcular las soluciones de la ecuación f ( x ) = 0. El método consiste en transformar dicha ecuación en una de punto ﬁjo g( x ) = x, de la cual al ir iterándola obtenemos la sucesión de iterados: xk+1 = g( xk ) ∀k ∈ N. En el método ya no partimos de un intervalo inicial [ a,b] para encontrar aproximaciones a la raíz, lo que hacemos es partir de un valor inicial x0 que esté lo más próximo posible a la raíz a calcular. Para determinar éste valor inicial la forma más sencilla es dibujar con Octave la gráﬁca de la función, aunque también podemos utilizar el Teorema de Bolzano para comprobar su existencia, esto lo mejor es hacerlo a mano; o utilizando previamente el método de la Bisección. Una cosa muy importante a tener en cuenta en el método de punto ﬁjo es que la función g no tiene porqué ser necesariamente única, basta tomar como ejemplo la sencilla función f ( x ) = x2 − 2x. Si en f ≡ 0 despejamos x lo podemos hacer de 3 formas diferentes: x2 x2 ⇒ g1 ( x ) = 2 2 √ 2 − 2x = 0 ⇒ x2 = 2x ⇒ g ( x ) = Forma 2 ⇒ x 2x 2 √ Forma 3 ⇒ x2 − 2x = 0 ⇒ x2 = 2x ⇒ g3 ( x ) = − 2x Forma 1 ⇒ x2 − 2x = 0 ⇒ x = A partir de las diferentes funciones de punto ﬁjo nos queda determinar cuál es la mejor para la raíz que deseamos aproximar, tanto en velocidad de converegencia como en cuál es la que está más cerca a la raíz a aproximar. 2.3.2. Convergencia Del Método De Punto Fijo Consideremos la función g : [ a,b] −→ R de clase C 1 ([ a,b]), tal que g( x ) ∈ [ a,b] ∀ x ∈ [ a,b]. Sea la sucesión de punto ﬁjo xk+1 = g( xk ) ∀k ∈ N. RESOLUCIÓN ECUACIONES
21. ECUACIONES CAPÍTULO 2 23 Supongamos también que ∃k ∈ R/| g ( x )| < K < 1 Se cumple que xk+1 −→ x ∗ para cualquier x0 ∈ [ a,b]. Además, el error de la iteración n-ésima por: k→∞ | xn − x ∗ | ≤ Kn | x0 − x1 | n ≥ 1 1−K También se tiene que si g es continua entonces: Si | g ( x ∗ )| > 1 ⇒Los iterados NO convergen a x ∗ Si | g ( x ∗ )| > 1 ⇒Los iterados Convergen Linealmente a x ∗ con tasa de convergencia c L = | g ( x ∗ )| Si g ( x ∗ ) = 0 ⇒Los iterados convergen cuadráticamente a x ∗ 2.3.3. Algoritmo Del Punto Fijo El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 function [] = pfijo(g,f,a,b,x1,tol,maxiter) filename = 'puntofijo.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 5 6 t = a:0.01:b; 7 8 9 10 fx = inline(f); gx = inline(g); func = fx(t); 11 12 13 14 15 16 17 18 19 20 21 22 iter = 1; x = x1; incr = tol+1; errores = incr; imagenes = [ ]; while incr > tol && iter <= maxiter x(iter+1) = feval(gx,x(iter)); incr = abs(x(iter+1)-x(iter)); errores =[errores,incr]; iter++; end 23 24 25 26 iter2 = 1:iter; fc1 = feval(fx,x); imagenes = [imagenes,fc1]; RESOLUCIÓN ECUACIONES
22. CAPÍTULO 2 24 27 28 29 30 31 e = abs(diff(x)); tcl = e(2:end) ./ e(1:end-1); %Tasa de convergencia lineal tcc = e(2:end) ./ e(1:end-1).^2; %Tasa de convergencia cuadratica %tccubica = e(3:end)./e(1:end-1).^3; %Tasa de convergencia cubica 32 33 34 35 36 XR = x(2:end); XL = x(1:end-1); format free; dg = diff(XR)./diff(XL); 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 if incr > tol sol = [ ]; disp('Insuficientes iteraciones'); else sol = x(end); for k=1:iter fprintf(fid,' %.0ftt %.15ftt %.15fn',iter2(k),x(k),errores(k)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %e',sol,iter,errores(end)); %fprintf(fid,'nnTasa De Convergencia Cubican %.15f',tccubica'); fclose(fid); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n', sol, iter); end clf(); figure(1), subplot(1,3,1); fplot(fx,[a-5,b+5]) set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({f},'location', 'southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,x,imagenes,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); RESOLUCIÓN ECUACIONES
23. ECUACIONES 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 CAPÍTULO 2 grid; title('Metodo Punto Fijo'); legend({f,'Aproximaciones'},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter2,errores,iter2,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en el Punto Fijo'); legend({'Errores','Puntos Errores'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); set (get (gca, 'ylabel'), 'rotation', 360); 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 figure(2); subplot(1,2,1); plot(iter2(1:length(tcl)),tcl) set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tasa Convergencia Lineal Punto Fijo'); legend({'T.C.L.'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('tcl'); subplot(1,2,2); plot(iter2(1:length(tcc)),tcc) set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tasa Convergencia Cuadratica Newton'); legend({'T.C.C.'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('tcc'); end RESOLUCIÓN ECUACIONES 25
24. CAPÍTULO 2 26 pﬁjo.m 2.4 El Método De Wegstein 2.4.1. Introducción También conocido como el método del promotor de convergencia de Wegstein. Éste método es una mejora sustancial del método del punto ﬁjo. Consideremos la ecuación a resolver f ( x ) = 0, y su respectiva función de punto ﬁjo g( x ); siendo x0 el valor o estimador inicial. Los pasos a seguir en el método son los siguientes: 1. Al comienzo del paso k-1-ésimo se conoce xk−1 2. Se actualiza xk−1 con el punto ﬁjo: xk = g ( xk−1 ) 3. Se calculan: A= g ( xk ) − xk x k − x k −1 α= 1 1−A 4. Se corrige xk con la sucesión de iterados: x k +1 = x k + α · ( g ( x k ) − x k ) ; ∀ k ∈ N 2.4.2. Algoritmo De Wegstein El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 5 6 7 8 9 function[]=wegstein(f,x0,epsilon,max1) fx = inline(f); for k=1:max1 x1=feval(fx,x0); y1=feval(fx,x1); delta=(y1-x1)/(x1-x0); alpha=1/(1-delta); x0=x1+alpha*(y1-x1); if abs(x1-x0)<epsilon RESOLUCIÓN ECUACIONES
25. ECUACIONES 10 11 12 13 14 15 16 17 CAPÍTULO 2 27 break; endif endfor err=abs(x1-x0)/2; printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n', x0, k); printf('Con un error de %.15en',err); end wegstein.m 2.5 Método De Newton 2.5.1. Introducción El método de Newton consiste en aproximar las soluciones de la ecuación f ( x ) = 0 mediante la recta tangente, por lo que como primera premisa se debe cumplir que la función f sea derivable de orden 1 en un entorno, lo suﬁcientemente grande, de la raíz a aproximar. Partimos siempre de una aproximación inicial de la raíz, x0 , la cual podemos estimar fácilmente dibujando la gráﬁca de la función f . Una vez determinado x0 , en la primera iteración debemos calcular la recta tangente al punto ( x0 , f ( x0 )). Dicha recta tangente se puede calcular a partir de la ecuación de una recta en la forma punto-pendiente: y − y0 = m · ( x − x0 ) Siendo y0 = f ( x0 ). Como se sabe de Bachiller, la pendiente es m = f ( x0 ). Tomamos un punto x1 que es la intersección de la recta tangente con el eje OX, luego y = 0, sustituimos y despejamos x1 : y = f ( x0 ) + f ( x0 ) · ( x − x0 ) ⇒ 0 = f ( x0 ) + f ( x0 ) · ( x1 − x0 ) ⇒ − f ( x0 ) = f ( x0 ) · ( x1 − x0 ) ⇒ − f ( x0 ) f ( x0 ) x1 − x0 = ⇒ x1 = x0 − f ( x0 ) f ( x0 ) Si iteramos éste procedimiento obtenemos la sucesión de iterados del método de Newton: x k +1 = x k − RESOLUCIÓN ECUACIONES f ( xk ) ; ∀k ∈ N f ( xk )
26. CAPÍTULO 2 28 Al ver la expresión de la sucesión del método observamos que también se debe cumplir en cada iteración que f ( xk ) = 0 ∀ k ∈ N o próximo a cero en un entorno pequeño, ya que en caso contrario dividiríamos por 0, provocando error. Podemos garantizar la convergencia cuadrática del método mediante el uso del siguiente teorema: Teorema 2.5.1 Convergencia De Newton Sea f ∈ C 2 ([ a, b]), consideremos x ∗ ∈ [ a, b] una raíz de la ecuación f ( x ) = 0, tal que f ( x ∗ ) = 0 (i.e. ; es una raíz simple) ⇒ xk −→ x ∗ para cualquier aproximación inicial x0 , con velocidad de convergencia cuak→∞ drática, O(n2 ). Si f ( x ∗ ) = 0 ⇒ El método de Newton converge linealmente, O(n). 2.5.2. Ejemplo Gráﬁco Del Método De Newton √ Consideremos la función f ( x ) = x2 − 2 con ecuación x2 − 2 = 0, sus soluciones son ± 2. Apro√ ximemos de las dos raíces x ∗ = 2 con el método de Newton. Obviamente f ∈ C 2 ([1,1.5]), con √ √ f ( 2) = 2 2 = 0, luego se cumplen las condiciones del teorema anterior, así que la convergencia es cuadrática. Tomemos como valor inicial x0 = 1.25. En la ﬁgura adjunta puedes ver las primeras 2 iteraciones para éste ejemplo, cabe destacar la rápidez del método, con tan sólo 2 iteraciones obtenemos una buena aproximación de la raíz, ya que las 4 primeras cifras decimales son exactas. 2.5.3. Algoritmo De Newton El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 function [] = newton(f,df,x1,a,b,tol,maxiter) 2 3 4 5 filename = 'newtonraphson.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 6 7 8 9 10 tic; fx = inline(f); dfx = inline(df); t = x1-5:0.01:x1+5; RESOLUCIÓN ECUACIONES
27. ECUACIONES CAPÍTULO 2 Método Newton f ( x ) = x2 − 2 29 y = x2 − 2 4 y = 2.5x − 3.5625 y = 2.85x − 4.030625 2 −4 −2 2 4 −2 −4 f ( x ) = x2 − 2 f ( x ) = 2x y = yk + f ( xk )( x − xk ) Iteración 1 x0 = 1.25 y0 = f (1.25) = −0.4375 f (1.25) = 2.5 y = 2.5x − 3.5625 Si y = 0 → x1 = 1.425 Iteración 2 x1 = 1.425 y1 = f (1.425) = 0.030625 f (1.425) = 2.85 y = 2.85x − 4.030625 Si y = 0 → x2 = 1.414254386 Con sólo 2 iteraciones obtenemos una aprox. de 4 cifras decimales exactas. Figura 1: Método Newton f ( x ) = x2 − 2 para aproximar 11 12 func = fx(t); 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 k = 1; x = x1; incr = tol+1; sol = [ ]; errores = incr; imagenes = [ ]; while incr > tol && k <= maxiter delta = -feval(fx,x(k))/feval(dfx,x(k)); x(k+1) = x(k) + delta; incr = abs(delta); errores =[errores,incr]; sol = [sol,x]; k = k+1; end 28 29 30 31 32 33 display('El algoritmo tiene un coste en tiempo de: ') toc; iter2 = 1:k; fc1 = feval(fx,sol); imagenes = [imagenes,fc1]; 34 35 e = abs(diff(x)); RESOLUCIÓN ECUACIONES √ 2
28. CAPÍTULO 2 36 37 30 tcl = e(2:end) ./ e(1:end-1); %Tasa de convergencia lineal tcc = e(2:end) ./ e(1:end-1).^2; %Tasa de convergencia cuadratica 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 if incr > tol sol = [ ]; disp('Diverge o insuficientes iteraciones'); else sol = x(end); for m=1:k fprintf(fid,' %.0ftt %.15ftt %.15fn',iter2(m),x(m),errores(m)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %enn',sol(end),k,errores(end)); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n', sol,k); end fprintf(fid,'Tasa de convergencia lineal = %f n',tcl); fprintf(fid,'Tasa de convergencia cuadratica = %f n',tcc); fclose(fid); 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 clf(); figure(1); subplot(1,3,1); fplot(fx,[a,b]); set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({f},'location', 'southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,sol,imagenes,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('Metodo Newton Raphson'); legend({f,'Aproximaciones'},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); RESOLUCIÓN ECUACIONES
29. ECUACIONES CAPÍTULO 2 set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter2,errores,iter2,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en Newton'); legend({'Errores', 'Puntos Errores'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); 80 81 82 83 84 85 86 87 88 89 90 91 92 figure(2); subplot(1,2,1); plot(iter2(1:length(tcl)),tcl) set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tasa Convergencia Lineal Newton'); legend({'T.C.L.'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('tcl'); subplot(1,2,2); plot(iter2(1:length(tcc)),tcc) set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tasa Convergencia Cuadratica Newton'); legend({'T.C.C.'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('tcc'); 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 end newton.m 2.6 Método De La Secante RESOLUCIÓN ECUACIONES 31
30. CAPÍTULO 2 32 2.6.1. Introducción El método de la secante tiene una velocidad de convergencia mayor que el de la bisección o de la regula, pero menor que el de Newton. Es un buen método para cuando el de Newton o el del punto ﬁjo fallan, o cuando el cálculo de la primera derivada es muy costoso. Para aplicarlo se toman dos estimaciones iniciales x1 y x2 de la raíz o solución de la ecuación f ( x ) = 0, siendo f una curva cualquiera, a poder ser continua en un entorno de la solución o soluciones de la ecuación. En el método lo primero es calcular una recta secante a la curva/función f , es decir; una recta que pase por los puntos ( x1 , f ( x1 )) y ( x2 , f ( x2 )), luego se toma como siguiente estimación un x3 que es la intersección de la secante con el eje OX. Si iteramos el proceso obtenemos como sucesión de iterados del método la expresión: x k +1 = x k − f ( xk ) ∀k ∈ N f ( x k ) − f ( x k −1 ) x k − x k −1 También puedes encontrar la sucesión expresada de forma equivalente como: x k +1 = x k − x k − x k −1 f ( xk ) ∀ k ∈ N f ( x k ) − f ( x k −1 ) √ 1+ 5 La velocidad de convergencia del método es φ = , es decir; su velocidad de convergencia 2 es el número áureo, es una convergencia superlineal. 2.6.2. Ejemplo Gráﬁco Método Secante √ Consideremos la función f ( x ) = x2 − 2 con ecuación x2 − 2 = 0, sus soluciones son ± 2. Aproxi√ memos de las dos raíces x ∗ = 2 con el método de la Secante. Tomamos como puntos de partida x1 = 1 y x2 = 2. Para calcular las expresiones de las rectas secantes en cada iteración vamos a utilizar la expresión: y = f ( x1 ) + x − x1 [ f ( x2 ) − f ( x1 )] x2 − x1 Las sucesivas aproximaciones a la raíz las obtendremos sustituyendo en la expresión anterior y = 0 y despejando x. x−1 3 ⇒ y = 3x − 4 =⇒ x3 = 4 = 1.333333333 . . . 3 1 4 2 x − 3 20 10 14 x2 = 2 =⇒ y = − + ⇒ y = x− =⇒ x4 = 7 = 1.4 5 2 9 9 3 3 3 x1 = 1 x2 = 2 =⇒ y = −1 + x1 = 4 3 RESOLUCIÓN ECUACIONES
31. ECUACIONES CAPÍTULO 2 33 7 x − 51 1 17 24 5 x2 = 2 =⇒ y = − + ⇒ y = x− =⇒ x5 = 3 25 25 5 5 5 x1 = 7 5 x1 = 24 17 24 x− 2 17 580 ⇒ y = 58 x − 82 =⇒ x = x2 = 2 =⇒ y = − + 6 10 2899 289 17 17 17 41 29 41 x− 1 29 1683 ⇒ y = 99 x − 140 =⇒ x = x2 = 2 =⇒ y = − + 7 17 841 841 29 29 29 x1 = 24 17 = 1.411764706 41 29 = 1.413793103 140 99 = 1.41414141 . . . Luego después de 5 iteraciones obtenemos la aproximación x ∗ = 1.4141414141 . . . , es decir, obtenemos 3 cifras decimales exactas. En la ﬁgura adjunta puedes ver las primeras 5 iteraciones para éste ejemplo. y = x2 − 2 Método Secante f ( x ) = x2 −2 y = 3x − 4 y = 10 x − 14 3 3 y = 17 x − 24 5 5 y = 58 x − 82 17 17 y = 99 x − 140 29 29 4 2 0.5 1 2 1.5 −2 Figura 2: Método Secante f ( x ) = x2 − 2 para aproximar RESOLUCIÓN ECUACIONES √ 2
32. CAPÍTULO 2 34 2.6.3. Algoritmo Secante. El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 function [] = secante(f,x1,x2,tol,maxiter) filename = 'secante.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 5 6 7 8 9 10 t = x1:0.01:x2; fx = inline(f); x = [x1,x2]; y = feval(fx,x); func = fx(t); 11 12 k = 2; 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 sol = [ ]; errores = [ ]; incr = tol + 1; delta = x2 - x1; while incr > tol && k < maxiter m = (y(k) - y(k-1))/delta; delta = -y(k)/m; x(k+1) = x(k) + delta; y(k+1) = feval(fx,x(k+1)); sol = [sol,x]; incr = abs(delta); errores =[errores,incr]; k++; end 28 29 iter = 1:k-2; 30 31 32 33 34 35 36 37 38 39 40 if incr > tol sol = [ ]; disp('Diverge o insuficientes iteraciones.') else sol = x(end); for m=1:k-2 fprintf(fid,' %.0ftt %.15ftt %.15fn',iter(m),x(m),errores(m)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %e',sol(end),k-2,errores(end)); fclose(fid); RESOLUCIÓN ECUACIONES
33. ECUACIONES 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 CAPÍTULO 2 printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n', sol,k-2); end clf(); subplot(1,3,1); fplot(fx,[x1-5,x1+5]) set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({f},'location', 'southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,x,y,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('Metodo Secante'); legend({f,'Aproximaciones'},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter,errores,iter,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en Secante'); legend('Errores','location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); set (get (gca, 'ylabel'), 'rotation', 360); end secante.m RESOLUCIÓN ECUACIONES 35
34. CAPÍTULO 2 2.7 36 Métodos Illinois Y Pegasus Ambos métodos son generalizaciones del método de la Regula Falsi. Utilizan un paso de secante modiﬁcado para acelerar la convergencia lineal de dicho método. Su algoritmo asociado comparte los dos primeros pasos con el de la Regula Falsi. 1) Se parte de un intervalo ] a, b[ que contiene a la raíz x ∗ . 2) Se calculan yb = f (b) y yi = f ( xi ), tal que yb · yi < 0 xi · y b − b · yi abscisa del punto de intersección con el eje OX de la recta secante y b − yi que pasa por los puntos ( xi , yi ) y (b, yb ). 3) Se deﬁne xk = 4) Se evalúa P := yb · yk  b = x , k Si P > 0 =⇒ y = y , b k  x = x , i k Si P < 0 =⇒ y = y , i k xi = x i yi = α · yi b=b yb = α · yb En donde se tiene: 1 Método Illinois =⇒ α = 2  yb α = ,  yb + yk Método Pegasus=⇒ yi α =  , yi + y k P>0 P<0 Ambos métodos tienen convergencia superlineal α ≈ 1.642. 2.7.1. Algoritmo De Illinois El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 5 6 function []=illinois(f,a,b,delta,epsilon,max1) %a,b extremos del intervalo que enmarca la raiz %a= %b= %delta es la tolerancia para la raiz y epsilon para el valor de la funcion %delta= RESOLUCIÓN ECUACIONES
35. ECUACIONES 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 CAPÍTULO 2 %epsilon= %max1= % max1 es el numero maximo de tolerancias fx = inline(f); ya=feval(fx,a); yb=feval(fx,b); %Evalua la funcion en los extremos del intervalo donde presuntamente se %encuentra la raiz. if ya*yb>0 disp('¡Ojo!: f(a)*f(b)>0') break, end %Aqui no hay raiz en el intervalo en estudio. for k=1:max1 dx=yb*(b-a)/(yb-ya); c=b-dx; %Se define la abcisa del punto c interseccion de la secante por %(a,ya),(b,yb) con el eje x ac=c-a; yc=feval(fx,c); %y se evalua su ordenada sobre la curva. if yc==0 break; %ya tenemos la solucion elseif yb*yc>0 %una vez que sabemos que no es la raiz debemos saber en que intervalo (a,c) %o (c,b) esta la raiz buscada, luego estudiamos los signos de yb*yc y de %ya*yc b=c; ya=(0.5)*ya; yb=yc; %en este caso el nuevo intervalo (a,c) y multiplicamos la ya por 0.5 else a=c; ya=yc; yb=(0.5)*yb; end %aqui el nuevo intervalo es (c,b) dx=min(abs(dx),ac); if abs(dx)<delta break; end if abs(yc)<epsilon break; end RESOLUCIÓN ECUACIONES 37
36. CAPÍTULO 2 52 53 54 55 56 57 58 59 38 end err=abs(b-a)/2; yc=feval(fx,c); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n',c, k); printf ('Con un error de %en',err); printf('Con f(raiz)= %en',yc); clear; end illinois.m 2.7.2. Algoritmo De Pegasus El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 function []=pegasus(f,a,b,delta,epsilon,max1) fx = inline(f); ya=feval(fx,a); yb=feval(fx,b); %Evalua la funcion en los extremos del intervalo donde presuntamente se %encuentra la raiz. if ya*yb>0 disp('¡Ojo!: f(a)*f(b)>0') break, end %Aqui no hay raiz en el intervalo en estudio. for k=1:max1 dx=yb*(b-a)/(yb-ya); c=b-dx; alpha=b/(b+c); printf ('alpha = %.15fn',alpha); %Se define la abcisa del punto c interseccion de la secante por %(a,ya),(b,yb) con el eje x ac=c-a; yc=feval(fx,c); %y se evalua su ordenada sobre la curva. if yc==0, break; %ya tenemos la solucion elseif yb*yc>0 %una vez que sabemos que no es la raiz debemos saber en que intervalo (a,c) %o (c,b) esta la raiz buscada, luego estudiamos los signos de yb*yc y de %ya*yc b=c; ya=(alpha)*ya; RESOLUCIÓN ECUACIONES
37. ECUACIONES 30 31 32 33 34 35 36 37 38 39 40 41 CAPÍTULO 2 39 yb=yc; %en este caso el nuevo intervalo (a,c) y multiplicamos la ya por 0.5 else a=c; ya=yc; yb=(alpha)*yb; end %aqui el nuevo intervalo es (c,b) dx=min(abs(dx),ac); if abs(dx)<delta,break,end if abs(yc)<epsilon,break,end end 42 43 44 45 46 47 48 49 err=abs(b-a)/2; yc=feval(fx,c); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n',c, k); printf ('Con un error de %en',err); printf('Con f(raiz)= %en',yc); clear; end pegasus.m 2.8 Método De Steffensen 2.8.1. Introducción El método de Steffensen es un método con velocidad de convergencia cuadrática O(n2 ), al igual que en Newton. Se basa en la aceleración de Aitken y el método de punto ﬁjo. El método lo iniciamos a partir de un valor inicial x0 que esté lo suﬁcientemente cerca de la raíz a aproximar, si no es así lo que puede suceder es que encontremos la aproximación de la raíz muy lentamente, perdiendo su velocidad de convergencia cuadrática; o que el método diverja. Tiene la ventaja de que la función sólo debe ser continua en un entorno de la raíz, ya que no precisa de calcular ninguna derivada de la función. Por contra, tiene el inconveniente de que en cada iteración debemos calcular la imagen de la función 2 veces, lo cual puede aumentar su costo si la expresión algebraica de la función es complicada. La sucesión de iteraciones del método es: RESOLUCIÓN ECUACIONES
38. CAPÍTULO 2 40 ( f ( x n ) − x n )2 ∀n ∈ N f ( f ( xn )) − 2 f ( xn ) + xn x n +1 = x n − Cuando decimos que debemos calcular dos veces el valor de la función nos referimos cuando calculamos f ( f ( xn )) = ( f ◦ f )( xn ) = f 2 ( xn ) ∀ n ∈ N Para asegurar que el método tiene convergencia cuadrática podemos hacerlo comprobando que la raíz a aproximar tenga multiplicidad simple, es decir, que no se anule en su primera derivada de la función. Esto se puede traducir en el siguiente teorema: Teorema 2.8.1 Convergencia Cuadrática Steffensen Sea x = g( x ) / g( p) = 0, con g ( p) = 1. Si ∃δ > 0 / g ∈ C 3 ([ p − δ, p + δ]) ⇒ El método de Steffensen es de orden O(n2 ) ∀ p0 ∈ [ p − δ, p + δ] 2.8.2. Combinando el método de 2 −Aitken y Punto Fijo Para Steffensen ˆ n= Para obtener Steffensen combinamos ambos métodos. El de 2 −Aitken nos dice que si { xn }∞ 1 es una sucesión que converge más rápidamente que una sucesión { xn }∞ 1 , se tiene que: n= ˆ xn = xn − Deﬁnicion 2.8.1 ( x n +1 − x n )2 ∀n ∈ N xn+2 − 2xn+1 + xn Diferencia Progresiva Dadada la sucesión { xn }∞ 1 se deﬁne la Diferencia Progresiva n= xn , mediante: x n = x n +1 − x n ∀ n ≥ 0 Las potencias de orden k se deﬁnen recursivamente como: k xn = k −1 ( xn ) ∀n ≥ 0 ∀k ≥ 2 Así pues, por la deﬁnición anterior tenemos lo siguiente: RESOLUCIÓN ECUACIONES
39. ECUACIONES CAPÍTULO 2 2 41 = ( x n +1 − x n ) = xn x n +1 − xn = ( x n +2 − x n +1 ) − ( x n +1 − x n ) = xn+2 − 2xn+1 + xn Luego la sucesión del método 2 −Aitken se puede reescribir como: ˆ xn = xn − Sobre el método siguiente modo: 2 −Aitken ( x n )2 ∀n ∈ N 2x n aplicamos el método de punto ﬁjo con función de punto ﬁjo g del ˆ x0 , x1 = g ( x0 ), x2 = g ( x1 ), x0 = { 2 } x0 , x3 = g ( x2 ), ˆ 1 = { x 2 } x1 , . . . , De esa forma es como actúa el algoritmo de Steffensen, cada 3 pasos en vez de aplicar el punto ﬁjo aplica la aceleración de Aitken para acelerar la convergencia. 2.8.3. Algoritmo De Steffensen El algoritmo que se ha realizado bajo el programa Octave para éste documento es el siguiente: 1 2 3 4 function []=steffensen(fx,x0,a,b,tolx,nmax) filename = 'steffensen.txt'; fid = fopen(filename,'w'); fprintf(fid,'kttx (k)tttttttErrorn'); 5 6 7 8 9 10 11 12 13 14 15 f = inline(fx); t = a:0.01:b; func = f(t); err=tolx+1; x=x0; phi=0; iter = 0; errores = [ ]; sol = [ ]; imagenes = [ ]; 16 17 while iter < nmax && err > tolx RESOLUCIÓN ECUACIONES
40. CAPÍTULO 2 xx=x; fxk=feval(f,x); tolf=tolx*abs(phi); if abs(fxk)<=tolf break end fxk2=feval(f,x+fxk); phi=(fxk2-fxk)/fxk; x=xx-fxk/phi; err=abs(x-xx); errores =[errores,err]; sol = [sol,x]; iter++; 18 19 20 21 22 23 24 25 26 27 28 29 30 31 42 end 32 33 34 35 iter2 = 1:1:iter; fc1 = feval(f,sol); imagenes = [imagenes,fc1]; 36 37 38 39 40 41 42 43 for k=1:iter fprintf(fid,' %.0ftt %.15ftt %.15fn',iter2(k),sol(k),errores(k)); end fprintf (fid,'nLa ósolucin es %.15f.nSe ha alcanzado en %d iteraciones.nCon un error de %e',sol(end),iter,errores(end)); fclose(fid); printf ('La solucion es %.15f, y se ha alcanzado en %d iteraciones.n', sol(end),iter); 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 clf(); subplot(1,3,1); fplot(f,[a,b]) set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); set(gca, 'box', 'off'); grid; title('La Funcion'); legend({fx},'location', 'southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,2); plot(t,func,x,imagenes,'*r') set (gca, 'xaxislocation', 'zero'); set (gca, 'yaxislocation', 'zero'); RESOLUCIÓN ECUACIONES
41. ECUACIONES set(gca, 'box', 'off'); grid; title('Metodo Steffensen'); legend({fx,'Aproximaciones'},'location','southoutside'); legend boxon; xlabel('y'); ylabel('x'); set (get (gca, 'ylabel'), 'rotation', 360); subplot(1,3,3); plot(iter2,errores,iter2,errores,'*r') set (gca, 'xaxislocation', 'bottom'); set (gca, 'yaxislocation', 'right'); set(gca, 'box', 'off'); grid; title('Tendencia de errores en Steffensen'); legend({'Errores', 'Puntos Errores'},'location','southoutside'); legend boxon; xlabel('Numero Iteraciones'); ylabel('Errores'); set (get (gca, 'ylabel'), 'rotation', 360); 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 CAPÍTULO 2 end steffensen.m RESOLUCIÓN ECUACIONES 43
42. 3 Ejercicio 3.0.1 EJERCICIOS RESUELTOS MÉTODOS ITERATIVOS Polinomio Bisección Determinar las raíces del polinomio p( x ) = 2x3 − 3x − 1 mediante el método de bisección mostrando el proceso paso a paso. Resolución Lo primero es determinar aproximadamente qué puntos son los que cortan al eje de abcisas. Una forma sencilla es dibujando la gráﬁca del polinomio. En Octave lo podemos hacer con los siguientes comandos: 1 2 3 4 5 x = linspace(-2,2); p = [2 0 -3 -1]; y = polyval(p,x); plot(x,y) grid on La gráﬁca de la función es la siguiente: En dicha gráﬁca podemos apreciar que el polinomio posee 3 raíces simples, una es claramente -1, otra esta cerca de -0.5, y la tercera se encuentra cerca de 1.5. De hecho podemos determinar algebraicamente dichas soluciones. Aplicando Rufﬁni con raíz -1 obtenemos el polinomio de grado 2, q( x ) = 2x2 − 2x − 1, si lo resolvemos por la fórmula general de segundo grado obtenemos las √ raíces: x = 1±2 3 45
43. CAPÍTULO 3 46 10 5 0 −5 −10 −15 −2 −1 0 1 2 Figura 1: Gráﬁca polinomio p( x ) = 2x3 − 3x − 1 0 −3 −1 −2 2 1 2 −2 −1 0 2 −1 Ambas raíces son números irracionales (inﬁnitas cifras decimales), eso lo sabemos porque 3 es primo y la raíz cuadrada de un número primo es siempre irracional. Luego hemos obtenido de forma algebraica las tres raíces del polinomio, con una simple calculadora podemos obtener una aproximación numérica de las dos soluciones irracionales. x1 = √ 1− 3 2 = −0.3660254038 x2 = √ 1+ 3 2 = 1.366025404 Vamos a aplicar el método de la bisección 1 de forma desglosada para la obtención de la raíz -1. Lo primero es establecer el intervalo de partida, en nuestro caso va a ser [ a,b] = [−1.5, − 0.5] Calculamos p( a) · p(b) y comprobamos si es negativo (i.e.; si hay cambio de signo, Teorema de Bolzano) como primer paso. Después calculamos el punto medio de [ a,b]. p( a) · p(b) = p(−1.5) · p(−0.5) = −3.25 × 0.25 = −0.8125 < 0 →Hay cambio de signo, podemos proseguir. m1 = a+b 2 = −1.5−0.5 2 = −1 Con esto establecemos 2 intervalos que son: [ a,m1 ] = [−1.5, − 1] [m1 ,b] = [−1, − 0.5] El siguiente paso es ver en cuál de dichos subintervalos hay un cambio de signo (la raíz está en el intervalo susodicho), nos quedaremos con el que tenga el cambio de signo: 1 El método se compone de la aplicación de 2 teoremas, el de Bolzano para comprobar si hay cambio de signo en cada subintervalo. Y el teorema de los Intervalos Encajados de Cantor, con el cual vamos dividiendo cada subintervalo en su mitad, creando así una sucesión de intervalos encajados monótona que converge. RESOLUCIÓN ECUACIONES
44. ECUACIONES CAPÍTULO 3 47 p( a) · p(m1 ) = p(−1.5) · p(−1) = −3.25 × 0 = 0 Pero aquí nos detenemos ya porque hemos encontrado la raíz buscada, nos da como solución √ -1; 1+ 3 exacta y con una iteración solamente. Volvamos ha hacer lo mismo pero para la raíz x2 = 2 . Tomamos el intervalo de partida [ a,b] = [1,1.5] y procedemos con el algoritmo de la bisección hasta que obtengamos 2 cifras decimales exactas respecto a nuestra solución exacta: p( a) · p(b) = p(1) · p(1.5) = −2 × 1.25 = −2.5 < 0 → Hay cambio de signo. m1 = a+b 2 = 1+1.5 2 = 1.25 → → [ a,m1 ] = [1,1.25] [m1 ,b] = [1.25,1.5] p( a) · p(m1 ) = p(1) · p(1.25) = −2 × (−0.84375) = 1.6875 > 0 →Aquí no se encuentra la raíz. p(m1 ) · p(b) = p(1.25) · p(1.5) = −0.84375 × 1.25 = −1.0546875 < 0 →Hay cambio de signo m2 = m1 + b 2 = 1.25+1.5 2 = 1.375 → → [m1 ,m2 ] = [1.25,1.375] [m2 ,b] = [1.375,1.5] p(m1 ) · p(m2 ) = p(1.25) · p(1.375) = −0.84375 × 0.07421875 = −0.06262207031 < 0 →Hay cambio de signo. m3 = m1 + m2 2 = 1.25+1.375 2 = 1.3125 → → [m1 ,m3 ] = [1.25,1.3125] [m3 ,m2 ] = [1.3125,1.375] p(m1 ) · p(m3 ) = p(1.25) · p(1.3125) = −0.84375 × (−0.4155273438) > 0 →Aquí no se encuentra la raíz. p(m3 ) · p(m2 ) = p(1.3125) · p(1.375) = −0.4155273438 × 0.07421875 < 0 →Hay cambio de signo. m4 = m3 + m2 2 = 1.3125+1.375 2 = 1.34375 → RESOLUCIÓN ECUACIONES
45. CAPÍTULO 3 → [m3 ,m4 ] = [1.3125,1.34375] 48 [m4 ,m2 ] = [1.34375,1.375] p(m3 ) · p(m4 ) = p(1.3125) · p(1.34375) = −0.4155273438 × (−0.178527832) > 0 →Aquí no se encuentra la raíz. p(m4 ) · p(m2 ) = p(1.34375) · p(1.375) = −0.178527832 × 0.07421875 < 0 →Hay cambio de signo. m5 = m4 + m2 2 = 1.34375+1.375 2 = 1.359375 → → [m4 ,m5 ] = [1.34375,1.359375] [m5 ,m2 ] = [1.359375,1.375] p(m4 ) · p(m5 ) = p(1.34375) · p(1.359375) = −0.178527832 × (−0.05414581302) > 0 →Aquí no hay raíz. p(m5 ) · p(m2 ) = p(1.359375) · p(1.375) = −0.05414581302 × 0.07421875 < 0 →Hay cambio de signo. m6 = m5 + m2 2 = 1.359375+1.375 2 = 1.3671875 √ Aquí nos detenemos porque si comparamos m6 = 1.3671875 con x2 = 1+2 3 = 1.366025404 observamos que las 2 primeras cifras decimales coinciden, como queríamos. Lo que hemos de tener en cuenta es que aunque el método de la bisección siempre converge 2 , lo hace generalmente muy lento. Hemos necesitado 6 iteraciones para tan sólo conseguir 2 cifras decimales exactas respecto a la solución algebraica exacta. Luego nuestra solución con el método de la bisección después de 6 iteraciones es x2 ≈ m6 = 1.3671875. Si queremos una mejor aproximación deberíamos seguir iterando el método de la bisección. √ De la misma forma que se ha procedido con x2 lo puedes hacer con x1 = 1−2 3 , te aconsejo que lo hagas porque es una buena forma de entender el método de la bisección. Vamos a realizarlo seguidamente con nuestro algoritmo de la bisección realizado con Octave y veamos qué soluciones obtenemos para los diferentes casos. Veamos qué pasa para la raíz -1 con intervalo inicial [−1.5, − 0.5], con una tolerancia de 0.001 y un número máximo de iteraciones de 20. Vemos que hemos obtenido el mismo resultado que cuando lo hemos hecho a mano, cosa que era de esperar. biseccion('2*power(x,3)-3*x-1',-1.5,-0.5,0.001,20) La solucion es -1.000000000000000. Se ha alcanzado en 1 iteraciones. 2 El método siempre converge porque el intervalo de partida siempre tiene un cambio de signo RESOLUCIÓN ECUACIONES
46. ECUACIONES CAPÍTULO 3 Ahora lo vamos a realizar para x2 = iteraciones como máximo: √ 1+ 3 2 49 en el intervalo [1,1.5], una tolerancia de 0.001 y 20 biseccion('2*power(x,3)-3*x-1',1,1.5,0.001,20) La solucion es 1.366210937500000. Se ha alcanzado en 9 iteraciones. Podemos apreciar que hemos necesitado 9 iteraciones para obtener 3 cifras decimales exactas. Si cambiamos la tolerancia a 0.000001 apreciamos que con 19 iteraciones obtenemos que las 5 primeras cifras decimales son exactas. biseccion('2*power(x,3)-3*x-1',1,1.5,0.000001,20) La solucion es 1.366024971008301. Se ha alcanzado en 19 iteraciones. Además lo podemos apreciar de forma gráﬁca con la siguiente ﬁgura: La Funcion Tendencia de errores en Biseccion Metodo Biseccion 0.5 1 0.4 1.1 1.2 1.3 1.4 y 1.5 Errores 200 x x 400 0.3 0.2 −1 0.1 0 −4 −2 0 2 4 y 2*power(x,3)-3*x-1 6 8 −2 2*power(x,3)-3*x-1 Aproximaciones 5 10 15 NumeroErrores Iteraciones Puntos Errores 20 Figura 2: Método Bisección p( x ) = 2x3 − 3x − 1 en [1,1.5] Los asteriscos rojos nos indican las diferentes aproximaciones que vamos obteniendo con el método en cada paso, en nuestro caso las aproximaciones son los puntos medios de cada subintervalo en el cual hay un cambio de signo. En la gráﬁca de la izquierda podemos apreciar la rapidez en la que converge nuestro método en cada caso particular. √ Para terminar vamos a ver qué ocurre con x1 = 1−2 3 , le vamos a aplicar una tolerancia de 0.000001 con un máximo de 20 iteraciones, y todo ello sobre el intervalo [−0.75, − 0.25] RESOLUCIÓN ECUACIONES
47. CAPÍTULO 3 50 biseccion('2*power(x,3)-3*x-1',-0.75,-0.25,0.000001,20) La solucion es -0.366024971008301. Se ha alcanzado en 19 iteraciones. Podemos apreciar que para alcanzar las 5 cifras decimales exactas hemos necesitado como antes 19 iteraciones. Las gráﬁcas de resolución son muy similares también a las anteriores. Metodo Biseccion −2 4 0.2 6 0.4 0.3 x 2 x −4 0.5 0.4 100 −6 Tendencia de errores en Biseccion 0.6 y −100 −0.8 −0.6 0.2 −0.4 Errores La Funcion 200 y −200 −0.2 −300 −0.4 2*power(x,3)-3*x-1 2*power(x,3)-3*x-1 Aproximaciones 0.1 0 5 10 15 NumeroErrores Iteraciones 0 20 Puntos Errores Figura 3: Método Bisección p( x ) = 2x3 − 3x − 1 en [−0.75, − 0.25] Ejercicio 3.0.2 Coseno Regula-Falsi 2 Halla la menor raíz positiva de la ecuación e− x − cos x = 0 con una tolerancia de 10−5 , utilizando el método de regula falsi y mostrando los pasos en las primeras iteraciones. Resolución 2 En ésta ecuación hay una solución obvia que es x = 0. Tomamos como función f ( x ) = e− x − cos x y dibujamos su gráﬁca para hacernos una idea de dónde se encuentran las raíces/soluciones. 1 2 3 x=-3:0.1:3; y=exp(-x.^2)-cos(x); plot(x,y), grid on Viendo la ﬁgura apreciamos que tenemos 2 raíces más cerca de -1 y 1. Vamos a calcular sólo la que está cerca de 1 porque la función es simétrica respecto el eje OY, cumpliéndose pues que la otra raíz es la misma que la que esta cerca de 1 pero cambiada de signo (su raíz opuesta). RESOLUCIÓN ECUACIONES
48. ECUACIONES CAPÍTULO 3 51 1 0.8 0.6 0.4 0.2 0 −0.2 −2 0 2 2 Figura 4: Función f ( x ) = e− x 2 2 f (− x ) = e−(− x) − cos(− x ) = e− x − cos x = f ( x ) ; ∀ x ∈ R Para aplicar el método de la Regula Falsi vamos a considerar el intervalo de partida [ a1 ,b1 ] = [1,2] f ( a1 ) = f (1) = −0.172423 f (b1 ) = f (2) = 0.434462 Siguiendo el método debemos calcular el punto c1 donde la recta que une los puntos ( a1 , f ( a1 )) y (b1 , f (b1 )) corta al eje OX (y = 0): c1 = a1 · f (b1 )−b1 · f ( a1 ) f (b1 )− f ( a1 ) = 1· f (2)−2· f (1) f (2)− f (1) = 0.434462+2·0.172423 0.434462+0.172423 = 1.284111 f (c1 ) = f (1.284111) = −0.090521 Ahora hay que ver en qué subintervalo se encuentra nuestra raíz aplicando Bolzano, los subintervalos son [ a1 ,c1 ] y [c1 ,b1 ]. Tenemos que f (c1 ) · f (b1 ) = −0.090521 × 0.434462 = −0.00393279347 < 0 Luego la raíz se encuentra en éste subintervalo. Si seguimos nuestro algoritmo de la Regula Falsi nombramos a2 = c1 , b2 = b1 , teniendo pues el subintervalo [ a2 ,b2 ] = [1.284111,2]. Calculamos c2 : RESOLUCIÓN ECUACIONES
49. CAPÍTULO 3 c2 = a2 · f (b2 )−b2 · f ( a2 ) f (b2 )− f ( a2 ) = 1.284111· f (2)−2· f (1.284111) f (2)− f (1.284111) = 52 0.5578974333+0.181042 0.434462+0.090521 = 1.40754926 f (c2 ) = f (1.40754926) = −0.02461859979 Realizamos un paso más: [ a2 ,c2 ] = [1.284111,1.40754926] [c2 ,b2 ] = [1.40754926,2] f (b2 ) · f (c2 ) = f (2) · f (1.40754926) = 0.434462 × (−0.02461859979) = −0.0106958461 < 0 Luego la raíz se encuentra en [c2 ,b2 ] := [ a3 ,b3 ] = [1.40754926,2]. Calculamos c3 c3 = a3 · f (b3 )−b3 · f ( a3 ) f (b3 )− f ( a3 ) = 1.40754926· f (2)−2· f (1.40754926) f (2)− f (1.40754926) = 0.6115266666+0.04923719958 0.434462+0.02461859979 = 1.439319951 f (c3 ) = −0.005119 Así pues, con 3 iteraciones del Método de la Regula Falsi obtenemos como raíz aproximada 1.439319951 con 2 cifras decimales exactas. Vamos a resolverlo con nuestro algoritmo, tomamos como intervalo de partida el [1,2], con una tolerancia de 0.0000001 y un número máximo de iteraciones de 20: regulafalsi('exp(-x.^2)-cos(x)',1,2,0.000001,20) La solucion es 1.447414195298149, y se ha alcanzado en 10 iteraciones. Con un error de 3.220810e-07 Como vemos hemos obtenido 6 cifras decimales exactas con 10 iteraciones, lo cual esta muy bien. En la siguiente tabla se puede apreciar iteración tras iteración el cálculo de la raíz con su respectivo error cometido. RESOLUCIÓN ECUACIONES
50. ECUACIONES CAPÍTULO 3 53 Iteración k Raíz xk Error ek 1 2 3 4 5 6 7 8 9 10 1.000000000000000 1.284111105170264 1.407549190540402 1.439319902661950 1.445849349137007 1.447114778447193 1.447357067800570 1.447403349491422 1.447412186157907 1.447413873217170 0.284111105170264 0.123438085370137 0.031770712121548 0.006529446475057 0.001265429310186 0.000242289353378 0.000046281690852 0.000008836666484 0.000001687059263 0.000000322080979 Esta tabla se genera cada vez que ejecutamos el método y se almacena en el ﬁchero regulafalsi.txt. Como sabemos que la función es simétrica respecto al eje OY podemos concluir que la otra raíz que nos falta por calcular es la opuesta a la que hemos calculado, es decir; la raíz es aproximadamente −1.447413873217170, sin tener que aplicar ningún método. Ejercicio 3.0.3 Coseno Punto Fijo a) Halla las raíces de la ecuación x = 2 · cos x por el método del punto ﬁjo, ilustrando el método paso a paso. b) Estudia el valor de la derivada de la función elegida g( x ) cerca de la raíz y justiﬁca el comportamiento del método. Resolución a) La función a la que le queremos hallar sus raíces es f ( x ) = x − 2 · cos x. Tomamos como función de punto ﬁjo la obvia, g1 ( x ) = 2 · cos x. Vamos a ver que esta función de punto ﬁjo es inestable, por lo que deberemos buscar otra como alternativa. Representamos la función g1 ( x ) = 2 · cos x junto a la bisectriz del primer cuadrante para obtener la estimación inicial. Para ello podemos introducir en Octave las siguientes órdenes: 1 2 3 x = linspace(-3,3); y = 2*cos(x); plot(x,y,x,x), grid on En el dibujo de ambas gráﬁcas podemos ver que la bisectriz corta a la función en el punto con abcisa 1, por lo que vamos a tomar como estimación inicial x1 = 1, y aplicamos el algoritmo de RESOLUCIÓN ECUACIONES
51. CAPÍTULO 3 54 2 0 −2 −2 0 2 Figura 5: Función f ( x ) = 2 cos( x ) recta y = x punto ﬁjo: x2 = g1 ( x1 ) = g1 (1) = 2 · cos(1) = 1.080604612 x3 = g1 ( x2 ) = g1 (1.080604612) = 2 · cos(1.080604612) = 0.941590072 x4 = g1 ( x3 ) = g1 (0.941590072) = 2 · cos(0.941590072) = 1.177006409 x5 = g1 ( x4 ) = g1 (1.17700409) = 2 · cos(1.17700409) = 0.7673820112 x6 = g1 ( x5 ) = g1 (0.7673820112) = 2 · cos(0.7673820112) = 1.439461363 Y no continuamos porque se aprecia claramente que los resultados obtenidos no se estabilizan, por el contrario, son oscilantes. Teóricamente hay una forma de comprobarlo mediante la condición de convergencia, se trata de ver si | g ( x )| < 1. Veámoslo en nuestro caso: g1 ( x ) = −2 · sin x → | g1 ( x )| = | − 2 · sin x | < 2 · | sin x | < 2 > 1; ∀ x ∈ R Por tanto, es claro que no se cumple la condición de convergencia para la función de punto ﬁjo g1 ( x ) = 2 · cos x que hemos considerado3 . Así pues, debemos buscar otra función de punto ﬁjo, 3 Dicha función no cumple el criterio de convergencia ni en un entorno alrededor de un punto, ya que el módulo de su derivada es mayor que 1 en todo el conjunto de números reales, por lo que nos vemos obligados a buscar otra función de punto ﬁjo. Si no pudiésemos encontrar otra deberíamos aplicar otro método iterativo como alternativa a buscar la solución. RESOLUCIÓN ECUACIONES
52. ECUACIONES CAPÍTULO 3 55 para ello basta hacer unos cuantos pasos algebraicos a partir de nuestra ecuación de partida: x = 2 · cos x → x 2 x x = cos x → arc cos( 2 ) = x → g2 ( x ) = arc cos( 2 ) Veamos si dicha función cumple el criterio de convergencia para un entorno de 1, | x | < 14 : g2 ( x ) = 2 −1 x 1−( 2 )2 → | g2 ( x )| = 2 −1 x 1−( 2 )2 ≤ 1 2 1 4− x 2 4 ≤ 1 4 √ 1 4− x 2 < 1; | x | < 1 Luego el criterio de convergencia se cumple en un entorno de 1, lo cual es suﬁciente para poder tomar como valor inicial x1 = 1, ya que pertenece a dicho entorno. Aplicamos pues el método de punto ﬁjo: x2 = g2 ( x1 ) = g2 (1) = arc cos( 1 ) = 1.047197551 2 x3 = g2 ( x2 ) = g2 (1.047197551) = arc cos( 1.047197551 ) = 1.019726744 2 x4 = g2 ( x3 ) = g2 (1.019726744) = arc cos( 1.019726744 ) = 1.035770367 2 x5 = g2 ( x4 ) = g2 (1.035770367) = arc cos( 1.035770367 ) = 1.026419399 2 x6 = g2 ( x5 ) = g2 (1.026419399) = arc cos( 1.026419399 ) = 1.031875939 2 Como vemos, tras 6 iteraciones el método se mantiene estable para esta función de punto ﬁjo. Apliquemos nuestro algoritmo de punto ﬁjo para una tolerancia de 0.0001 y 20 iteraciones como tope: pﬁjo('acos(x/2)','x-2*cos(x)',0,1.5,1,0.0001,20) La solucion es 1.029893453594281, y se ha alcanzado en 14 iteraciones. Con un error de 7.308525e-05 Vamos a dibujar la función f ( x ) = x − 2 · cos x en el intervalo [−5,5] con las órdenes: 1 2 3 clf() fplot('x-2*cos(x)',[-5,5]) grid on 4 Notar que fuera de dicho entorno esta segunda función hará que nuestro método sea inestable. RESOLUCIÓN ECUACIONES
53. CAPÍTULO 3 56 Podemos apreciar que dicha función sólo corta al eje OX en la solución que hemos obtenido, por lo que el problema queda terminado. Para asegurarnos deberíamos hacer un razonamiento más teórico sobre su crecimiento y decrecimiento de la función, pero eso te lo dejo a ti lector como deberes por ser un sencillo ejercicio de segundo bachiller. 6 x-2.*cos(x) 4 2 0 −2 −4 −6 −6 −4 −2 0 2 4 6 Figura 6: Función f ( x ) = x − 2 cos( x ) Ejercicio 3.0.4 Raíces Quintas a) Utiliza el método de Newton para aproximar la raíz quinta de 2 con tres cifras decimales exactas, ilustrando el método paso a paso. b) Comprueba la convergencia cuadrática del método de Newton. Resolución a) Muy bien, lo primero es establecer la función a la que le vamos a aplicar el método. Nos dicen √ de aproximar la raíz quinta de 2, lo cual es lo mismo que decir x = 5 2, vamos a manipularlo algebraicamente hasta obtener el polinomio de quinto grado con solución la raíz quinta de 2: x= √ 5 2 → x5 = 2 → x5 − 2 = 0 Ya tenemos el polinomio o función a la que aplicar el método de Newton, f ( x ) = x5 − 2. Seguif (x ) damente debemos aplicar el método, que de forma general es: xn = xn−1 − f ( xn−1 ) , partiendo de n −1 RESOLUCIÓN ECUACIONES
54. ECUACIONES CAPÍTULO 3 57 un valor inicial de partida x0 . Pero primero calculemos la derivada, f ( x ) = 5x4 . √ En nuestro caso ya sabemos que la solución exacta es 5 2, pero lo que queremos es aproximarla con 3 cifras decimales exactas (esta solución es irracional), así que va a ser fácil determinar el valor inicial, vamos a tomar x0 = 1, porque si calculamos con una calculadora tenemos que √ 5 2 = 1.148698355 vemos que 1 está cerca de la raíz.5 Perfecto, iteremos: x1 = 1 − f (1) 1 = 1 + = 1.2 f (1) 5 x2 = 1.2 − f (1.2) 0.48832 = 1.2 − = 1.152901235 f (1.2) 10.368 x3 = 1.152901235 − f (1.152901235) f (1.152901235) x4 = 1.148728887 − f (1.148728887) 0.000265811 = 1.148728887 − = 1.148698357 f (1.148728887) 8.706431225 = 1.152901235 − 0.036856916 8.83361408 = 1.148728887 Lo que debemos apreciar es que de todos los métodos que llevamos aplicados hasta ahora en este documento, y salvando la diferencia importante de que todas las funciones han sido diferentes; es que éste método ha sido el más rápido con diferencia. Si lo realizamos con nuestro algoritmo de Octave con 10−7 obtenemos lo siguiente: newton('x.^5-2','5*x.^4',1,-1,2,1e-7,20) La solucion es 1.148698354997035, y se ha alcanzado en 6 iteraciones. b) Y ahora veamos que la convergencia es cuadrática. En nuestro algoritmo de Newton ya lo hemos añadido todo para ver bajo gráﬁcas qué tipo de convergencia alcanzamos. Viendo en la ﬁgura adjunta apreciamos que no se cumple la lineal porque la gráﬁca tiende a 0 en vez de tender a 1. En la gráﬁca de la cuadrática apreciamos que tiende cerca de 2, por lo que concluimos que la convergencia es 2. Además, el algoritmo también nos crea el ﬁchero newtonraphson.txt en el cual podemos ver los valores de ambas tasas de convergencia. En nuestro caso: Tasa de convergencia lineal = 0.235494 Tasa de convergencia lineal = 0.088587 Tasa de convergencia lineal = 0.007317 Tasa de convergencia lineal = 0.000053 5 No podemos tomar como valor inicial 0 o en un entorno pequeño de 0 porque f (0) = 0, y ya sabemos que el método de Newton falla en valores que hacen que su primera derivada se anule. RESOLUCIÓN ECUACIONES
55. CAPÍTULO 3 58 Tasa de convergencia cuadratica = 1.177469 Tasa de convergencia cuadratica = 1.880882 Tasa de convergencia cuadratica = 1.753737 Tasa de convergencia cuadratica = 1.741193 De ésta segunda forma podemos apreciar muy bien que la Tasa de Convergencia Lineal tiende claramente a 0 y la cuadrática aproximadamente a 1.7, que es cerca de 2. Tasa Convergencia Lineal Newton Tasa Convergencia Cuadratica Newton 0.25 2 0.2 1.8 0.15 tcl tcc 1.6 0.1 5 · 10−2 2 3 Numero Iteraciones 1 4 1.4 1.2 0 1 T.C.L. 2 3 Numero Iteraciones 4 1 T.C.C. Figura 7: Tasas De Convergencia Ejercicio 3.0.5 Función Billar Para jugar al billar en una mesa circular, hemos de golpear la bola Q con la bola P, tras un impacto I en la banda. Conocido el radio R de la mesa y las posiciones en coordenadas cartesianas (cuyo origen es el centro de la mesa) de los puntos P( x p , y p ) y Q( xq , yq ), el punto de impacto viene deﬁnido por el ángulo central t. Tras un análisis geométrico del problema, se prueba que los valores de t que proporcionan los puntos de impacto posibles son los ceros de la función: f (t) = x p sin t − y p cos t ( R cos t − x p )2 + ( R sin t − y p )2 + xq sin t − yq cos t ( R cos t − xq )2 + ( R sin t − yq )2 Consideramos un billar de radio unidad, P = (0.6, − 0.3) y Q = (−0.6, − 0.3) . Determina los puntos de impacto utilizando los métodos de bisección, regula falsi, Newton y secante. Analiza la convergencia de los distintos métodos y compara los resultados obtenidos. RESOLUCIÓN ECUACIONES
56. ECUACIONES CAPÍTULO 3 59 ESQUEMA I t x O Q P Figura 8: Billar Circular Resolución a) Método de bisección: En los códigos verás un ﬁchero llamado fbillar.m en el cual está introducida la función impacto I(t) para poder llamarla cada vez que ejecutamos uno de nuestros algoritmos iterativos. Hagámonos una idea de la función impacto visualizando su gráﬁca, que puedes apreciar en la ﬁgura. 1 % Para que los metodos iterativos funcionen es biseccion('fbillar(t)',1,2,1e-4,100) 2 3 4 5 6 7 8 9 10 11 12 function [y] = fbillar(t) P=[0.6,-0.3]; Q = [-0.6,-0.3]; y1 = P(1)*sin(t)-P(2)*cos(t); y2 = Q(1)*sin(t)-Q(2)*cos(t); den1 = sqrt((cos(t)-P(1)).^2+(sin(t)-P(2)).^2); den2 = sqrt((cos(t)-Q(1)).^2+(sin(t)-Q(2)).^2); y1 = y1./den1; y2 = y2./den2; y = y1 + y2; end fbillar.m Al ser una función trigonométrica sabemos que es una función periódica 6 Los trozos de función se repiten cada cierto periodo T. RESOLUCIÓN ECUACIONES 6 de periodo 2π, por
57. CAPÍTULO 3 60 La Funcion Billar 1 0.5 −6 −4 −2 2 4 6 −0.5 −1 Figura 9: Función Billar lo que nos vamos a ﬁjar en el intervalo [0,2π ]. Sabiendo que las funciones seno y coseno son de periodo 2π con lo que no es difícil de comprobar que nuestra función impacto es de periodo 2π, I (t) = I (t + 2π ); ∀t ∈ R. Como el método de bisección siempre converge queda ver cuáles son los intervalos a utilizar, en [0,2π ] vemos que hay 4 raíces reales en las que podemos situarlas en los 4 respectivos intervalos: [1,2], [3,4], [4,5] [5,6]. Ejecutamos el algoritmo de bisección para el intervalo [1,2] con el comando: biseccion('fbillar(t)',1,2,1e-4,100) Y obtenemos la solución: La solucion es 1.570739746093750. Se ha alcanzado en 14 iteraciones. Con un error de 1.220703e-04 Antes que nada te aconsejo que mires con detenimiento la gráﬁca de la función en este intervalo y la solución que hemos obtenido, igual así no lo recuerdas pero la solución viene en radianes y es una aproximación de π que son 90◦ . Luego aquí hemos hecho un poco el primo ya que no es 2 necesario aplicar ningún método iterativo para saber que la solución exacta es un ángulo de 90◦ , no te desesperes si no has caido, es normal y como mal menor piensa que el juego del billar nos da una forma de aproximar π , aunque sea un poco burda. 2 Por tanto tenemos que las raíces distintas en [0,2π ] son: RESOLUCIÓN ECUACIONES
58. ECUACIONES CAPÍTULO 3 61 Intervalo Solución En Radianes Solución En Grados Iteraciones [1,2] [3,4] [4,5] [5,6] 1.570739746093750 3.871276855468750 4.712341308593750 5.553405761718750 89.996758164622406 221.807825145036162 269.997268607577325 318.186712070118517 14 14 14 14 En la tabla observamos que al tener los 4 intervalos la misma longitud nos da que el método de la bisección con una tolerancia de 1e-4 siempre nos hace el mismo número de iteraciones, 14 en éste caso. Quizá te estés preguntando por qué no calculamos las raíces negativas, pues no las calculamos porque son las opuestas que las positivas por ser nuestra función periódica, por lo que ya sabemos cuáles son. Además, al ser negativas lo que debemos saber es que los ángulos son negativos, por lo que las trayectorias de las bolas se hacen en sentido positivo; mientras que las raíces positivas nos indican que las trayectorias son en sentido negativo.7 b) Método Regula Falsi: Lo aplicamos con una tolerancia de 1e-12 y un máximo de 50 iteraciones, obteniendo los siguientes resultados: Intervalo [1,2] [3,4] [4,5] [5,6] Solución En Radianes Solución En Grados Errores Iteraciones 1.570796326794897 3.871320309816781 4.712388980384690 5.553457650952138 90.000000000000000 221.810314895779868 270.000000000000000 318.189685104193813 2.220446e-16 2.304823e-13 8.881784e-16 6.670220e-13 5 13 6 35 Se aprecia que aun habiendo aumentado considerablemente la tolerancia el método es mucho más rápido que el de la bisección, incluso obtenemos en algunos una tolerancia mayor que la de 1e-12. Además, observamos que la primera solución ya nos la marca como 90 grados, al igual que la tercera que nos dice que son 270 8 c) Método de Newton: Lo primero va a ser crear un ﬁchero con código para implementar la primera derivada de la función del problema, al ﬁchero le he llamado fbillarderivada.m y el código es el siguiente: 1 2 3 4 function [dy] = fbillarderivada(t) P=[0.6,-0.3]; Q = [-0.6,-0.3]; x = P(1); y = P(2); 7 El sentido positivo es el sentido contrario a las agujas del reloj, y el negativo es en el sentido de las agujas del reloj. depende del ordenador, de la cantidad de cifras decimales que indiquemos en el programa que aparezcan por pantalla y del tipo de redondeo de Octave, entre otros muchos factores. 8 Esto RESOLUCIÓN ECUACIONES
59. CAPÍTULO 3 5 6 7 8 x1 = y1 = den1 den2 62 Q(1); Q(2); = sqrt(1-2*(x*cos(t)+y*sin(t))+x^2+y^2); = sqrt((1-2*(x*cos(t)+y*sin(t))+x^2+y^2).^3); 9 10 11 num1 = x*cos(t)+y*sin(t); num2 = (x*sin(t)-y*cos(t)).^2; 12 13 14 dy1 = num1/den1; dy2 = num2/den2; 15 16 17 den11 = sqrt(1-2*(x1*cos(t)+y1.*sin(t))+x^2+y^2); den21 = sqrt((1-2*(x1*cos(t)+y1*sin(t))+x^2+y^2).^3); 18 19 20 num11 = x1*cos(t)+y1*sin(t); num21 = (x1*sin(t)-y1*cos(t)).^2; 21 22 23 dy3 = num11/den11; dy4 = num21/den21; 24 25 26 dy = dy1 - dy2 + dy3 - dy4; end fbillarderivada.m Para calcular la derivada de la función lo puedes hacer a mano con mucha paciencia o usar el programa de matemáticas simbólicas wxMaxima de código libre. Una vez que lo tenemos todo, y para no hacer el problema muy largo vamos a calcular con Newton sólo la primera raíz y comparar con los métodos anteriores, así que en Octave introducimos y obtenemos: newton('fbillar(t)','fbillarderivada(t)',1,0,1.5,1e-6,20) La solucion es 1.570796326794897, y se ha alcanzado en 5 iteraciones. Si comparamos con el método de la bisección es claro que obtenemos una mayor precisión en Newton y en un número sensiblemente menor de iteraciones, tan sólo 5 iteraciones. También podemos comparar si el método de Newton tiene convergencia lineal o cuadrática mediante las siguientes gráﬁcas. Pero éste problema tiene un handicap muy a tener en cuenta, y es que calcular la derivada de nuestra función y sustituir los valores en cada iteración tiene un coste muy elevado, aunque hoy en día con los ordenadores que tenemos casi no se nota, a no ser que tengamos grandes cálculos. En nuestro caso le he aplicado las funciones tic y toc a nuestro algoritmo de newton con ésta función y no tarda más allá de 0.11 segundos. Aunque sí es cierto que lo mejor es no aplicar RESOLUCIÓN ECUACIONES
60. ECUACIONES CAPÍTULO 3 Tasa Convergencia Lineal Newton ·10−2 5 63 Tasa Convergencia Cuadratica Newton ·10−2 8 4 6 4 tcc tcl 3 2 2 1 1 2 2.5 1.5 Numero Iteraciones T.C.L. 3 0 1 2 2.5 1.5 Numero Iteraciones 3 0 T.C.C. Figura 10: Tasas Convergencia Billar Newton Newton y tirar por otro lado, porque el tiempo que invertimos es bastante calculando la primera derivada de nuestra función a mano o con ordenador; ése es el verdadero motivo por el que en éste prolema el método de Newton no es recomendable. d) Método de la Secante: En éste método no tenemos que calcular ninguna derivada, por lo que si obtenemos unos resultados similares a los obtenidos en Newton sería nuestra verdadera alternativa para calcular las raíces: secante('fbillar(t)',1,2,1e-6,10) La solucion es 1.570796326794897, y se ha alcanzado en 4 iteraciones. Efectivamente, incluso con una iteración menos que con Newton obtenemos idéntico resultado y sin el costo adicional de tener que calcular la primera derivada. Por tanto, el método de la secante es la mejor alternativa de las 4 vistas aquí para nuestro peoblema.9 9 Aunque no lo hemos visto aquí, es obvio que el método del punto ﬁjo en éste problema es prácticamente inviable, no sólo por el elevado costo que supone encontrar la función de punto ﬁjo, sino también comprobar si el módulo de su primera derivada cumple o no el criterio de convergencia. Otra alternativa muy interesante es aplicar el algoritmo de Steffensen. RESOLUCIÓN ECUACIONES
61. CAPÍTULO 3 Ejercicio 3.0.6 64 Problema Catenaria Se denomina catenaria a la forma que adopta una cadena o un cable suspendido de sus 2 extremos. Eligiendo adecuadamente el sistema de referencia, la ecuación de la catenaria es: x y = λ cosh( ), donde el origen de coordenadas está situado en la vertical del punto más λ bajo de la cadena, a distancia λ del mismo. Sea d la diferencia de altura entre el punto de la cadena de abcisa igual a 1 y el de altura mínima: d = y(1) − y(0). Conociendo d, podemos 1 determinar λ resolviendo la ecuación λ + d = λ cosh( ). Supongamos d = 0.1. λ a) Determina con precisión de milésimas el valor de λ utilizando la ecuación de punto ﬁjo 1 λ = λ cosh( ) − d λ . b) Comprueba que la función de punto ﬁjo 1 λ= arg cosh λ+d λ es también equivalente a la ecuación 1 λ + d = λ cosh( ) λ . Itera con esta ecuación y la misma estimación inicial para hallar λ. c) ¿Qué relación existe entre la derivada respecto a de las funciones de iteración consideradas en los apartados anteriores y el proceso de convergencia? Resolución 1 a) Tenemos que nuestra función de punto ﬁjo es g(λ) = λ cosh λ − d. Vamos a dibujar la función junto con la bisectriz ppal. del plano para ver dónde cortan y qué valor inicial podemos tomar. 10 1 2 3 x=1:0.1:6; plot(x,x.*cosh(1./x)-0.1,x,x), grid on; legend({'g(x)','Bisectriz'},'location','southoutside'); 10 Ten cuidado, no dibujes la función incluyendo valores de λ en ]0,1[ porque sus imágenes se disparan y no verás nada. El que sus imágenes se disparen es debido a que la función cosh es unacombinación de la función exponencial, cosh x = 1 x −x ) 2 (e + e RESOLUCIÓN ECUACIONES
62. ECUACIONES CAPÍTULO 3 65 En la ﬁgura se aprecia que ambas funciones se cortan alrededor de 5, por lo que tomaremos como valor inicial 4 en nuestro método de punto ﬁjo. 6 5 4 3 2 1 1 2 3 g(x) Bisectriz 4 5 6 Figura 11: Catenaria x cosh( 1 ) − 0.1 x Lo vamos a aplicar con una tolerancia de 1e-4: pﬁjo('x.*cosh(1./x)-0.1','x-x.*cosh(1./x)-0.1',1,5,4,1e-4,300) La solucion es 5.011784191306107, y se ha alcanzado en 255 iteraciones. Y vemos que sí obtenemos la solución pero el número de iteraciones necesarias es excesivo, así que debemos encontrar otra función de punto ﬁjo, la que vamos a usar en b). b) Vamos a demostrar cómo llegamos a esa función de punto ﬁjo. λ + d = λ cosh 1 λ → λ+d = cosh λ 1 λ → 1 = arg cosh λ λ+d λ →λ= 1 arg cosh λ+d λ Y ahora que ya lo hemos demostrado lo que vamos a hacer es dibujarla junto con la bisectriz: 1 2 3 x=4:0.1:6; plot(x,1./(acosh((x+0.1)./x)),x,x), grid on; legend({'g(x)','Bisectriz'},'location','southoutside'); RESOLUCIÓN ECUACIONES
63. CAPÍTULO 3 66 6 5.5 5 4.5 4 4 4.5 5 g(x) Bisectriz Figura 12: Catenaria 5.5 6 1 acosh( x+0.1 ) x En la ﬁgura observamos que obtenemos la misma solución que con la anterior función de punto ﬁjo, apliquemos el método y a ver qué pasa: pﬁjo('1./(acosh((x+0.1)./x))','x-1./(acosh((x+0.1)./x))',1,5,4,1e-4,20) La solucion es 5.016512483824275, y se ha alcanzado en 15 iteraciones. Con un error de 6.658831e-05 ¡Bingo!, con ésta función de punto ﬁjo obtenemos la aproximación deseada y además mucho más rápidamente, sólo con 15 iteraciones, luego ésta segunda función de punto ﬁjo es la óptima para nuestro problema de la catenaria.11 c) Para realizar el apartado lo que haremos es ver la tendencia para ambas funciones su Tasa de Convergencia Lineal y compararlas. En la primera función tiende a 1, por eso la aproximación es lenta y mala. En la segunda función la tasa tiende a 0.5<1, por esa razón el método converge rápidamente y con una aproximación muy buena. 11 Atención, ésta es la óptima para d = 0.1, pero puede que no lo sea para otros valores diferentes de d, te animo a que experimentes por tu cuenta, es un buen ejercicio de aprendizajes, ¡ánimooo ! ¡No muerde ! RESOLUCIÓN ECUACIONES
64. ECUACIONES CAPÍTULO 3 67 Tasa Convergencia Lineal Punto Fijo 0.98 0.98 0.97 tcl 0.98 0.97 0.97 0 200 100 Numero Iteraciones 0.97 300 T.C.L. Figura 13: TCL Catenaria x cosh( 1 ) − 0.1 x Tasa Convergencia Lineal Punto Fijo 0.52 0.5 0 5 10 Numero Iteraciones T.C.L. Figura 14: TCL Catenaria RESOLUCIÓN ECUACIONES 1 acosh( x+0.1 ) x tcl 0.54
65. CAPÍTULO 3 Ejercicio 3.0.7 68 Raíces Décimas Compara diferentes métodos para resolver la ecuación p( x ) = x10 − 1 = 0, con una tolerancia de 10−6 . Resolución Bueno la ecuación tiene 10 raíces entre complejas y reales, de las reales tenemos 2, que son 1 y -1. Nos vamos a centrar en 1. Apliquemos el método de la bisección: biseccion('x.^10-1',0.5,1.5,1e-6,30) La solucion es 1.000000000000000. Se ha alcanzado en 1 iteraciones. biseccion('x.^10-1',0,1.5,1e-6,30) La solucion es 0.999999761581421. Se ha alcanzado en 21 iteraciones. Bien, lo primero que observamos es que hay que tener en cuenta qué intervalo escogemos, en el primer caso obtenemos la solución exacta y en sólo 1 iteración, eso se debe a que el punto medio del intervalo coincide exactamente con la solución. Sin embargo, en el segundo caso su punto medio no es la solución, además no obtenemos la solución exacta y le cuesta hacerlo en 21 iteraciones. Por tanto, en éste caso no hace falta que sigamos buscando métodos porque hemos obtenido la solución ideal. Obviamente esto no ocurrirá en la mayoría de los casos, porque los métodos están hechos para aproximar raíces o para buscar raíces que no podemos (o es difícil) calcularlas mediante operaciones algebraicas. Veamos qué pasa con el método de la Regula Falsi: regulafalsi('x.^10-1',0.5,1.5,1e-6,200) La solucion es 0.999989952440931, y se ha alcanzado en 135 iteraciones. Con un error de 9.723436e-07 regulafalsi('x.^10-1',0,1.5,1e-6,200) La solucion es 0.999989784169578, y se ha alcanzado en 158 iteraciones. Con un error de 9.886274e-07 Con éste método apreciamos que no obtenemos la solución exacta y además es muy lento, esto se debe a que en los intervalos elegidos la función casi no varía y esta cerca de cero, como puedes apreciar en la ﬁgura adjunta, puedes ver en rojo que las iteraciones se acumulan todas en el intervalo elegido con una distancia entre ellas muy pequeña. RESOLUCIÓN ECUACIONES
66. ECUACIONES CAPÍTULO 3 Tendencia de errores en Regula Falsi ·10−2 2 Metodo Regula Falsi La Funcion ·108 1.5 60 1 1.5 0.5 Errores 40 x x 69 20 −5 1 0.5 5 y 0.6 −0.5 0.8 1 y 1.2 1.4 1.6 50 100 Numero Iteraciones x.^10-1 Aproximaciones x.^10-1 Errores Puntos Errores Figura 15: Regula Falsi p( x ) = x10 − 1 en [0.5,1.5] Pasemos al método de Newton a ver qué pasa. Lo primero que debes tener en cuenta es que no podemos tomar como punto inicial el cero o puntos alrededor de un entorno pequeño del 0, porque su derivada se anula. Apliquemos el método: newton('x.^10-1','10*x.^9',0.5,0,1,1e-6,50) El algoritmo tiene un coste en tiempo de: Elapsed time is 0.030031 seconds. La solucion es 1.000000000000000, y se ha alcanzado en 44 iteraciones. Aquí obtenemos la solución exacta pero con muchas iteraciones, con lo cual perdemos la convergencia cuadrática del método, debido a que las rectas tangentes que se calculan en cada iteración varían muy poco. Eso se puede apreciar muy bien en las gráﬁcas de tasas de convergencia en la ﬁgura adjunta. Y ahora veamos qué ocurre con el método de la secante. secante('x.^10-1',0,1.5,1e-6,200) La solucion es 0.051574707031250, y se ha alcanzado en 5 iteraciones. Vaya, aquí el método claramente falla de forma estrepitosa, obtenemos una solución incorrecta aun no teniendo problemas de divergencia. Para solucionarlo debemos tener en cuenta el método de la bisección, ya que además de converger siempre, se tiene en consideración en cada iteración un subintervalo que contiene a la raíz mediante la partición en mitades iguales los subintervalos. Si en el método de la bisección elegimos el intervalo [0.9375,1.1255] vemos que el método RESOLUCIÓN ECUACIONES
67. CAPÍTULO 3 Tasa Convergencia Lineal Newton 70 Tasa Convergencia Cuadratica Newton 1 0.6 tcl 4 tcc 6 0.8 0.4 2 0.2 0 20 30 10 40 Numero Iteraciones 0 50 T.C.L. 0 20 30 10 40 Numero Iteraciones 0 50 T.C.C. Figura 16: TC Newton p( x ) = x10 − 1 en [0.5,1.5] sólo utiliza 3 iteraciones. Si lo aplicamos ahora en el método de la secante se obtienen resultados coherentes. secante('x.^10-1',0.9375,1.1255,1e-6,20) La solucion es 0.999999999999997, y se ha alcanzado en 7 iteraciones. Puedes probar éste mismo intervalo en el de la Regula Falsi o en el de Newton eligiendo un valor inicial contenido en dicho intervalo. newton('x.^10-1','10*x.^9',1.1255,0,1,1e-6,50) El algoritmo tiene un coste en tiempo de: Elapsed time is 0.002492 seconds. La solucion es 1.000000000000979, y se ha alcanzado en 6 iteraciones. Ejercicio 3.0.8 Ecuación Logística La ecuación logística x = ax (1 − x ), a > 0 es una ecuación de punto ﬁjo que modela una población, cuyo crecimiento esta limitado. Es muy conocida porque los iterados pueden presentar un comportamiento caótico. RESOLUCIÓN ECUACIONES
Newton And Neville Interpolation

References: Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 resolución 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 RESOLUCIÓN