Artículos

Diseño y análisis de la convergencia y estabilidad de métodos iterativos para la resolución de ecuaciones no lineales

Design and analysis of convergence and stability of iterative methods for solving nonlinear equations

Armando Gabriel Solís Zúñiga
Instituto Tecnológico de Costa Rica, Costa Rica
Alicia Cordero Barbero
Universitat Politècnica de València, España
Juan Ramón Torregrosa Sánchez
Universitat Politècnica de València, España
Juan Pablo Soto Quirós
Instituto Tecnológico de Costa Rica, Costa Rica

Diseño y análisis de la convergencia y estabilidad de métodos iterativos para la resolución de ecuaciones no lineales

Revista Digital: Matemática, Educación e Internet, vol. 21, núm. 2, pp. 1-27, 2021

Instituto Tecnológico de Costa Rica

Derechos Reservados © 2021 Revista digital Matemática, Educación e Internet

Recepción: 25 Noviembre 2019

Aprobación: 30 Agosto 2020

Resumen: Una corriente en métodos numéricos es la elaboración de nuevos métodos iterativos para la resolución de ecuaciones no lineales; se busca que estos métodos sean óptimos en contraste con su orden de convergencia y el número de evaluaciones funcionales en comparación con métodos usualmente conocidos. El presente trabajo muestra el diseño de una nueva familia paramétrica de métodos iterativos basada en la familia de métodos de Chun que contiene el esquema de Ostrowski como ca- so particular. Se realiza además un análisis por medio de dinámica compleja para visualizar planos de parámetros y dinámicos para seleccionar los parámetros que presentan un comportamiento más estable para el esquema en estudio y así formar un método iterativo más eficiente.

Palabras clave: métodos iterativos, ecuaciones no lineales, métodos óptimos, dinámica compleja, análisis de convergencia.

Abstract: One stream in numerical analysis is the creation of new iterative methods for the resolution of non-linear equations; optimal processes are sought in contrast with their order of convergence and the number of functional evaluations compared to usual known methods. This article shows a design of a new parametric family of iterative methods based on the Chun’s family of methods which contains the particular case of the Ostrowski’s scheme. Through an analysis with complex dynamic it is intended to visualize dynamic planes and parameters’ planes to choose the best parameter who brings more stable behavior for the scheme under study and make it more efficient.

Keywords: iterative methods, non-linear equations, optimal methods, complex dynamic, analysis of convergence.

1. Introducción

La resolución de ecuaciones es un aspecto fundamental en el área de matemática. Tal como explica Kung y Traub (1974) en [11], algunas de las aplicaciones de ecuaciones matemáticas son: “estimación de la distancia en redes de sensores inalámbricos, desarrollo de osciladores caóticos en circuitos computacionales, entre otros”. Es por ello la importancia de hallar una solución a una ecuación. Un caso aplicado es encontrar la solución a la ecuación a x c o s h ( x ) = 0 , donde a es una constante. Esta ecuación se utiliza en el análisis e implementación de matrices de puertas programables (FPGA por sus siglas en inglés) en computación como se puede apreciar en [13].

Habitualmente, se introduce este concepto con ecuaciones lineales, en las cuales se usan algoritmos sencillos para despejar su incógnita involucrada; sin embargo, existen ecuaciones en las que no se puede determinar sus raíces por métodos algebraicos, es decir, no se pueden obtener sus soluciones de manera exacta. A raíz de esto, ha surgido la necesidad de crear métodos alternativos para la resolución de dichas ecuaciones, tales como los métodos iterativos.

Diversos investigadores han propuesto varios métodos. Estos métodos generan una sucesión de números reales que converge a una solución de la ecuación, si es que ésta existe. Sin embargo, algunos de estos métodos no son del todo eficientes. Por ejemplo: la velocidad de convergencia no es rápida, ya que implica el cálculo adicional de derivadas. Otro problema presentado para ciertos métodos es si la función posee ceros múltiples, como la función f ( x ) = x 3 donde posee una solución, el cero, con multiplicidad 3; lo cual genera que el esquema iterativo se indefina y la sucesión no converja a la solución. Dados estos problemas, surge el interés de investigar sobre la creación de nuevos métodos, perfeccionando las técnicas de algún procedimiento específico para dar mayores eficiencias al método de solución.

Como se menciona en Chapra y Canale (2017) [3]: “un modelo matemático se define, de manera general, como una formulación o una ecuación que expresa las características esenciales de un sistema físico o de un proceso en términos matemáticos”. Por ejemplo, en secundaria se realiza el estudio de ecuaciones lineales de la forma a x + b = 0 , donde a y b son números reales con a 0 , y cuya solución es: x = - b ± b 2 4 a c 2 a .

Sin embargo, no todas las ecuaciones tienen una forma explícita de solución. Por ejemplo, para la ecuación e x x 2 = 0 no existe una forma algebraica para encontrar la solución exacta; pero al graficar la función f ( x ) = e x x 2 se puede observar que ésta interseca al eje X en dos puntos (ver Figura1), lo que significa que la ecuación sí tiene, al menos, una solución en los reales.

Gráfica de la función f(x)=ex-x-2
Figura 1:
Gráfica de la función f(x)=ex-x-2

En este escrito se hace mención de métodos iterativos que aproximan, bajo determinadas condiciones, una solución α ∈R de la ecuación no lineal f ( x ) = 0 donde f : I R R está definida en un intervalo abierto I . Existen muchos métodos iterativos como el método de Newton, el esquema de Halley, el método de Chebyshev y otros más. El método de Newton es probablemente el más conocido y el que más se utiliza para resolver la ecuación f ( x ) = 0 . Su algoritmo es el siguiente:

x k + 1 = x k f ( x k ) f ' ( x k ) , k = 0 , 1 , . . . (1)

con f ' ( x k ) 0 para todo k 0 y x 0 algún valor inicial dado.

Sin embargo, existen otros métodos más eficientes para la resolución de ecuaciones no lineales. Basado en el esquema de una familia de métodos de Chun, se propone una familia de métodos parametrizada y, por medio de sistemas dinámicos, se estudia la estabilidad de los mismos.

Finalmente se realizan pruebas numéricas para comparar los métodos propuestos con algunos esquemas que se utilizan generalmente para la resolución de ecuaciones no lineales.

2. Métodos iterativos para ecuaciones no lineales

Para poder establecer una clasificación de cómo un esquema para resolución de ecuaciones no lineales es, o no, estable, se deben tener algunos conceptos claros.

En primera instancia, los conceptos de métodos iterativos: no solo el orden de convergencia de un método es relevante, si no también conceptos como eficiencia e índice de eficiencia. Se mencionará también la conjetura de un método óptimo para establecer una relación entre el orden de convergencia y número de evaluaciones funcionales por iteración.

2.1. Conceptos básicos

Definición 1 (Multiplicidad de la raíz de una función)

Sea f : I R R y α un cero de f ; o sea f ( α ) = 0 . El valor α se le llama cero de multiplicidad k de f si existe un número real c ; c 0 tal que:

lim x α | f ( x ) | | x α | k = c .

Si k = 1 , entonces α recibe el nombre de raíz simple.

Traub en [14] puntualiza las Definiciones2-7:

Definición 2 (Convergencia)

Sea α R , x n R , n = 0 , 1 , 2 , . . . entonces, se dice que la sucesión { x n } converge a α si

lim n | x n α | = 0

Definición 3 (Orden de convergencia)

Si, además de cumplir las condiciones de convergencia, existen: una constante 0 < c 1 , un natural n 0 0 y un entero p 1 tales que, para todo n > n 0 :

| x n + 1 α | c | x n α | p

o bien:

lim n | x n + 1 α | | x n α | p = c ;

entonces, se dice que la sucesión { x n } converge con al menos orden p .

Si p = 2 o 3 , entonces se dice que es de orden cuadrático o cúbico, respectivamente. Por ejemplo, dado el método de Newton presentado en (1), se dice que es de orden cuadrático para algún vecindario de α , es decir, que existen una constante c tal que lim k | x k + 1 α | | x k α | 2 = c .

Definición 4 (Orden de convergencia por medio de la ecuación del error)

Si e n = x n α denota el error en la n-ésima iteración, la relación:

e n + 1 = c e n p + a ( e n p + 1 ) , c R

se conoce como ecuación del error. Al valor p obtenido se le conoce como el orden del método. Esta definición es equivalente a la Definición 3.

