SIMULACIONES NUMÉRICAS DE SISTEMAS

DINÁMICOS CAÓTICOS OSCILATORIOS QUE

CONSERVAN SU CAOS

Por

Jessica Zaqueros Martı́nez

Tesis sometida como requisito parcial para

obtener el grado de

MAESTRÍA EN CIENCIAS EN EL ÁREA DE CIENCIAS

COMPUTACIONALES

En el

Instituto Nacional de Astrofı́sica, Óptica y Electrónica

Noviembre, 2018
Tonantzintla, Puebla, México

Supervisada por

Dr. Gustavo Rodrı́guez Gómez, INAOE

Dr. Felipe Orihuela Espina, INAOE
Dr. Esteban Tlelo Cuautle, INAOE

©INAOE 2018
Derechos Reservados
El autor otorga al INAOE el permiso de reproducir
y distribuir copias de esta tesis en su
totalidad o en partes mencionando la fuente.
A mi familia por su amor incondicional.
 Agradecimientos

Las palabras no son suficientes para expresar mi agradecimiento por todo lo apren-
dido durante estos dos años de estancia en el Instituto Nacional de Astrofı́sica, Óptica y
Electrónica, fue más de lo que imaginé.

Agradezco a cada uno de mis profesores por contribuir en mi formación, de manera

especial a mis asesores los doctores Gustavo, Felipe y Esteban por su invaluable apoyo,
esfuerzo y dedicación hacia mı́. Asimismo a mi comité revisor los doctores Alejandro,
Alfonso y la doctora Claudia por las sugerencias y consejos brindados que permitieron
mejorar este escrito.

Agradezco infinitamente a mi familia por ser fortaleza para mı́ en todo momento y
porque me inspiran a siempre dar lo mejor, y a mis amigos por la vivencias inolvidables
que permitieron conocernos mejor.

Finalmente, agradezco a mi paı́s por darme la oportunidad de realizar este posgrado

a través de la beca No. 609491 otorgada por el Consejo Nacional de Ciencia y Tecnologı́a
(CONACYT).

V
 Índice general

Agradecimientos v

Índice general vi

Índice de figuras xi

Índice de tablas xii

Acrónimos xv

Resumen xvii

Abstract xix

1. Introducción 1
1.1. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2. Problemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3. Justificación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4. Preguntas de investigación . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5. Hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.6. Objetivo General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.6.1. Objetivos particulares . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7. Racional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.8. Contribución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.9. Organización de la tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2. Marco Teórico 9
2.1. Ecuaciones diferenciales . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2. Caracterización de sistemas dinámicos . . . . . . . . . . . . . . . . . . . 10

VII
 2.3. Espacios métricos y lı́mites . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4. Definiendo al caos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.1. Condiciones necesarias para sistemas caóticos . . . . . . . . . . . 16
2.4.2. Exponentes de Lyapunov . . . . . . . . . . . . . . . . . . . . . . 17
2.4.3. Espacio de fase . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4.4. Periodograma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.5. Diagramas de bifurcación . . . . . . . . . . . . . . . . . . . . . . 20
2.4.6. Condiciones suficientes para sistemas caóticos . . . . . . . . . . 21
2.4.7. Espacio de proyecciones de retraso . . . . . . . . . . . . . . . . . 22
2.5. Soluciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.1. Ensombrecimiento . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.2. Métodos Multipasos Lineales . . . . . . . . . . . . . . . . . . . . 25
2.5.3. Métodos Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . 28
2.5.4. Métodos de extrapolación . . . . . . . . . . . . . . . . . . . . . . 29
2.5.5. Métodos especiales . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5.6. Errores relativo y absoluto . . . . . . . . . . . . . . . . . . . . . . 33
2.6. Sistemas caóticos elegidos . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.6.1. Sistema de Lorenz . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.6.2. Atractor de Rossler . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.6.3. Sistema de Chen . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.6.4. Sistema de Chua . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6.5. Ecuaciones de Rabinovich–Fabrikant . . . . . . . . . . . . . . . . 36
2.6.6. Ecuaciones de Lorenz–Stenflo . . . . . . . . . . . . . . . . . . . . 37
2.7. Resumen de capı́tulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3. Trabajo relacionado 39
3.1. Aplicaciones de los sistemas caóticos . . . . . . . . . . . . . . . . . . . . 39
3.2. Problemas en simulaciones caóticas . . . . . . . . . . . . . . . . . . . . . 41
3.2.1. Alta sensibilidad a las condiciones iniciales . . . . . . . . . . . . 41
3.2.2. Errores de discretización: truncamiento y redondeo . . . . . . . 42
3.2.3. Supresión de caos . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
 3.2.4. Caos falso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2.5. Duración de las simulaciones . . . . . . . . . . . . . . . . . . . . 45
3.3. Métodos numéricos empleados para resolver sistemas caóticos . . . . . . 45
3.3.1. Orientados a las aplicaciones . . . . . . . . . . . . . . . . . . . . 46
3.3.2. Analizando los métodos numéricos . . . . . . . . . . . . . . . . . 47
3.3.3. Métodos numéricos para problemas oscilatorios . . . . . . . . . . 49
3.4. Resumen de capı́tulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4. Planteamiento del problema y metodologı́a propuesta 53

4.1. Importancia del problema dentro de las ciencias computacionales . . . . 53
4.2. Metodologı́a propuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2.1. Selección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2.2. Simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2.3. Evaluación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2.4. Exponentes de Lyapunov . . . . . . . . . . . . . . . . . . . . . . 58
4.2.5. Entropı́a de Kolmogorov–Sinai . . . . . . . . . . . . . . . . . . . 60
4.3. Resumen de capı́tulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5. Experimentos y resultados 65
5.1. Primer experimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.2. Segundo experimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.3. Tercer experimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3.1. Evaluación y comparación . . . . . . . . . . . . . . . . . . . . . . 75
5.3.2. Espacio de fase . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.3.3. Periodograma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.3.4. Número de evaluaciones . . . . . . . . . . . . . . . . . . . . . . . 83
5.3.5. Exponentes de Lyapunov . . . . . . . . . . . . . . . . . . . . . . 84
5.3.6. Entropı́a de Kolmogorov–Sinaı́ . . . . . . . . . . . . . . . . . . . 86
5.3.7. Conclusiones del tercer experimento . . . . . . . . . . . . . . . . 90
5.4. Resumen de capı́tulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6. Conclusiones y trabajo futuro 93

 6.1. Hallazgos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.2. Trabajo futuro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

A. Condiciones de simulación 97

B. Algoritmo polinomio trigonométrico Gautschi 99

Bibliografı́a 103
 Índice de figuras

2.1. Dependencia sensible a las condiciones iniciales . . . . . . . . . . . . . . 17

2.2. Ejemplo de atractor extraño . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3. Ejemplos de periodogramas . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4. Diagrama de bifurcación para el mapeo logı́stico . . . . . . . . . . . . . . 21
2.5. Espacio de fase y de proyecciones de restrasos . . . . . . . . . . . . . . . 23

4.1. Diagrama de la metodologı́a propuesta . . . . . . . . . . . . . . . . . . . 54

4.2. Espacio de fase de un sistema dinámico discreto . . . . . . . . . . . . . . 59
4.3. Estimación del máximo exponente de Lyapunov . . . . . . . . . . . . . . 60
4.4. Ejemplos de entropı́as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5. Aproximación de la entropı́a de la función seno(x) + ruido. . . . . . . . 63

5.1. Supresión de caos en el sistema caótico de Lorenz . . . . . . . . . . . . . 67

5.2. Supresión de caos en las ecuaciones de Rabinovich–Fabrikant . . . . . . 68
5.3. Supresión de caos en las ecuaciones de Lorenz–Stenflo . . . . . . . . . . 69
5.4. Caos falso en el sistema no caótico de Lorenz . . . . . . . . . . . . . . . . 71
5.5. Caos falso en las ecuaciones de Rabinovich–Fabrikant . . . . . . . . . . . 72
5.6. Caos falso en las ecuaciones de Lorenz–Stenflo . . . . . . . . . . . . . . . 72
5.7. Espacio de fase del sistema de Lorenz . . . . . . . . . . . . . . . . . . . . 79
5.8. Espacio de fase del atractor de Rossler . . . . . . . . . . . . . . . . . . . . 80
5.9. Espacio de fase del sistema de Chen . . . . . . . . . . . . . . . . . . . . . 81
5.10. Espacio de fase del sistema de Chua . . . . . . . . . . . . . . . . . . . . . 82

XI
 Índice de Tablas

5.1. Frecuencias aproximadas de los sistemas caóticos seleccionados . . . . . 74

5.2. Pasos de integración utilizados por el B–E . . . . . . . . . . . . . . . . . 76
5.3. Resultados numéricos de los métodos de sofisticados con Er = 10−2 . . . 77
5.4. Resultados numéricos de los métodos de sofisticados con Er = 10−8 . . . 78
5.5. Pasos de integración utilizados por el B–E . . . . . . . . . . . . . . . . . 78
5.6. Periodogramas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.7. Número de evaluaciones realizadas por los métodos numéricos . . . . . . 87
5.8. Estimación del máximo exponente de Lyapunov . . . . . . . . . . . . . . 88
5.9. Estimación del espectro de los exponentes de Lyapunov . . . . . . . . . 89
5.10. Estimación de la entropı́a de K–S . . . . . . . . . . . . . . . . . . . . . . . 91

XIII
 Acrónimos

A continuación mostramos los acrónimos utilizados en esta tesis.

AB–AM: Adams Bashforth – Adams Moulton

ADM: Método de Descomposición de Adomian (siglas en inglés)

B–E: Euler hacia atrás (siglas en inglés)

CNS: Simulación Numérica Limpia (siglas en inglés)

CSA: Recocido Simulado de Caos (siglas en inglés)

DQ: Cuadratura Diferencial (siglas en inglés)

EDO: Ecuación Diferencial Ordinaria

F–E: Euler hacia adelante (siglas en inglés)

FIFO: First In, First Out

K–S: Kolmogorov–Sinaı́

MML: Método Multipaso Lineal

MSRM: Método de Relajación Espectral Multietapa (siglas en inglés)

ODE: Ordinary Diferential Equation

ODEX: Método de extrapolación

PI: Path Integration

PVI: Problema de Valores Iniciales

R–K: Runge–Kutta

XV
R–K F: Runge–Kutta Felbergh

TSP: Problema del viajero (siglas en inglés)

 Resumen

El uso de simulaciones numéricas para estudiar la dinámica de los sistemas caóticos

tiene algunas deficiencias asociadas que incluyen la supresión del caos para los sistemas
caóticos verdaderos, o la inducción del caos en sistemas que de otra manera no serı́an
caóticos. En consecuencia, los resultados de las simulaciones numéricas pueden apartar-
se sustancialmente de la verdadera solución. Aunque muchas simulaciones numéricas de
la familia Euler o Runge-Kutta emplean estrategias numéricas para contener el error, a
menudo no contemplan las caracterı́sticas especiales de los sistemas caóticos. Una de es-
tas particularidades ocurre cuando los sistemas caóticos son oscilatorios y, por lo tanto,
sus soluciones abarcan muchas frecuencias. En este trabajo, proponemos la aplicación de
métodos numéricos basados en polinomios trigonométricos sobre las elecciones tradicio-
nales porque son capaces de capturar comportamientos oscilatorios. Nuestra hipótesis
es que tal elección resultará en la preservación del verdadero comportamiento caótico
durante simulaciones más largas. Especı́ficamente, mostramos que las simulaciones de
sistemas caóticos basados en el método de Gautschi preservan el caos para simulacio-
nes de al menos 50,000 unidades de tiempo, superando comparativamente el desempeño
de otras opciones tradicionales. Las familias de simulación numérica consideradas aquı́
para comparación incluyen las estrategias de Euler, Runge-Kutta, Adams-Bashforth y
Adams-Moulton. Tres experimentos se llevaron a cabo. Los dos primeros proporciona-
ron evidencia empı́rica de la supresión del caos en los sistemas caóticos con el método
de Euler hacia atrás, y la generación del caos en los sistemas no caóticos con el método
de Euler hacia adelante, justificando la necesidad de considerar las particularidades de
los sistemas caóticos. El tercer experimento mostró que el método de Gautschi simula
satisfactoriamente los sistemas caóticos de prueba. Durante la evaluación, los resultados
numéricos se evaluaron cualitativamente a través de su espacio de fase y periodograma y
se evaluaron cuantitativamente en términos de los exponentes de Lyapunov, la entropı́a

XVII
de Kolmogorov - Sinaı́ y el número de evaluaciones realizadas por el método numérico
en cada paso al modelo que representa el sistema caótico como aproximación del coste
computacional. El método de Gautschi exhibió espacios de fase y espectros de frecuen-
cia con un comportamiento similar a los métodos alternativos, exponentes de Lyapunov
positivos estadı́sticamente similares y una entropı́a de K-S más alta, requiriendo una can-
tidad de evaluaciones hasta cuarenta y cinco veces menor que las alternativas probadas.
Nuestros resultados sugieren que las particularidades del sistema caótico no se pueden
ignorar y que una elección adecuada del sistema numérico es fundamental para garan-
tizar la búsqueda de buenas soluciones. Además, apoyan nuestra hipótesis para el uso
de métodos polinomios trigonométricos en sistemas caóticos oscilatorios. Los resultados
obtenidos tienen repercusiones en una amplia variedad de dominios como la climato-
logı́a o las telecomunicaciones, donde tener simulaciones confiables implica predicciones
climáticas más precisas o comunicaciones más seguras.

Palabras clave: Sistemas dinámicos, caos, supresión de caos, caos falso, algoritmos
numéricos, método de Gautschi.
 Abstract

The use of numerical simulations for studying the dynamic of chaotic systems has
some associated shortcomings including chaos suppression for true chaotic systems, or
the induction of chaos in otherwise non chaotic systems. In consequence, numerical si-
mulations results may substantially depart from the true solution. Although many nume-
rical simulations from Euler or Runge-Kutta family employ numerical strategies to con-
tain error, they often do not contemplate the special characteristics of chaotic systems.
One of these particularities occurs when chaotic systems are oscillatory and hence, their
solutions span over many frequencies. In this work, we propose the application of nume-
rical methods based on trigonometric polynomials over traditional choices because they
are able to capture oscillatory behaviours. We hypothesize that such choice shall result
in the preservation of true chaotic behaviour during longer simulations. Specifically, we
show that simulations of chaotic systems based on Gautschi method preserve the chaos
for simulations of at least 50,000 time units, comparatively surpassing the performance
of other traditional choices. Numerical simulation families considered here for compari-
son include Euler, Runge-Kutta, Adams-Bashforth and Adams-Moulton strategies. Three
experiments were carried out. The first two experiments provided empirical evidence of
chaos suppression in chaotic systems with the backward Euler method, and chaos gene-
ration in non-chaotic systems with forward Euler method, justifying the need of consider
chaotic systems particularities. The third experiment showed that the Gautschi method
satisfactorily simulates the test chaotic systems. During evaluation, numerical results
were assessed qualitatively through their phase space and periodogram and quantitati-
vely evaluated in terms of the Lyapunov exponents, the Kolmogorov-Sinai entropy and
the number of evaluations performed by the numerical method at each step to the model
representing the chaotic system as a proxy of computational cost. The Gautschi method
exhibited phase spaces and frequency spectrums with a behaviour similar to alternative

XIX
methods, and positive Lyapunov exponents statistically similar and higher KS entropy
whilst requiring a number of evaluations up to forty-five times smaller than tested al-
ternatives. Our results suggest that the particularities of the chaotic system cannot be
ignored and that an appropriate choice of numerical system is critical to ensure finding
good solutions. Further, they support our hypothesis for using trigonometric polyno-
mials methods on oscillatory chaotic systems. The obtained results have repercussions
in a wide variety of domains such as climatology or telecommunications, where having
reliable simulations imply more accurate climate predictions or safer communications.

Keywords: Chaotic systems, chaos, chaos suppression, false chaos, numerical al-
gorithms, Gautschi method.
 1
Introducción

En el lenguaje cotidiano el caos implica la existencia de comportamientos imprede-

cibles. La palabra usualmente lleva una connotación negativa que involucra desorganiza-
ción o confusión indeseables. En el ámbito cientı́fico este comportamiento no es necesa-
riamente indeseable y se entiende como una caracterı́stica inherente de ciertos sistemas
naturales. El caos sólo es un desorden en apariencia, tiene muy poco que ver con el azar.
Aunque los sistemas caóticos parecen evolucionar de forma aleatoria y errática, éstos
tienen en realidad un orden interno subyacente. Los sistemas caóticos son deterministas
pero difı́ciles de predecir.

No existe una definición universal del caos en matemáticas. No obstante, los siste-
mas caóticos se caracterizan por tres principios importantes (Bradley, 2009; Smith, 2007):

La sensibilidad a los parámetros y a las condiciones iniciales.

Son deterministas.
Son no lineales.

La alta sensibilidad a las condiciones iniciales se conoce a menudo como efecto mariposa.
Las diferencias minúsculas en las condiciones iniciales pueden aparecer debido al ruido o
la precisión de la máquina utilizada para aproximar las soluciones, y generan resultados
completamente divergentes en los sistemas dinámicos caóticos. En general, en sistemas
caóticos la predicción a largo plazo resulta imposible (Shafique et al., 2018).
 1. Introducción

Históricamente, el precursor de la teorı́a del caos es el matemático francés Henri

Poincaré, quién en 1890 durante el estudio del problema especial de los tres cuerpos en
mecánica celeste concluı́a que pequeñas diferencias en las condiciones iniciales podrı́an
resultar en soluciones diferentes después de un largo tiempo (Sprott & Sprott, 2003).
Sin embargo, no fue hasta 1960 que, con la llegada de las computadoras, los cientı́ficos
pudieron solucionar, visualizar y estudiar a los problemas caóticos. Para entonces, ya se
empezaba a reconocer el caos en otros campos. Por ejemplo, Lorenz (1963) descubrió
fenómenos caóticos en el modelado del clima. Desde entonces la teorı́a del caos ha sido
explotada en muchos otros campos de la ciencia .

La mayorı́a de los sistemas dinámicos caóticos no cuentan con soluciones analı́ti-

cas. Para resolverlos se recurre a soluciones numéricas. Las simulaciones numéricas pue-
den proporcionar resultados erróneos porque sólo se obtienen aproximaciones y no so-
luciones precisas. En consecuencia, es importante determinar cuándo y en qué circuns-
tancias se tienen soluciones confiables (Varsakelis & Anagnostidis, 2016; Corless, 1994).

Las simulaciones numéricas son una herramienta útil para entender los sistemas
dinámicos no lineales y explorar el terreno de los sistemas caóticos, pero tiene sus li-
mitaciones. Las computadoras tienen precisión finita y en consecuencia es inevitable
generar errores al evaluar expresiones de punto flotante. Además, las computadoras se
encuentran en el contexto de tiempo discreto e inevitablemente existen errores cuando
son utilizadas para simular sistemas de tiempo continuo (Parker & Chua, 1989). Antes de
continuar, es importante mencionar que en este trabajo cuando nos referimos al tiempo,
estamos hablando de la variable independiente del sistema en cuestión sin ser necesaria-
mente una magnitud fı́sica.

Al aproximar las soluciones de los sistemas dinámicos caóticos, en general, éstas

son resueltas numéricamente empleando estrategias como el Euler o la familia de los
Runge-Kutta (Owolabi & Atangana, 2017; Noorani et al., 2007; Varsakelis & Anagnostidis,
2016) porque se encuentran disponibles en diversos paquetes de software como Matlab
(MathWorks USA) o Mathematica (Wolfram) y en consecuencia son fáciles de utilizar. A

2
1. Introducción Simulaciones numéricas

menudo, no se escogen estrategias numéricas que capturen las caracterı́sticas especiales

de este tipo de sistemas, las cuales podrı́an dar mejores aproximaciones que los métodos
comúnmente utilizados (Bardhi et al., 2010).

Actualmente, tampoco existen mecanismos que permitan identificar las buenas si-
mulaciones numéricas de sistemas dinámicos caóticos ni tampoco discernir cuándo el
sistema dinámico simulado realmente posee caos (Corless, 1994; Varsakelis & Anagnos-
tidis, 2016).

1.1. Simulaciones numéricas

Las simulaciones numéricas se definen de forma intuitiva como: Amplias coleccio-

nes de métodos y aplicaciones que imitan el comportamiento real de los sistemas por
medio de una computadora con el software apropiado.

Las simulaciones numéricas utilizan algoritmos que aproximan numéricamente la

solución de procesos fı́sicos representados por modelos matemáticos. Son utilizadas co-
mo una herramienta del análisis numérico, y en conjunto con gráficas y cálculos ma-
temáticos simbólicos establecen, resuelven e interpretan modelos matemáticos compli-
cados (Atkinson, 2017).

Los resultados de las simulaciones numéricas pueden estar bastante alejados de la

solución verdadera del proceso fı́sico en cuestión. Esto puede deberse a distintos factores:

Algoritmos numéricos utilizados.

El modelo no representa apropiadamente al proceso fı́sico.
La cantidad de ecuaciones que tiene el modelo matemático.
Errores de redondeo dependientes de la computadora utilizada para realizar la si-
mulación.
El problema es mal condicionado: al efectuar pequeños cambios en los parámetros

3
 Problemática 1. Introducción

de entrada se originan grandes cambios en los parámetros de salida.

Un ejemplo de una simulación numérica es el siguiente. Supongamos que queremos

resolver el siguiente sistema y ′ (t) = cos(t) sujeto a la condición inicial y(0) = 0 en el
intervalo [0, π]. Dado que se trata de la función cos(t), la solución analı́tica del sistema es
y(t) = sen(t). Sin embargo, podemos resolver tal sistema con un método del tipo Euler
y la solución aproximada serı́a de la siguiente forma

yn+1 = yn + hcos(tn ), con y0 = y(0). (1.1)

donde h > 0 y tn = 0 + nh, n = 1, 2, 3, . . . .

El método Euler discretiza la solución del sistema, por lo que se tiene la solución
aproximada solo en un conjunto de puntos y no en el intervalo continuo como sucede
con la solución analı́tica.

1.2. Problemática

Cuando tenemos un sistema dinámico caótico oscilatorio, sus condiciones iniciales

y un intervalo de tiempo grande en el que pretendemos resolverlo, es necesario contar
con el o los métodos numéricos que nos permitan simular digitalmente el comportamien-
to del sistema durante el intervalo dado, asegurando que conserven el caos inherente al
sistema caótico. En este trabajo entendemos por grande, simulaciones de al menos 50, 000
unidades de tiempo por el contexto de las aplicaciones de los sistemas caóticos (sección
1.3).

Además, después de realizar las simulaciones numéricas debemos disponer de métri-

cas que permitan cuantificar el caos e identificar la calidad de los resultados numéricos.
En consecuencia, es conveniente tener una metodologı́a que permita eligir al método
numérico que cumpla con las condiciones exigidas.

4
1. Introducción Justificación

En este trabajo para examinar las simulaciones numéricas se eligieron tres tipos
de métodos numéricos: estándares (utilizan paso de integración fijo, p.e. Euler o R–K),
sofisticados (paso de integración y orden variable con o sin control de error) y especiales
(con paso de integración fijo pero basados en polinomios trigonométricos).

Concisamente es importante garantizar que al simular sistemas dinámicos caóticos

los resultados numéricos obtenidos son confiables. Mas aún, es imprescindible contar con
mecanismos que nos permitan elegir adecuadamente estrategias numéricas para sistemas
dinámicos caóticos y conocer cuándo realmente los resultados numéricos del sistema
simulado son o no caóticos porque puede suceder que el sistema dinámico es caótico,
pero los resultados indican lo contrario (Varsakelis & Anagnostidis, 2016).

1.3. Justificación

El caos afecta las predicciones meteorológicas, que nos afectan directamente a

través del clima, e indirectamente a través de las consecuencias económicas tanto del
clima como de los propios pronósticos. El caos también juega un papel importante en las
cuestiones del cambio climático y nuestra capacidad de prever la fuerza y los impactos
del calentamiento global (Smith, 2007).

Hoy la teorı́a del caos se utiliza en la criptografı́a (Soriano-Sánchez et al., 2016;

Garcı́a-Martı́nez et al., 2015; Behnia et al., 2008; Wang et al., 2016; Belazi et al., 2016;
Zhou et al., 2015), en la creación de modelos de población en la biologı́a (Pham et al.,
2015); estudiando la turbulencia en la mecánica de fluidos (Shafique et al., 2018). En eco-
nomı́a se utiliza para predecir el comportamiento del mercado de valores (Smith, 2007).
En astronomı́a el caos se sigue utilizado para describir el movimiento de muchos cuer-
pos planetarios, para predecir mejor los caminos de asteroides y si pueden entrar o no en
contacto con la Tierra (Shafique et al., 2018). En años más recientes se ha aplicado en la
predicción y/o control de la dinámica del cerebro humano (Shafique et al., 2018). Aunque
esta lista de aplicaciones no es exhaustiva, proporciona una idea de cómo el caos está

5
 Preguntas de investigación 1. Introducción

presente en diversos campos de las ciencias directa o indirectamente relacionados con el

área de Ciencias Computacionales.

1.4. Preguntas de investigación

¿Cuáles pruebas se pueden aplicar a los resultados numéricos de las simulaciones

de sistemas dinámicos de tal modo que si son positivas impliquen que el sistema simulado
es caótico?

¿Cuáles métodos numéricos nos permiten obtener simulaciones de sistemas caóti-

cos durante intervalos grandes de tiempo con bajo costo en tiempos de ejecución con-
servando el caos?

¿Cuáles métricas permiten cuantificar el caos en los resultados numéricos de las

simulaciones de sistemas que se saben a priori que son caóticos?

1.5. Hipótesis

En sistemas dinámicos caóticos oscilatorios, el uso de métodos numéricos basados

en polinomios trigonométricos permite realizar simulaciones que conservan el caos en
intervalos de tiempo grandes.

1.6. Objetivo General

Conservar el caos de sistemas dinámicos caóticos oscilatorios cuando se simulan

numéricamente y mantenerlo durante periodos de tiempo grandes mediante el uso de

6
1. Introducción Racional

estrategias numéricas basadas en polinomios trigonométricos.

1.6.1. Objetivos particulares

1. Establecer casos en que los métodos numéricos suprimen el caos y evaluarlos me-
diante el espacio de fase.

2. Establecer casos en que los métodos numéricos generan caos falso y evaluarlos con
el espacio de fase.

3. Conservar el caos de un sistema dinámico caótico cuando se simula numéricamen-

te, mantenerlo durante grandes periodos de tiempo y evaluarlo con métricas como
el espacio de fase, el periodograma, los exponentes de Lyapunov y la entropı́a de
Kolmogorov-Sinaı́ (K–S).

4. Evaluar el beneficio de métodos especiales con respecto a los métodos estándares

y sofisticados comparando los resultados de las métricas calculadas en el objetivo
anterior.

1.7. Racional

Parker & Chua (1989) y Gautschi (1961) mencionan que al momento de solucionar
problemas es importante tomar en cuenta las caracterı́sticas especiales que éstos poseen.
De modo que si no se contemplan y se aplican métodos numéricos genéricos puede con-
ducirnos a resultados erróneos. Por ejemplo, el método de Runge-Kutta puede generar
soluciones estables que son el resultado de las discretizaciones y no una caracterı́stica
de la ecuación diferencial subyacente, por lo que podrı́an conducir a resultados compu-
tacionales erróneos (Vadillo, 1997).

Ahora bien, existen métodos numéricos especiales que fueron pensados para apro-

7
 Contribución 1. Introducción

vechar caracterı́sticas especı́ficas de los sistemas a resolver. Los métodos basados en po-
linomios trigonométricos son de esta clase, porque se benefician de las propiedades os-
cilatorias de los problemas a resolver (Gautschi, 1961; Lambert, 1973).

Los sistemas caóticos oscilatorios tienen como caracterı́stica natural las oscilacio-
nes en sus soluciones. En consecuencia, al utilizar métodos especiales basados en po-
linomios trigonométricos para solucionar sistemas caóticos oscilatorios se espera tener
mejores resultados que con los métodos de propósito general, porque los especiales apro-
vechan la caracterı́stica oscilatoria inherente a los sistemas caóticos oscilatiorios mien-
tras que los de propósito general no.

1.8. Contribución

La contribución de este trabajo de investigación es una metodologı́a que permite

seleccionar un método numérico especial (basado en polinomios trigonométricos) para
aproximar la solución numérica de un sistema caótico preservando su caos en periodos
de tiempo arbitrarios.

1.9. Organización de la tesis

La organización del presente escrito es la siguiente. En el capı́tulo 2 exponemos de

manera sucinta la teorı́a necesaria para entender la tesis. En el capı́tulo 3 presentamos
una recopilación de trabajos relevantes relacionados al tema que ocupa a esta investiga-
ción. En el capı́tulo 4 describimos detalladamente la problemática y la solución propuesta.
En el capı́tulo 5 presentamos las pruebas realizadas y los resultados de las mismas. Final-
mente, en el capı́tulo 6 proporcionamos las conclusiones generales de la tesis, ası́ como
el posible trabajo a futuro.

8
 2
Marco Teórico

En este capı́tulo proporcionamos las nociones básicas para la comprensión del tra-
bajo propuesto. Primero, damos la caracterización de un sistema caótico representado por
ecuaciones diferenciales ordinarias. Posteriormente describimos los métodos numéricos
que utilizamos en este trabajo, y finalmente proporcionamos los sistemas caóticos utili-
zados para realizar los experimentos en la presente investigación.

2.1. Ecuaciones diferenciales

El caos es una área de estudio en el campo de la dinámica no lineal, la cual forma

parte de un campo de estudio más amplio: los sistemas dinámicos. Un sistema dinámico es
aquél que evoluciona en el tiempo. Los sistemas dinámicos pueden ser estocásticos, los
cuales evolucionan de acuerdo con algún proceso aleatorio como lanzar una moneda, o
deterministas, en los que el futuro es determinado únicamente por el pasado de acuerdo
con alguna regla o fórmula matemática (Sprott & Sprott, 2003).

Además, los sistemas dinámicos pueden ser disipativos o conservativos. En los pri-
meros actúan fuerzas que disipan o pierden energı́a mecánica (como la fricción), mientras
que en los conservativos sólo actúan fuerzas que conservan la energı́a mecánica (Boyce
et al., 1969; Drazin, 1992).
 Caracterización de sistemas dinámicos 2. Marco Teórico

En esta tesis estudiaremos sistemas dinámicos caóticos que son no lineales, de-
terministas y conservativos. Las ecuaciones que rigen a estos sistemas son ecuaciones
diferenciales ordinarias (EDOs) no lineales. Una EDO es una igualdad que involucra una
función y sus derivadas (Boyce et al., 1969). Un sistema es lineal si la suma de solucio-
nes también es una solución: permite la superposición de soluciones (Smith, 2007). La
propiedad anterior no es cierta para sistemas no lineales.

La notación matemática que utilizaremos en lo que resta del capı́tulo es la siguiente:

1. Negritas para vectores.

dq x
2. Superı́ndices para q-ésimas derivadas : x(q) (tk ) = .
dtq t=tk

3. El superı́ndice T indica la transpuesta del vector.

2.2. Caracterización de sistemas dinámicos

Un sistema autónomo de EDOs no depende explı́citamente de la variable inde-

pendiente (usualmente t). Por ejemplo, x′ = x2 es autonómo, mientras que x′ = x2 + t
no lo es (Boyce et al., 1969).

El orden de un sistema de EDOs es la derivada de orden mayor que aparece en el

sistema (Robinson, 2004). Por ejemplo, la segunda ley del movimiento de Newton
d2 x
m 2 = F, (2.1)
dt
es de segundo orden.
Definición 1. Un sistema dinámico continuo autónomo de orden n-ésimo se define por
la ecuación de estado
dn x(t)
= g(x), x(t0 ) = x0 (2.2)
dtn
donde x(t) es el estado en el tiempo t.

10
2. Marco Teórico Caracterización de sistemas dinámicos

La ecuación (2.2) es equivalente a un sistema de EDOs de primer orden mediante

la siguiente transformación. Sea x ∈ Rn con x = (x1 , . . . , xn ) tal que

x1 = x (2.3)
x2 = x(1) (2.4)
x3 = x(2) (2.5)
..
. (2.6)
xn = x(n−1) (2.7)

Entonces, la ecuación (2.2) se puede reescribir de la siguiente manera

dx
= f(x) x(t0 ) = x0 (2.8)
dt
donde dx/dt = [x(1) , x(2) , . . . , x(n−1) ]T y f : Rn → Rn tal que f(x) = [x2 , x3 , . . . , xn , g(x)]T .
La función f es llamada campo vectorial en Rn (Parker & Chua, 1989). Puesto que el
campo vectorial no depende del tiempo, el tiempo inicial siempre puede ser tomado como
t0 = 0. El sistema dinámico (2.8) es no lineal si el campo vectorial f es no lineal.

La solución de la ecuación (2.8) es x(t) tal que x(t0 ) = x0 . Podemos considerar

la solución como el camino que sigue una partı́cula que comienza en un punto dado x0
y se mueve con velocidad f en Rn para t > 0. Para describir esta ruta es útil definir la
función uniparamétrica φt : Rn → Rn . Para mostrar explı́citamente la dependencia de
la condición inicial, a menudo la solución de (2.8) se escribe como φt (x0 ).
Definición 2. La familia de funciones de un solo parámetro φt : Rn → Rn que satisface
las condiciones dadas por las ecuaciones (2.9) y (2.10) es conocido como flujo.

φt1 +t2 = φt1 ◦ φt2 (2.9)

φ0 (x) = x (2.10)

Definición 3. El conjunto de puntos {φt (x0 ) : −∞ < t < ∞} es nombrado trayectoria

a través de x0 .

11
 Caracterización de sistemas dinámicos 2. Marco Teórico

Mediante un ejemplo se explican las dos definiciones anteriores. Considere el sis-

tema dinámico de primer orden dado por la expresión (2.11). La solución al sistema está
dada por la ecuación (2.12).
dx
= −5x x(t0 ), (2.11)
dt

x(t) = x0 e−5t . (2.12)

Sea φt (x0 ) = x(t), entonces

φt1 +t2 (x0 ) = (φt1 ◦ φt2 )(x0 ) (2.13)

= φt2 (φt1 (x0 )) (2.14)
= φt2 (x0 e−5t1 ) (2.15)
= x0 e−5t1 e−5t2 (2.16)
= x0 e−5(t1 +t2 ) . (2.17)

Además

φt (x0 ) = x0 exp−5(t) (2.18)

Los desarrollos anteriores implican que φt (x) = xe−5t es el flujo porque cumple
con las condiciones dadas por las ecuaciones (2.9) y (2.10). Mientras que φt (x0 ) = x0 e−5t
es una trayectoria. La diferencia entre flujo y trayectoria, es que el primero es una familia
de funciones donde el valor de x es variable y, la trayectoria es una sola función en la
que x es un número fijo ( x = x0 ).
Definición 4. Un punto de equilibrio xeq de un sistema autónomo es una solución que
satisface
xeq = φt (xeq ) para todo t. (2.19)

Un punto de equilibrio xeq de un sistema autónomo es una solución constante de

(2.8) (Parker & Chua, 1989). En un punto de equilibrio, el campo vectorial f se desvanece

12
2. Marco Teórico Caracterización de sistemas dinámicos

y, salvo casos patológicos (Parker & Chua, 1989), f(x) = 0 implica que x es un punto de
equilibrio.
Definición 5. φt (x∗ ) es una solución periódica de un sistema autónomo si, para toda t,

φt (x∗ ) = φt+T (x∗ ) (2.20)

para algún periodo mı́nimo T > 0 (Parker & Chua, 1989).

Una solución periódica es aislada si posee una vecindad (ver definición (2.3)) que
no contenga otra solución periódica. En este caso, tal solución se denomina ciclo lı́mite.
Definición 6. El estado estacionario se refiere al comportamiento asintótico de un siste-
ma cuando t → ∞ (Parker & Chua, 1989).

Para que tenga sentido, se requiere que el estado estacionario sea acotado. A con-
tinuación, se presenta la definición de acotado. Algunas veces, un punto de equilibrio (o
una solución periódica) es llamado estado estacionario (Boyce et al., 1969).
Definición 7. Una función ρ : R → R es acotada si existen M y m tales que m ≤ ρ(x) ≤
M para todo ρ(x) ∈ R.
Definición 8. Una función cuasi-periódica es aquella que puede expresarse como una
suma numerable de funciones periódicas (Parker & Chua, 1989)
X
x(t) = hi (t) (2.21)
i

donde hi tiene un periodo mı́nimo Ti y frecuencia ωi = 1/Ti . Además, debe existir un

conjunto finito de frecuencias de base{ω̂1 , . . . , ω̂p } que cumple las siguientes propiedades:

1. Es linealmente independiente, esto es, no existe un conjunto no nulo de enteros

{k1 , . . . , kp } tal que k1 ω̂1 + · · · + kp ω̂p = 0

2. Forma una base integral finita para el ωi , es decir, para cada i, ωi = |k1 ω̂1 + · · · +
kp ω̂p | para algunos enteros {k1 , . . . , kp }.

Una solución cuasi-periódica con p frecuencias base es llamada p-periódica.

13
 Espacios métricos y lı́mites 2. Marco Teórico

2.3. Espacios métricos y lı́mites

A continuación se presentan los conceptos necesarios para entender una de las

definiciones del caos que se describe en la sección 2.4.
Definición 9. Un conjunto M se dice que es un espacio métrico si para cualesquiera dos
elementos x, y de M hay asociado un número real d(x, y) tal que

1. d(x, y) > 0 si x 6= y; d(x, x) = 0,

2. d(x, y) = d(y, x),

3. d(x, y) ≤ d(x, r) + d(r, y), para cualquier r ∈ M .

Cualquier función que cumpla con las tres propiedades anteriores es llamada función dis-
tancia o métrica (Rudin, 1976).

Para las siguientes definiciones, sea M un espacio métrico, x, y en M y E subcon-

junto de M .
Definición 10. Una vecindad del elemento x es un conjunto Vr (x) formado por todos los
puntos y tales que d(x, y) < r. Al número r se le llama radio de Vr (x).
Definición 11. x es un punto interior de E si existe una vecindad V de x tal que V ⊂ E.
Definición 12. Se dice que E es abierto si todos sus puntos son interiores.
Definición 13. Sean A y Ω dos conjuntos. Suponga que para cada elemento α ∈ A hay
asociado un conjunto de Ω, representado por Eα . El conjunto cuyos elementos son los
conjuntos Eα es representado por {Eα } y es llamado colección de conjuntos o familia de
conjuntos (Rudin, 1976).
Definición 14. Se llama cubierta abierta del conjunto E en un espacio métrico M a la
colección {Gα } de subconjuntos abiertos de M , tales que E ⊂ ∪α Gα
Definición 15. Un subconjunto K de un espacio métrico M es compacto si toda cubierta
abierta de K contiene una subcubierta finita.

La definición (15) quiere decir que si {Gα } es una cubierta abierta de K, entonces

14
2. Marco Teórico Definiendo al caos

existe un número finito de ı́ndices α1 , α2 , . . . , αn tales que

K ⊂ Gα1 ∪ . . . ∪ G αn . (2.22)

También, es necesario definir los siguientes lı́mites. Sea {an } una sucesión de núme-
ros reales, el lı́mite inferior está dado por la ecuación (2.23) y el lı́mite superior por
la ecuación (2.24) (Rudin, 1976).

lı́m inf an = sup ı́nf ak = sup{ı́nf{ak : k ≥ n} : n ≥ 0}, (2.23)

n→∞ n≥0 k≥n

lı́m sup an = ı́nf sup ak = ı́nf{sup{ak : k ≥ n} : n ≥ 0} (2.24)

n→∞ ≥0 k≥n

donde sup denota la menor de las cotas inferiores e ı́nf la mayor de las cotas inferiores.

Otra definición importante es la de conjunto desordenado (scrambled, en inglés)

dada por Li & Yorke (1975), la cual se encuentra en la definición 16.
Definición 16. Sea f : M → M , con M un espacio métrico compacto. Sea J un subcon-
junto infinito de M . Se dice que J es un conjunto desordenado si para cada par x, y ∈ J
con x 6= y se cumplen las siguientes condiciones

lı́m sup |f n (x) − f n (y)| > 0 (2.25)

n→∞

lı́m inf |f n (x) − f n (y)| = 0 (2.26)

n→∞

2.4. Definiendo al caos

Como se mencionó antes, no existe una definición ampliamente aceptada de caos.

Desde un punto de vista práctico, Parker & Chua (1989) lo definen como se presenta a
continuación.

15
 Definiendo al caos 2. Marco Teórico

Definición 17. El caos se define como un comportamiento de estado estacionario acotado

que no es un punto de equilibrio, es no periódico y tampoco cuasi periódico.

Otra definición de caos que es comúnmente utilizada es la propuesta por Li & Yorke
(1975), la cual está dada por la definición 18.
Definición 18. El caos se define como la presencia de un conjunto desordenado.

La definición dada por Li & Yorke (1975) fue motivada por las estructuras complejas
que presentan las soluciones de los sistemas caóticos. La condición dada por la inecuación
(2.25) indica que las soluciones de x a y son cercanas un número infinito de veces, y por
la condición dada en la ecuación (2.26) la distancia entre las soluciones de x a y también
excede un valor fijo positivo un número infinito de veces.

La mayorı́a de los sistemas caóticos muestran una dependencia sensible a las con-
diciones iniciales (Parker & Chua, 1989), como se ilustra en la Figura 2.1.
Definición 19. Se define a la dependencia sensible a las condiciones iniciales como las
trayectorias cercanas que divergen y eventualmente no se correlacionan (Parker & Chua,
1989).

Una de las ventajas de la definición dada por Li & Yorke (1975) es que excluye
algunos casos que presentan dependencia sensible a las condiciones iniciales pero que
usualmente son considerados como no caóticos. Sin embargo, también incluye algunos
casos que generalmente son considerados como no caóticos (Hunt & Ott, 2015). Como
el ámbito general de esta investigación no es definir el caos, al lector interesado en una
discusión más amplia sobre este tema se le recomienda consultar Hunt & Ott (2015).

2.4.1. Condiciones necesarias para sistemas caóticos

Aún cuando no se tiene una definición universal del caos, la comunidad matemática
coincide en la necesidad de que se cumplan ciertas condiciones para la consideración del
caos, hasta cierto punto, “definiendo” el caos por exclusión.

16
2. Marco Teórico Definiendo al caos

4
x(0) = 0.7000
x(0) = 0.7001
2
x(t)

-2

-4
0 5 10 15 20 25 30
t[s]

Figura 2.1: La dependencia sensible a las condiciones iniciales se ilustra con dos solucio-
nes de la variable x(t) del sistema caótico de Chua. Las condiciones iniciales solo difieren
al 0.01 %, puesto que una es 0.7000 y la otra 0.7001, sin embargo, después de 15s las
soluciones empiezan a divergir (Parker & Chua, 1989).

A continuación, reportamos las condiciones necesarias que presentan los sistemas

caóticos, no se trata de un informe exhaustivo pero si comprende las popularmente refe-
renciadas. Comenzamos con los exponentes de Lyapunov que son los más utilizados al
momento de caracterizar sistemas caóticos.

2.4.2. Exponentes de Lyapunov

Definición 20. Sea f una función tal que f : R → R y

|f (x0 + ǫ) − f (x0 )| ∼ expλ . (2.27)

Sea f n (x) = f (f (f (. . . f (x) . . .))), entonces

ǫ|Dx [f n (x0 )]| ∼ ǫ expnλ n → ∞. (2.28)

17
 Definiendo al caos 2. Marco Teórico

El exponente de Lyapunov λ es (Drazin, 1992)

( N −1 )
1 X
λ = lı́m ln |f ′ (xn )| . (2.29)
N →∞ N n=1

El exponente de Lyapunov es una medida cuantitativa de la separación de las solu-

ciones del sistema caótico a pequeñas perturbaciones en las condiciones iniciales. Si un
sistema es caótico entonces al menos un exponente de Lyapunov es positivo (Parker &
Chua, 1989). Cabe notar que si al menos un exponente es positivo no implica caos. Por
lo tanto, los exponentes de Lyapunov son una herramienta necesaria pero no suficiente
para caracterizar el caos. El siguiente ejemplo, tomado de Drazin (1992), es muestra de
ello. Considere la función dada por la ecuación (2.30).

g(a, x) = an x. (2.30)

Si a > 0 entonces λ = ln(a) > 0 pero la función no es caótica.

Definición 21. El espectro es el conjunto de todos los exponentes de Lyapunov. En un
sistema existen tantos exponentes de Lyapunov como variables dependientes (Sprott &
Sprott, 2003).

2.4.3. Espacio de fase

Otra herramienta necesaria pero no suficiente para caracterizar el caos es el espacio

de fase.
Definición 22. Un estado es un punto en la trayectoria (ver definición (3)) que define
completamente la condición al tiempo t del sistema (Smith, 2007).
Definición 23. El espacio de fase es el conjunto de todos los posibles estados de un
sistema dinámico (Robinson, 2004).

El atractor de un sistema caótico es un objeto que es complicado geométricamente

en el sentido que posee dimensión fraccional. El objeto geométrico en el espacio de fase

18
2. Marco Teórico Definiendo al caos

2
z(t)

-2

-0.2 -2
0 0
0.2
2
y(t) x(t)

Figura 2.2: Atractor extraño del sistema caótico de Chua con dos enrollamientos.