Definición 5

Se le llama método multipaso a aquel que tiene la siguiente estructura;

y n ( 1 ) = G 1 ( x n )

y n ( 2 ) = G 2 ( x n )

x n + 1 = G ( x n , y n ( 1 ) , y n ( 2 ) , . . . )

donde G 1 , G 2 , G 3 , . . . , G son operadores.

El primero en introducir este concepto fue Traub en [14] con su método multipaso:

(2)

Más adelante en la Subsección2.3, se profundizará la importancia de los métodos multipasos para el trabajo.

Definición 6 (Evaluaciones funcionales)

Al número de veces en el que la función f y sus derivadas están siendo evaluadas por cada iteración se le llama número de evaluaciones funcionales y se denota por d .

Por ejemplo en el esquema de Newton presentado en (1) posee dos evaluaciones funcionales: f ( x k ) , y f ' ( x k ) ; por lo que d = 2 .

En el esquema de Traub presentado en (2) posee tres evaluaciones funcionales: f ( x k ) , f ' ( x k ) y f ( y k ) ; por lo que d = 3 .

Definición 7 (Eficiencia de un método)

Sea M un método iterativo. Entonces la eficiencia de M se define como:

I = p d ,

donde p es el orden de convergencia y d es el número de evaluaciones funcionales por iteración.

Ostrowski en [12] precisa la siguiente definición:

Definición 8 (índice de eficiencia de un método)

Sea M un método iterativo. Entonces el índice de eficiencia de M se define como:

EI = p 1 / d

donde p es el orden de convergencia y d es el número de evaluaciones funcionales por iteración.

2.1.1. Conjetura de Kung y Traub

Kung y Traub en [11] establecen la siguiente conjetura: “El orden de convergencia de cualquier método multipaso no puede ser mayor a 2 d 1 (llamado orden óptimo), donde d es el número de evaluaciones funcionales por iteración”.

Bajo esta conjetura, se estudiarán métodos los cuales cumplan que su orden de convergencia p sea igual a 2 d 1 ; mientras se cumpla esta condición se dirá que el método es óptimo.

2.1.2. Orden de convergencia computacional y orden de convergencia computacional apro- ximado

Cuando se prueban nuevos métodos, ya sea para verificar el orden de convergencia o para estimar qué tanto difiere del orden teórico en la implementación práctica, se puede calcular el orden de con- vergencia computacional (COC), definido por Weerakoon y Fernando en [15] como:

p l n | ( x k + 1 α ) / ( x k α ) | l n | ( x k α ) / ( x k 1 α ) | , k = 1 , 2 , . . .

donde x k + 1 , x k y x k 1 son las últimas tres aproximaciones sucesivas de α obtenidas en el proceso iterativo. Dado que el valor de α es desconocido en la práctica, se utiliza el orden de convergencia computacional aproximado (ACOC) definido por Cordero y Torregrosa en [6] como:

p ACOC = p = l n | ( x k + 1 x k ) / ( x k x k 1 ) | l n | ( x k x k 1 ) / ( x k 1 x k 2 ) | , k = 2 , 3 , . . .

2.2. Herramientas dinámicas

Las herramientas dinámicas son conceptos que se utilizarán para el estudio del análisis de estabilidad de un método ya que así se crearán mapeos en para visualizar el comportamiento de un método para converger a una solución en comparación con métodos más conocidos como Newton de la ecuación (1), el cual es bastante estable.

Blanchard en [2] establece las Definiciones 9-17:

Definición 9 (Punto fijo)

Dado un operador G : ̂ ̂ ; ̂ = { ± } definido sobre la esfera de Riemann, un punto fijo z ̂ es aquel que se mantiene invariante al operador G , o sea:

G ( z ̂ ) = z ̂

Por ejemplo, para una ecuación cualquiera f ( z ) = 0 , z 0 una de sus soluciones y G un operador para aproximar las soluciones de f ; es claro que G ( z 0 ) = z 0 . Por lo tanto se puede concluir que todas las raíces de una ecuación serán, eventualmente, puntos fijos del operador G .

Sin embargo, al resolver la ecuación G ( z ) = z podrían aparecer puntos fijos que no correspondan a ninguna raíz de la función f . Estos puntos se llaman puntos fijos extraños. Los puntos fijos extraños no son deseables desde el punto de vista numérico, ya que si se toma como estimación inicial un valor que esté cerca de alguno de estos, existe la posibilidad de que el método numérico converja a él, es decir, a un punto que no es una solución.

Definición 10 (Órbita de un punto)

Dada una función racional : ̂ ̂ definida sobre la esfera de Riemann, la órbita de un punto z 0 es la sucesión de puntos:

{ z 0 , k ( z 0 ) , R 2 ( z 0 ) , . . . , R n ( z 0 ) , . . . }

del plano complejo generada por el método iterativo.

Definición 11 (Punto k - periódico)

Dada una función racional : ̂ ̂ definida sobre la esfera Riemann, Un punto es k -periódico si se cumple que:

R k ( z ̂ ) = z ̂ y R p ( z ̂ ) z ̂ para p = 1 , 2 , . . . , k 1

Definición 12 (Estabilidad de puntos fijos)

Los puntos fijos se pueden clasificar según el comportamiento del operador derivada sobre ellos, así, un punto fijo z * de R puede ser:

Atractor, si | R ' ( z * ) | < 1

Superatractor, si | R ' ( z * ) | = 0

Repulsor, si | R ' ( z * ) | > 1

Parabólico o indiferente, si | R ' ( z * ) | = 1

Definición 13 (Estabilidad de puntos periódicos)

La estabilidad de los puntos k -periódicos se puede determinar calculando cada | R ' ( z i * ) | con i = 1 , . . . , k . Luego, se dice que la órbita periódica de periodo k es atractora si se cumple que:

| R ' ( z 0 * ) R ' ( z 1 * ) · · · R ' ( z k * ) | < 1

Análogamente, si la expresión es igual a cero se le llama superatractora, si es mayor que uno se llama repulsora y si es igual a uno se le llama indiferente.

Definición 14 (Cuenca de atracción)

La cuenca de atracción de un atractor z es el conjunto de preimágenes de cualquier orden:

A ( z ) = { z 0 ̂ : R n ( z 0 ) z , n }

La cuenca de atracción es el conjunto de todas las órbitas posibles de un atractor, por tanto indica la región de convergencia de un método iterativo.

Definición 15 (Conjunto de Fatou y de Julia)

Al conjunto de todas las cuencas de atracción se llama conjunto de Fatou, .

A la frontera entre las componentes conexas de las cuencas de atracción, se le llama conjunto de Julia, .

Estos conjuntos son complementarios, o sea: . = ̂

Definición 16 (Cuenca de atracción inmediata)

La cuenca de atracción inmediata de un punto fijo z * es la componente conexa de A ( z * ) que contiene a z * .

Definición 17 (Punto Crítico)

Un punto crítico de un operador R es un punto z * para el que la derivada de R en un punto fijo se anula, o sea:

R ' ( z * ) = 0

Eventualmente los puntos fijos que sean superatractores van a ser críticos, sin embargo podrían existir puntos que no sean fijos pero sí cumplan la anterior igualdad; a dichos puntos se le llaman puntos críticos libres.

Teorema 1 (Teorema de Fatou-Julia)

La cuenca de atracción inmediata de un punto k -periódico contiene, al menos, un crítico.

García en [8] explica que: “a la hora de realizar el análisis dinámico de un nuevo método es necesario aplicarlo sobre un problema concreto, es decir, sobre una ecuación no lineal específica; por eso, una manera natural de abordar el estudio dinámico de un método iterativo en cuestión es aplicándolo sobre polinomios de bajo grado”.

En este caso se hará el estudio sobre un polinomio de grado dos genérico, ya que es la función no lineal más sencilla.

Artidiello en [1] y García en [8] aseveran que, en la práctica, cuando se modela el método sobre ecuaciones no lineales sencillas, si el estudio dinámico sobre el método iterativo en el polinomio cuadrático se comporta de forma inestable o necesita muchas condiciones para que converja, entonces sobre funciones no lineales más complicadas también; por lo que, dependerá del comportamiento del método sobre dicho polinomio para definir cómo es el método y sacar conclusiones sobre su desempeño en general. Si bien es cierto, no se puede generalizar partiendo del estudio sobre funciones cuadráticas, en la práctica se puede verificar con las pruebas numéricas pues los resultados aplicados en funciones no lineales más complejas coinciden con la teoría estudiada en polinomios de grado dos.