al que se atraen las trayectorias caóticas se denomina atractor extraño (Parker & Chua,
1989). Sin embargo, hasta el momento no existe una definición generalmente aceptada
de atractor extraño (Drazin, 1992).

Una forma de distinguir un atractor extraño de uno no caótico es mediante la exis-

tencia de un exponente de Lyapunov positivo. Un atractor extraño posee al menos tres
exponentes de Lyapunov. Entonces, el caos no se puede presentar en sistemas continuos
de primer o segundo orden (Parker & Chua, 1989). En la Figura 2.2 se muestra el atractor
extraño del sistema caótico de Chua con dos enrollamientos.

Una condición necesaria para aceptar que un comportamiento es caótico es la si-

P
guiente: λ1 > 0 y λi < 0 (Parker & Chua, 1989).

2.4.4. Periodograma

Para observar las soluciones de un sistema caótico en el dominio de la frecuencia

una representación común es el periodograma. Existen diferentes versiones del perio-

19
 Definiendo al caos 2. Marco Teórico

dograma, y aquı́ se utilizó la magnitud al cuadrado de la transformada discreta de Fourirer

de la solución en cuestión. El periodograma de los sistemas caóticos presenta una amplia
gama de frecuencias y es parecido al ruido (Parker & Chua, 1989). Además, es común que
el periodograma tenga picos en la frecuencias dominantes de la solución.

50 50

0 0
FFT(x)2

FFT(x)2
-50 -50

-100 -100
0 1 2 3 4 0 1 2 3 4
Frequencia (Hz) Frequencia (Hz)

(a) No caótico (b) Caótico

Figura 2.3: Peridogramas de la ecuación de Duffing. (a) Solución no caótica y (b) solución
caótica. En (b) se presenta una amplia gama de frecuencias y es parecido al ruido.

2.4.5. Diagramas de bifurcación

Otra condición necesaria pero no suficiente para el caos es la presencia de una

secuencia infinita de duplicaciones de periodo en los diagramas de bifurcación. Para ex-
plicar estos diagramas, primero consideremos la ecuación (2.31) llamada mapeo logı́stico.
Éste surge a partir de un modelo de población: Sea una especie de insectos en una isla
con fertilidad r, es decir, cada insecto produce un promedio r de descendencia en un ciclo
de vida. Hay una cantidad fija de alimentos en la isla, de modo que, si la población de
insectos aumenta, la fertilidad media disminuye.

xn+1 = rxn (1 − xn ) (2.31)

El mapeo logı́stico muestra una variedad de comportamientos y tiene transiciones entre

éstos a medida que cambiamos el parámetro r. Tales transiciones en sistemas dinámicos
se llaman bifurcaciones y son mostradas en la gráfica de la Figura 2.4.

20
2. Marco Teórico Definiendo al caos

1.0

(1 - ) 0.8
0.6
)=