Con estos conceptos se generará el estudio de estabilidad de los métodos a proponer.

2.3. Métodos multipasos

En esta sección se hará referencia de cómo se generaron algunos métodos multipasos, algunas gene- ralizaciones que realizaron diversos autores y deducciones para construir un nuevo método; mismas que se utilizarán para realizar la propuesta de nuevos esquemas iterativos.

Traub en [14] introduce la idea de métodos multipasos, propone un método multipaso utilizando el esquema de Newton sin cambiar de paso la derivada, con el fin de realizar una menor cantidad de evaluaciones funcionales para aumentar el orden de convergencia del método.

Método de Traub de orden 3 en [14]:

y k = x k f ( x k ) f ' ( x k )

x k + 1 = y k f ( y k ) f ' ( x k )

Como se mencionó antes, en este caso d = 3; luego 2d−1 = 4; siguiendo la conjetura de Kung y Traub en [11], este método no es óptimo ya que es de orden 3 y podría buscar un método de orden 4 con la misma cantidad de evaluaciones funcionales. A partir de dicha conjetura, Jarratt en [9] propone el siguiente método:

y k = x k 2 3 f ( x k ) f ' ( x k )

x k + 1 = y k 1 2 ( 3 f ' ( y k ) + f ' ( x k ) 3 f ' ( y k ) f ' ( x k ) ) f ( x k ) f ' ( x k ) (3)

Este método es óptimo ya que el número de evaluaciones funcionales es 3 y su orden de convergencia p es 4, cumpliendo así la conjetura de Kung y Traub ya que p=2d-1.

Artidiello en [1], llama función peso a aquella función que, bajo ciertas condiciones iniciales, el orden de convergencia puede aumentar gracias a multiplicar esta función. Chun en [5] propone el esquema presentado en (4) con una función peso H ( t k ) donde t k = f ( y k ) f ( x k ) .

y k = x k f ( x k ) f ' ( x k ) ,

x k + 1 = y k H ( t k ) f ( y k ) f ' ( x k ) (4)

Particularmente si la función H cumple ciertas condiciones, entonces el método es de orden 4, haciendo el método óptimo ya que utiliza únicamente tres evaluaciones funcionales.

El siguiente teorema es un caso particular al Teorema 3.3.1 demostrado por Artidiello en [1].

Teorema 2

Sea f : I una función suficientemente diferenciable en cada punto del intervalo abierto I tal que α I es una raíz simple de la ecuación no lineal f ( x ) = 0 . Sea H : cualquier función suficientemente diferenciable tal que H ( 0 ) = 1 , H ' ( 0 ) = 1 y | H ' ' ( 0 ) | < . Entonces, si la estimación inicial x 0 está suficientemente cerca de la solución α , la sucesión { x k } k 0 obtenida usando la expresión con función peso de Chun presentada en (4) converge a α y el orden de convergencia local es al menos 4 y su ecuación del error es:

e k + 1 = ( ( 5 H ' ' ( 0 ) 2 ) C 2 3 C 2 C 3 ) e k 4 + O ( e k 5 )

donde C j = 1 j ! f ( j ) ( α ) f ' ( α ) , j = 2 , 3 , . . . y e k = x k α , k .

La prueba de este teorema (elaboración propia) puede ser encontrada en la sección Anexos en la subsección A.1.

Algunos métodos iterativos de esta forma son:

El método de Chun en [4]:

y k = x k f ( x k ) f ' ( x k )

x k + 1 = y k f ( x k ) + 2 f ( y k ) f ( x k ) f ( x k ) f ' ( x k ) (5)

El método de Ostrowski en [14], que es un caso particular para β = 2 de la familia de King en [10]:

y k = x k f ( x k ) f ' ( x k ) ,

x k + 1 = y k f ( x k ) + ( 2 + β ) f ( y k ) f ( x k ) + β f ( y k ) f ( x k ) f ' ( x k ) ; β (6)

El método de Kung y Traub en [11]:

y k = x k f ( x k ) f ' ( x k ) ,

x k + 1 = y k f ( x k ) 2 ( f ( x k ) f ( y k ) ) 2 f ( x k ) f ' ( x k ) (7)

Método de Zhao et al. en [16]:

y k = x k f ( x k ) f ' ( x k )

x k + 1 = y k 1 + 2 t k + t k 2 1 4 t k 2 f ( x k ) f ' ( x k ) ; t k = f ( y k ) f ( x k ) (8)

Todos los métodos anteriores tienen la estructura de la familia de Chun presentado en (4), ya que poseen una función peso de la forma H ( t k ) con t k = f ( y k ) f ( x k ) .

En el caso del esquema de Chun presentado en (5), la función peso es:

H ( t k ) = 1 + 2 t k ; t k = f ( y k ) f ( x k )

En el caso del esquema de King presentado en (6), la función peso es:

H ( t k ) = 1 + ( 2 + β ) t k 1 + β t k ; t k = f ( y k ) f ( x k )

En el caso del esquema de Kung y Traub presentado en (7), la función peso es:

H ( t k ) = 1 ( 1 t k ) 2 ; t k = f ( y k ) f ( x k )

En el caso del esquema de Zhao presentado en (8), la función peso es:

H ( t k ) = 1 + 2 t k + t k 2 1 4 t k 2 ; t k = f ( y k ) f ( x k )

Basado en los esquemas anteriores, para generar un método iterativo con la estructura de la familia de Chun de (4), se debe tener el mismo esquema pero con una función peso que cumpla las condiciones del teorema 2. Para ello el método a proponer sigue la estructura anterior con funciones peso racionales pero el denominador posee una función cuadrática:

Familia de métodos MA1:

y k = x k f ( x k ) f ' ( x k )

x k + 1 = y k 1 1 2 t k + γ t k 2 f ( y k ) f ' ( x k ) ; t k = f ( y k ) f ( x k ) (9)

donde γ es un parámetro real. Esta es una familia de métodos ya que γ varía.

En general, ya existen escritos a cerca del estudio dinámico de la familia de Chun con función peso H ( t ) = α t + β γ t + δ (con α , β , y δ son parámetros reales), es decir, donde su numerador y su denominador sean lineales como el estudio generado en [1] por Artidiello.

Se propone realizar el estudio de la dinámica de una familia de métodos con una función peso racional con alguna de sus partes cuadrática. Es importante aclarar que si se escogiera que tanto el numerador como el denominador sean cuadráticas, habrían más parámetros a considerar, por lo que el estudio dinámico se extiende.

Se presenta la siguiente función peso H ( t ) = α 1 + β t + γ t 2 con α , β , γ parámetros reales. Sin embargo, otro detalle a considerar es que se debe cumplir ciertas condiciones sobre H descritas en el teorema 2 para que el método sea óptimo; específicamente que: H ( 0 ) = 1 , H ' ( 0 ) = 2 y | H ' ' ( 0 ) | < . Entonces:

H ( 0 ) = 1 α = 1

H ' ( 0 ) = 2 β = 2

Luego, la función peso propuesta queda de la forma: H ( t ) = 1 1 2 t + γ t 2 .

En la sección 3 del escrito se realiza el estudio de estabilidad del método MA1 para determinar cuáles valores de γ hacen que el método sea más eficiente.

Cabe destacar que este método es nuevo y no ha sido encontrado en literaturas anteriores.

2.3.1. Clases de conjugación

Artidiello en [1] explica que: Las clases de conjugación permiten simplificar los cálculos, trabajando simultáneamente con todos los polinomios cuadráticos.”

Definición 18 (Clases de conjugación)

Sean f y g dos funciones analíticas en la esfera de Riemann. Una conjugación analítica entre f y g es un difeomorfismo h en la esfera de Riemann tal que h f h - 1 = g .

2.3.2. Teorema del escalado y transformación afín

Las clases conjugadas y el Teorema del escalado, si se cumple para el método en estudio, permite generalizar el estudio dinámico sobre la función no lineal a la que se le hace el estudio, en este caso, a un polinomio cuadrático cualquiera; lo anterior quiere decir que gracias a este aspecto se puede garantizar que el estudio dinámico realizado contempla cualquier polinomio de grado dos.

Teorema 3 (Teorema del Escalado para la familia de Chun)