0.4
0.2
(

0.0
0 1 2 3 4

Figura 2.4: Diagrama de bifurcación para el mapeo logı́stico. Las bifurcaciones son tran-
siciones entre diversos comportamientos que ocurren a medida que cambia el parámetro
r.

2.4.6. Condiciones suficientes para sistemas caóticos

Hasta el momento solo hemos reportado condiciones necesarias para la presencia

de caos en los resultados numéricos de sistemas caóticos. La siguiente condición implica
caos en el sentido de Li & Yorke (sección 2.4). Es la única condición de suficiencia de la
que tenemos constancia.
Definición 24. Sea un sistema dinámico con un espacio de fase M, el cual es d-dimensional.
Considere que M es partido en cajas con volumen δ. Sea L ∈ N fijo, τ > 0. La expresión
p(i1 , . . . , iL ) denota la probabilidad que y(τ ), y(2τ ), . . . , y(Lτ ) se encuentre dentro de la
caja i1 , i2 , . . . , iL , respectivamente. Entonces la K entropı́a de K–S está dada por

1
K = lı́m lı́m lı́m (KL+1 − KL ) (2.32)
τ →0 δ→∞ L→∞ Lτ

donde
X
KL ≡ − p(i1 , . . . , iL ) ln(p(i1 , . . . , iL )). (2.33)
i1 ,i2 ,...,iL

21
 Definiendo al caos 2. Marco Teórico

La entropı́a de K-S describe la velocidad a la que un sistema dinámico crea o pierde

temporalmente información. Los sistemas caóticos, los sistemas estocásticos y los ci-
clos lı́mite están asociados con entropı́a de K–S positiva, infinita y cero, respectivamente
(Varsakelis & Anagnostidis, 2016).

Existe un lı́mite inferior para la entropı́a de K–S, generalmente denominada en-

tropı́a K2 , el cual es posible calcular mediante el algoritmo de Grassberger & Procaccia
(1983). La positividad de K2 (K2 > 0), constituye una condición suficiente para el com-
portamiento caótico (Grassberger & Procaccia, 1983).

2.4.7. Espacio de proyecciones de retraso

A continuación, se muestra la justificación teórica que permite estimar los expo-

nentes de Lyapunov y la entropı́a de K–S a través del espacio de proyecciones de retraso.

Considere la trayectoria que pasa a través de x0 del sistema (2.8). Sea τ > 0 y
constante llamado retraso, y avance sobre la trayectoria dando pasos de tamaño τ . Sea
m ∈ N fijo, y considere el nuevo espacio m-dimensional con los m vectores observados
cada τ tiempo. A este espacio se le conoce como el espacio de proyecciones de retraso.

El siguiente ejemplo pretende ilustrar el concepto de espacio de proyecciones de

retraso. Sean las ecuaciones (2.34) y (2.35) las soluciones a un sistema como el del péndulo
no amortiguado.

x = sen(x) (2.34)
y = cos(x) (2.35)

El espacio de fase de tal sistema se muestra en la Figura 2.5(a) , el cual es un cı́rculo

que revela la dependencia exacta de cada estado después de un ciclo. Con el espacio de
proyecciones de retraso, se puede reconstruir este retrato del espacio de fase usando sólo
una observación. Por ejemplo, con la variable y, y un retraso τ = 5 se reconstruye el

22
2. Marco Teórico Definiendo al caos

espacio de fase, tal como se observa en la Figura 2.5(b).

1 1

0.5 0.5

y(t+5)
0 0
y

-0.5 -0.5

-1 -1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
x y(t)

(a) Espacio de fase (b) Espacio de proyecciones de retraso

Figura 2.5: El espacio de proyecciones de retraso reconstruye el espacio de fase: con la

variable y, y un retraso τ = 5.

Una caracterı́stica notable de los sistemas caóticos es que en principio no importa

con cual variable se caracteriza. La mayorı́a de las propiedades como el máximo expo-
nente de Lyapunov o la dimensión del atractor, se encuentran en (casi) cualquier variable
y sus retrasos. Más aún, a veces no es necesario o incluso no es deseable reconstruir el
espacio de fase completo de la variable estudiada, porque con frecuencia la dimensión
del atractor será más pequeña que la dimensión de tal espacio (Sprott & Sprott, 2003).

Es suficiente con reconstruir un nuevo espacio en el cual un atractor equivalente

pueda ser embebido. Este procedimiento es llamado reconstrucción de espacio de pro-
yecciones de retraso (Sprott & Sprott, 2003) y fue propuesto por Packard et al. (1980) y
Takens (1981). El Teorema de Whitney 1 proporciona el fundamento teórico que permite
realizar tal reconstrucción (Parker & Chua, 1989).
Teorema 1. Teorema de Whitney. Si M es una variedad N -dimensional existe una función
f tal que f : M → R2 (embebe a M en R2 ).

Una variedad diferenciable es un objeto n dimensional en el que se puede hacer

análisis (derivar e integrar). Por ejemplo el cı́rculo es una variedad unidimensional, el
cual puede ser embebido en R2 pero no en R.

23
 Soluciones numéricas 2. Marco Teórico

2.5. Soluciones numéricas

En la presente sección se proporcionan las soluciones numéricas utilizadas para

aproximar a los sistemas caóticos. Para la mayorı́a de los sistemas caóticos no se tienen
o no se conocen sus soluciones analı́ticas por lo que se recurre a los métodos numéricos
para encontrar soluciones numéricas o aproximadas, los cuales son ejecutados compu-
tacionalmente. Sin embargo, las computadoras representan a las variables con precisión
finita y además discretizan el tiempo continuo. Encontrar la solución numérica de una
ecuación diferencial implica encontrar un mapeo que aproxime el flujo (ver definición
(2)) (Sprott & Sprott, 2003). Recordemos que un mapeo es una función que preserva es-
tructuras.

En el área de ecuaciones diferenciales, la ecuación (2.36) es conocida como el pro-

blema de valores iniciales (PVIs) para una ecuación diferencial ordinaria de primer
orden
y ′ = f (t, y), y(a) = η. (2.36)

Una propiedad esencial de la mayorı́a de los métodos numéricos para resolver (2.36)
es la discretización, la cual busca una solución aproximada, no en el intervalo continuo
a ≤ t ≤ b, sino en el conjunto discreto de puntos {tn | n = 0, 1, . . . , (b − a)/h} donde
la sucesión de puntos {tn } está definida por tn = a + nh, n = 0, 1, 2, . . ., con h > 0
llamado paso de integración. Es importante mencionar que a es la condición inicial y b es
el tiempo final en que se encontrará la solución.

2.5.1. Ensombrecimiento

La trayectoria de un sistema caótico es sensible a las condiciones iniciales y además

existen pequeños errores numéricos en el cálculo de la misma. En este sentido, ¿cómo
estamos seguros de que el atractor extraño no es el resultado de falsas observaciones

24
2. Marco Teórico Soluciones numéricas

numéricas?.

El lema del ensombrecimiento (shadowing, en inglés) establece que aunque la tra-

yectoria generada por la computadora no sea la teórica para el conjunto de condiciones
iniciales elegido, esta trayectoria ruidosa será parecida a otra trayectoria teórica con un
conjunto de condiciones iniciales ligeramente diferente al elegido. Se dice que la tra-
yectoria teórica ensombrece a la trayectoria calculada a partir de un modelo aproximado
(Grebogi et al., 1990). Si los errores arrojan a la trayectoria fuera del atractor, ésta rápi-
damente regresará y seguirá una de las infinitas trayectorias del atractor hasta que los
errores otra vez la expulsen del atractor (Sprott & Sprott, 2003).

En conclusión, el lema de ensombrecimiento permite obtener soluciones que per-

tenecen al atractor extraño y aunque no son las soluciones teóricas son muy parecidas a
éstas, además permite tener ciertas caracterı́sticas como el rango de los posibles valores
(Grebogi et al., 1990; Sprott & Sprott, 2003).

2.5.2. Métodos Multipasos Lineales

En los siguientes apartados proporcionamos los métodos numéricos que utilizamos

en el desarrollo de este trabajo.

Sea yn una aproximación para la solución teórica en tn , y(tn ), y sea fn ≡ f (tn , yn ).

Si el método numérico para determinar la sucesión de puntos {yn } toma la forma de una
relación lineal entre yn+j , fn+j , j = 0, 1, . . . , k, es llamado un método multipaso lineal
(MML) de k-pasos. El MML general puede escribirse acorde a la ecuación (2.37).

k
X k
X
αj yn+j = h βj fn+j , (2.37)
j=0 j=0

donde αj y βj son constantes, αk = 1, α02 + β02 6= 0.

25
 Soluciones numéricas 2. Marco Teórico

El problema de determinar la solución y(t) de la ecuación (2.36) es sustituido por el

de encontrar una sucesión yn que satisfaga la ecuación en diferencias (2.37). En general, y
en esta investigación, la función f (t, y) es no lineal por lo que la ecuación (2.37) será una
ecuación en diferencias no lineal. La ventaja de utilizar la ecuación (2.37) es que permite
calcular la sucesión yn numéricamente. Para encontrar la solución de la ecuación (2.37)
es necesario proporcionar un conjunto de valores iniciales y0 , y1 , . . . , yk−1 . Hairer et al.
(1993) y Lambert (1973) exhiben diversos métodos para tal propósito.

El método numérico representado por la ecuación (2.37) se dice que es explı́cito si

βk = 0 e implı́cito si βk 6= 0 (Lambert, 1973). Para un método explı́cito la ecuación (2.37)
genera el valor actual yn+k directamente de los términos yn+j , fn+j , j = 0, 1, 2, . . . , k −1,
cuyos valores en la etapa actual del cálculo ya han sido calculados. Por el contrario, si
un método es implı́cito no se realiza directamente, puesto que evalúa en cada etapa del
cálculo la ecuación (2.38) para obtener la solución de la ecuación (2.37).

yn+k = hβk f (tn+k , yn+k ) + g (2.38)

donde g es una función conocida de los valores previamente calculados yn+j , fn+j , j =
0, 1, 2, . . . , k − 1

Aún está presente el problema de determinar los coeficientes αj y βj que apare-

cen en la ecuación (2.37). Existen diversas técnicas para obtener tales coeficientes, por
ejemplo, utilizar la expansión de polinomios de Taylor. Los autores Hairer et al. (1993) y
Lambert (1973) presentan ésta y otras técnicas de manera detallada.

A continuación, ejemplificamos cómo mediante la expansión de los polinomios de

Taylor se obtienen los coeficientes de un método numérico en particular. La ecuación
(2.39) representa la expansión de Taylor para y(tn + h) en tn .

h2 (2)
y(tn + h) = y(tn ) + hy (1) (tn ) + y (tn ) + . . . (2.39)
2!
26
2. Marco Teórico Soluciones numéricas

Si se trunca la expansión en los primeros dos términos y se sustituye a y (1) (tn ) por
la ecuación diferencial (2.36) se obtiene una aproximación a la solución de la ecuación
(2.36) dada por la expresión (2.40).

y(tn + h) ≃ y(tn ) + hf (tn , y(tn )). (2.40)

El error en la aproximación está dado por la ecuación (2.41).

h2 (2) h3
y (tn ) + y (3) (tn ) + . . . (2.41)
2! 3!

Ahora, si en la ecuación (2.40) reemplazamos y(tn ), y(tn + h) por yn ,y yn + 1 se

obtiene la igualdad dada por la ecuación (2.42).

yn+1 = yn + hfn . (2.42)

La ecuación (2.42) representa un método de un paso lineal, conocido popularmente

como el método de Euler, el cual es el más simple de todos los métodos multipaso.

Métodos Adams Bashforth y Adams Moulton

Los MML para los que el lado derecho de la ecuación (2.37) es igual a α0 yn +
α1 yn+1 + α2 yn+2 , donde α0 = 0, α1 = −1 y α2 = 1 son llamados métodos Adams.
Los métodos Adams que son explı́citos son llamados Adams-Bashforth mientras que
los implı́citos son Adams-Moulton (Hairer et al., 1993).

La ecuación (2.43) representa un método Adams Moulton de 2 pasos.

h
yn+2 − yn+1 = (5fn+2 + 8fn+1 − fn ) (2.43)
12
27
 Soluciones numéricas 2. Marco Teórico

2.5.3. Métodos Runge-Kutta

La filosofı́a de Runge, después desarrollada por Kutta y Heun, es lograr métodos

numéricos de orden alto sacrificando la linealidad pero preservando sólo un paso a dife-
rencia de los MML (Lambert, 1973).

Considere el siguiente método numérico

1 1
yn+1 − yn = hf (tn + h, yn + hf (tn , yn )). (2.44)
2 2

Observe que a diferencia de los métodos antes mencionados, en la ecuación (2.44)

para encontrar a yn+1 es necesario evaluar de manera recursiva a la función f . Este tipo de
métodos es conocido como métodos Runge–Kutta (R–K). La ecuación (2.44) representa
un R–K de dos etapas.

La formalización matemática de tales métodos es la siguiente. Los métodos R–K

de R etapas explı́citos están definidos por

yn+1 − yn = hφ(tn , yn , h),

R
X
φ(t, y, h) = cr kr ,
r=1

k1 = f (t, y), (2.45)

r−1
!
X
kr = f t + har , y + h brs ks , r = 2, 3, . . . , R,
s=1
r−1
X
ar = brs , r = 2, 3, . . . , R.
s=1

donde brs y cr son constantes.

Además, se debe satisfacer la siguiente condición

R
X
cr = 1 (2.46)
r=1

28
2. Marco Teórico Soluciones numéricas

para que el método encuentre aproximaciones cercanas a la solución (Hairer et al., 1993).

Un método R–K de R etapas implica R evaluaciones funcionales por etapa. Cada

una de las funciones kr (t, y, h), r = 2, 3, . . . , R, puede ser interpretada como una apro-
ximación a la derivada y ′ (t), y la función φ(t, y, h) como una media ponderada de estas
aproximaciones, lo que significa que, en general, son métodos costosos.

2.5.4. Métodos de extrapolación

Sea y(t; h) una aproximación a la solución teórica y(t) del PVI representado por la
ecuación (2.36). Suponga que para un N fijo, y(t; h) tiene una expansión de la forma

y(t; h) = y(t) + A1 h + A2 h2 + . . . + AN hN + RN (h) (2.47)

RN (h) = O(hN +1 ) cuando h → 0 (2.48)

donde los coeficientes y(t), A1 , . . . , AN son independientes de h. De aquı́ en adelante se

resume esta afirmación escribiendo

y(t; h) ∼ y(t) + A1 h + A1 h + A2 h2 + . . . (2.49)

O(g(n)) denota el conjunto de todas las funciones p : N → R≥0 tales que p(n) ≤
cg(n) para todo n ≥ n0 con n0 ∈ N y c ∈ R≥0 .

Las ecuaciones (2.50) y (2.51) permiten calcular y(t; h0 ) y y(t; 12 h0 ).

y(t; h0 ) = y(t) + A1 h0 + A2 h20 + . . . (2.50)

 
1 1 1
y t; h0 = y(t) + A1 h0 + A2 h20 + . . . (2.51)
2 2 4

Note que y(t; h0 ) = y(t) + O(h0 ) y que y(t; 12 h0 ) = y(t) + O( 12 h0 ) cuando h → 0,

entonces  
1 1
2y t; h0 − y(t; h0 ) = y(t) − A2 h2 + . . . = O(h20 ). (2.52)
2 2

29
 Soluciones numéricas 2. Marco Teórico

Por lo tanto, el término izquierdo de la ecuación (2.52) es una mejor aproximación a y(t)
que y(t; h0 ) o y(t; 12 h0 ). Ésta es la idea básica de la extrapolación de Richardson (Lambert,
1973).

En general, se puede encontrar una combinación lineal con la propiedad

S
X
cs,S y(t; hs ) = y(t) + O(hS+1
s ), h → 0. (2.53)
s=0

La extrapolación polinomial se realiza para aproximar y(t) en los puntos básicos

a + jH, j = 0, 1, 2, . . ., hasta que a + jH ≥ b. H es llamado el paso de integración básico.
A continuación, se explica cómo se realiza este proceso.

Primero se elige un paso de integración h0 = H/N0 , donde N0 es un entero positi-

vo. Se aplica el método numérico N0 empezando en t = a para obtener una aproximación
y(a + H; h0 ) a la solución teórica y(a + H). Un segundo paso de integración h1 = H/N1 ,
con N1 un entero positivo mayor que N0 , es elegido y el método es aplicado N1 veces,
empezando otra vez en a, para obtener otra aproximación y(a + H; h1) a la solución
teórica y(a + H).

Procediendo de esta forma, para la sucesión de pasos de integración {hs } donde

hs = H/Ns , Ns | s = 0, 1, . . . , S, se obtiene una sucesión de aproximaciones {y(a +
H; hs ) | s = 0, 1, . . . , S} a la solución teórica y(a + H).

(0)
Tomando en cuenta la ecuación (2.49) y haciendo as = y(a + H : hs ), se calculan
los elementos del siguiente arreglo

(0)
h 0 a0
(0) (1)
h 1 a1 a0
(0) (1) (2)
h 2 a2 a1 a0
(0) (1) (2) (3)
h 3 a3 a2 a1 a0
.. .. .. .. ..
. . . . . ,

30
2. Marco Teórico Soluciones numéricas

donde
m−1
am−1 m−1
(s+1) − a(s)
am
(s) = a(s+1) + , m = 1, 2, . . . , s = 0, 1, . . . (2.54)
(hs /hm+s )2 − 1

Luego, se elige la última entrada de la diagonal principal del arreglo como la apro-
ximación final a y(a + H) y se denota por y ∗ (a + H; H).

Para obtener la solución numérica en el siguiente punto básico a + 2H, se aplica

el procedimiento descrito anteriormente al nuevo PVI y ′ = f (t, y) sujeto a y(a + H) =
y ∗ (a + H; H). Este proceso se repite hasta que a + jH ≥ b.

2.5.5. Métodos especiales

Si la ecuación (2.36) tiene alguna propiedad especial entonces se pueden diseñar

métodos numéricos especiales que aprovechen esta propiedad. Cuando las soluciones
son periódicas u oscilan con una frecuencia conocida (o una estimación de la misma). Se
han desarrollado una clase de métodos basados en polinomios trigonométricos propuesta
por Gautschi (1961) que es particularmente apropiada.

Sea T el periodo, conocido o estimado de la solución y ω = 2π/T la frecuencia.

Entonces el método basado en polinomios trigonométricos está dado por la ecuación
(2.55) (Gautschi, 1961).
k
X k
X
αj yn+j = h βj (v)fn+j v = ωh, αk = +1. (2.55)
j=0 j=0

Observe que la diferencia de la ecuación (2.55) con los MML (descritos en la sección
2.5.2) es que los coeficientes βj (v) se encuentran en función de h. El método basado en
polinomios trigonométricos que fue empleado en el capı́tulo 5 es el siguiente.

α0 yn + α1 yn+1 + α2 yn+2 = h(β0 (v)fn + β1 (v)fn+1 + β2 (v)fn+2 ), (2.56)

31
 Soluciones numéricas 2. Marco Teórico

donde α0 = 0, α1 = −1, α2 = 1, β2 (v) = 0,

cos(v) − 1 tan( v2 )
β0 (v) = =− , (2.57)
v sin(v) v
− cos(2v) + cos(v)
β1 (v) = , (2.58)
v sin(v)
quedando de la siguiente manera

yn+2 − yn+1 = h(β1 fn+1 + β0 fn ). (2.59)

Las ecuaciones (2.57) y (2.58) son las soluciones analı́ticas para β0 (v) y β1 (v). Sin
embargo expandiendo β0 y β1 por series de Taylor, consideremos las siguientes aproxi-
maciones para su implementación computacional
1 1 1 4
β0 (v) = − (1 + v 2 + v ) (2.60)
2 12 120
3 1 1 4
β1 (v) = (1 − v 2 + v ) (2.61)
2 4 120
Las aproximaciones a β0 y β1 podrı́an tomar más términos, pero para los fines de esta
investigación con ese número de términos es suficiente.

De aquı́ en adelante al método numérico dado por la ecuación (2.59) será nombrado
método de Gautschi. La implementación computacional del método de Gautschi no se
encuentra disponible en repositorios públicos en la web debido a que requiere el manejo
de estructuras de datos como colas FIFO, lo cual tiene su grado de complejidad.

La implementación que realizamos de este método se encuentra en el Apéndice B.

Se trata de un método multipaso donde la aproximación de y en el tiempo tn+2 requiere
las aproximaciones de y en los tiempos tn+1 y tn . Se sabe que y0 es igual a las condiciones
iniciales, y en la primera iteración (n = 0) de la implementación realizada se utilizó un
R–F 2 para encontrar el valor de y1 , porque de esta forma se preserva el orden del método.
Una vez calculado y1 , los valores y0 y y1 se almacenaron en una estructura de datos del

32
2. Marco Teórico Sistemas caóticos elegidos

tipo cola FIFO (en inglés First In, First Out) y el tiempo fue actualizado (n = n + 1).
A partir de la segunda iteración el método de Gautschi fue utilizado para aproximar
los valores de yn+2 , porque con la estructura de datos se guardaban sólo las dos últimas
aproximaciones de y en cada iteración del método necesarias para calcular y en el tiempo
actual.

2.5.6. Errores relativo y absoluto

Al hablar de aproximaciones es importante medir el error en los resultados. Existen

diferentes tipos de errores, pero los que interesan en los métodos numéricos son los
siguientes.

Sea y ∗ una aproximación a la solución y de la ecuación (2.36), el error absoluto

Ea está dado por la ecuación (2.62).

Ea =| y − y ∗ | . (2.62)

El error relativo Er viene dado por la ecuación (2.63).

Ea
Er = , con y 6= 0. (2.63)
y

2.6. Sistemas caóticos elegidos

En esta sección presentamos los modelos matemáticos de los sistemas caóticos que
utilizamos en el capı́tulo 5. Los sistemas caóticos los elegimos porque son frecuentemente
mencionados en la literatura y tienen caracterı́sticas afines a otros sistemas caóticos.

33
 Sistemas caóticos elegidos 2. Marco Teórico

2.6.1. Sistema de Lorenz

El sistema de Lorenz fue propuesto como un modelo simplificado de la dinámica

de la turbulencia inducida por la convección térmica en forma de anillos que ocurre en
la atmósfera (Sprott & Sprott, 2003). El sistema describe el comportamiento futuro que
caracteriza la corriente de convección en un depósito de fluido que se calienta por debajo
(Lorenz, 1963). En estas condiciones, el fluido es arrastrado, formando una especie de
carrete de fotos. Esta corriente se estabiliza si la temperatura de la fuente de calor no es
muy alta. En otro caso, la velocidad del flujo fluctúa, y puede invertirse, mostrando un
comportamiento caótico.

Las ecuaciones diferenciales que representan al sistema de Lorenz son

dx
= σ(y − x) (2.64)
dt
dy
= Rx − y − xz (2.65)
dt
dz
= xy + γz (2.66)
dt

donde σ, R, γ representan parámetros.

El sistema de Lorenz describe la turbulencia con tres estados (dos distribuciones

de temperatura y una de velocidad). La variable x es proporcional a la intensidad del
movimiento convectivo, y es proporcional a la diferencia de temperatura entre los ele-
mentos ascendente y descendente del fluido y z es proporcional a las desviaciones del
perfil vertical de temperatura respecto de su valor de equilibrio (Lorenz, 1963).

2.6.2. Atractor de Rossler

En 1970, Otto Rossler, un médico no practicante, se interesó en el caos y concibió

un sistema autónomo de tres dimensiones más simple que el de Lorenz. El atractor de

34
2. Marco Teórico Sistemas caóticos elegidos

Rossler. El cual está dado por las ecuaciones (2.68) a (2.70).

dx
= −y − z (2.67)
dt
dy
= x + ay (2.68)
dt
dz
= b + xz − cz, (2.69)
dt
(2.70)

donde a, b, c son parámetros reales.

El atractor de Rossler tiene el mismo número de términos que el sistema de Lorenz

pero sólo un término no lineal: xz (Sprott & Sprott, 2003; Kontorovich et al., 2009).

2.6.3. Sistema de Chen

El sistema de Chen viene dado por las ecuaciones (2.72) a (2.73).

dx
= a(y − x) (2.71)
dt
dy
= (c − a)x − xz + cy (2.72)
dt
dz
= xy − bz. (2.73)
dt

Los sistemas de Lorenz y Chen tienen términos no lineales en forma de productos

de dos variables dependientes. Aunque ambos sistemas son parecidos, análisis teóricos y
simulaciones numéricas sugieren que el atractor de Chen tiene estructuras topológicas
diferentes a las del sistema de Lorenz. También, es interesante notar que el exponente
de Lyapunov (positivo) del sistema de Chen es cercano a 2.0272, mientras que el corres-
pondiente exponente para el sistema de Lorenz es cercano a 0.9056. En consecuencia,
el sistema de Chen es más sensible a las condiciones iniciales que el sistema de Lorenz
(Eftekhari & Jafari, 2012).

35
 Sistemas caóticos elegidos 2. Marco Teórico

2.6.4. Sistema de Chua

El caos se puede observar en circuitos electrónicos. Un ejemplo popular es el cir-

cuito de Chua, el cual contiene un inductor y dos capacitores acoplados a un par de
amplificadores que proporcionan resistencia negativa no lineal (Sprott & Sprott, 2003).
Aunque existen diversas versiones del sistema de Chua, elegimos trabajar con el sistema
de Chua con tres enrollamientos.

dx
= −αx + αy − αG(x) (2.74)
dt
dy
= x−y+z (2.75)
dt
dz
= −βy (2.76)
dt
(2.77)

m1 (x + b3 ) − d3 x < −b3




m (x + b2 ) − d2 −b3 ≤ x < −b2


 2


G(x) = m1 x −b2 ≤ x < b2


m2 (x − b2 ) + d2 b2 ≤ x < b3






m (x − b ) + d

x > b3
1 3 3

donde α, β, m1 , m2 , m1 , b2 y b3 son constantes reales, ademas d2 = m1 b2 y d3 = d2 +

m2 (b3 − b2 ).

2.6.5. Ecuaciones de Rabinovich–Fabrikant

Las ecuaciones de Rabinovich–Fabrikant son un sistema de tres EDOs no lineales

acoplado, dado por las ecuaciones (2.78), (2.79) y (2.80). Este sistema fue propuesto por
Rabinovich y Fabrikant para describir ondas en un medio no equilibrado (Rabinovich &
Fabrikant, 1979).

dx
= y(z − 1 + x2 ) + γx (2.78)
dt
36
2. Marco Teórico Resumen de capı́tulo

dy
= x(3z + 1 − x2 ) + γy (2.79)
dt
dz
= −2z(α + xy) (2.80)
dt

donde α y γ son constantes positivas.

2.6.6. Ecuaciones de Lorenz–Stenflo

Las ecuaciones de Lorenz–Stenflo fueron introducidas por Stenflo para describir

las ondas atmosféricas (Stenflo, 1996). Están representadas por el siguiente sistema de
cuatro EDOs no lineales acopladas

dx
= σ(y − x) + sw (2.81)
dt
dy
= rx − y − xz (2.82)
dt
dz
= xy − βz (2.83)
dt
dw
= −x − sw (2.84)
dt

donde σ, r, β y s son constantes positivas.

2.7. Resumen de capı́tulo

En este capı́tulo revisamos conceptos claves que conducen a dos definiciones es-
pecı́ficas del caos puesto que no existe una definición universal del mismo.

También listamos las condiciones necesarias para la presencia de caos, de manera

no exhaustiva pero si las más aceptadas por la comunidad cientı́fica. Reportamos una
condición suficiente para la presencia de caos en los resultados numéricos de las simula-
ciones de sistemas dinámicos.

37
 Resumen de capı́tulo 2. Marco Teórico

Explicamos que las soluciones numéricas son aproximaciones a la solución teórica

del problema a resolver.

Ası́ mismo, proporcionamos las bases para entender los métodos numéricos que
utilizamos en el capı́tulo 5 y finalmente presentamos el conjunto de sistemas caóticos
elegidos para realizar las pruebas numéricas.

38
 3
Trabajo relacionado

En este capı́tulo presentamos una revisión de trabajos relevantes al propósito de

esta investigación. Primero abordamos las aplicaciones de los sistemas caóticos, los pro-
blemas que se presentan al realizar simulaciones haciendo énfasis en sistemas caóticos y
posteriormente presentamos los métodos numéricos que se han empleado para la solu-
ción numérica de sistemas los caóticos.

3.1. Aplicaciones de los sistemas caóticos

Los sistemas caóticos tienen aplicación en muchos dominios. Sin ser exhaustivos
repasamos algunos con la idea de dejar ver su versatilidad.

El comportamiento caótico se presenta en neuronas y en toda la actividad cerebral.

En el ámbito de redes neuronales, se han propuesto algoritmos que utilizan la dinámica
caótica de los osciladores de Stuart-Landau como neuronas artificiales para reconocer
patrones fractales generados aleatoriamente y se realizan aplicaciones en las que aparece
la geometrı́a fractal, como la identificación de áreas u objetos en la superficie de la Tierra
en imágenes satelitales o la clasificación de tumores cancerı́genos en imágenes médicas
(da Silva & Zhao, 2016).

En criptografı́a, se han propuesto esquemas de cifrado para las comunicaciones pri-

 Aplicaciones de los sistemas caóticos 3. Trabajo relacionado

vadas, como el esquema de enmascaramiento caótico (Belazi et al., 2016; Soriano-Sánchez

et al., 2016; Garcı́a-Martı́nez et al., 2015; Behnia et al., 2008). Recientemente, se han for-
mulado nuevas técnicas de cifrado y esquemas de comunicación. Algunas generan se-
cuencias pseudo-aleatorias iterando mapas caóticos con alta sensibilidad a los valores
iniciales y aleatoriedad (Wu et al., 2017). Algunas otras se centran en la robustez y la
seguridad (Silva-Garcı́a et al., 2018).

Las señales caóticas generadas por un sistema caótico representan una atractiva
clase de señales para modelar fenómenos fı́sicos. La detección, estimación, análisis y
caracterización de este tipo de señales representa un desafı́o significativo en el ámbito
del procesamiento de señales. Las señales caóticas han sido observadas en el ruido blanco,
y una forma de estimarlas es mediante la estimación bayesiana (Pantaleón et al., 2003).
Otros trabajos separan las señales caóticas del ruido no lineal (Hasler et al., 2002).

En optimización, se investiga la técnica de optimización del problema del viajero

(TSP, por sus siglas en inglés) mediante redes neuronales que utilizan dinámica caótica
(Tokuda et al., 1997). Otros trabajos introducen sistemas caóticos al recocido simulado
que reemplazan a las distribuciones normales, proponiendo un algoritmo de optimización
denominado recocido simulado de caos (CSA, por sus siglas en inglés). Los resultados de
simulación de la optimización de funciones complejas tı́picas muestran que CSA mejora
la convergencia y es eficiente, aplicable y fácil de implementar (Mingjun & Huanwen,
2004).

En biologı́a, los sistemas caóticos se utilizan para modelar el comportamiento de

poblaciones (Nee, 2018).

En quı́mica existen reacciones que presentan comportamiento caótico las cuales

son modeladas y analizadas por los investigadores (Olabodé et al., 2018).

En fı́sica, procesos fı́sicos de la dinámica de fluidos y termodinámica son modelados

mediante sistemas caóticos (Lorenz, 1963; Deutsch, 2018).

40
3. Trabajo relacionado Problemas en simulaciones caóticas

3.2. Problemas en simulaciones caóticas

Al simular sistemas dinámicos caóticos se presentan diversos problemas debidos

a la abstracción. La abstracción se da en cuatro niveles: la realidad fı́sica del problema
estudiado, el modelo continuo de esa realidad fı́sica, la discretización numérica de ese
modelo matemático y la simulación en punto flotante de tal discretización (Corless, 1994).
En cada una de las etapas se presentan diversos problemas. A continuación, se describen
los problemas comúnmente encontrados y los trabajos que han tratado de solucionarlos
o aliviarlos.

3.2.1. Alta sensibilidad a las condiciones iniciales

Ya se introdujo en la sección 2.4 la caracterı́stica inherente a los sistemas caóti-

cos de ser altamente sensibles a las condiciones iniciales. Una pequeña diferencia en las
condiciones iniciales podrı́a conducir a una gran diferencia en las simulaciones genera-
das por computadora después de un largo perı́odo de tiempo (Liao, 2017). Este problema
es popularmente conocido como efecto mariposa (Lorenz, 1963). Diversos autores han
reportado este problema: Góra & Boyarsky (1988); Corless et al. (1991); Corless (1994);
Varsakelis & Anagnostidis (2016).

Actualmente, existe una técnica llamada Simulación numérica limpia (CNS por
sus siglas en inglés) que trata de subsanar este problema. La CNS está basada en series
de Taylor con al menos 50 términos de aproximación y utiliza datos con 300 dı́gitos de
precisión. La CNS evalúa los resultados con otra serie de Taylor con un orden mayor al
utilizado (mayor a 50). La CNS fue implementada en supercomputadoras (Liao, 2017).

Un reto que aún queda abierto es aliviar el problema de alta sensibilidad con me-
nos requerimientos de hardware que el utilizado por las supercomputadoras. Una de las
razones es que, por ejemplo, las aplicaciones criptográficas que se mencionaron en la
sección 3.1 no requieren tanto poder de cómputo como el de una supercomputadora.

41
 Problemas en simulaciones caóticas 3. Trabajo relacionado

3.2.2. Errores de discretización: truncamiento y redondeo

Existen al menos dos tipos de discretización: en tiempo y en espacio. La prime-

ra consiste en la transición de las ecuaciones diferenciales a ecuaciones en diferencias
(también es conocida como error de truncamiento). Mientras que la discretización en es-
pacio es la transición de un espacio de fase continuo a uno discreto (un conjunto finito
de puntos) (Blank, 1994).

Cuando se resuelven sistemas dinámicos mediante una computadora se presentan

perturbaciones causadas por discretizaciones en tiempo y espacio. En tiempo porque se
utilizan métodos numéricos para resolverlos. Y en espacio porque la computadora utiliza
la aritmética de punto flotante, la cual pretende representar a la aritmética de los números
reales en un subconjunto finito de los números racionales. Un tipo de error particular que
ocasiona la aritmética de punto flotante es el error de redondeo (Fryska & Zohdy, 1992;
Corless, 1994; Polhill et al., 2006).

Las consecuencias de los errores de redondeo pueden ser drásticas. Por ejemplo,
una solución periódica puede ser reemplazada por otra periódica con un múltiplo del
periodo de la solución original o una solución caótica puede ser sustituida por una pe-
riódica. Este último fenómeno es conocido como supresión de caos (Blank, 1994). También
puede pasar que la solución a un sistema que no es caótico sea reemplazada por una so-
lución caótica, efecto conocido como caos falso o caos computacional (Yamaguti & Ushiki,
1981; Corless, 1994).

Los errores de redondeo son inevitables para cualquier simulación numérica (Blank,
1994; Fryska & Zohdy, 1992; Polhill et al., 2006), y para aliviar este tipo de errores no es
suficiente con variar la precisión utilizada (Corless, 1994; Liao, 2017). Algunas de las
técnicas más utilizadas para tratar de reducir los efectos de los errores de redondeo son
las siguientes. Reescribir las expresiones equivalentes en la aritmética real pero no ne-
cesariamente en la aritmética de punto flotante para que sean implementados en código
computacional. Escribir los números de punto flotante en cadenas (strings) con un núme-

42
3. Trabajo relacionado Problemas en simulaciones caóticas

ro especı́fico de decimales o utilizar precisión extendida (Polhill et al., 2006), entre otras.
Pero como se mencionó, estas técnicas sólo reducen los efectos de redondeo, pero no los
erradican por completo.

Debido a que la supresión de caos y el caos falso son problemas que pueden apa-
recer en las simulaciones numéricas, las siguientes subsecciones son dedicadas a estos
problemas desde la perspectiva del error de truncamiento.

3.2.3. Supresión de caos

Al simular sistemas caóticos los métodos numéricos pueden suprimir el caos real.
Por ejemplo, el método implı́cito de Euler causa este problema cuando se aplica al sistema
de Rossler (Corless et al., 1991) o el oscilador de Duffing (De Markus, 2001).

En un estudio teórico se discute cómo los métodos de integración implı́citos tipo

Euler pueden producir resultados discretos no caóticos de sistemas caóticos represen-
tados por ecuaciones diferenciales, y mediante simulaciones numéricas del sistema de
Rossler con el método de Euler hacia atrás se proporciona evidencia empı́rica que reafir-
ma el estudio teórico (Grote & Meyer-Spasche, 1998).

Recientemente, otro análisis teórico realizado con un sistema caótico genérico y

el método de Euler hacia atrás concluye que para pasos de integración suficientemente
grandes se suprime el caos (Varsakelis & Anagnostidis, 2016).

A pesar de que este problema está presente cuando se simulan sistemas caóticos,
el problema continúa estando poco estudiado.

43
 Problemas en simulaciones caóticas 3. Trabajo relacionado

3.2.4. Caos falso

El caos falso es la introducción de comportamiento caótico en las soluciones de un

sistema no caótico.

En una investigación se muestra empı́ricamente con dos sistemas no caóticos que el

uso de pasos de integración incorrectos usando los esquemas numéricos de Euler, Runge-
Kutta y Taylor puede conducir al caos falso únicamente debido al método numérico uti-
lizado (Lorenz, 1989).

En un estudio teórico-empı́rico se manifiesta que se puede presentar un compor-

tamiento caótico en los resultados numéricos utilizando un método Runge-Kutta de una
ecuación diferencial no lineal, aunque no exista caos en la ecuación diferencial original,
debido al tamaño de paso utilizado (Hirai & Adachi, 1994).

En otros trabajos teóricos se demuestra que los métodos numéricos basado en es-
quemas de diferencias como el Euler hacia adelante pueden producir soluciones falsas
y tienden a volverse inestables para pasos de integración mayores a un umbral teórico
que depende del sistema y del método numérico utilizado (Grote & Meyer-Spasche, 1998;
Varsakelis & Anagnostidis, 2016). También se reporta que tales soluciones son estables
solo para pasos de integración menores a un umbral teórico (Varsakelis & Anagnostidis,
2016).

En otra investigación teórica-empı́rica se presenta un ejemplo en el que el método

de Runge-Kutta-Felberg (paso variable) no integra correctamente un acoplamiento de dos
osciladores simples, porque presenta un comportamiento caótico cuando se sabe que el
sistema no lo es. En este caso, el Runge-Kutta-Felberg falla porque permite utilizar pasos
de integración mayores a 6 que provocan que no se capture la dinámica real del sistema
simulado (Skufca, 2004).

El problema del caos falso podrı́a presentarse en métodos más elaborados que los
del tipo Euler y Runge–Kutta, por lo que este problema aún está siendo estudiado.

44
3. Trabajo relacionado Métodos numéricos empleados para resolver sistemas caóticos

3.2.5. Duración de las simulaciones

La duración de las simulaciones numéricas de sistemas caóticos es importante por

el contexto de las aplicaciones mencionadas en la sección 3.1, pero a menudo los méto-
dos numéricos colapsan antes de poder lograr simulaciones largas debido a los errores
de truncamiento cometidos por el método o los errores de redondeo causados por la
computadora.

El problema de las simulaciones largas ha sido abordado en una investigación

empı́rica que trata de responder a la pregunta de cuánto tiempo duran las simulaciones
numéricas caóticas. Mediante un ejemplo de un sistema mecánico (doble rotor pateado)
se muestra que las soluciones largas generadas por computadora no se correspondı́an
con las soluciones verdaderas, este fenómeno era causado por un exponente de Lyapu-
nov que fluctuaba alrededor de cero. La longitud de la simulación está relacionada con
los exponentes de Lyapunov: entre más cercanos a cero se encuentren, más cortas son
las simulaciones (Sauer et al., 1997).

La técnica CNS, explicada en la sección 3.2.1, también pretende ser una solución al
problema de las simulaciones largas (Liao, 2017).

Los retos que aún quedan abiertos acerca de este problema son: extender la inves-
tigación de los exponentes de Lyapunov mediante análisis teóricos, ası́ como plantear
soluciones que no requieran tanto poder de cómputo como el de las supercomputadoras.

3.3. Métodos numéricos empleados para resolver siste-

mas caóticos

A continuación se presenta una recopilación de los métodos numéricos utilizados

para simular sistemas caóticos con el propósito de observar cuáles son los más utiliza-

45
 Métodos numéricos empleados para resolver sistemas caóticos 3. Trabajo relacionado

dos y por qué. Primero se proporciona una colección de trabajos orientados a aplicar los
resultados de las simulaciones de sistemas caóticos. La siguiente sección agrupa a los tra-
bajos que proponen métodos numéricos novedosos y analizan el rendimiento al resolver
sistemas caóticos. Finalmente, persiguiendo el objetivo principal de la tesis, indagamos
sobre los métodos numéricos utilizados para integrar sistemas oscilatorios.

3.3.1. Orientados a las aplicaciones

Los trabajos presentados a continuación tienen como objetivo principal utilizar los
resultados numéricos de sistemas caóticos en aplicaciones especı́ficas dando por hecho
que los resultados numéricos son confiables.

En aplicaciones criptográficas se utilizan los métodos de integración Euler (Silva-

Garcı́a et al., 2018) y Runge–Kutta 4 (Zhi-Hong et al., 2005; Garcı́a-Martı́nez et al., 2015),
además otro trabajo reporta que utiliza Matlab pero sin especificar el método numérico
(Wu et al., 2017).

En electrónica, estudiando circuitos eléctricos se utilizó el Runge–Kutta 4 (Sprott,

2000).

En fı́sica, para estudiar ecuaciones de reacción-difusión se utilizó el Runge–Kutta

4 (Owolabi & Atangana, 2017), y para estudiar las caracterı́sticas de vibración no lineal
de un rotor de fricción se empleó métodos del tipo Runge–Kutta (Chu & Zhang, 1998).

En meteorologı́a, para un estudio del clima se utilizaron los métodos del tipo Euler
y Runge-Kutta (Lorenz, 1989).

Aunque esta revisión no es exhaustiva, con ella se pretende dar una idea de los
métodos numéricos comúnmente utilizados cuando se trata de realizar aplicaciones uti-
lizando sistemas caóticos.

46
3. Trabajo relacionado Métodos numéricos empleados para resolver sistemas caóticos

3.3.2. Analizando los métodos numéricos

El objetivo principal de los siguientes trabajos es encontrar métodos numéricos

novedosos que permitan solucionar sistemas caóticos.

Los métodos Runge–Kutta 4 y el método de descomposición de Adomian (ADM,

por sus siglas en inglés) fueron comparados al integrar el sistema de Chen, concluyendo
que el ADM funciona bien para el sistema de Chen, el cual es altamente caótico. El inter-
valo de integración en el que realizaron las simulaciones fue: [0, 7] (Noorani et al., 2007).
Es conveniente probar el ADM en un mayor número de sistemas caóticos para observar
su eficiencia.

Otra propuesta utiliza el método de cuadratura diferencial (DQ por sus siglas en
inglés) para resolver los sistemas dinámicos caóticos: Lorenz, Chen, Genesio y Rössler.
En esta investigación se comparan las soluciones del método de integración DQ y las del
Runge–Kutta 4, y se reporta que el método DQ produce una precisión mejor que el R–K
4 usando pasos de integración mayores. De los resultados obtenidos, el estudio concluye
que el esquema de integración de tiempo de DQ es fiable y computacionalmente eficiente.
Las simulaciones reportadas son de hasta 30s. La investigación enfatiza el hecho de tener
cuidado al aplicar el método DQ a los sistemas caóticos porque el esquema de integración
DQ puede producir resultados inexactos para sistemas caóticos con un paso de tiempo
demasiado grande (Eftekhari & Jafari, 2012).

En otro trabajo se propone emplear métodos de elemento espectral para sistemas

dinámicos temporales no lineales que presentan comportamiento caótico. Este enfoque
requiere un alto costo con respecto al tiempo CPU empleado (Bar-Yoseph et al., 1996).

Al estudiar los movimientos de un sistema de un rotor se empleó la técnica de

integración numérica de tipo implı́cito Newmark-β con el método de Newton Raphson
para resolver las ecuaciones diferenciales no lineales de tal sistema de forma iterativa.
Encontraron que dicho sistema presenta movimientos cuasi periódicos, subarmónicos y
caóticos (Harsha et al., 2003). Sobre esta técnica numérica se ha mostrado que el método

47
 Métodos numéricos empleados para resolver sistemas caóticos 3. Trabajo relacionado

de aceleración promedio de Newmark produce resultados falsos en la simulación de un

sistema de un solo grado de libertad con un resorte bilineal y es sensible tanto al paso de
tiempo de integración como al número de pasos de tiempo por ciclo de forzado Koh &
Liaw (1991).

Path Integration (PI) es un método numérico que emplea ecuaciones diferenciales

estocásticas para investigar la dinámica de osciladores no lineales que pueden exhibir
una respuesta caótica a la excitación armónica. En un estudio empı́rico se valida PI con
el oscilador de Duffing y se concluye que el método PI tiene el potencial de ser una herra-
mienta efectiva para investigar las dinámicas de los osciladores no lineales que exhiben
respuesta caótica. Sin embargo, PI añade cálculo porque utiliza un R–K de orden cuatro
junto con operaciones de probabilidad (Naess & Skaug, 2002).

Otro método numérico empleado para resolver el sistema de Chen es el método

de colocación de Bessel (Yüzbas, 2012). Dicho método numérico consiste en reducir el
problema a un sistema de ecuaciones algebraicas no lineales al expandir las soluciones
aproximadas por medio de polinomios de Bessel con coeficientes desconocidos. Con la
ayuda de los puntos de colocación y las operaciones matriciales de las derivadas, se cal-
culan los coeficientes desconocidos de los polinomios de Bessel. El método requiere hacer
algunos cálculos matriciales antes y al momento de su implementación y se requieren
polinomios hasta de orden 15, lo cual implica un costo computacional alto.

El método de relajación espectral multietapa (MSRM, por sus siglas en inglés) es

propuesto para resolver PVI no lineales con propiedades caóticas. Este método resuelve
sistemas caóticos como Lorenz, Chen, Liu, Rikitake, Rössler, Genesio-Tesi y Arneodo-
Coullet. La precisión y la validez del método propuesto se prueba con los métodos ba-
sados en Runge-Kutta y Adams-Bashforth-Moulton, que se encuentran en Matlab. Los
resultados numéricos indican que el MSRM es un método preciso, eficiente y confiable
para resolver problemas de valor inicial con un comportamiento caótico. Las simulacio-
nes realizadas son de hasta 100s (Motsa et al., 2013).

El método Bubnov-Galerkin basado en aproximaciones spline, que tienen como

48
3. Trabajo relacionado Métodos numéricos empleados para resolver sistemas caóticos

base polinomios trigonométricos, fue utilizado para analizar numéricamente el compor-

tamiento cualitativo de las soluciones del problema de Kolmogorov (determinar el flujo
de un fluido viscoso incompresible). Tal problema está formado por ecuaciones parciales
y presentan un comportamiento caótico. En la investigación justifican que se utilizó el
método Bubnov-Galerkin porque presenta alta precisión y buenas aproximaciones (Evs-
tigneev et al., 2015). Este trabajo es mencionado porque ofrece la posibilidad de obtener
resultados confiables al aplicar métodos basados en polinomios trigonométricos.

3.3.3. Métodos numéricos para problemas oscilatorios

Muchos sistemas caóticos son oscilatorios, por lo que utilizar métodos numéricos
para integrar sistemas caóticos podrı́a resultar útil. A continuación se presenta una re-
copilación de trabajos haciendo énfasis en los métodos numéricos utilizados.

Una caracterı́stica deseable de un método numérico para resolver EDOs es que

aproveche las propiedades especı́ficas que la solución pueda tener (Paternoster, 1998).
Siguiendo esta idea, en la investigación desarrollada por Gautschi (1961) se propone
utilizar métodos numéricos basados en polinomios trigonométricos que requieren con
antelación una estimación de la frecuencia para resolver sistemas oscilatorios. Tiempo
después, otro trabajo explotó esta idea para derivar métodos de tipo Nyström y Milne-
Simpson. Los resultados fueron comprados con los métodos numéricos propuestos por
Gautschi (1961): tipo Adams y Stormer, al simular sistemas oscilatorios. El estudio con-
cluye que a diferencia de los métodos propuesto por Gautschi, los esquemas numéricos
propuestos son insensibles a la sub-estimación y sobre-estimación de la frecuencia del
sistema simulado (Neta & Ford, 1984).

Un estudio compara dos métodos sofisticados para resolver sistemas oscilatorios.

Los métodos se basan en los algoritmos Crank-Nicolson y Gear, ambos son semi-implı́ci-
tos y el de Gear es de orden variable y paso variable. Se simula numéricamente sistemas
de reacción-difusión. Con ambos métodos, la solución tiene la misma forma cualitativa,

49
 Métodos numéricos empleados para resolver sistemas caóticos 3. Trabajo relacionado

aunque existe diferencia en la forma en que convergen. Además, el sistema simulado tie-
ne un comportamiento que es temporalmente caótico, por lo que una solución numérica
no puede representar más que su forma cualitativa. El trabajo concluye que la diferencia
de convergencia entre los dos esquemas numéricos es beneficiosa, ya que los métodos so-
lo coinciden cuando ambos han convergido (Sherratt, 1997). Ésto apoya al punto de vista
de Parker & Chua (1989) acerca de utilizar al menos dos métodos para simular numéri-
camente sistemas caóticos. En general, la comparación de las soluciones dadas por dos
métodos proporciona una prueba simple para la convergencia.

En otra investigación se utilizan diversas modificaciones al esquema de integración

R–K para integrar numéricamente sistemas oscilatorios, pero sus implementaciones son
costosas computacionalmente, y sólo integran intervalos de a lo más 100s (Franco, 2004).

Otra propuesta utiliza funciones trigonométricas (en lugar de algebraicas como en

integradores clásicos). Desarrolla una nueva clase de métodos de Runge-Kutta multi-
etapa, de un solo paso, tamaño de paso variable y coeficientes variables implı́citos para
resolver problemas EDOs oscilatorios (Nguyen et al., 2007).

En otro estudio se propone emplear un método numérico basado en funciones lla-

madas Φ, para integrar problemas altamente oscilatorios con un tamaño de paso grande.
Los resultados se comparan con métodos sofisticados como el ODE (Ordinary Diferential
Equation) de Gear obteniendo buenos resultados (Garcı́a-Alonso & Reyes, 2009).

El método Haar wavelet fue propuesto para obtener la solución numérica aproxi-
mada de la ecuación de Van der Pol (la cual es oscilatoria con cierto conjunto de paráme-
tros). Los resultados se comparan contra un método explı́cito y se obtienen resultados
más precisos. Simulan periodos de tiempo de a lo más 4 segundos (Saha Ray & Patra,
2013).

Existe todo un abanico de métodos para resolver problemas oscilatorios, pero de

acuerdo a los propósitos de esta tesis los métodos basados en polinomios trigonométri-
cos parecen ser los más adecuados porque capturan las propiedades oscilatorias de los

50
3. Trabajo relacionado Resumen de capı́tulo

sistemas caóticos oscilatorios.

3.4. Resumen de capı́tulo

Parker & Chua (1989) indican que las simulaciones de sistemas caóticos requieren
una interpretación cuidadosa y siempre deben verificarse. Si es posible, se deben usar
dos o más rutinas diferentes para integrar el mismo sistema. Como las rutinas introdu-
cen errores diferentes y los errores se propagan de manera diferente, ésta es una buena
verificación de la validez de la integración. La intuición y el razonamiento fı́sico también
deben usarse para verificar la razonabilidad de un resultado de integración.

Con el propósito de reducir los efectos de los errores de redondeo, en esta inves-
tigación, las implementaciones de los métodos numéricos utilizadas fueron tomadas de
repositorios como Netlib UTK & ORNL (2005) y ACM (2005), los cuales proporcionan
código validado y verificado.

Para solucionar el problema de la supresión de caos o del caos falso es importante

apoyarse de los análisis teóricos de estabilidad como los realizados en las investigaciones
de Grote & Meyer-Spasche (1998) y Varsakelis & Anagnostidis (2016).

De acuerdo al trabajo reportado por Evstigneev et al. (2015) parece ser convenien-
te utilizar los métodos numéricos basados en polinomios trigonométricos para simular
sistemas caótico oscilatorios.

51
 4
Planteamiento del problema y
metodologı́a propuesta

En este capı́tulo se aborda la importancia del problema planteado y la solución

propuesta de forma detallada. También, se explica la forma en que aproximamos los ex-
ponentes de Lyapunov y la entropı́a de K–S.

4.1. Importancia del problema dentro de las ciencias compu-

tacionales

Los sistemas caóticos tienen aplicaciones en criptografı́a. Una propuesta plantea

un algoritmo de cifrado de imágenes basado en caos. La idea principal del algoritmo es la
siguiente: intercambiar las posiciones de los pı́xeles de la imagen en el dominio espacial y
simultáneamente cambiar los valores en escala de grises de la imagen. Los valores de los
pı́xeles son cambiados mediante el sistema caótico de Chen, el cual es integrado numéri-
camente con el esquema Runge-Kutta. La alta sensibilidad a las condiciones iniciales del
sistema de Chen es aprovechada, puesto que las condiciones iniciales del sistema de Chen
forma parte de la llave privada para cifrar y descifrar la imagen (Zhi-Hong et al., 2005).
4. Planteamiento del problema y metodologı́a propuesta Metodologı́a propuesta

Sistemas dinámicos continuos

Los sistemas dinámicos continuos que se eligieron están representados por EDOs
y se dividen en dos grupos, los cuales se describen a continuación.

Primer grupo Comprende un conjunto de sistemas caóticos, que son comúnmente re-
portados en la literatura y con caracterı́sticas afines a otros sistemas caóticos: el
sistema de Lorenz, el atractor de Rössler, el sistema de Chen y el sistema de Chua
con 3 enrollamientos.

Segundo grupo Está conformado por sistemas dinámicos continuos que presentan com-
portamiento caótico o no caótico dependiendo de los valores que se asignen a sus
parámetros. Tales sistemas son: el sistema de Lorenz, las ecuaciones de Rabinovich-
Fabrikant y las ecuaciones de Lorenz-Stenflo.

Métodos de integración

El racional para elegir los métodos de integración fue elegir los más populares para
sistemas caóticos, ası́ como los comúnmente aplicados para resolver sistemas dinámicos
representados por EDOs que no presentan comportamiento caótico.

Los métodos de integración que seleccionamos se dividen en tres grandes grupos:

estándares, sofisticados y especiales, los cuales son descritos a continuación.

Métodos estándares En este trabajo, llamamos métodos estándares a los esquemas numéri-
cos que utilizan paso de integración fijo. Los algoritmos que elegimos fueron: un
Euler hacia atrás con dos correcciones (B–E, por sus siglas en inglés) y un Euler
hacia adelante (F–E, por sus siglas en inglés ). Se escogieron estos esquemas porque
son de los más sencillos en cuanto a número de operaciones y la literatura confirma
que el primero suprime el caos de los sistemas caóticos mientras que el segundo

55
 Metodologı́a propuesta 4. Planteamiento del problema y metodologı́a propuesta

genera caos falso en sistemas no caóticos (Corless et al., 1991; Corless, 1994; Var-
sakelis & Anagnostidis, 2016). Estos algoritmos se tomaron del repositorio Netlib
(UTK & ORNL, 2005).

Métodos sofisticados En este trabajo, nombramos métodos sofisticados a aquellos esque-

mas que tienen paso de integración y orden variable con o sin control de error. Las
estrategias numéricas que se eligieron son las siguientes: un Runge–Kutta Felbergh
45 (R–K F 45), un método de extrapolación (ODEX) y un método Adams Bashforth -
Adams Moulton (AB-AM). Los tres esquemas numéricos cuentan con paso variable
y control de error, y el AB-AM es de orden variable. Tales algoritmos son utiliza-
dos para integrar numéricamente sistemas de ecuaciones diferenciales de primer
orden (Hairer et al., 1993). Estas estrategias fueron tomadas del repositorio Netlib
(UTK & ORNL, 2005) y de (ACM, 2005).

Métodos especiales En este trabajo, los métodos especiales son aquellos que tienen paso
de integración fijo, pero están basados en polinomios trigonométricos. Se eligió el
método de Gautschi de segundo orden (ver sección 2.5.5) porque fue pensado para
sistemas oscilatorios (Gautschi, 1961), y parte de la hipótesis de esta investigación
es confirmar que este tipo de métodos capturan la propiedad oscilatoria de los
sistemas caóticos oscilatorios.

Es importante mencionar que en esta investigación, se utiliza código computacio-

nal de repositorios como Netlib (UTK & ORNL, 2005) y (ACM, 2005) porque es código
validado y verificado, por lo que ayuda a reducir los efectos de los errores de redondeo
(explicados en la sección 3.2.2).

4.2.2. Simulación

En esta etapa se simularon los sistemas dinámicos continuos elegidos con los méto-
dos numéricos propuestos. Se realizaron tres experimentos. A continuación se explica el
propósito de cada uno.

56
4. Planteamiento del problema y metodologı́a propuesta Metodologı́a propuesta

Primer experimento El objetivo es establecer casos en que los métodos numéricos su-
primen el caos en los resultados numéricos de los sistemas dinámicos del segundo
grupo de prueba que presentan comportamiento caótico.

Segundo experimento El objetivo es establecer casos en que los métodos numéricos ge-
neran caos falso en los resultados numéricos de sistemas dinámicos del segundo
grupo de prueba que presentan comportamiento no caótico.

Tercer experimento El objetivo es conservar el caos del primer grupo de los sistemas
caóticos de prueba cuando se simulan numéricamente con los métodos estándares
(sólo el Euler hacia atrás), sofisticados y especiales. Y mantenerlos durante grandes
periodos de tiempo.

4.2.3. Evaluación

En esta sección se explica la forma en que evaluamos las simulaciones numéricas.

Antes de continuar, es importante resaltar que interpretar los resultados numéricos de
un sistema caótico no es un trabajo meramente algorı́tmico, debemos contar con cierta
intuición y conocimiento para ello, tal como lo indican las investigaciones de Hegger
et al. (2007); Kantz & Schreiber (2003). Con tal premisa, el primer y segundo experimento
únicamente fueron evaluados cualitativamente mediante el espacio de fase.

Para evaluar el tercer experimento, primero se seleccionaron los resultados numéri-

cos que lograron terminar todo el intervalo de simulación. Posteriormente, estos resul-
tados fueron evaluados con las siguientes métricas (explicadas a detalle en el capı́tulo
2).

1. Cualitativas

a) Espacio de fase: Sugiere el caos presente en los resultados numéricos. Per-

mite observar los atractores extraños gráficamente.

57
 Metodologı́a propuesta 4. Planteamiento del problema y metodologı́a propuesta

b) Periodograma: Sugiere el caos presente en los resultados numéricos. Cuan-

do el sistema es caótico, el periodograma presenta una amplia gama de fre-
cuencias y posee picos en las más altas.

2. Cuantitativas

a) Número de evaluaciones: En esta investigación se calculó una aproxima-

ción del costo computacional mediante el número de evaluaciones realizadas
por el método numérico al modelo que representa al sistema caótico.

b) Exponentes de Lyapunov: Es una medida cuantitativa necesaria pero no

suficiente (ver capı́tulo 2). Se calculó con el programa TISEAN (Hegger et al.,
2007).

c) Entropı́a de Kolmogorov-Sinai: Medida cuantitativa suficiente pero no ne-

cesaria (ver capı́tulo 2). Se calculó con el programa TISEAN (Hegger et al.,
2007).

4.2.4. Exponentes de Lyapunov

El cálculo de los exponentes de Lyapunov en los resultados de las simulaciones

numéricas es una tarea que presenta obstáculos debido a que se deben interpretar co-
rrectamente las aproximaciones. A pesar de que existe una amplia gama de técnicas para
obtenerlos, en este trabajo se eligieron las presentadas por TISEAN (Hegger et al., 2007)
por estar ampliamente documentadas.

En esta sección se explica con un ejemplo la forma en que se aproximaron los

exponentes de Lyapunov y el espectro de los mismos (ver sección 2.4.2 para su definición)
utilizando las rutinas lyap k y lyap spec de TISEAN (Hegger et al., 2007).

Con ayuda del siguiente ejemplo se explica el uso de las rutinas lyap k y lyap spec.
Considere el sistema dinámico discreto dado por las ecuaciones (4.1-4.3) (Kantz & Schrei-

58
4. Planteamiento del problema y metodologı́a propuesta Metodologı́a propuesta
5000

1 4000

3000

z(t)
0

2000
-1

1000

1
0 0
-1
-1.5 -1 -0.5 0 0.5 1 1.5
y(t) x(t)

Figura 4.2: Espacio de fase del sistema dinámico discreto de las ecuaciones (??,2,3).

ber, 2003).

xn+1 = 1.76 − yn2 − 0.1zn (4.1)

yn+1 = xn (4.2)
zn+1 = yn , (4.3)

sujeto a las condiciones iniciales : xn = 0.1, yn = 0.1 y zn = 0.1. El valor teórico del
máximo exponente de Lyapunov es λ ≈ 0.225 (Kantz & Schreiber, 2003).

Antes de calcular los exponentes de Lyapunov y el espectro fue necesario iterar

5000 veces el sistema dinámico con ayuda de Matlab v R2017a (MathWorks USA) para
generar un conjunto de 5000 puntos de la forma (x, y, z). En la Figura 4.2 se muestra el
espacio de fase generado por estos 5000 puntos.

Después, se ejecutó la rutina lyap k con los datos obtenidos. La Figura 4.3 mues-
tra la gráfica de los resultados que se obtuvieron. La rutina proporciona el logaritmo del
factor de estiramiento (Hegger et al., 2007). El máximo exponente de Lyapunov es la pen-
diente de las lı́neas rectas aproximadas (lı́nea continua) a las funciones obtenidas (lı́nea
punteada). En este caso la pendiente aproximada es de 0.1917, con un error relativo de
0.148 con respecto al valor teórico.

Finalmente, con la rutina lyap spec y los datos obtenidos de la iteración del
sistema dinámico discreto se calculó el espectro de Lyapunov. Los resultados son 0.2400

59
 Metodologı́a propuesta 4. Planteamiento del problema y metodologı́a propuesta

-2
Log(S)

-4

-6

-8
0 10 20 30 40 50
n

Figura 4.3: La estimación del máximo exponente de Lyapunov se obtiene a partir del valor
de las pendientes de las lı́neas rectas aproximadas (lı́nea continua) al logaritmo del factor
de estiramiento (lı́nea punteada).

y 0.1766. El error relativo del primer exponente es 0.066 con respecto a 0.225 (valor
teórico).

La rutina lyap spec al calcular el espectro de los exponentes de Lyapunov tam-

bién calcula el máximo exponente. En los resultados obtenidos podemos observar que
estas estimaciones no son iguales, lo cual sucede porque las rutinas que los aproximan
emplean diferentes algoritmos, para mayor detalle consultar Hegger et al. (2007).

4.2.5. Entropı́a de Kolmogorov–Sinai

Otra métrica que no es fácil de calcular es la entropı́a de de K–S, de la misma forma

que en la sección anterior se eligió TISEAN (Hegger et al., 2007) para aproximarla. En
esta sección se ejemplifica el cálculo de la entropı́a de K–S, explicada en la sección 2.4.6.
Para ello se utiliza la función seno y la rutina d2 de TISEAN (Hegger et al., 2007).

La rutina d2 utiliza los siguientes parámetros.

60
4. Planteamiento del problema y metodologı́a propuesta Metodologı́a propuesta

1. Los datos a los cuales les calcularemos su entropı́a.

2. La dimensión del espacio de proyecciones de retraso (sección 2.4.7).

3. El retraso considerado en tal espacio.

Lo primero que se realizó fue evaluar la función seno(x) con x ∈ (0, 50) de 0.1 en
0.1 para generar una serie temporal en Matlab vR2017a (MathWorks USA). Lo anterior
se repitió con la función seno(x) + seno(5x) y seno(x) + ruido. El ruido se generó con
la instrucción rand(5001,1) de Matlab vR2017a (MathWorks USA).

Después, siguiendo el teorema de Whitney (ver sección 2.4.7), la dimensión del

espacio de proyecciones se establece igual a 3. Y de acuerdo con Kantz & Schreiber (2003)
el retardo lo podemos elegir como un cuarto del periodo de la función en cuestión. En
este caso, la función seno tiene periodo 2π. Por lo que el retardo serı́a aproximadamente
6/4 = 3/2. Lo hemos tomado como 1.

La rutina d2 genera varios archivos de salida. Estamos interesados en el archivo

con extensión .h2 porque proporciona aproximaciones a la entropı́a de K–S (Hegger
et al., 2007).

La función seno(x) es determinista, teóricamente la entropı́a de tal función debe

ser cero. En la Figura 4.4(a) se muestran las aproximaciones de la entropı́a de K–S para
la función seno(x). De acuerdo con Kantz & Schreiber (2003) la entropı́a de K–S es la
meseta observada en la gráfica. Con la dimensión del espacio de proyecciones igual a 3
la aproximación de la entropı́a es 0.0689, con un error estándar de 0.059 y un intervalo
de confianza al 95 % (0.0592, 0.0782).

Con la función seno(x) + seno(5x) la entropı́a teórica de K–S también deberı́a