Sea f ( z ) una función analítica en la esfera de Riemann y T ( z ) = α z + β ¸ una aplicación afín con α 0 . Si g ( z ) = λ f T ( z ) , entonces ( T G g , γ T - 1 ) ( z ) = G f , γ ( z ) , es decir, G f , γ , es analíticamente conjugada a G g , γ , mediante T . Donde G f , γ = z f ( z ) f ' ( z ) H ( t f ( z ) , γ ) f ( z f ( z ) f ' ( z ) ) f ' ( z ) y G g , γ = z g ( z ) g ' ( z ) H ( g f ( z ) , γ ) g ( z g ( z ) g ' ( z ) ) g ' ( z ) son los operadores correspondientes a la familia de métodos de Chun (4) sobre f y g respectivamente, con f ( t f ( z ) , γ ) y H ( t g ( z ) , γ ) funciones suficientemente diferenciables tales que ambas funciones cumplen que H ( 0 , γ ) = 1 , H ' ( 0 , γ ) = 2 , | H ' ' ( 0 , γ ) | < y γ es un parámetro; con t f ( z ) = f ( z f ( z ) f ' ( z ) ) f ( z ) y t g ( z ) = g ( z g ( z ) g ' ( z ) ) g ( z ) . La prueba de este teorema (elaboración propia) puede ser encontrada en la sección Anexos en la subsección A.2.

Definición 19

Sean R 1 y R 2 funciones racionales en la esfera de Riemann. Se dice que R 1 y R 2 son conjugadas si existe una transformación de Möbius en la esfera de Riemann tal que: M R 1 M - 1 = R 2

Luego del cálculo de conceptos teóricos de herramientas dinámicas como: puntos fijos, puntos críticos, puntos periódicos y cuencas de atracción para el método MA1 descrito en (9); el teorema del escalado ayudará a realizar mapeos de planos dinámicos y zonas de estabilidad que se puedan interpretar mejor.

3. Análisis de Estabilidad de MA1

Para comparar el método MA1 presentado en (9) con algunos métodos conocidos, se realizará un estudio de estabilidad sobre el plano complejo para ecuaciones no lineales específicas y así visualizar qué tan eficiente es el método.

Se analizará el método sobre un polinomio genérico de grados dos p ( z ) = ( z a ) ( z b ) con la función peso H ( t ) = 1 1 2 t + γ t 2 , con γ un parámetro.

y ( z ) = z p ( z ) p ' ( z )

O f ( z ) = y ( z ) H ( t ) p ( y ( z ) ) p ' ( z )

= a b z 2 a + b 2 z + ( a z ) 2 ( b z ) 2 ( a + b 2 z ) 3 ( 1 + ( a z ) ( b + z ) ( 2 ( a + b 2 z ) 2 + ( a z ) ( b + z ) γ ) ( a + b 2 z ) 4 ) (10)

con t = p ( y ( z ) ) p ( z )

Dicho resultado es el operador de la función p que depende de a y de b.

Ahora se realizará el conjugado del operador anterior con la Transformación de Möbius, sobre el operador para simplificarlo ya que manda las soluciones del polinomio a y b a 0 e para que el resultado ya no dependa de a ni de b.

M ( z ) = z a z b

O p ( z ) = M ( O f ( M - 1 ( z ) ) ) = z 4 ( ( 1 + z ) 2 + γ ) 1 + z ( 2 + z + z γ )

Gracias a la Definición18 a cerca de las clases de conjugación, el Teorema 3 del escalado y la Definición 19, se puede asegurar que O p y O f son analíticamente conjugados por lo que su dinámica asociada es equivalente, permitiendo así realizar el estudio en casos más sencillos. Se calcularán los valores de γ para los que O p se pueda simplificar.

Igualando el numerador del operador O p a cero se obtiene:

s n 1 ( γ ) = 1 - γ

s n 2 ( γ ) = 1 + - γ

Luego igualando el denominador del operador O p a cero se obtiene:

s d 1 ( γ ) = 1 1 + γ

s d 2 ( γ ) = 1 1 + γ

Se iguala cada término s n i ( γ ) con cada término s d i ( γ ) ; con i { 1 , 2 } , para saber en qué valores de γ dichos términos son iguales y el operador O p se pueda simplificar. Haciendo todas las igualdades posibles se obtiene que O p se puede simplificar si γ = 0 , γ = ± 2 i , γ = 1 o si γ = 2 . Estos valores de γ son los que se utilizarán para la variación del número de puntos fijos extraños y para las pruebas numéricas.

3.1. Puntos fijos

Como se ha mencionado anteriormente, los puntos fijos del operador O p ( z ) son los puntos tales que O p ( z ) = z ; eventualmente lo cumplen z = 0 y z = , pues son las raíces del polinomio p ( z ) luego de la transformación de Möbius.

Para verificar, se resuelve la expresión O p ( z ) = z y una de sus raíces es 0. Para corroborar que es un punto fijo se calcula el operador inverso de la siguiente manera:

z 1 = 1 O p ( 1 z ) = z 4 ( 1 + 2 z + z 2 + γ ) 1 + 2 z + z 2 ( 1 + γ )

y evaluando en z = 0 su imagen es 0 por lo que también es fijo.

Sin embargo, al resolver la ecuación O p ( z ) = z , aparecen otras 5 soluciones las cuales son puntos fijos extraños:

1. p f 0 ( γ ) = 1

2. p f 1 ( γ ) = 1 4 ( 3 1 4 γ 6 4 γ + 6 1 4 γ )

3. p f 2 ( γ ) = 1 4 ( 3 1 4 γ + 6 4 γ + 6 1 4 γ )

4. p f 3 ( γ ) = 1 4 ( 3 + 1 4 γ 6 4 γ 6 1 4 γ )

5. p f 4 ( γ ) = 1 4 ( 3 + 1 4 γ + 6 4 γ 6 1 4 γ )

Como se puede observar, los puntos fijos extraños dependen de γ , por lo que su estabilidad, así como los puntos críticos, también. Se analizará la variación de números de puntos fijos extraños, la estabilidad y puntos críticos independientes para valores del parámetro γ específicos.

Lema 1

El número de puntos fijos extraños del operador O p ( z ) = z 4 ( ( 1 + z ) 2 + γ ) 1 + z ( 2 + z + z γ ) es 5 excepto en el siguiente caso: si γ = 0 entonces hay 3 puntos fijos extraños. La prueba de este lema (elaboración propia) puede ser encontrada en la sección Anexos en la subsección A.5.

Existen otros valores de γ que hacen que el operador O p se simplifique pero mantiene la misma cantidad de puntos fijos extraños:

- Si γ = ± 2 i , entonces el operador es O p ( z ) = z 4 ( ± 2 i + ( 1 + z ) 2 ) 1 + z ( 2 + ( 1 ± 2 i ) z ) y al resolver O p ( z ) = z se obtienen los siguientes puntos fijos extraños:

z = 1

z = 1 4 ( 3 1 8 i ( 6 8 i ) + 6 1 8 i )

z = 1 4 ( 3 1 8 i + ( 6 8 i ) + 6 1 8 i )

z = 1 4 ( 3 + 1 8 i ( 6 8 i ) 6 1 8 i )

z = 1 4 ( 3 + 1 8 i + ( 6 8 i ) 6 1 8 i )

Si γ = 1 , entonces el operador es O p ( z ) = z 5 ( 2 + z ) 1 + 2 z y al resolver O p ( z ) = z se obtienen los siguientes puntos fijos extraños

z = 1

z = 1 4 ( 3 5 2 + 6 5 )

z = 1 4 ( 3 5 + 2 + 6 5 )

z = 1 4 ( 3 5 2 6 5 )

z = 1 4 ( 3 5 + 2 6 5 )

Si γ = 2 , entonces el operador es O p ( z ) = z 4 ( 2 + ( 1 + z ) 2 ) 1 ( 2 + z ) z ) y al resolver O p ( z ) = z se obtienen los siguientes puntos fijos extraños

z = 1

z = 1 2 ( 3 5 )

z = 1 2 ( 3 + 5 )

z = i

z = i

3.2. Estabilidad de los puntos fijos

La importancia al obtener los puntos fijos es poder hacer un estudio de su estabilidad para saber cuándo son valores atractores, repulsores o indiferentes dependiendo del valor del parámetro γ . Para calcular la estabilidad de los puntos fijos se procede a evaluar en la derivada del operador, o sea, en O p ' ( z ) .

Teorema 4