ser cero porque se trata de una función determinista. En la gráfica de la Figura 4.4(b)
se muestran las aproximaciones a la misma para las dimensiones del espacio de proyec-
ciones m = 2, 3. Con m = 3, la entropı́a calculada es 0.1636, con un error estándar de
0.1253 y un intervalo de confianza al 95 % (−0.0432, 0.3704).

61
 5
Experimentos y resultados

A continuación, se presentan detalladamente los experimentos y resultados que se

realizaron en la presente investigación. Todos los experimentos fueron realizados utili-
zando lenguaje Fortran, ver apéndice (A) para las especificaciones del hardware y soft-
ware utilizados.

5.1. Primer experimento

El objetivo de este experimento es establecer casos en que el Euler hacia atrás su-
prime el caos.

La hipótesis es que el método Euler hacia atrás suprime caos cuando se realiza una
elección inadecuada del paso de integración.

En este experimento se fijan los parámetros de los sistemas dinámicos continuos

del segundo grupo de prueba (ver sección 4.2.1) con los cuales presentan un comporta-
miento caótico. Tales valores son derivados de estudios teóricos de las ecuaciones que
representan a los sistemas. También, se establece el conjunto de condiciones iniciales y
el intervalo de simulación.

Para utilizar el método de Euler hacia atrás, es necesario proporcionar el tamaño

 Primer experimento 5. Experimentos y resultados

del paso. Basados en los resultados de un análisis teórico hecho por Varsakelis & Anag-
nostidis (2016) se realizó lo siguiente. Dada una cota hc1 para el paso de integración, si el
paso de integración elegido hs1 para simular el sistema caótico satisface que hs1 > hc1 ,
entonces los resultados numéricos suprimen al caos. En caso contrario, si hs1 ≤ hc1 ,
el caos se conserva. Tomando en cuenta la cota hc1 , la forma de seleccionar el hs1 que
suprime o mantiene el caos fue mediante prueba y error.

La métrica de evaluación de este experimento es cualitativa: el espacio de fase.

Los valores de los parámetros, las condiciones iniciales, intervalos de simulación y

cotas hc1 que se utilizaron en este experimento se presentan a continuación.

1. Sistema de Lorenz: Los valores de los parámetros son σ = 10, R = 28, γ = − 83 ,

las condiciones iniciales son (1, 1, 1), el intervalo de simulación es [0, 50] y la cota
es hc1 = 0.0018.

2. Ecuaciones de Rabinovich-Fabrikant: Los valores de los parámetros son α =

1.1, γ = 0.87, las condiciones iniciales son (−1.0824, 0.1141, 0.2873), el intervalo
de simulación es [0, 50] y la cota es hc1 = 0.007.

3. Ecuaciones de Lorenz–Stenflo: Los valores de los parámetros son σ = 1, r = 26,

β = 0.7, s = 1.5, las condiciones iniciales son (1, 1, 1, 1), el intervalo de simulación
es [0, 50] y la cota es hc1 = 0.0503.

En la figura 5.1 se muestra que con un paso de integración h = 0.01 el méto-

do numérico Euler hacia atrás suprime el caos de los resultados numéricos del sistema
caótico de Lorenz. Un caso similar se presenta en la Figura 5.2 con un paso de integración
h = 0.04 para las ecuaciones de Rovinovich-Fabrikant. Finalmente, en la Figura 5.3 se
observa el mismo comportamiento con h = 0.07 para las ecuaciones de Lorenz-Stenflo.

Con los resultados de los experimentos numéricos concluimos que una elección
incorrecta del paso de integración conduce a la supresión de caos en sistemas caóticos,
con lo que la hipótesis del experimento no es rechazada.

66
5. Experimentos y resultados Primer experimento

(a) B–E h = 0.001 (b) B–E h = 0.01

h=0.001
15 h=0.01
10

5
x(t)

-5

-10

-15

0 10 20 30 40 50
t

(c) Solución x(t)

Figura 5.1: Supresión de caos en el sistema caótico de Lorenz utilizando el método

numérico Euler hacia atrás. Las gráficas a) y b) muestran el espacio de fase obtenido
con h = 0.001 y h = 0.01 respectivamente. Mientras que la gráfica c) muestra la solución
x(t) para cada paso de integración. Con h = 0.01 los resultados numéricos presentan un
comportamiento estable: el caos es suprimido.

67
 Primer experimento 5. Experimentos y resultados

(a) B–E h = 0.005 (b) B–E h = 0.04

-0.6 h=0.005
h=0.04
-0.8

-1
x(t)

-1.2

-1.4

-1.6

-1.8
0 10 20 30 40 50
t

(c) Solución x(t)

Figura 5.2: Supresión de caos en las ecuaciones de Rabinovich–Fabrikant con el método

numérico Euler hacia atrás. Las gráficas a) y b) muestran el espacio de fase obtenido con
h = 0.005 y h = 0.04 respectivamente. Mientras que la gráfica c) muestra la solución x(t)
para cada paso de integración. Con h = 0.04 los resultados numéricos suprimen el caos.

68
5. Experimentos y resultados Primer experimento

(a) B–E h = 0.001 (b) B–E h = 0.07

h=0.001
6 h=0.07
4

2
x(t)

-2

-4

-6
0 10 20 30 40 50
t

(c) Solución x(t)

Figura 5.3: Supresión de caos en las ecuaciones de Lorenz–Stenflo con el método numéri-
co Euler hacia atrás. Las gráficas a) y b) muestran el espacio de fase obtenido con
h = 0.001 y h = 0.07 respectivamente. Mientras que la gráfica c) muestra la solución
x(t) para cada paso de integración. Con h = 0.07 los resultados numéricos suprimen el
caos.

69
 Segundo experimento 5. Experimentos y resultados

5.2. Segundo experimento

El objetivo de este experimento es establecer casos en que Euler hacia adelante

genera caos falso en los resultados numéricos de los sistemas dinámicos continuos no
caóticos.

La hipótesis de este experimento es que con una elección inadecuada del paso de
integración se genera caos falso en los resultados numéricos.

Los pasos seguidos para realizar este experimento son análogos a los del primer
experimento. Derivado de estudios teóricos (Varsakelis & Anagnostidis, 2016) se propor-
cionan los parámetros de los sistemas dinámicos continuos de prueba del grupo 2 (ver
sección 4.2.1) que dan pie a un comportamiento no caótico. Se precisan las condiciones
iniciales y el intervalo de simulación.

Los resultados del estudio teórico proporcionan una hc2 , tal que si el paso de inte-
gración elegido hs2 satisface que hs2 > hc2 , entonces los resultados numéricos presentan
un comportamiento caótico. En otro caso (hs2 ≤ hc2 ) el comportamiento no caótico se
mantiene (Varsakelis & Anagnostidis, 2016).

La manera de encontrar el hs2 que genera o no un comportamiento caótico falso

fue mediante prueba y error.

La métrica de evaluación de este experimento es el espacio de fase.

Los valores de los parámetros, de las condiciones iniciales, de los intervalos de

simulación y de las cotas hc2 que se utilizaron son los siguientes.

1. Sistema de Lorenz: Los valores de los parámetros son σ = 10, R = 20, γ = − 83 ,

las condiciones iniciales son (1, 1, 1), el intervalo de simulación es [0, 50] y la cota
es hc2 = 0.0028.

70
5. Experimentos y resultados Segundo experimento

(a) F–E h = 0.001 (b) F–E h = 0.01

Figura 5.4: Caos falso en el sistema no caótico de Lorenz con el método numérico Euler
hacia adelante. Las gráficas a) y b) muestran el espacio de fase obtenido con h = 0.001
y h = 0.01 respectivamente. Con un paso de integración h = 0.01 se obtienen resultados
numéricos que presentan un comportamiento caótico: surge caos falso.

2. Ecuaciones de Rabinovich–Fabrikant: Los valores de los parámetros son α =

0.1, γ = 0.9, las condiciones iniciales son (1, 1, 1), el intervalo de simulación es
[0, 100] y la cota es hc2 = 0.009.

3. Ecuaciones de Lorenz–Stenflo: Los valores de los parámetros son σ = 1, r = 26,

β = 0.7, s = 0.1, las condiciones iniciales son (1, 1, 1, 1), el intervalo de simulación
es [0, 50] y la cota es hc2 = 0.03.

En la figura 5.4 se muestra que con un paso de integración h = 0.01 el méto-

do numérico Euler hacia adelante arroja resultados numéricos que presentan caos falso
para el sistema de Lorenz. Un caso similar se presenta en la Figura 5.5 con un paso de
integración h = 0.05 para las ecuaciones de Rovinovich-Fabrikant. En la Figura 5.6 se
presenta el mismo comportamiento con h = 0.035 para las ecuaciones de Lorenz-Stenflo.

Con los resultados obtenidos en este experimento se concluye que una selección
inadecuada del paso de integración en el método Euler hacia adelante se presenta caos
falso. En consecuencia, la hipótesis planteada para este experimento no es rechazada.

71
 Segundo experimento 5. Experimentos y resultados

(a) F–E h = 0.005 (b) F–E h = 0.05

Figura 5.5: Caos falso en las ecuaciones de Rabinovich–Fabrikant con el método numérico
Euler hacia adelante. Las gráficas a) y b) muestran el espacio de fase obtenido con h =
0.005 y h = 0.05 respectivamente. Con un paso de integración h = 0.05 los resultados
numéricos presentan un comportamiento caótico.

(a) F–E h = 0.01 (b) F–E h = 0.035

Figura 5.6: Caos falso en las ecuaciones de Lorenz–Stenflo con el método numérico Euler
hacia adelante. Las gráficas a) y b) muestran el espacio de fase obtenido con h = 0.01
y h = 0.035 respectivamente. Con un paso de integración h = 0.035 los resultados
numéricos presentan un comportamiento caótico.

72
5. Experimentos y resultados Tercer experimento

5.3. Tercer experimento

El objetivo de este experimento es encontrar a los métodos numéricos, entre los

métodos estándares, sofisticados y especiales, que conserven el caos de los sistemas caóti-
cos de prueba del grupo 1 (ver sección 4.2.1) y que lo mantengan durante grandes inter-
valos de simulación.

La hipótesis de este experimento es que con el método especial de Gautschi se

conserva y mantiene el caos en intervalos de tiempo grandes y con menor número de
evaluaciones que los métodos sofisticados.

El intervalo de simulación se fijó en [0, 50000] debido al contexto de las aplica-

ciones (ver sección 3.1). Los métodos numéricos utilizados en este experimento son los
siguientes.

Métodos estándares Sólo se utilizó el Euler hacia atrás. Para utilizar este método se debe
proporcionar el paso de integración. La estrategia que se siguió para encontrar
un paso h que permitiera integrar todo el intervalo de simulación fue prueba y
error. Comenzando con un h = 0.1 y disminuyéndolo hasta encontrar un h que
permitiera simular el intervalo completo.

Métodos sofisticados Se utilizaron el AB-AM, ODEX y R–K F 45. Los parámetros que se
deben dar para utilizar estos métodos son el error relativo y el absoluto (ver sección
2.5.6). En este trabajo se decidió que los errores absoluto y relativo serı́an iguales
y se trabajó con dos magnitudes de error: 10−2 y 10−8 . Además, como este tipo de
métodos son de paso variable, otro parámetro que se debe fijar es la frecuencia de
salida fs . Ésta es una cota superior para el paso de integración puesto que h ≤ fs .
Para definir el valor de este parámetro se utilizó prueba y error comenzando con
un valor de 10 y disminuyéndolo hasta que el método numérico lograba integrar
el sistema en todo el intervalo de simulación.

Métodos especiales Se utilizó el método de Gautschi (ver sección 2.5.5). En este tipo de

73
 Tercer experimento 5. Experimentos y resultados

Sistema caótico ω aprox.

Sistema de Lorenz 24.0

Atractor de Rossler 1.5
Sistema de Chen 24.0
Sistema de Chua 3.0

Tabla 5.1: Frecuencias aproximadas de los sistemas caóticos seleccionados.

métodos también se debe proporcionar el tamaño de paso de integración. El cual

depende de la frecuencia (ω) del sistema a simular. En consecuencia, lo primero que
se realizó fue aproximar numéricamente las frecuencias de los sistemas caóticos. En
la tabla 5.1 mostramos las frecuencias obtenidas correspondientes a cada sistema
caótico. Contando con esta información, también se utilizó prueba y error para
hallar un paso de integración v = ωh adecuado, iniciando con h = 2−1 = 1/2. Por
ejemplo, si ω = 1.5 entonces comenzamos probando con v = 1.5 × 1/2 = 1.5/2 y
disminuyendo la potencia de 2 hasta obtener simulaciones en todo el intervalo de
prueba. La implementación computacional de este método numérico se encuentra
en el apéndice B.

El primer grupo de sistemas dinámicos continuos (sección 4.2.1) se utilizó para

realizar este experimento, es importante resaltar que estos sistemas presentan un com-
portamiento caótico. Los parámetros de los sistemas caóticos utilizados y las condiciones
iniciales se presentan a continuación.

1. Sistema de Lorenz: Los valores de los parámetros son σ = 10, R = 28, γ = −8/3
y las condiciones iniciales son (−15.8, −17.48, 35.64).

2. Atractor de Rossler: Los valores de los parámetros son a = b = 0.2, c = 5.7 y

las condiciones iniciales son (−9, 0, 0).

74
5. Experimentos y resultados Tercer experimento

3. Sistema de Chen: Los valores de los parámetros son a = 35, b = 3, c = 28 con

condiciones iniciales (−10, 0, 37).

4. Sistema de Chua: Los valores de los parámetros son α = 10, β = 11.5, m1 =

−0.276, m2 = −3.036, m1 = −0.276, b2 = 0.8, b3 = 1.4, d2 = m1 b2 y d3 =
d2 + m2 (b3 − b2 ) y las condiciones iniciales son (0.1, 0, 0).

En la Tabla 5.2 se muestran los pasos de integración encontrados que permitieron

integrar numéricamente los sistemas caóticos de prueba durante el intervalo simulación
completo utilizando el método Euler hacia atrás. En todos los casos el B–E integra el
intervalo completo de simulación.

En la Tabla 5.3 se muestran los tiempos alcanzados por las simulaciones numéri-
cas y la frecuencia de salida utilizada por cada método integrador sofisticado para cada
sistema caótico, tomando un error relativo de 10−2 . De forma análoga, en la Tabla 5.4
se presentan los resultados obtenidos utilizando un error relativo de 10−8 . En todos los
casos sólo el AB–AM integra completamente el intervalo de simulación.

Los pasos de integración que se utilizaron para integrar los sistemas caóticos de
prueba con el método de Gautschi se presentan en la Tabla 5.5. En todos los casos de
prueba el método de Gatuschi logró simular el intervalo completo de simulación.

Cabe aclarar que no se realizaron repeticiones en ninguna de las simulaciones de-

bido a que los resultados son deterministas (establecidos el conjunto de parámetros y las
condiciones de simulación siempre devuelven los mismos resultados).

5.3.1. Evaluación y comparación

Las métricas utilizadas para evaluar los resultados numéricos de este experimento
son las siguientes: el espacio de fase y el periodograma de forma cualitativa; el número de
evaluaciones, los exponentes de Lyapunov y la entropı́a de K–S de manera cuantitativa.

75
 Tercer experimento 5. Experimentos y resultados

Sistema caótico Paso integración h

Sistema de Lorenz 0.002

Atractor de Rossler 0.01
Sistema de Chen 0.01
Sistema de Chua 0.01

Tabla 5.2: Se muestran los pasos de integración obtenidos que permiten simular los siste-
mas caóticos durante el intervalo completo mediante el método B–E.

5.3.2. Espacio de fase

La primera comparación que se llevó a cabo fue cualitativa: se analizó el atractor

extraño de los resultados numéricos mediante las gráficas del espacio de fase. Lo que
se esperaba observar es que los atractores de cada sistema caótico lucieran semejantes
cualitativamente cuando son integrados con distintos métodos numéricos.

En las Figuras 5.7 a 5.10 se presentan los espacios de fase de los sistemas de Lorenz,
Rossler, Chen y Chua respectivamente, obtenidos con los métodos numéricos que logra-
ron simular los 50000s. En cada gráfica el tiempo de simulación se codifica mediante el
color: colores más claros indican resultados más tardı́os.

Para el sistema de Lorenz, las estrategias numéricas muestran un atractor extraño

con forma bastante similar excepto el de la Figura 5.7(b) que fue obtenido utilizando el
método de AB–AM 10−2 porque se observa un ciclo en el atractor de Lorenz.

En el caso del atractor de Rossler, mostrado en la Figura 5.8, los espacios de fase
son parecidos en su forma. Para el sistema de Chen, los espacios de fase mostrados en
las Figuras 5.9(b) y (c) generados con AB–AM 10−2 y B–E respectivamente, presentan
diferencias frente a los de las Figuras 5.9(a) y (d) integrados con AB–AM 10−8 y Gauts-
chi respectivemnte. Estos últimos muestran un atractor extraño más parecido al que se

76
5. Experimentos y resultados Tercer experimento

Sistema Método Frecuencia

Caos Tiempo s
caótico numérico salida

Sistema AB-AM Si 50,000 2

de ODEX No 1,720 10
Lorenz R-K F 45 No 23 10
Atractor AB-AM Si 50,000 10
de ODEX No 5,180 10
Rossler R-K F 45 No 170 10
Sistema AB-AM Si 50,000 10
de ODEX No 1,000 10
Chen R-K F 45 No 18 2
Sistema AB-AM Si 50,000 0.001
de ODEX No 15 0.0001
Chua R-K F 45 No 150 10

Tabla 5.3: Resultados numéricos de los métodos sofisticados utilizando un error relativo
de 10−2 . En la primera columna están los sistemas caóticos de prueba, posteriormente los
métodos sofisticados. En la siguiente columna se reporta si los métodos numéricos logra-
ron completar los 50,000s de simulación, en caso afirmativo se escribe un “Si”, en caso
contrario un “No”. En la cuarta columna se proporciona el máximo tiempo de simulación
alcanzado. Finalmente se reporta la frecuencia de salida utilizada.

77
 Tercer experimento 5. Experimentos y resultados

Sistema Método Frecuencia

Caos Tiempo s
dinámico numérico salida

Sistema AB-AM Si 50,000 2

de ODEX No 732 2
Lorenz R-K F 45 No 22 1
Atractor AB-AM Si 50,000 10
de ODEX No 4,320 10
Rossler R-K F 45 No 20 10
Sistema AB-AM Si 50,000 1
de ODEX No 730 10
Chen R-K F 45 No 2 2
Sistema AB-AM Si 50,000 1
de ODEX No 1,044 0.1
Chua R-K F 45 No 19 0.1

Tabla 5.4: Resultados numéricos de los métodos sofisticados utilizando un error relativo
de 10−8 . La columna “Caos” reporta si el método numérico integró todo el intervalo de
prueba.

Sistema caótico Paso integración h

Sistema de Lorenz 0.333

Atractor de Rossler 0.167
Sistema de Chen 0.25
Sistema de Chua 0.273

Tabla 5.5: Se muestran los pasos de integración obtenidos que permiten simular todo el
intervalo de prueba utilizando el método de Gautschi.

78
5. Experimentos y resultados Tercer experimento

104 104
5 5
40
40
4 4
30 30
3 3
z(t)

z(t)
20 20
2 2
10 10
1 1

20 0 20 0
0 10 0 10
-20 -10 0 -20 -10 0
y(t) x(t) y(t) x(t)

(a) AB–AM, Er = 10−8 (b) AB–AM, Er = 10−2

104 104
5 5
40 40
4 4
30 30
3 3
z(t)

z(t)
20 20
2 2
10 10
1 1

20 0 20 0
0 0 10 0 0 10
-20 -10 -20 -10
y(t) x(t) y(t) x(t)

(c) B–E h = 0.007 (d) Gautschi h = 0.333

Figura 5.7: Espacio de fase del sistema de Lorenz integrado con los métodos numéricos
que lograron simular los 50000s. El tiempo de simulación se codifica mediante el color:
colores más claros indican muestras más tardı́as. Los espacios de fase son parecidos en su
forma excepto el (b), donde se utilizó el método de AB–AM con Er = 10−2 que presenta
un ciclo en el atractor de Lorenz.

encuentra en el artı́culo de Eftekhari & Jafari (2012).

Con el sistema de Chua, observado en la Figura 5.10, las estrategias numéricas

muestran un atractor extraño con forma bastante similar puesto que se presentan tres
enrollamientos. Es importante mencionar que, en el caso del espacio de fase generado
con el AB–AM 10−2 solo se graficaron 500s debido a que el periodo de salida es pequeño
0.001, provocando un archivo de salida bastante grande (alrededor de 4 gigabytes) para
ser procesado por completo por el software que grafica los espacios de fase.

79
 Tercer experimento 5. Experimentos y resultados

104 104
5 5

4 4
20 20
3 3
z(t)

z(t)

10 10
2 2

1 1
0 0
5 10 5 10
0 5 0 5
-5 0 0 -5 0 0
-10 -5 -10 -5
y(t) x(t) y(t) x(t)

(a) AB–AM, Er = 10−8 (b) AB–AM, Er = 10−2

104 104
5 5

4 4
20
20
3 3
z(t)

z(t)

10
2 10 2

1 1
0 0
5 10 5 10
0 5 0 5
-5 0 0 -5 0 0
-10 -5 -10 -5
y(t) x(t) y(t) x(t)

(c) B–E h = 0.01 (d) Gautschi h = 0.167

Figura 5.8: Análogamente a la Figura 5.7 se muestran los espacios de fase del atractor de
Rossler. Los espacios de fase son parecidos en su forma.

80
5. Experimentos y resultados Tercer experimento

104 104
50 5 5
40
40 4 4
30
30
z(t)

z(t)
3 3
20 20
2 2
10 10
1 1
20 20
0 0 0 0
-20 -20 20
-20 -10 0 10 20 -20 -10 0 10
y(t) x(t) y(t) x(t)

(a) AB–AM, Er = 10−8 (b) AB–AM, Er = 10−2

104 104
5 60
30 4 4
40
z(t)

z(t)

20 3 3

2 20 2
10
1 1
20 20
0 0 0 0
-20 10 20 -20 10 20
-20 -10 0 -20 -10 0
y(t) x(t) y(t) x(t)

(c) B–E h = 0.01 (d) Gautschi h = 0.25

Figura 5.9: Análogamente a la Figura 5.7 se muestran los espacios de fase del sistema
de Chen. Los espacios de fase (b) y (c), generados con AB–AM (Er = 10−2 ) y B–E
respectivamente, presentan diferencias frente a los espacios de fase (a) y (d), generados
con AB–AM (Er = 10−8 ) y Gautschi respectivamente. Estos últimos muestran un atractor
extraño más parecido al que se encuentra en el artı́culo de Eftekhari & Jafari (2012).

81
 Tercer experimento 5. Experimentos y resultados

104
5 5

4
z(t)

0 3

-5 1
0.5 0
-0.5 0
-2 0 2
y(t) x(t)