Los puntos fijos z = 0 y z = son superatractores independientemente del valor de γ , y para z = 1 se obtienen los siguientes resultados:

1. Atractor si | γ + 2 8 3 | < 8 en particular si γ = 8 .

2. Indiferente si | γ + 2 8 3 | = 8

3. Repulsor si | γ + 2 8 3 | > 8

La estabilidad de z = 1 se puede observar en el plano complejo en la Figura2.

Estabilidad del punto fijo extraño pf0(γ)
Figura 2:
Estabilidad del punto fijo extraño pf0(γ)

Visualmente, si γ está dentro de la circunferencia del cono, entonces es atractor. Si γ está sobre la circunferencia, es indiferente y si está fuera de la circunferencia entonces es repulsor.

Teorema 5

Los puntos fijos extraños p f 1 ( γ ) y p f 2 ( γ ) son:

1. Atractores si 0 . 2 3 8 2 2 8 7 3 < γ < 0 . 2 5 ; en particular superatractor si γ 0 . 2 4 6 2 1 1 1 9

2. Indiferentes si γ 0 . 2 3 8 2 2 8 7 3 , γ 0 . 2 5

3. Repulsores para valores de γ diferentes de los anteriores

Mientras que los puntos fijos extraños p f 3 ( γ ) y p f 4 ( γ ) son:

1. Atractores si 2 1 < γ < 1 2 ; en particular superatractor si γ 1 6 . 2 4 6 5

2. Indiferentes si γ 2 1 , γ 1 2

3. Repulsores para valores de γ diferentes de los anteriores.

No se puede comprobar de forma analítica el teorema ya que al sustituir estos puntos fijos extraños en la funciones de estabilidad, aparecen raíces de polinomios de alto grado. El análisis numérico de este teorema (elaboración propia) puede ser encontrado en la sección Anexos en la subsección A.4.

La estabilidad de z = pf.(γ) y de z = pf.(γ) se puede visualizar en el plano complejo en la Figura 3a.

La estabilidad de z = pf.(γ) y de z = pf.(γ) se puede visualizar en el plano complejo en la Figura 3b.

Los valores para los que Op. γ cambia su comportamiento, es decir, pasa de ser atractor a repulsor, se han obtenido numéricamente por lo que son valores aproximados. En las figuras anteriores se pueden apreciar las regiones donde el parámetro varía su estabilidad, si γ está dentro de esas cuencas entonces el método MA1 será inestable.

3.3. Puntos críticos

Los puntos críticos, como se ha mencionado antes, se obtienen de resolver la ecuación O p ' ( z ) = 0 . La importancia de los puntos críticos del operador es que, además de los puntos fijos que son superatractores, se pueden obtener puntos que no son ni las raíces ni puntos fijos extraños, a estos se les llaman puntos críticos libres.

García [8] menciona que estos puntos van a tener interés ya que cada cuenca de atracción tiene al menos un punto crítico, por lo que los críticos libres podrían estar en una cuenca de atracción de alguna de las soluciones de la ecuación, o bien estar en la cuenca de algún punto fijo extraño u órbita periódica atractora.

Los puntos críticos del operador O p ( z ) son z = 0 y z = (corresponden a la transformación de Möbius de las raíces del polinomio) ya que son puntos fijos superatractores y los críticos libres:

1. c r 1 ( γ ) = 8 + 3 γ 8 ( 1 + γ ) 1 7 γ 2 8 γ 3 8 1 + 2 γ + γ 2 1 2 ( 8 + 3 γ ) 2 8 ( 1 + γ ) 2 4 + 4 γ 2 ( 1 + γ ) 1 2 + 2 γ + γ 2 2 ( 1 + γ ) λ

2. c r 2 ( γ ) = 8 + 3 γ 8 ( 1 + γ ) 1 7 γ 2 8 γ 3 8 1 + 2 γ + γ 2 + 1 2 ( 8 + 3 γ ) 2 8 ( 1 + γ ) 2 4 + 4 γ 2 ( 1 + γ ) 1 2 + 2 γ + γ 2 2 ( 1 + γ ) λ

3. c r 3 ( γ ) = 8 + 3 γ 8 ( 1 + γ ) + 1 7 γ 2 8 γ 3 8 1 + 2 γ + γ 2 1 2 ( 8 + 3 γ ) 2 8 ( 1 + γ ) 2 4 + 4 γ 2 ( 1 + γ ) 1 2 + 2 γ + γ 2 2 ( 1 + γ ) + λ

4. c r 4 ( γ ) = 8 + 3 γ 8 ( 1 + γ ) + 1 7 γ 2 8 γ 3 8 1 + 2 γ + γ 2 + 1 2 ( 8 + 3 γ ) 2 8 ( 1 + γ ) 2 4 + 4 γ 2 ( 1 + γ ) 1 2 + 2 γ + γ 2 2 ( 1 + γ ) + λ

con λ = 1 + 2 γ + γ 2 ( 4 ( 8 + 3 γ ) 1 + γ ( 8 + 3 γ ) 3 8 ( 1 + y ) 3 + ( 8 + 3 γ ) ( 1 2 + 2 γ + γ 2 ) ( 1 + γ ) 2 ) 1 7 γ 2 8 γ 3

Particularmente c r 2 ( γ ) = 1 c r 1 ( γ ) y c r 4 ( γ ) = 1 c r 3 ( γ )

Lema 2

El número de puntos críticos extraños del operador es 4, excepto en los siguientes casos:

1. Si γ = 0 , entonces no hay puntos críticos libres.

2. Si γ = 1 7 8 , entonces hay 2 puntos críticos libres, ambos de multiplicidad dos.

3. Si γ = 8 , entonces hay 3 puntos críticos libres, uno de ellos de multiplicidad dos.

4. Si γ = 4 , entonces hay 2 puntos críticos libres.

La prueba de este lema (elaboración propia) puede ser encontrada en la sección Anexos en la subsección A.5.

Para determinar qué valores de γ hacen que varíe el número de puntos críticos libres, se procede a resolver todas las ecuaciones posibles: c r i ( γ ) = c r j ( γ ) ; con i , j { 1 , 2 , 3 , 4 } y i j . Resolviendo todas las ecuaciones posibles, se obtiene que el operador O p ' se simplifica si γ = 0 , γ = 1 7 8 , γ = 8 o si y = 4 . Algunos de estos valores también serán utilizados para las pruebas numéricas.

3.4. Planos de parámetros

Los planos de parámetros son gráficos de los valores críticos independientes del método que permiten determinar de manera visual qué valores de γ harán que el método sea estable o inestable.

El plano de parámetros se genera partiendo con un punto crítico independiente como valor inicial para el operador O p y se aplica el método iterativo MA1 para todos los valores del parámetro definidos en una malla sobre el plano complejo. Si en un punto, donde γ es un valor específico, el método converge en menos de 50 iteraciones a alguna de las raíces del polinomio z = 0 o z = entonces a dicho punto se le asigna el color rojo; en otro caso entonces al punto se le asigna el color negro.

Particularmente para el método MA1 se obtienen dos planos paramétricos: uno en donde se visualiza el comportamiento del método con alguno de los puntos críticos c r 1 ( γ ) o c r 2 ( γ ) y otro con alguno de los puntos críticos c r 3 ( γ ) o c r 4 ( γ ) esto debido a que c r 2 ( γ ) = 1 c r 1 ( γ ) y c r 4 ( γ ) = 1 c r 3 ( γ ) por lo que solo hay 2 puntos críticos independientes.

Plano de parámetros del método MA1
Figura 4
Plano de parámetros del método MA1

Los puntos del plano complejo que aparecen en negro en las Figuras 4a y 4b corresponden a valores del parámetro γ en donde el método iterativo MA1 no converge a las raíces del polinomio cuadrático, tomando como estimación inicial alguno de los puntos críticos libres c r i ( γ ) con i { 1 , 2 , 3 , 4 } . Significa entonces que existen, al menos, tres cuencas de convergencia, la del 0, la del y otra (o más) de algún punto fijo extraño atractor o alguna órbita periódica atractora.

Las regiones de los círculos negros que se pueden apreciar a la izquierda de la figura4bcorrespon- den con los diagramas de estabilidad de los puntos fijos extraños p f i ( γ ) con i { 1 , 2 , 3 , 4 } que se pueden apreciar en las Figuras2-3b, donde dichos puntos son atractores. Las otras regiones en negro corresponden a diagramas de estabilidad de órbitas periódicas atractoras.

3.5. Planos dinámicos

Los planos dinámicos son otra forma de ampliar la información obtenida con los planos paramétricos. Con los planos dinámicos se pueden visualizar las cuencas de atracción de puntos fijos o periódicos del método MA1 con algún valor del parámetro γ en particular.

Para construir un plano dinámico se define un mallado de puntos en el plano complejo, cada uno de los cuales se toma como estimación inicial del método iterativo, si partiendo desde un punto específico del mallado el método iterativo converge a algún cero de la función, entonces se le asignará un determinado color, particularmente de color naranja si converge a z = 0 y color azul si converge a z = ; en caso de que el punto tomado como estimación inicial no converja a ninguna de las raíces de la función en un máximo de iteraciones establecido, se le asignará el color negro.

Todas las gráficas son del método MA1 sobre un polinomio cuadrático con la Transformación de Möbius.

En las Figuras 5c-5g se puede observar los planos dinámicos de valores de γ que se encuentran en la zona de estabilidad. Esto se puede verificar visualmente ya que en los gráficos el método solo converge a las soluciones, las cuales en este caso luego de aplicar la Transformación de Möbius son 0 e .

En las Figuras 5h-5l se puede observar los planos dinámicos de valores de γ que se encuentran fuera de la zona de estabilidad. Esto se puede verificar visualmente ya que en los gráficos el método converge a otros valores además de la solución, dando como resultado zonas negras.

En la Figura 5a se puede observar el plano dinámico del método de Ostrowski presentado en (6) aplicado a un polinomio cuadrático con la Transformación de Möbius, y en la Figura 5b se puede observar el plano dinámico del método de Kung y Traub de la ecuación 7, que se obtiene de MA1 si γ = 1 . Estos valores de γ , como es de esperar, están en la zona de estabilidad; es esperable ya que generan métodos conocidos que se sabe de antemano que son esquemas estables.

En la Figura 5l se pueden observar las cuencas de atracción correspondientes a γ = 4 , 7 + 5 , 4 i ; este valor es importante de resaltar ya que en el plano dinámico que se obtiene aparece una órbita periódica de periodo 3; con este resultado, gracias al Teorema de Sarkovski (ver [7]), se puede probar que existen órbitas periódicas de cualquier periodo.

3.6. Pruebas numéricas

Las pruebas numéricas que aparecen en el siguiente apartado han sido calculadas utilizando aritmé- tica de precisión variable, con 7000 dígitos de mantisa y una tolerancia máxima de 1 0 - 3 0 0 en Matlab R2015a en una computadora portátil con procesador Intel(R) Core(TM) i5-7200U CPU@2.50GHz 2.70GHz con una memoria RAM de 8GB.


Plano dinámico del método MA1 para distintos valores de γ.
Figura 5
Plano dinámico del método MA1 para distintos valores de γ.

En este apartado se presentarán resultados obtenidos al utilizar el método familia MA1 con algunos valores γ específicos para estimar ceros de las siguientes funciones:

f 1 ( x ) = s e n ( x ) x 2 + 1 , = 1 . 4 0 9 6 2 4

f 2 ( x ) = ( x 1 ) 3 1 , = 2

f 3 ( x ) = a r c t a n ( x ) , = 0

f 4 ( x ) = ( s e n ( x ) x 2 ) 2 , = 0 , doble

Se comparará el método familia MA1 con otros métodos iterativos conocidos: el esquema de Newton presentado en (1), el esquema de Jarratt presentado en (3), el esquema de Chun presentado en (5), el esquema de Ostrowski presentado en (6), el esquema de Kung y Traub presentado en (7) y el esquema de Zhao presentado en (8).

En las siguientes tablas se mostrarán los resultados obtenidos de funciones anteriormente mencionadas, teniendo en cuenta que el criterio para que el código se detenga es, en todos los métodos, | x k + 1 x k | + | f ( x k + 1 ) | < 1 0 - 3 0 0 así se puede verificar si la sucesión { x n } converge y su límite es la solución de la ecuación no lineal f ( x ) = 0 .

Para cada una de las funciones se mostrará la estimación inicial utilizada x 0 , la raíz la que tiende el proceso, las estimaciones del error cometido en la última iteración: | x k + 1 x k | y | f ( x k + 1 ) | , el número de iteraciones que ha necesitado para converger y ρ como el orden de convergencia computacional aproximado ACOC.

Tabla 1
Resultados obtenidos por los métodos iterativos sobre la función f1
Resultados obtenidos por los métodos iterativos sobre la función f1

En las Tablas 1 y 2 se puede observar que en general, la mayoría de los métodos convergen a la solución de la ecuación; el número de iteraciones es en general el mismo y el orden de convergencia computacional aproximado coincide con el orden teórico.

En la Tabla 2 se puede observar que el método MA1 diverge cuando γ = 8 y cuando γ = 2 2 . Este resultado es esperable ya que en el análisis realizado dichos valores de γ se encuentran fuera de la zona de estabilidad de MA1.

En la Tabla 3 se puede observar que en general, la mayoría de métodos divergen, esto se debe a que, como todos son una composición del método de Newton; dicho método tiene problemas en el cálculo de la derivada cerca del cero de la función arco-tangente. Sin embargo, aquellas funciones que convergen, lo hacen con un número bajo de iteraciones y su ACOC es mayor al teórico.

Tabla 2
Resultados obtenidos por los métodos iterativos sobre la función f2
Resultados obtenidos por los métodos iterativos sobre la función f2

Tabla 3
Resultados obtenidos por los métodos iterativos sobre la función f3
Resultados obtenidos por los métodos iterativos sobre la función f3

Algo particular es que en la Tabla 4, la gran mayoría de métodos convergen lentamente: esto se sabe ya que su ACOC es igual a 1, en otros casos el método diverge; esto es esperable ya que la función posee ceros con multiplicidad lo cual genera problemas para los métodos iterativos usuales. Sin embargo, el método MA1 cuando γ = 4 realizó un número considerablemente menor de iteraciones a los demás, con un orden de convergencia mayor.

De repente se puede pensar que dicho método funciona para funciones con multiplicidad; sin embargo, luego de realizar pequeñas pruebas con ecuaciones como x 2 = 0 , el método diverge; una corriente de trabajo podría someterse al método para γ = 4 para funciones con multiplicidad para determinar en qué condiciones converge.

Tabla 4
Resultados obtenidos por los métodos iterativos sobre la función f4
Resultados obtenidos por los métodos iterativos sobre la función f4

4. Conclusiones

En el escrito se realizó un análisis para el diseño de métodos iterativos multipasos para la resolución de ecuaciones no lineales con el fin de encontrar métodos más eficientes en comparación con esquemas usualmente utilizados.

El orden de convergencia no es lo único que define cuán eficiente y estable es un método iterativo; la eficiencia, el índice de eficiencia y el orden de convergencia computacional generan una mejor visión de la eficacia del método en comparación con otros, pues entre mayor sean estos valores, hay mayor certeza que el método converja en un menor número de iteraciones y en un menor tiempo de ejecución.

Además un orden de convergencia mayor no significa un mejor método, para esto entra el papel de sistemas dinámicos. Por ejemplo, si se tiene un método con un orden de convergencia alto, entonces, garantiza que converja en menos iteraciones que otro método de un orden menor. Sin embargo, cabe la posibilidad de que este método tenga un radio de convergencia muy pequeño; o sea, que el valor inicial debe estar muy cercano a la solución de la ecuación para garantizar convergencia, caso contrario el método diverge o converge a valores que no sean soluciones; en otras palabras, prácticamente se debe conocer la solución de la ecuación y esto no es útil pues en la práctica la solución no se conoce y es lo que se busca. Con esto, se concluye la necesidad de planos dinámicos para concluir también si el método es viable o no.

En contraste con la teoría sobre dinámica compleja, la práctica numérica dio resultados favorables sobre los métodos en estudio: algunos valores de γ que se obtuvieron de la zona de convergencia (zona roja) de los planos de parámetros hicieron que el método converja, mientras que algunos valores de γ de la zona de divergencia (zona negra) hicieron que el método fallara.

Se logró encontrar valores de γ para los que el método MA1 es mejor que los métodos iterativos conocidos. Particularmente γ = 1 fue un valor bueno en el sentido de que requirió de un menor número de iteraciones en comparación con otros valores de γ y métodos conocidos.

Referencias