(a) AB–AM, Er = 10−8 (b) ∗AB–AM, Er = 10−2

104 104
5 5 5
4
4 4
2
z(t)

z(t)

0 3 0 3

-2 2 2
-4
1 -5 1
0.5 0.5 0
0 0 0
-0.5 -0.5
-2 0 2 -2 0 2
y(t) x(t) y(t) x(t)

(c) B–E h = 0.01 (d) Gautschi h = 0.273

Figura 5.10: De forma análoga a la Figura 5.7 se muestran los espacios de fase del sistema
de Chua. Las estrategias numéricas muestran un atractor extraño con forma bastante simi-
lar porque presentan tres enrollamientos. Para el espacio de fase (d) se grafican solo 500s
debido a que el periodo de salida es pequeño 0.001, y el archivo resultante es bastante
grande para ser procesado por el software de graficación.

82
5. Experimentos y resultados Tercer experimento

5.3.3. Periodograma

Tal como se menciona en la sección del marco teórico 2.4.1, una condición necesa-
ria pero no suficiente para la presencia de caos es que el espectro de frecuencias de los
sistemas caóticos sea parecido al ruido y que presente picos en las frecuencias dominan-
tes.

En la Tabla 5.6 grafican los periodogramas correspondientes a la variable de estado

x de cada sistema caótico integrado numéricamente con los distinto métodos. Se mues-
tra sólo una variable porque las demás son redundantes. Cada uno de los periodogramas
muestran un comportamiento parecido al ruido y bastantes picos en las frecuencias do-
minantes. En el caso del sistema de Chua estamos mostrando el periodograma de los
primero 500s por lo cual luce diferente al resto porque el archivo completo no pudo ser
procesado por Matlab vR2017 (MathWorks USA) para calcular el periodograma.

El periodograma del sistema de Lorenz integrado con AB–AM 10−2 presenta un

comportamiento diferente a los resueltos con las otras estrategias numéricas porque pre-
senta dos picos altos (cercanos al 60) en las frecuencias 0 y 0.05, y los picos restantes solo
alcanzan el valor 20. Mientras que los picos en los periodogramas restantes del sistema
de Lorenz sólo alcanzan valores cercanos al 30. Tal fenómeno está en sintonı́a con el
hecho de que el espacio de fase de los resultados del sistema de Lorenz integrado con
AB–AM 10−2 presenta un comportamiento distinto a los restantes (para mayor detalle
ver sección 5.3.2).

5.3.4. Número de evaluaciones

En la Tabla 5.7 se reporta el número de evaluaciones realizada por cada método in-
tegrador a cada sistema caótico como una medida del costo computacional. En todos los
casos los métodos más costosos (en cuanto al número de evaluaciones) son los sofistica-
dos y estándares. Por ejemplo, en el atractor de Rossler el método B–E resulta ser el más

83
 Tercer experimento 5. Experimentos y resultados

costoso, contradiciendo al hecho de que los métodos tipo Euler son los menos costosos.
La estrategia numérica menos costosa en todos los sistemas de prueba fue la de Gautschi,
acorde a nuestra hipótesis en este experimento que plantea que los métodos especiales
debieran ser ventajosos para estos casos.

5.3.5. Exponentes de Lyapunov

El máximo exponente de Lyapunov y el conjunto de todos los exponentes de Lya-

punov (espectro) se calculó para cada una de las simulaciones numéricas. En la Tabla 5.8
se exhibe la estimación al máximo exponente de Lyapunov para cada sistema caótico in-
tegrado con los distintos métodos numéricos. En cada caso, se reporta el valor de interés,
el error estándar y el intervalo de confianza al 95 %.

84
 5. Experimentos y resultados
AB–AM 10−2 AB–AM 10−8 B–E Gautschi

60 40 40 40

40 20
20 20
FFT(x)2

FFT(x)2

FFT(x)2

FFT(x)2
Sistema de 20
0
0 0 0
-20

Lorenz -20

0 0.05 0.1 0.15 0.2

-40
0 0.05 0.1 0.15 0.2
-20
0 0.02 0.04 0.06 0.08
-20
0 0.02 0.04 0.06 0.08 0.1
Frequency (Hz) Frequency (Hz) Frequency (Hz) Frequency (Hz)

60 60 60 60

40 40 40 40
FFT(x)2

FFT(x)2

FFT(x)2

FFT(x)2
Atractor de 20 20 20 20

0 0 0 0

Rossler -20
0 0.01 0.02 0.03 0.04
-20
0 0.01 0.02 0.03 0.04
-20
0 0.01 0.02 0.03 0.04
-20
0 0.01 0.02 0.03 0.04 0.05
Frequency (Hz) Frequency (Hz) Frequency (Hz) Frequency (Hz)
85

40 40 40 40

30
20
20 20
FFT(x)2

FFT(x)2

FFT(x)2

FFT(x)2
Sistema de 20

10
0

0 0
-20
0

Chen -10
0 0.01 0.02 0.03 0.04
-40
0 0.1 0.2 0.3 0.4
-20
0 0.02 0.04 0.06 0.08 0.1 0.12
-20
0 0.05 0.1 0.15
Frequency (Hz) Frequency (Hz) Frequency (Hz) Frequency (Hz)

40 40 40 40

20

Tercer experimento
20
20 20
0
FFT(x)2

FFT(x)2

FFT(x)2
Sistema de 0
FFT(x)2

-20
0 0
-40 -20
-60 -20 -20
-40

Chua -80

-100
0 50 100 150 200 250 300 350 400 450
-60
0 0.1 0.2 0.3 0.4
-40
0 0.02 0.04 0.06 0.08
-40
0 0.02 0.04 0.06 0.08
Frequency (Hz)
Frequency (Hz) Frequency (Hz) Frequency (Hz)

Tabla 5.6: Gráficas de los periodogramas, cada una muestra un comportamiento parecido al ruido y picos en las frecuencias
dominantes lo que sugiere un comportamiento caótico. En el caso del sistema de Chua con AB– AM Er = 10−2 se muestra el
periodograma de los primeros 500s por lo cual luce diferente al resto.
 Tercer experimento 5. Experimentos y resultados

La estimación del máximo exponente de Lyapunov del sistema de Lorenz integra-

do con el método AB–AM 10−2 es distinto a los reportados para el mismo sistema con
los métodos numéricos restantes, puesto que tal valor es 0.2987 mientras que los otros
se encuentran entre 0.6 y 0.8. Lo que reafirma que se presentan un espacio de fase y
periodograma diferentes al resto (ver secciones 5.3.2 y 5.3.3).

El sistema de Chua integrado con la estrategia AB–AM 10−2 sólo toma en cuenta
los primeros 500s de simulación numérica porque procesar los resultados completos es
una tarea que la rutina lyap k no logró finalizar debido al gran tamaño del archivo.

En cuanto al espectro de los exponentes de Lyapunov recordemos que la condición

P
necesaria para que se presente caos es que λ1 > 0 y λi < 0 (ver sección 2.4.1). En la
Tabla 5.9, se reportan las aproximaciones obtenidas a través de la rutina lyap spec.
Todos los resultados numéricos de las simulaciones cumplen con la condición de positi-
vidad, excepto en el caso del sistema de Chua integrado con el método AB–AM 10−2 . Lo
cual puede deberse a que la rutina lyap spec no logra procesar los resultados com-
pletamente.

Las aproximaciones al máximo exponente de Lyapunov (λ1 ) de las Tablas 5.8 y

5.9 no son iguales, porque las rutinas lyap k y lyap spec que los calculan utilizan
algoritmos distintos para estimarlos. Para mayor detalle se refiere al lector a consultar
Hegger et al. (2007).

5.3.6. Entropı́a de Kolmogorov–Sinaı́

Por último, para cada uno de los resultados numéricos se estimó la entropı́a de K–S
(ver sección 4.2.5). En la tabla 5.10 se reporta la estimación de la entropı́a de K–S para cada
sistema caótico integrado con los distintos métodos numéricos. Todas las estimaciones
de la entropı́a son positivas, por lo que satisfacen la condición de suficiencia para afirmar
que todas las soluciones numéricas son caóticas (ver sección 2.4.6). Cabe resaltar que los

86
5. Experimentos y resultados Tercer experimento

Sistema Método
Evaluaciones
dinámico numérico

AB-AM 10−2 2,113,641

Sistema
AB-AM 10−8 10,944,766
de
B–E 25,000,000
Lorenz
Gautschi 300,000

AB-AM 10−2 626,576

Atractor
AB-AM 10−8 2,061,890
de
B–E 5,000,000
Rossler
Gautschi 450,000

AB-AM 10−2 4,059,489

Sistema
AB-AM 10−8 17,869,913
de
B–E 5,000,000
Chen
Gautschi 400,000

AB-AM 10−2 10,000,048

Sistema
AB-AM 10−8 4,429,370
de
B–E 5,000,000
Chua
Gautschi 183,334

Tabla 5.7: Comparación del número de evaluaciones que realiza el método integrador al
modelo de los sistemas caóticos. En cada sistema caótico resaltamos con azul el resultado
del método numérico con menor número de evaluaciones. En todos los casos resulta ser
el método de Gautschi.

87
 Tercer experimento 5. Experimentos y resultados

Sistema Método Máx. exp. Error Intervalo de confianza

caótico numérico Lyapunov estándar (95 %)
AB-AM 10−2 0.2987 0.0256 (0.2564, 0.3409)
Sistema
AB-AM 10−8 0.7998 0.0614 (0.6985, 0.9011)
de
B–E 0.6885 0.0428 (0.6179, 0.7591)
Lorenz
Gautschi 0.6324 0.0322 (0.5791, 0.6856)
AB-AM 10−2 0.3396 0.0440 (0.2671, 0.4121)
Atractor
AB-AM 10−8 0.2593 0.0582 (0.1632, 0.3554)
de
B–E 0.3047 0.0999 (0.1399, 0.4696)
Rossler
Gautschi 0.3597 0.0822 (0.2733, 0.4460)
AB-AM 10−2 0.6146 0.0216 (0.5789, 0.6502)
Sistema
AB-AM 10−8 0.6158 0.0761 (0.4902, 0.7414)
de
B–E 0.5975 0.0483 (0.5177, 0.6772)
Chen
Gautschi 0.5916 0.0559 (0.4995, 0.6838)
AB-AM 10−2 0.1934∗ 0.0666 (0.0835, 0.3033)
Sistema
AB-AM 10−8 0.7081 0.0646 (0.6014, 0.8148)
de
B–E 1.0766 0.01275 (0.8662, 1.2871)
Chua
Gautschi 1.1406 0.1097 (0.9596, 1.3215)

Tabla 5.8: Estimación del máximo de Lyapunov para cada sistema caótico integrado con
los distintos métodos numéricos. En cada caso, reportamos el valor de interés, el error
estándar y el intervalo de confianza al 95 %. El asterisco (∗) denota que sólo tomamos
en cuenta los primeros 500s de simulación numérica porque la rutina lyap k no logra
procesar los resultados completamente.

88
5. Experimentos y resultados Tercer experimento

Sistema Método
λ1 λ2 λ3
caótico numérico
AB-AM 10−2 0.2002 −0.4758 −1.4900
Sistema
AB-AM 10−8 0.9723 −0.0778 −1.3427
de
B–E 0.7123 −0.6061 −1.8019
Lorenz
Gautschi 0.7604 −0.5886 −1.7689
AB-AM 10−2 0.1630 −0.0860 −2.4630
Atractor
AB-AM 10−8 0.2295 −0.0946 −1.2476
de
B–E 0.2075 −0.0796 −1.3058
Rossler
Gautschi 0.2583 −0.0891 −1.9116
AB-AM 10−2 0.3837 −0.1481 −1.0297
Sistema
AB-AM 10−8 0.6523 0.0986 −0.8156
de
B–E 0.4669 −0.3203 −1.3576
Chen
Gautschi 0.6473 0.0241 −8.7107
AB-AM 10−2 0.2333∗ 0.03847 −0.2135
Sistema
AB-AM 10−8 0.6666 0.0719 −4.1726
de
B–E 1.7021 0.2252 −2.2349
Chua
Gautschi 1.1314 0.1261 −1.5135

Tabla 5.9: Estimación del espectro de los exponentes de Lyapunov. Casi todos los resulta-
P
dos numéricos cumplen con la condición: λ1 > 0 y λi < 0, excepto el sistema de Chua
(∗). Lo cual puede deberse a que solo se toma en cuenta los primeros 500s de simulación
numérica porque la rutina lyap spec no logra procesar los resultados completamente.

89
 Resumen de capı́tulo 5. Experimentos y resultados

resultados obtenidos con el método de Gautschi presentan el mayor valor de la entropı́a

K–S en todos los casos de estudio.

De forma análoga a las métricas anteriores, para el sistema de Chua sólo se toma-
ron en cuenta los primeros 500s de simulación numérica porque la rutina d2 no logra
procesar los resultados completamente.

5.3.7. Conclusiones del tercer experimento

De los resultados obtenidos en este experimento se concluye que el objetivo plan-

teado se cumplió, puesto que los métodos estándares (Euler hacia atrás), sofisticados
(AB-AM con error relativo 10−8 ) y especiales (Gautschi), conservaron el caos y lo man-
tuvieron durante grandes intervalos de simulación para los sistemas caóticos de prueba.
Sin embargo, se debe tener cuidado al momento de utilizar el Euler hacia atrás porque
puede suprimir el caos. Cabe resaltar que el método Gautschi conservó y mantuvo el
caos durante grandes intervalos de tiempo con el menor número de evaluaciones com-
parado con los otros métodos (estándares y sofisticados). Además, proporciona la mayor
entropı́a K–S en todos los casos de estudio. En consecuencia, la hipótesis planteada para
este experimento no es rechazada.

5.4. Resumen de capı́tulo

En el este capı́tulo se presentaron los tres experimentos realizados siguiendo la

metodologı́a propuesta en el capı́tulo 4. Con el primer experimento se obtuvieron casos
en los cuales el método Euler hacia atrás suprime el caos. Mientras que con el segun-
do experimento se presentaron casos en los que el método Euler hacia adelante genera
caos falso. Aunque se mostraron casos particulares, el estudio se puede extender a otros
métodos numéricos de paso fijo.

90
5. Experimentos y resultados Resumen de capı́tulo

Sistema Método Entropı́a Error Intervalo de confianza

caótico numérico K–S estándar (95 %)

AB-AM 10−2 0.1652 0.0241 (0.1652, 0.2051)

Sistema
AB-AM 10−8 0.7081 0.0289 (0.6603, 0.7558)
de
B–E 0.9910 0.0372 (0.9296, 1.0523)
Lorenz
Gautschi 1.3124 0.0528 (1.2252, 1.3996)

AB-AM 10−2 0.3477 0.0250 (0.3065, 0.3889)

Atractor
AB-AM 10−8 0.2965 0.0235 (0.2577, 0.3353)
de
B–E 0.5653 0.0337 (0.5096, 0.6210)
Rossler
Gautschi 0.5895 0.0334 (0.5344, 0.6446)

AB-AM 10−2 1.1553 0.0534 (1.0671, 1.2435)

Sistema
AB-AM 10−8 0.6580 0.0298 (0.6089, 0.7072)
de
B–E 1.2529 0.0553 (1.1616, 1.3442)
Chen
Gautschi 1.2904 0.0548 (1.2001, 1.3808)

AB-AM 10−2 0.3562∗ 0.0293 (0.3080, 0.4045)

Sistema
AB-AM 10−8 0.9055 0.0361 (0.8458, 0.9651)
de
B–E 1.0095 0.0429 (0.9388, 1.0802)
Chua
Gautschi 1.2344 0.0532 (1.1467, 1.322)

Tabla 5.10: Estimación de la entropı́a de K–S para cada sistema caótico integrado con
los distintos métodos numéricos. En cada caso, reportamos el valor de interés, el error
estándar y el intervalo de confianza al 95 %. Todas las estimaciones satisfacen la condición
de positividad con lo que se afirma que todas las soluciones numéricas son caóticas (ver
sección 2.4.6). El asterisco (∗) denota que solo se tomaron en cuenta los primeros 500s
de simulación numérica porque la rutina utilizada para estimarlos no logra procesarlos
completamente.

91
 Resumen de capı́tulo 5. Experimentos y resultados

Finalmente, en el tercer experimento se encontró que los métodos numéricos: Euler

hacia atrás, AB-AM con error 10−8 y Gautschi conservan el caos y lo mantienen durante
el intervalo de simulación [0, 50′ 000]. Además, con los resultados del costo computacio-
nal se afirma que el método de Gautschi es menos costoso que los métodos alternativos
porque realiza un menor número de evaluaciones al modelo que representa a los siste-
mas caóticos. También, es el que presenta mayor entropı́a K–S. Es importante notar que
el método de Gautschi es de paso fijo, y al contar con la frecuencia o una subestima-
ción de la misma, encontrar el paso adecuado de integración es una tarea que no ofrece
dificultad.

92
 6
Conclusiones y trabajo futuro

El problema planteado al principio de este trabajo es el siguiente. Al tener un sis-

tema caótico, sus condiciones iniciales y un intervalo de tiempo durante el cual preten-
demos obtener su solución, esperamos contar con el o los métodos numéricos que nos
permitan conservar el caos durante el intervalo de simulación completo.

Los sistemas caóticos oscilatorios cuentan con caracterı́sticas especiales: tienen alta
sensibilidad a las condiciones iniciales, son no lineales y son oscilatorios sin ser periódi-
cos, entre otras más. Además, la mayorı́a de ellos no cuentan con soluciones analı́ticas,
lo cual los vuelve difı́ciles de analizar.

Considerando que los sistemas caóticos oscilatorios presentan oscilaciones se eligió

al método de Gautschi para solucionar tal problema. Puesto que el método de Gautschi
fue pensado para resolver problemas oscilatorios. Con el método numérico elegido lo
que hacı́a falta era una forma determinar que tal método conservaba el caos durante
grandes intervalos de simulación del sistema caótico, pero no la encontramos al revisar
la literatura. En consecuencia, nos dimos a la tarea de realizar una metodologı́a para tal
fin.

La metodologı́a presentada consta de tres etapas: selección, simulación y evalua-

ción. En la primera etapa seleccionamos los sistemas caóticos, los métodos numéricos
de integración y las métricas de evaluación. En la segunda etapa aproximamos numéri-
camente los sistemas caóticos con los métodos numéricos seleccionados. En la tercera
 Hallazgos 6. Conclusiones y trabajo futuro

etapa evaluamos mediante las métricas los resultados numéricos.

Los sistemas dinámicos elegidos se dividieron en dos grupos. El primero fue con-
formado por sistemas dinámicos caóticos no lineales, deterministas y conservativos, que
son representativos de las áreas de estudio en que son empleados. Y el segundo grupo
consta de sistemas que dependiendo del valor que se asigne a sus parámetros presentan
o no un comportamiento caótico.

Entre los métodos numéricos que seleccionamos además del Gautschi se encontra-
ban: un Euler hacia atrás, un Euler hacia adelante, un Runge-Kutta Felberg, un ODEX y
un Adams Bashforth–Adams Moulton.

Las métricas fueron elegidas por su importancia teórica. Aunque la entropı́a de

Kolmogorov–Sinaı́ proporciona una condición suficiente para la presencia de caos se
emplearon otras métricas, como los exponentes de Lyapunov o el espacio de fase que son
condiciones necesarias para la presencia de caos, porque tenı́amos solo una aproximación
al valor teórico de la entropı́a, la cual podrı́a estar alejada del verdadero valor.

6.1. Hallazgos

Con los resultados de los primeros dos experimentos mostrados en el capı́tulo 5, se

cumplieron los objetivos particulares 1 y 2 de esta investigación, planteados en el capı́tulo
1, porque se lograron establecer casos en los que el método Euler hacia atrás suprime el
caos y el Euler hacia adelante genera caos falso. Con este tipo de métodos se debe cuidar
la selección del paso de integración.

Mediante los resultados del tercer experimento se alcanzaron los objetivos parti-
culares 3 y 4 de esta investigación, puesto que se conservó el caos y se mantuvo durante
grandes intervalos de simulación para los sistemas caóticos de prueba con los métodos:
Euler hacia atrás, AB–AM con un error relativo de 10−8 y Gautschi. Los resultados del

94
6. Conclusiones y trabajo futuro Trabajo futuro

tercer experimento también sugieren que el método de Gautschi es una opción recomen-
dable para simular sistemas caóticos puesto que es un método que no es costoso respecto
al número de evaluaciones realizadas al modelo comparado con las estrategias numéri-
cas B–E o AB–AM. Además, los espacios de fase y periodogramas que se obtuvieron son
similares a los generados con los métodos numéricos alternativos, los exponentes de Lya-
punov son positivos y estadı́sticamente similares y la entropı́a K–S presenta los valores
más altos en todos los casos de estudio. Con base en la evidencia empı́rica, la hipótesis
de esta investigación planteada en el capı́tulo 1 no es rechazada.

Además, tal como se explicó en la sección 2.5.5 debido a la implementación rea-

lizada del método de Gautschi, éste realiza sólo una evaluación después de la primera
iteración. En consecuencia, en cuanto al número de evaluaciones el método de Gautschi
es equiparable al método de Euler. Otra ventaja remarcable del método de Gautschi es
que cuando se tiene la frecuencia o una subestimación de la misma, seleccionar el paso
de integración adecuado es una tarea que no ofrece dificultad.

6.2. Trabajo futuro

Es necesario contar con más estudios teóricos que proporcionen cotas o valores
para los cuales el caos no se suprima al utilizar métodos de paso fijo, por ejemplo para
los Runge–Kutta o incluso para el método de Gautschi.

También, serı́a interesante extender la presente investigación a sistemas no caóti-

cos porque debido a la limitación del tiempo disponible no logramos realizar tal tarea.

Otra dirección posible para el trabajo futuro es evaluar el costo-beneficio de otros

métodos numéricos cuando son empleados para solucionar sistemas caóticos. Una pro-
puesta es el método de ajuste exponencial, el cual da mejores resultados al momento de
solucionar problemas oscilatorios que los métodos basados en polinomios trigonométri-
cos (Berghe & Daele, 2009).

95
 Trabajo futuro 6. Conclusiones y trabajo futuro

Una dirección más serı́a utilizar los resultados numéricos que se obtuvieron para
aplicaciones de criptografı́a.

96
 A
Condiciones de simulación

A continuación presentamos las caracterı́sticas de hardware y software que utiliza-

mos para llevar a cabo los experimentos en este trabajo.

Las pruebas las realizamos en una computadora de 1.8 GHz (2 procesadores) con
32.0 GB de RAM con un procesador Intel(R) Xeon(R) CPU E5-2603 y sistema de 64 bits.

El software que utilizamos fue Visual Studio 2015 Shell (Integrated) versión 14.0.2310
7.0 D14REL. El compilador fue PGI Visual Fortran versión 17.5. Además, la computado-
ra cuenta con una tarjeta NVIDIA y las banderas correspondientes fueron activadas en
Visual Studio al realizar las simulaciones numéricas.
 B
Algoritmo polinomio
trigonométrico Gautschi

En este apartado se encuentra el código utilizado para implementar el algoritmo

del polinomio trigonométrico de Gautschi, el cual es explı́cito, con orden algebraico 2 y
orden trigonométrico 1:

yn+2 = yn+1 + h(β1 fn+1 + β0 fn ) (B.1)

subroutine G a u t s h i O r d e n 2 ( dsedo , n u e c d f , b a n d e r a , b e t a 0 ,
b e t a 1 , y1 , dy1 , sy , sdy , t , w, y , h , dy )

i m p l i c i t none

external dsedo ! m o d e l o c o n l a s edo ’ s

! T i p o s de l a s v a r i a b l e
i n t e g e r ( kind = 4 ) : : i , b a n d e r a , n u e c d f
! bandera = 1 primera llamada
! n u e c d f = numero d e e d o s
 B. Algoritmo polinomio trigonométrico Gautschi

r e a l ( kind = 8 ) : : b e t a 1 , b e t a 0 , h , t , w, nu , nu2 , nu4

r e a l ( kind = 8 ) : : y ( n u e c d f ) , dy ( n u e c d f ) , y1 ( n u e c d f ) ,
dy1 ( n u e c d f )
! y ’ s v a r i a b l e s de e s t a d o
! dy ’ s d e r i v a d a s

r e a l ( kind = 8 ) : : sy ( n u e c d f ) , sdy ( n u e c d f )
! v a r i a b l e s a u x i l i a r e s de almacenamiento

i f ( b a n d e r a == 1 ) then
! Primera llamada al a l g o r i t m o
! Calculamos c o e f i c i e n t e s del algoritmo
nu = w ∗ h
nu2 = nu ∗ ∗ 2
nu4 = nu ∗ ∗ 4
! E s c a l a m o s b e t a 1 y b e t a 0 p o r un f a c t o r d e 1 / 1 2 0
beta1 = ( 1 . 0 d0 / 8 0 . 0 d0 ) ∗ ( 1 2 0 . 0 d0 − 3 0 . 0 d0 ∗ nu2 + nu4 )
beta0 = − ( 1 . 0 d0 / 2 4 0 . 0 d0 ) ∗ ( 1 2 0 . 0 d0 + 1 0 . 0 d0 ∗ nu2 + nu4 )

! Actualizamos la derivada fn al tiempo t = tn

c a l l d s e d o ( n u e c d f , t , y , dy )

! E s t i m a n o s yn +1 p o r m e di o d e un RK2 : y = yn + 1 , f n = dy
c a l l d r t r a p ( dsedo , n u e c d f , y1 , dy1 , t , y , h , dy )

! A c t u a l i z a m o s e l t i e m p o : t = t n +1

100
B. Algoritmo polinomio trigonométrico Gautschi

t = t + h
! S a l v a m o s l o s v a l o r e s d e yn +1 y f n
do i = 1 , n u e c d f
sy ( i ) = y ( i ) ! yn +1
sdy ( i ) = dy ( i ) ! fn
enddo

bandera = 2
else
! Llamadas p o s t e r i o r e s a l metodo
! S e c c i o n p a r a e s t i m a r yn +2
! C a l c u l a m o s d e r i v a d a f n + 1 , t = t n +1
c a l l d s e d o ( n u e c d f , t , y , dy ) ! dy = f n +1

! A l g o r i t m o de G a u t s h i e x p l i c i t o ,
! orden a l g e b r a i c o 2
do i = 1 , n u e c d f
! yn +2 = yn +1 + h ( b1 f n +1 + bo f n )
y ( i ) = sy ( i ) + h ∗ ( b e t a 1 ∗ dy ( i ) + b e t a 0 ∗ sdy ( i ) )
enddo

! A c t u a l i z a m o s t i e m p o t = t n +2
t = t + h
! S a l v a m o s f n +1 y yn +1
do i = 1 , n u e c d f
sy ( i ) = y ( i )
sdy ( i ) = dy ( i )
enddo

101
 B. Algoritmo polinomio trigonométrico Gautschi

endif

return
end subroutine G a u t s h i O r d e n 2

102
 Bibliografı́a

ACM (2005). Collected algorithms. [Online]. Last accessed 2018-06-06. URL:

http://calgo.acm.org/.

Atkinson, K. E. (2017). Numerical analysis. Encyclopædia Britannica, (pp. 1–8). Last

accessed 2018-08-20. URL: https://www.britannica.com/science/numerical-analysis.

Bar-Yoseph, P. Z., Fisher, D., & Gottlieb, O. (1996). Spectral element methods for non-
linear temporal dynamical systems. Computational Mechanics, 18(4), 302–313. DOI:
10.1007/BF00364145.

Bardhi, A., Cipo, P., & Braneshi, M. (2010). Study of a ferroresonant circuit using analytic
harmonic balance, numerical integration of nonlinear ODE and experimental methods.
Proceedings of EPE-PEMC 2010 - 14th International Power Electronics and Motion Control
Conference, (pp. 12–19). DOI: 10.1109/EPEPEMC.2010.5606848.

Behnia, S., Akhshani, A., Mahmodi, H., & Akhavan, A. (2008). Chaotic Cryptographic
Scheme Based on Composition Maps. International Journal of Bifurcation and Chaos,
18(01), 251–261. DOI: 10.1142/S0218127408020288.

Belazi, A., Abd El-Latif, A. A., & Belghith, S. (2016). A novel image encryption scheme
based on substitution-permutation network and chaos. Signal Processing, 128, 155–170.
DOI: 10.1016/j.sigpro.2016.03.021.

Berghe, G. V. & Daele, M. V. (2009). Trigonometric polynomial or exponential fitting

approach? Journal of Computational and Applied Mathematics, 233(4), 969–979. DOI:
10.1016/j.cam.2009.08.103.

Blank, M. (1994). Pathologies generated by round-off in dynamical systems. Physica D:

Nonlinear Phenomena, 78(1-2), 93–114. DOI: 10.1016/0167-2789(94)00103-0.
 BIBLIOGRAFÍA BIBLIOGRAFÍA

Boyce, W. E., DiPrima, R. C., & Haines, C. W. (1969). Elementary differential equations
and boundary value problems, volume 9. Wiley New York. ISBN: 978-0-470-45831-0.

Bradley, L. (2009). Chaos and fractals. [Online]. Last accessed 2018-06-06. URL:
http://www. stsci. edu/˜ lbradley/seminar/.

Chu, F. & Zhang, Z. (1998). Bifurcation and chaos in a rub-impact Jeffcott rotor system.
Journal of Sound and Vibration, 210(1), 1–18. DOI: 10.1006/jsvi.1997.1283.

Corless, R. (1994). What good are numerical simulations of chaotic dynamical systems?
Computers & Mathematics with Applications, 28(10-12), 107–121. DOI: 10.1016/0898-
1221(94)00188-X.

Corless, R. M., Essex, C., & Nerenberg, M. A. (1991). Numerical methods can suppress
chaos. Physics Letters A, 157(1), 27–36. DOI: 10.1016/0375-9601(91)90404-V.

da Silva, F. A. O. & Zhao, L. (2016). A Network of Neural Oscillators for Fractal Pattern
Recognition. Neural Processing Letters, 44(1), 149–159. DOI: 10.1007/s11063-015-9473-
y.

De Markus, A. S. (2001). Detection of the onset of numerical chaotic instabilities by

Lyapunov exponents. Discrete Dynamics in Nature and Society, 6, 121–128. DOI:
doi:10.1155/S1026022601000127.

Deutsch, J. M. (2018). Eigenstate Thermalization Hypothesis. DOI: 10.1088/1361-

6633/aac9f1.

Drazin, P. G. (1992). Nonlinear systems, volume 10. Cambridge University Press. ISBN:
9781139172455.

Eftekhari, S. A. & Jafari, A. A. (2012). Numerical simulation of chaotic dynamical systems

by the method of differential quadrature. Scientia Iranica, 19(5), 1299–1315. DOI:
10.1016/j.scient.2012.08.003.

Evstigneev, N. M., Magnitskii, N. A., & Silaev, D. A. (2015). Qualitative Analysis of Dy-
namics in Kolmogorov’s Problem on a Flow of a Viscous Incompressible Fluid. Diffe-
rential Equations, 51(10), 1292–1305. DOI: 10.1134/S0012266115100055.

104
BIBLIOGRAFÍA BIBLIOGRAFÍA

Franco, J. M. (2004). Runge-Kutta methods adapted to the numerical integration

of oscillatory problems. Applied Numerical Mathematics, 50(3-4), 427–443. DOI:
10.1016/j.apnum.2004.01.005.

Fryska, S. T. & Zohdy, M. A. (1992). Computer dynamics and shadowing of chaotic orbits.
Physics Letters A, 166(5-6), 340–346. DOI: 10.1016/0375-9601(92)90719-3.

Garcı́a-Alonso, F. & Reyes, J. A. (2009). A new method for exact integration of some per-
turbed stiff linear systems of oscillatory type. Applied Mathematics and Computation,
215(7), 2649–2662. DOI: 10.1016/j.amc.2009.09.005.

Garcı́a-Martı́nez, M., Ontañón-Garcı́a, L. J., Campos-Cantón, E., & Čelikovský, S. (2015).

Hyperchaotic encryption based on multi-scroll piecewise linear systems. Applied Mat-
hematics and Computation, 270, 413–424. DOI: 10.1016/j.amc.2015.08.037.

Gautschi, W. (1961). Numerical integration of ordinary differential equations based on

trigonometric polynomials. Numer.\ Math., 3, 381–397.

Góra, P. & Boyarsky, A. (1988). Why computers like Lebesgue measure. Computers and
Mathematics with Applications, 16(4), 321–329. DOI: 10.1016/0898-1221(88)90148-4.

Grassberger, P. & Procaccia, I. (1983). Estimation of the Kolmogorov entropy from a

chaotic signal. Phys. Rev. A, 28(4), 2591–2593. DOI: 10.1103/PhysRevA.28.2591.

Grebogi, C., Hammel, S. M., Yorke, J. A., & Sauer, T. (1990). Shadowing of physical tra-
jectories in chaotic dynamics: Containment and refinement. Phys. Rev. Lett., 65, 1527–
1530. DOI: 10.1103/PhysRevLett.65.1527.

Grote, K. & Meyer-Spasche, R. (1998). Euler-like discrete models of the logistic differen-
tial equation. Computers & Mathematics with Applications, 36(10-12), 211–225. DOI:
10.1016/S0898-1221(98)80022-9.

Hairer, E., Nørsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I
(2Nd Revised. Ed.): Nonstiff Problems. Berlin, Heidelberg: Springer-Verlag. ISBN: 978-
3-540-56670-0.

105
 BIBLIOGRAFÍA BIBLIOGRAFÍA

Harsha, S., Sandeep, K., & Prakash, R. (2003). Quasi-periodic, Subharmonic and Chaotic
Motions of a Rotor Bearing System. International Journal of Nonlinear Sciences and
Numerical Simulation, 4(4), 361–372. DOI: 10.1515/IJNSNS.2003.4.4.361.

Hasler, M., Mazzini, G., Ogorzalek, M., Rovatti, R., & Setti, G. (2002). Scanning the
special issue - special issue on applications of nonlinear dynamics to electronic and
information engineering. Proceedings of the IEEE, 90(5), 631–640. DOI: 10.1109/J-
PROC.2002.1014999.

Hegger, R., Kantz, H., Schreiber, T., & Olbrich, E. (2007). Tisean 3.0.1: Nonlinear time
series analysis. [Online]. Last accessed 2018-06-06. URL: https://www.pks.mpg.de/˜ti-
sean/Tisean 3.0.1/.

Hirai, K. & Adachi, T. (1994). Chaos and bifurcation in numerical computation by the
Runge-Kutta method. International Journal of Systems Science, 25(11), 1695–1706. DOI:
10.1080/00207729408949307.

Hunt, B. R. & Ott, E. (2015). Defining chaos. Chaos: An Interdisciplinary Journal of Non-
linear Science, 25(9), 097618. DOI: 10.1063/1.4922973.

Kantz, H. & Schreiber, T. (2003). Nonlinear Time Series Analysis. Cambridge University
Press, 2 edition. DOI: 10.1017/CBO9780511755798.

Koh, C. & Liaw, C. (1991). Effects of time step size on the response of a bilinear system,
I: Numerical study. Journal of Sound and Vibration, 144(1), 17–29. DOI: 10.1016/0022-
460X(91)90729-4.

Kontorovich, V., Beltrán, L., Aguilar, J., Lovtchikova, Z., & Tinsley, K. R. (2009). Cumulant
analysis of rössler attractor and its applications. Open Cybernetics & Systemics Journal,
3, 29–39. DOI: 10.2174/1874110X00903010029.

Lambert, J. D. (1973). Computational Methods in Ordinary Differential Equations. John

Wiley and Sons. ISBN: 0-471-51194-3.

Li, T.-Y. & Yorke, J. A. (1975). Period three implies chaos. The American Mathematical
Monthly, 82(10), 985–992. DOI: 10.2307/2318254.

106
BIBLIOGRAFÍA BIBLIOGRAFÍA

Liao, S.-j. (2017). On the clean numerical simulation (CNS) of chaotic dynamic systems.
Journal of Hydrodynamics, 29(5), 729–747. DOI: 10.1016/S1001-6058(16)60785-0.

Lorenz, E. N. (1963). Deterministic Nonperiodic Flow. DOI: 10.1175/1520-

0469(1963)020¡0130:DNF¿2.0.CO;2.

Lorenz, E. N. (1989). Computational chaos-a prelude to computational instability. Physica

D: Nonlinear Phenomena, 35(3), 299–317. DOI: 10.1016/0167-2789(89)90072-9.

Mingjun, J. & Huanwen, T. (2004). Application of chaos in simulated annealing. Chaos,

Solitons and Fractals, 21(4), 933–941. DOI: 10.1016/j.chaos.2003.12.032.

Motsa, S. S., Dlamini, P., & Khumalo, M. (2013). A new multistage spectral relaxation
method for solving chaotic initial value systems. Nonlinear Dynamics, 72(1-2), 265–
283. DOI: 10.1007/s11071-012-0712-8.

Naess, A. & Skaug, C. (2002). The study of chaotic attractors of nonlinear systems by path
integration. Proceedings of the IUTAM Symposium on Nonlinear Stochastic Dynamics.

Nee, S. (2018). Survival and weak chaos. ROYAL SOCIETY OPEN SCIENCE, 5, 1–15. DOI:
10.1098/rsos.172181.

Neta, B. & Ford, C. H. (1984). Families of methods for ordinary differential equations ba-
sed on trigonometric polynomials. Journal of Computational and Applied Mathematics,
10(1), 33–38. DOI: 10.1016/0377-0427(84)90066-9.

Nguyen, H. S., Sidje, R. B., & Huu Cong, N. (2007). Analysis of trigonometric implicit
Runge-Kutta methods. Journal of Computational and Applied Mathematics, 198(1), 187–
207. DOI: 10.1016/j.cam.2005.12.006.

Noorani, M. S. M., Hashim, I., Ahmad, R., Bakar, S. A., Ismail, E. S., & Zakaria, A. M.
(2007). Comparing numerical methods for the solutions of the Chen system. Chaos,
Solitons and Fractals, 32(4), 1296–1304. DOI: 10.1016/j.chaos.2005.12.036.

Olabodé, D. L., Miwadinou, C. H., Monwanou, A. V., & Orou, J. B. C. (2018). Horses-
hoes chaos and its passive control in dissipative nonlinear chemical dynamics. DOI:
10.1088/1402-4896/aacef0.

107
 BIBLIOGRAFÍA BIBLIOGRAFÍA

Owolabi, K. M. & Atangana, A. (2017). Numerical simulations of chaotic and complex

spatiotemporal patterns in fractional reaction–diffusion systems. Computational and
Applied Mathematics. DOI: 10.1007/s40314-017-0445-.

Packard, N. H., Crutchfield, J. P., Farmer, J. D., & Shaw, R. S. (1980). Geometry from a
time series. Phys. Rev. Lett., 45, 712–716. DOI: 10.1103/PhysRevLett.45.712.

Pantaleón, C., Vielva, L., Luengo, D., & Santamarı́a, I. (2003). Bayesian estimation of
chaotic signals generated by piecewise-linear maps. Signal Processing, 83(3), 659–664.
DOI: 10.1016/S0165-1684(02)00485-1.

Parker, T. S. & Chua, L. O. (1989). Practical Numerical Algoithms for Chaotic Systems.
Springer-Verlag. ISBN: 9781461281214.

Paternoster, B. (1998). Runge-Kutta(-Nyström) methods for ODEs with periodic solutions

based on trigonometric polynomials. Applied Numerical Mathematics, 28(2-4), 401–412.
DOI: 10.1016/S0168-9274(98)00056-7.

Pham, V.-T., Volos, C. K., & Vaidyanathan, S. (2015). Multi-scroll Chaotic Oscillator Ba-
sed on a First-Order Delay Differential Equation, (pp. 59–72). Springer International
Publishing.

Polhill, J. G., Izquierdo, L. R., & Gotts, N. M. (2006). What every agent-based modeller
should know about floating point arithmetic. Environmental Modelling and Software,
21(3), 283–309. DOI: 10.1016/j.envsoft.2004.10.011.

Rabinovich, M. I. & Fabrikant, A. (1979). Stochastic self-modulation of waves in none-

quilibrium media. Sov. Phys. JETP 50, (pp. 311).

Robinson, J. C. (2004). An Introduction to Ordinary Differential Equations. Cambridge

University Press. DOI: 0.1017/CBO9780511801204.

Rudin, W. (1976). Principles of mathematical analysis, Third Edition (International Student

Edition). Japan: McGraw-Hill Kogakusha LTD. ISBN: 0-07-054235-X.

108
BIBLIOGRAFÍA BIBLIOGRAFÍA

Saha Ray, S. & Patra, A. (2013). Haar wavelet operational methods for the numerical
solutions of fractional order nonlinear oscillatory Van der Pol system. Applied Mathe-
matics and Computation, 220, 659–667. DOI: 10.1016/j.amc.2013.07.036.

Sauer, T., Grebogi, C., & Yorke, J. A. (1997). How long do numerical chaotic solutions
remain valid? Physical Review Letters, 79(1), 59–62. DOI: 10.1103/PhysRevLett.79.59.

Shafique, A., Sayeed, M., & Tsakalis, K. (2018). Nonlinear Dynamical Systems with Chaos
and Big Data: A Case Study of Epileptic Seizure Prediction and Control, (pp. 329–369).
Springer International Publishing.

Sherratt, J. A. (1997). A comparison of two numerical methods for oscillatory reaction-

diffusion systems. Applied Mathematics Letters, 9659(2), 1–3. DOI: 10.1016/S0893-
9659(97)00001-3.

Silva-Garcı́a, V. M., Flores-Carapia, R., Renterı́a-Márquez, C., Luna-Benoso, B., & Aldape-
Pérez, M. (2018). Substitution box generation using Chaos: An image encry-
ption application. Applied Mathematics and Computation, 332, 123–135. DOI:
10.1016/j.amc.2018.03.019.

Skufca, J. D. (2004). Analysis Still Matters: A Surprising Instance of Failure

of Runge–Kutta–Felberg ODE Solvers. SIAM Review, 46(4), 729–737. DOI:
10.1137/S003614450342911X.

Smith, L. (2007). Chaos: a very short introduction, volume 159. Oxford University Press.
ISBN: 978-0-19-285378-3.

Soriano-Sánchez, A. G., Posadas-Castillo, C., Platas-Garza, M. A., & Elizondo-González,

C. (2016). Chaotic Synchronization of CNNs in Small-World Topology Applied to Data
Encryption, (pp. 337–362). Springer International Publishing.

Sprott, J. C. (2000). Simple chaotic systems and circuits. American Journal of Physics,
68(8), 758–763. DOI: 10.1119/1.19538.

Sprott, J. C. & Sprott, J. C. (2003). Chaos and time-series analysis, volume 69. Oxford
University Press. ISBN: 0-19-850840-9.

109
 BIBLIOGRAFÍA BIBLIOGRAFÍA

Stenflo, L. (1996). Generalized Lorenz equations for acoustic-gravity waves in the at-
mosphere. Physica Scripta, 53(1), 83.

Takens, F. (1981). Detecting strange attractors in turbulence. In D. Rand & L.-S. Young
(Eds.), Dynamical Systems and Turbulence, Warwick 1980 (pp. 366–381). Berlin, Heidel-
berg: Springer Berlin Heidelberg. ISBN: 978-3-540-38945-3.

Tokuda, I., Nagashima, T., & Aihara, K. (1997). Global bifurcation structure of chaotic
neural networks and its application to traveling salesman problems. Neural Networks,
10(9), 1673–1690. DOI: 10.1016/S0893-6080(97)00023-3.

UTK & ORNL (2005). Netlib repository. [Online]. Last accessed 2018-06-06. URL:
http://www.netlib.org/.

Vadillo, F. (1997). On spurious fixed points of Runge-Kutta methods. Journal of Compu-

tational Physics, 132(1), 78–90. DOI: 10.1006/jcph.1996.5615.

Varsakelis, C. & Anagnostidis, P. (2016). On the susceptibility of numerical methods

to computational chaos and superstability. Communications in Nonlinear Science and
Numerical Simulation, 33, 118–132. DOI: 10.1016/j.cnsns.2015.09.007.

Wang, B., Zhong, S. M., & Dong, X. C. (2016). On the novel chaotic secure communication
scheme design. Communications in Nonlinear Science and Numerical Simulation, 39,
108–117. DOI: 10.1016/j.cnsns.2016.02.035.

Wu, J., Liao, X., & Yang, B. (2017). Color image encryption based on chaotic sys-
tems and elliptic curve ElGamal scheme. Signal Processing, 141, 109–124. DOI:
10.1016/j.sigpro.2017.04.006.

Yamaguti, M. & Ushiki, S. (1981). Chaos in numerical analysis of ordinary differen-

tial equations. Physica D: Nonlinear Phenomena, 3(3), 618–626. DOI: 10.1016/0167-
2789(81)90044-0.

Yüzbas, S. (2012). A numerical scheme for solutions of the Chen system. Mathematical
Methods in the Applied Sciences, 35(8), 885–893. DOI: 10.1002/mma.1581.

110
BIBLIOGRAFÍA BIBLIOGRAFÍA

Zhi-Hong, G., Huang, F., & Guan, W. (2005). Chaos-based image encryption algorithm.
Physics Letters, Section A: General, Atomic and Solid State Physics, 346(1-3), 153–157.
DOI: 10.1016/j.physleta.2005.08.006.

Zhou, G., Zhang, D., Liu, Y., Yuan, Y., & Liu, Q. (2015). A novel image encryption
algorithm based on chaos and Line map. Neurocomputing, 169, 150–157. DOI:
10.1016/j.neucom.2014.11.095.

111