[1] S. Artidiello, Diseño, implementación y convergencia de métodos iterativos para resolver ecuaciones y sistemas no lineales utilizando funciones peso. Diss. Ph. D. thesis, Universitat Politècnica de València, 2014.

[2] P. Blanchard, The Dynamics of Newton’s method. Proceedings of Symposia in Applied Mathematics 49 (1994) 139-152.

[3] S. C. Chapra, y R. P. Canale, Métodos numéricos para ingenieros. Quinta Edición, McGraw-Hill, 2017.

[4] C. Chun, Construction of Newton-like iterative methods for solving nonlinear equations, Numerische Mathematik 104 (2006) 297-315.

[5] C. Chun, Some fourth-roder iterative methods for solving nonlinear equation, Applied Mathematica and Computation 196 (2008) 454-459.

[6] A. Cordero, J.R. Torregrosa, Variants of Newton’s method using fifth-order quadrature formulas, Appl. Math. Comput. 190 (2007) 686-698.

[7] R. L. Devaney The Mandelbrot Set, the Farey Tree and the Fibonacci sequence. Am.Math. Monthly, 106(4) (1999) 289-302.

[8] J. García, Análisis dinámico y numérico de familias de métodos iterativos para la resolución de ecuaciones no lineales y su extensión a espacios de Banach, Departamento de Matemática Aplicada, Universitat Politècnica de València, 2017.

[9] P. Jarratt, Some fourth order multipoint iterative methods for solving equations, Math.Comp. 20 (1966) 434-437.

[10] R. F. King, A family of fourth-order methods for nonlinear equations, SIAM Journal Numerical Analisis 10(5) (1973) 876-879.

[11] H.T. Kung, J.F. Traub, Optimal order of one-point and multipoint iteration, J. Assoc.Comput. Math. 21 (1974) 643-651.

[12] A.M. Ostrowski, Solution of equations and systems of equations, Academic Press, New York, 1960.

[13] K. Rajagopal, S. T. Kingni, G. F. Kuiate, V. K. Tamba, y V. Pham, Autonomous Jerk Osci- llator with Cosine Hyperbolic Nonlinearity: Analysis, FPGA Implementation, and Syn- chronization, Advances in Mathematical Physics, vol. 2018.

[14] J.F. Traub, Iterative methods for the solution of equations, Prentice-Hall, Englewood Cliffs, New Jersey, 1964.

[15] S. Weerakoon, T.G.I. Fernando, A variant of Newton’s method with accelerated third-order convergence, Appl. Math. Letters 13 (2000) 87-93.

[16] L. Zhao, X. Wang, W. Guo, New families of eighth-order methods with high efficiency index for solving nonlinear equations, Wseas Transactions on Mathematics 11 (2012) 283-293.

Anexos

A continuación se presentan las pruebas de algunos teoremas y lemas descritos en el trabajo. Cabe destacar que todas las demostraciones son de elaboración propia, sin embargo fueron hechas utilizando técnicas desarrolladas por los autores: Artidiello en [1] y García en [8].

A.1 Prueba del Teorema 2

Teorema:

Sea f : I una función suficientemente diferenciable en cada punto del intervalo abierto I tal que α I es una raíz simple de la ecuación no lineal f ( x ) = 0 . Sea H : cualquier función suficientemente diferenciable tal que H ( 0 ) = 1 , H ' ( 0 ) = 2 y | H ' ' ( 0 ) | < . Entonces, si la estimación inicial x 0 está suficientemente cerca de la solución α , la sucesión { x k } k 0 obtenida usando la expresión con función peso de Chun presentada en (4) converge a α y el orden de convergencia local es al menos 4 y su ecuación del error es:

e k + 1 = ( ( 5 H ' ' ( 0 ) 2 ) C 2 3 C 2 C 3 ) e k 4 + O ( e k 5 ) ,

donde C j = 1 j ! f ( j ) ( α ) f ' ( α ) , j = 2 , 3 , . . . y e k = x k α , k .

Demostración.

Sea α un cero simple de la función f (i.e., f ( α ) = 0 y f ' ( α ) 0 ); x k = α + e k . Para probar la convergencia local del método iterativo hacia la solución de f ( x ) = 0 se utiliza el desarrollo en serie de Taylor de la función alrededor de la solución:

f ( x k ) = f ( α + e k ) = f ( α ) + f ' ( α ) e k + 1 2 f ' ' ( α ) e k 2 + 1 3 ! f ' ' ' ( α ) e k 3 + 1 4 ! f ( 4 ) ( α ) e k 4 + O e k 5

= f ' ( α ) [ e k + 1 2 f ' ' ( α ) f ' ( α ) e k 2 + 1 3 ! f ' ' ' ( α ) f ' ( α ) e k 3 + 1 4 ! f ( 4 ) ( α ) f ' ( α ) e k 4 + O ( e k 5 ) ]

= f ' ( α ) [ e k + C 2 e k 2 + C 3 e k 3 + C 4 e k 4 + O ( e k 5 ) ]

Procediendo igual para determinar f ' se tiene que:

f ' ( x k ) = f ' ( α ) [ 1 + 2 C 2 e k + 3 C 3 e k 2 + 4 C 4 e k 3 + O ( e k 4 ) ]

Luego:

f ( x k ) f ' ( x k ) = f ' ( α ) [ e k + C 2 e k 2 + C 3 e k 3 + C 4 e k 4 + O ( e k 5 ) ] f ' ( α ) [ 1 + 2 C 2 e k + 3 C 3 e k 2 + 4 C 4 e k 3 + O ( e k 4 ) ]

= e k C 2 e k 2 + ( 2 C 2 2 2 C 3 ) e k 3 + [ 8 C 2 3 + C 2 ( 4 C 2 2 3 C 3 ) + 1 0 C 2 C 3 3 C 4 ] e 4 + O ( e k 5 ) 1 1 )

luego, como y k = x k f ( x k ) f ' ( x k ) ; entonces por la ecuación (11) y como x k = α + e k se sigue que:

y k = x k f ( x k ) f ' ( x k ) = α + C 2 e k 2 + ( 2 C 2 2 + 2 C 3 ) e k 3 + ( 4 C 2 3 7 C 2 C 3 + 3 C 4 ) e k 4 + O ( e k 5 )

f ( y k ) = f ' ( α ) [ C 2 e k 2 2 ( C 2 2 C 3 ) e k 3 + ( 5 C 2 3 7 C 2 C 3 + 3 C 4 ) e k 4 + O ( e k 5 )

Dado t k = f ( y k ) f ( x k ) ; como f ( y k ) tiende a cero más rápido que f ( x k ) , se puede asegurar que t k 0 cuando k .

Entonces el desarrollo de serie Taylor alrededor de 0 para una función H que depende de t k sería:

H ( t k ) = H ( 0 ) + H ' ( 0 ) t k + 1 2 H ' ' ( 0 ) t k 2

luego la ecuación del error es la siguiente

e k + 1 = y k H ( t k ) · f ( y k ) f ' ( x k )

= ( 1 + H ( 0 ) ) C 2 e k 2 + ( ( 2 + 4 H ( 0 ) H ' ( 0 ) ) C 2 2 2 ( 1 + H ( 0 ) ) C 3 ) e k 3

+ [ ( 4 1 3 H ( 0 ) + 7 H ' ( 0 ) H ' ' ( 0 ) 2 ) C 2 3 + ( 7 + 1 4 H ( 0 ) 4 H ' ( 0 ) ) C 2 C 3 3 ( 1 + H ( 0 ) ) C 4 ] e k 4 + O ( e k 5 )

Pero como H ( 0 ) = 1 y H ' ( 0 ) = 2 , entonces la ecuación del error se reduce a lo siguiente:

e k + 1 = ( ( 5 H ' ' ( 0 ) 2 ) C 2 3 C 2 C 3 ) e k 4 + O ( e k 5 )

Y como | H ' ' ( 0 ) | < , entonces el método es de, al menos, orden 4 y óptimo ya que d = 3 tiene evaluaciones iterativas y p = 2 d 1 = 4 .

A.2. Prueba del Teorema 3

Teorema:

Sea f ( z ) una función analítica en la esfera de Riemann y T ( z ) = α z + β ¸ una aplicación afín con α 0 . Si g ( z ) = λ f T ( z ) , entonces ( T G g , γ T - 1 ) ( z ) = G f , γ ( z ) , es decir, G f , γ , es analíticamente conjugada a G g , γ ,mediante T . Donde G f , γ = z f ( z ) f ' ( z ) H ( t f ( z ) , γ ) f ( z f ( z ) f ' ( z ) ) f ' ( z ) y son los operadores correspondientes a la familia de métodos de Chun (4) sobre f y g respectivamente, con H ( t f ( z ) , γ ) y H ( t g ( z ) , γ ) funciones suficientemente diferenciables tales que ambas funciones cumplen que H ( 0 , γ ) = 1 , H ' ( 0 , γ ) = 2 , | H ' ' ( 0 , γ ) | < y γ es un parámetro; con t f ( z ) = f ( z f ( z ) f ' ( z ) ) f ( z ) y t g ( z ) = g ( z g ( z ) g ' ( z ) ) g ( z ) .

Demostración:

Teniendo en cuenta que T ( x y ) = T ( x ) T ( y ) + β y g ' ( z ) = α λ f ' ( T ( z ) ) entonces:

g ( T - 1 ( z ) ) = λ f ( T ( T - 1 ( z ) ) ) = λ f ( z )

g ' ( T - 1 ( z ) ) = α λ f ( T ( T - 1 ( z ) ) ) = α λ f ( z )

T ( T - 1 ( z ) f ( z ) α f ' ( z ) ) = T ( T - 1 ( z ) ) T ( f ( z ) α f ' ( z ) ) + β

= z α f z α f ' z β + β

= z f ( z ) f ' ( z )

t g ( T - 1 ( z ) ) = g ( T - 1 ( z ) g ( T - 1 ( z ) ) g ' ( T - 1 ( z ) ) ) g ( T - 1 ( z ) )

= λ f ( T ( T - 1 ( z ) f ( z ) α f ' ( z ) ) λ f ( z )

= f ( z f ( z ) f ' ( z ) ) f ( z )

= t f ( z )

Haciendo uso de los resultados anteriores entonces:

( T G g , γ T - 1 ) = T ( G g , γ ( T - 1 ( z ) ) )

= T ( T - 1 ( z ) g ( T - 1 ( z ) ) g ' ( T - 1 ( z ) ) H ( t f ( T - 1 ( z ) ) , γ ) g ( T - 1 ( z ) g ( T - 1 ( z ) ) g ' ( T - 1 ( z ) ) ) g ' ( T - 1 ( z ) )

= T ( T - 1 ( z ) λ f ( z ) α λ f ' ( z ) H ( t f ( z ) , γ ) λ f ( z f ( z ) f ' ( z ) ) α λ f ' ( z )

= T ( T - 1 ( z ) f ( z ) α f ' ( z ) T ( H ( t f ( z ) , γ ) f ( z f ( z ) f ' ( z ) ) α f ' ( z ) ) + β

= z f ( z ) α f ' ( z ) α H ( t f ( z ) , γ ) f ( z f ( z ) f ' ( z ) ) α f ' ( z ) β + β

= G f , γ

Por lo tanto, ( T G g , γ T 1 ) ( z ) = G f , γ ( z ) , es decir, G f , γ y G g , γ son analíticamente conjugadas mediante T ( z ) .

A.3. Prueba del Teorema 4

Teorema:

Los puntos fijos z = 0 y z = son superatractores independientemente del valor de γ , y para z = 1 se obtienen los siguientes resultados:

1. Atractor si | γ + 2 8 3 | < 8 en particular superatractor si γ = 8 .

2. Indiferente si | γ + 2 8 3 | = 8

3. Repulsor si | γ + 2 8 3 | > 8

Demostración:

Evaluando en 0 se obtiene como imagen 0, por lo que 0 es superatractor. Para el caso del se derivala inversa del operador O p , o sea en ( 1 O p ( 1 z ) ) ' y se evalúa en z = 0 ; el resultado es 0, por lo que el es superatractor.

Ahora se evaluará en O p ' ( z ) el punto fijo extraños p f 0 = 1 para poder clasificar la estabilidad del mismo.

O p ' ( 1 ) = 2 γ + 1 6 4 + γ

Entonces si se tiene | 2 γ + 1 6 4 + γ | < 1 implica que | 2 γ + 1 6 | < | 4 + γ | , como γ , entonces dados a , b arbitrarios tales que γ = a + b i , entonces:

| 2 a + 1 6 + 2 b i | < | 4 + a + b i |

( 2 a + 1 6 ) 2 + ( 2 b ) 2 < ( 4 + a ) 2 + b 2

3 a 2 + 5 6 a + 2 4 0 + 3 b 2 < 0

( a + 2 8 3 ) 2 + b 2 < 6 4

Por lo tanto O p ' ( 1 ) < 1 , si y solo si | γ + 2 8 3 | < 8

A.4 Análisis numérico del Teorema 5

Teorema:

Los puntos fijos extraños p f 1 ( γ ) y p f 2 ( γ ) son:

1. Atractores si 0 . 2 3 8 2 2 8 7 3 < γ < 0 . 2 5 ; en particular superatractor si γ 0 . 2 4 6 2 1 1 1 9 .

2. Indiferentes si γ 0 . 2 3 8 2 2 8 7 3 , γ 0 . 2 5 .

3. Repulsores para valores de γ diferentes de los anteriores.

Mientras que los puntos fijos extraños p f 3 ( γ ) y p f 4 ( γ ) son:

1. Atractores si 2 1 < γ < 1 2 ; en particular superatractor si γ 1 6 . 2 4 6 5 .

2. Indiferentes si γ 2 1 , γ 1 2 .

3. Repulsores para valores de γ diferentes de los anteriores.

Análisis numérico:

Para comprobar la estabilidad de los puntos fijos extraños que dependen de γ se hará de forma numérica usando diagramas de estabilidad, que son representaciones tridimensionales en donde el punto es atractor, indiferente o repulsor; esto debido a que no se puede comprobar de forma analítica ya que al sustituir estos puntos fijos extraños en la funciones de estabilidad, aparecen raíces de polinomios de alto grado.

Al realizar O p ' ( p f 1 ) O p ' ( p f 2 ) se obtiene como resultado 0. Quiere decir que la estabilidad de p f 1 ( γ ) y p f 2 ( γ ) es la misma, por lo que se procede a analizar solamente una de ellas. Análogamente se procede a realizar con p f 3 γ y p f 4 ( γ ) pues O p ' ( p f 3 ) O p ' ( p f 4 ) = 0

A.5. Prueba del Lema 1

Si γ = 0 , entonces el operador es O p ( z ) = z 4 , cumpliendo el test de Cayley, y los puntos fijos son la solución de z 4 = z , o sea, las soluciones de z ( z 1 ) ( 1 + z + z 2 ) = 0 ; por lo tanto hay tres puntos fijos extraños: z = 1 , z = 1 3 y z = ( 1 ) 2 / 3 . Además, el método que se obtiene es el método de Ostrowski.

A.6. Prueba del Lema 2

1. Si γ = 0 , entonces el operador es O p ' ( z ) = 4 z 3 y los puntos críticos son la solución de 4 z 3 = 0 ,o sea, las solución es z = 0 con multiplicidad tres, y como el cero es una raíz, no hay puntos críticos libres.

2. Si γ = 1 7 8 , entonces el operador es O p ' ( z ) = 2 z 3 ( 2 0 + 2 3 z + 2 0 z 2 ) ( 8 + 1 6 z + 2 5 z 2 ) 2 y los puntos críticos son las raíces de O p ' ( z ) = 0 . Las soluciones son z = 0 con multiplicidad tres y z = 1 4 0 ( 2 3 ± 3 i 1 1 9 ) , ambas con multiplicidad dos.

3. Si γ = 8 , entonces el operador es O p ' z = 4 ( 1 + z ) 2 z 3 ( 7 + 2 2 z + 7 z 2 ) ( 1 + 2 z 7 z 2 ) 2 y los puntos críticos son las raíces de O p ' ( z ) = 0 . Las soluciones son z = 0 con multiplicidad tres, z = 1 con multiplicidad dos y z = 1 7 ( 1 1 ± 6 2 ) .

4. Si γ = 4 , entonces el operador es O p ' ( z ) = 4 z 3 ( 3 + 8 z + 3 z 2 ) ( 1 + 3 z ) 2 y los puntos críticos son las raíces de O p ' ( z ) = 0 . Las soluciones son z = 0 con multiplicidad tres y z = 1 3 ( 4 ± 7 ) .

HTML generado a partir de XML-JATS4R por