Cálculo de los exponentes característicos de Lyapunov. Aplicación a la detección de caos

Traballo Fin de Grao Cálculo de los exponentes característicos de Lyapunov. Aplicación a la detección de caos. Áurea Mene Hevia 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Cálculo de los exponentes característicos de Lyapunov. Aplicación a la detección de caos. Áurea Mene Hevia 07/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Dedicado a: Jerónimo Rodríguez García por su atención y esfuerzo en la realización de este trabajo y por enseñarme grandes cosas. Mi familia y mi pareja por su gran apoyo durante estos años de sacrificio. I II Trabajo propuesto Área de Coñecemento: Matemática aplicada Título: Cálculo de los exponentes característicos de Lyapunov. Aplicación a la detección de caos. Breve descrición do contido Descripción de una herramienta que nos permite detectar órbitas caóticas, que será el cálculo de los exponentes característicos de Lyapunov. Aplicaciones de los mismos. Recomendacións Se recomienda haber superado las materias del Grado en Matemáticas relacionadas con las ecuaciones diferenciales ordinarias y con su resolución numérica. Outras observacións III Índice general Resumen VIII Introducción XI 1. Resultados básicos para EDOs 1 1.1. Definición del problema de valor inicial . . . . . . . . . . . . . . . . . . . . . 1 1.2. Existencia y unicidad de solución . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. Dependencia con respecto a las condiciones iniciales . . . . . . . . . . . . . . 2 2. Integradores para problemas de valor inicial 5 2.1. Método de Euler explícito . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Runge-Kutta explícitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3. Paso de tiempo adaptativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1. Runge-Kutta encajados . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.2. ODE45.Dormand-Prince . . . . . . . . . . . . . . . . . . . . . . . . . 13 3. Cálculo de los exponentes de Lyapunov 15 3.1. Motivación e introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2. Exponentes característicos de Lyapunov . . . . . . . . . . . . . . . . . . . . 18 3.2.1. CálculodelosLCE............................ 22 3.3. Exponentes de Lyapunov de orden superior . . . . . . . . . . . . . . . . . . 24 4. Aplicaciones 29 4.1. AtractordeLorenz................................ 29 4.2. Rössler....................................... 36 4.3. HyperChaos.................................... 38 Apéndices 45 Álgebra exterior 47 V XII INTRODUCCIÓN cia de caos para algunos de los valores de los parámetros; nos referimos a los sistemas de Lorenz [3], Rössler [4] e Hyperchaos [5]. Para alguno de ellos se ha utilizado estos exponentes para observar los valores de los parámetros para los cuales se produce un cambio de comportamiento (diagramas de bifurcación). El documento finaliza con un anexo en el que se incluyen los códigos generados, algún concepto técnico de utilidad sobre el álgebra exterior y una sección de referencias. Capítulo 1 Resultados básicos para EDOs En este trabajo vamos a centrarnos en el cálculo de los exponentes característicos de Lyapunov, que son un indicador cualitativo de la sensibilidad frente al cambio de las condiciones iniciales. Antes de introducir estas magnitudes, que se harán en el capítulo 3, introduciremos algunos resultados básicos sobre las ecuaciones diferenciales ordinarias que relacionaremos con los exponentes. Especial mención deben de recibir los resultados relativos a la continuidad y diferenciabilidad con respecto a las condiciones iniciales y en especial, y relacionado con este último concepto, cabe mencionar la ecuación variacional, es decir, el sistema linealizado, que es de especial importancia para el cáluclo de los exponentes característicos de Lyapunov. En este capítulo seguiremos el libro [1]. 1.1. Definición del problema de valor inicial Definición 1.1. Un problema de valor inicial es una ecuación diferencial x0(t) = f(x(t)),(1.1) con f:E→Rndonde Ees un subconjunto abierto de Rn, junto con un punto (t0,x0)∈E llamada condición inicial. Definición 1.2. Sea el problema de valor inicial    x0(t) = f(x(t)), x(t0) = x0. (1.2) Una solución para dicho problema será una función x(t)diferenciable en un intervalo I ∀t∈ I, x(t)∈E, que es una solución de la ecuación diferencial (1.1) que satisface x(t0) = x0 con t0∈I. 1 2CAPÍTULO 1. RESULTADOS BÁSICOS PARA EDOS 1.2. Existencia y unicidad de solución Para garantizar que el problema de valor inicial (1.2) tenga una única solución vamos a dar unas hipótesis previas sobre la función fque denotaremos por Hf. Hipótesis Hf: 1. La función f:E→Rndebe de ser continua en ambas variables, es decir, f∈ C(E×Rn). 2. ftiene que ser Lipschitziana globalmente, es decir, ∃L≥0constante, tal que kf(x)−f(¯x)k ≤ Lkx−¯xk ∀ x, ¯x∈Rn. Bajo estas condiciones tendremos una única solución x(t)∈ C1(I). Teorema 1.3. Teorema Fundamental de Existencia y Unicidad Sea E un subconjunto abierto de Rnque contiene a x0y supongamos que f verifica las hipótesis Hf. Entonces el problema de valor inicial    x0(t) = f(x(t)), x(t0) = x0, (1.3) tiene una única solución x(t)en todo R. Observación: Nótese que la hipótesis sobre fde ser globalmente lipschitziana es más fuerte que la hipótesis habitual de ser localmente lipschitziana. Se ha considerado esta hipótesis más fuerte por los siguientes motivos: Para simplificar el discurso, para justificar el buen funcionamiento de los métodos numéricos y por último para los casos en los que fsea localmente lipschitziana se puede generar una función ¯ fque en un entorno de las condiciones iniciales coincide con la foriginal y es globalmente lipschitziana, esto proporcionará un resultado de existencia y unicidad local para la foriginal. 1.3. Dependencia con respecto a las condiciones iniciales Dado que en este trabajo nos interesamos por la evolución de dos trayectorias con condiciones iniciales similares comenzamos enunciando algún resultado básico de la teoría cualitativa de la dependencia respecto a las condiciones iniciales. Teorema 1.4. Sea Eun subconjunto abierto de Rnque contiene a x0y sea fverificando las hipótesis Hf. Supongamos que el problema de valor inicial (1.2) tenga solución x(t, x0) definida en el intervalo cerrado I. Entonces existe δ > 0y una constante positiva Ltal que para todo y∈Nδ(x0)el problema de valor inicial    x0(t) = f(x(t)), x(0) = y, 1.3. DEPENDENCIA CON RESPECTO A LAS CONDICIONES INICIALES 3 tiene una única solución x(t, y)definida en I, la cual satisface kx(t, y)−x(t, x0)k⩽eLtky−x0k y l´ım y→x0 x(t, y) = x(t, x0) uniformemente para todo t∈I. Nótese que L es la constante de Lipschitz de las hipótesis Hf. Este teorema nos proporciona el primer resultado cualitativo que nos aporta continuidad y que la dependencia temporal depende de una constante Kque está relacionada con la constante de Lipschitz de la función f. Nuestro objetivo será calcular una magnitud que nos permita calcular una estimación similar a esta de tipo cualitativo y que veremos en el capítulo 3, para ello utilizaremos el corolario 1.6 que jugará un papel muy importante. Observación. Mencionar que el anterior teorema no suele ser óptimo en la práctica. Véase con el siguiente ejemplo: Tengamos la siguiente ecuación diferencial x0= 2x, (1.4) Entonces kf(x)−f(¯x)k= 2kx−¯xkcuya constante de Lipschitz es 2y el comportamiento de kx(t, (t0, x0)) −x(t, (t0,¯x0))kviene dado por e2t(x0−¯x0)con lo que su estimación es bastante precisa. Sin embargo, si tomamos esta otra ecuación diferencial x0=−2x, (1.5) la constante de Lipschitz es nuevamente 2pero el comportamiento de kx(t, (t0, x0)) − x(t, (t0,¯x0))kviene dado por e−2t(x0−¯x0)con lo cual la estimación resulta ser muy pobre por lo tanto necesitaremos buscar un resultado cualitativo que nos solvente estos errores. Teorema 1.5. Dependencia con respecto a las condiciones iniciales Sea Eun subconjunto abierto de Rnque contiene a x0y supongamos que fverifica las hipótesis Hf y f∈ C1(E). Entonces existe un δ > 0tal que para todo y∈Nδ(x0)el problema de valor inicial    x0(t) = f(x(t)), x(t0) = y, (1.6) tiene una única solución x(t, y)con x∈ C1(G)donde G=R×Nδ(x0)⊂Rn+1. Además para cada y∈Nδ(x0),x(t, y)es dos veces continuamente diferenciable respecto de t, con t∈R. 4CAPÍTULO 1. RESULTADOS BÁSICOS PARA EDOS La demostración del anterior teorema se puede ver en [1] de la cual se deduce el siguiente corolario. Corolario 1.6. Bajo las hipótesis del anterior teorema, Φ(t, y) = ∂x ∂y (t, y), con t∈I e y∈Nδ(x0)si y solo si Φ(t, y)es la matriz solución de    Φ0(t, y) = Df[x(t, y)]Φ, Φ(t0, y) = Id, (1.7) con t∈I, y∈Nδ(x0)y donde Id denota la matriz identidad y Df denota la matriz jacobiana evaluada en la solución de (1.6). Observación. El sistema que se plantea en (1.7) a pesar de que es un sistema lineal, no es autónomo cuya dependencia temporal se produce a través de la solución del sistema no lineal por lo que su resolución va a ser difícil, ya que vamos a necesitar resolver a la par el sistema linealizado con el no lineal. Esto nos va a resultar muy costoso computacionalmente. Capítulo 2 Integradores para problemas de valor inicial Para el cálculo de los exponentes característicos de Lyapunov se requerirá resolver la ecuación diferencial ordinaria no lineal, así como el sistema linealizado alrededor de la órbita. Esta tarea en general no puede realizarse con métodos analíticos con lo cual procederemos a realizarlos con métodos numéricos. Para ello contaremos con métodos de integración para problemas de valor inicial, centrándonos en métodos explícitos por su sencillez. Para este capítulo seguiremos el libro [2] y los apuntes de la asignatura Métodos Numéricos en Optimización y Ecuaciones Diferenciales del Grado en Matemáticas de la Universidad de Santiago de Compostela. Dado el siguiente problema de valor inicial    x0(t) = f(x(t)), x(t0) = x0, (2.1) nuestro objetivo es buscar una solución x(t)a dicho problema. Además bajo las hipótesis Hf mencionadas en el capítulo anterior tenemos que la solución será única y que x(t)∈ C1(I). Los métodos numéricos que resuelven (2.1) no nos proporcionan la solución exacta del problema sino que nos dan una aproximación de esta. El funcionamiento de estos métodos numéricos es el siguiente: En primer lugar se realiza una malla del intervalo I= [a, b]con a, b ∈Ry posteriormente se calculan las aproximaciones de la solución exacta en los diferentes nodos equispaciados por un paso h. tn=a+ (n−1)h, n = 1, ..., N 5 6CAPÍTULO 2. INTEGRADORES PARA PROBLEMAS DE VALOR INICIAL t1t2 ab = = t3tN } h h=b−a N−1, que se conoce como stepsize. Es decir, lo que se pretende es calcular aproximaciones de la solución exacta en los diferentes nodos, aproximaciones de x(t1), x(t2), ..., x(tN). Con lo cual tendremos lo siguiente, xn≈x(tn), n = 1, ..., N. 2.1. Método de Euler explícito El método de Euler explícito tiene la siguiente estructura,    x1=x0, xn+1 =xn+hf(xn), n = 1, ..., N −1, (2.2) Una de las interpretaciones de este método se basa en la definición de la derivada pues, x0(tn)≈x(tn+1)−x(tn) hsi hes pequeño. Entonces, f(x(tn)) ≈x(tn+h)−x(tn) h, que es la definición de la primera derivada. De lo que se deduce lo siguiente, x(tn+1)≈x(tn) + hf(x(tn)) si hes pequeño. Por otra parte el orden de convergencia de este método es 1. Se verifica l´ım h→0+( m´ax n∈{1,...,N}kxn−x(tn)k) = 0, donde m´ax n∈{1,...,N}kxn−x(tn)k= Θ(h), es decir, el método es convergente. El error se comporta de forma lineal con respecto a h, como una constante real positiva Cque multiplica ah. Para finalizar esta sección vamos a realizar un análisis de la convergencia del Euler explícito bajo las hipótesis Hf. 2.1. MÉTODO DE EULER EXPLÍCITO 7 Definición 2.1. Definimos el error local de discretización en tn+1 como Tn+1 =hτn+1, n = 1, ..., N −1, donde τn+1 =x(tn+1)−x(tn) h−f(x(tn)), n = 1, ..., N −1, siendo x(t)una solución de la EDO x0=f(t, x). Proposición 2.2. Supongamos las hipótesis Hf. 1. l´ım h→0( m´ax 2⩽n⩽Nkτnk) = 0 para cualquier solución x(t)de la EDO x0=f(t, x). Que un método verifique este primer punto nos indica que el método es consistente. 2. Si x∈ C2([a, b]), entonces m´ax 2⩽n⩽Nkτnk= Θ(h). Cabe preguntarse si en la estructura del método de Euler explícito realizamos dos pequeñas perturbaciones del original si obtendremos soluciones parecidas. Proposición 2.3. Supongamos las hipótesis sobre Hf. Sea el primer sistema perturbado:    z1=η zn+1−zn h=f(zn) + δn+1, n = 1, ..., N −1. (2.3) El segundo sistema perturbado:    ¯z1= ¯η ¯zn+1−¯zn h=f(¯zn) + ¯ δn+1, n = 1, ..., N −1. (2.4) Entonces se tiene que m´ax 1⩽n⩽Nkzn−¯znk⩽eL(b−a)h N X n=2 kδn−¯ δnk. Siendo Lla constante de Lipschitz de las hipótesis Hf. Por lo tanto el método de Euler explícito es estable, es decir, una pequeña perturbación en el problema no afecta notablemente a la solución del mismo, es decir ambas soluciones se van a mantener cerca. Proposición 2.4. Si un método es consistente y estable entonces el método será convergente. Teorema 2.5. Convergencia del Euler explícito. Bajo las hipótesis de Hf. 1. l´ım h→0( m´ax 1⩽n⩽Nkxn−x(tn)k)=0 2. Si x∈ C2([a, b]), entonces m´ax 1⩽n⩽Nkxn−x(tn)k= Θ(h). (Orden de convergencia 1). 8CAPÍTULO 2. INTEGRADORES PARA PROBLEMAS DE VALOR INICIAL 2.2. Runge-Kutta explícitos Como vimos en la anterior sección el método de Euler explícito es un método con orden de convergencia 1. Si nos interesamos por métodos de alto orden podemos pensar en los métodos Runge-Kutta (RK). Los métodos de tipo Runge-Kutta son métodos de un solo paso,y son aquellos que se pueden escribir de la siguiente manera:    x1=x0 xn+1 =xn+hPs i=1 biki, n = 1, ..., N −1, s ∈N, (2.5) donde, ki=f(xn+h s X j=1 aijkj), Lo que caracteriza el orden de un método Runge-Kutta es el error del método, que tiene la siguiente forma Error ≡Chk, donde Ces una constante real positiva y kes el orden del método. Por otra parte sdenota el número de etapas del Runge-Kutta, es decir, el número de veces que se evalúa la función en cada paso. Este concepto es importante porque la evaluación de la función requiere un coste computacional, por esta razón son preferidos los métodos con un número mínimo de etapas como sea posible. Una característica interesante de los métodos Runge-Kutta es que no es necesario calcular derivadas de fpara avanzar, a cambio el precio a pagar consiste en evaluar más veces la función con el consiguiente coste computacional. Es habitual presentar este método con las tablas de Butcher. c1 c2 ⋮ cs a11 a1s as1ass ⋮ ⋮ ⋮ ⋮ a21 a2s b1bs ⋮ Figura 2.1: Estructura de una tabla de Butcher Observación: Nótese que al considerar sistemas autónomos, el vector cde la figura 2.1 no juega ningún papel. Consideramos de todos modos métodos Runge-Kutta que cumplan 2.2. RUNGE-KUTTA EXPLÍCITOS 9 la llamada condición de fila, es decir, ci=X j aij. Ver [2] para las implicaciones correspondientes. Definición 2.6. Un Runge-Kutta se dice explícito si la matriz Ade la tabla de Butcher es una matriz triangular inferior (aij = 0 si i⩽j). Las fórmulas explícitas serían las siguientes:                k1=f(xn) k2=f(xn+ha21k1) . . . ks=f(xn+h[as1k1+... +as,s−1ks−1]) (2.6) Proposición 2.7. Una condición necesaria y suficiente para que el método sea convergente es que Ps i=1 bi= 1. Algunos ejemplos de los métodos Runge-Kutta más importantes son: 0 0 1 Figura 2.2: Tabla de Butcher Euler explícito 0 0 0 00 0 0 0 0 0 0 0 0 1 1/2 1/2 1/2 1/2 0 11/6 2/6 1/62/6 Figura 2.3: Tabla de Butcher Runge-Kutta clásico Este último es un Runge-Kuta clásico de orden 4, es decir, Error(h) = ch4. 16 CAPÍTULO 3. CÁLCULO DE LOS EXPONENTES DE LYAPUNOV A continuación consideramos una perturbación infinitesimal de la nueva condición inicial (¯ t, ¯x)en la dirección unitaria ¯v. De este modo tenemos una nueva solución x(t, (¯ t, ¯x+δ¯v)).(3.3) Nos planteamos a qué velocidad evoluciona la distancia entre ambas órbitas en tiempos cercanos a ¯ t, comparada a la distancia inicial en ¯ tdada por δ, es decir, 1 δkx(t, (¯ t, ¯x+δ¯v)) −x(t, (¯ t, ¯x))k, t ∈(¯ t−, ¯ t+).(3.4) Dado que supondremos los tiempos cercanos a ¯ ty las perturbaciones son pequeñas sabemos que x(t, (¯ t, ¯x+δ¯v)) = x(t, (¯ t, ¯x)) + δΦ(t, (¯ t, Id))¯v+ Θ(δ2), con lo cual la magnitud a estudiar se comporta como kΦ(t, (¯ t, Id))¯vk, t ∈(¯ t−, ¯ t+).(3.5) En las anteriores expresiones la función matricial Φ(t, (¯ t, ¯ Φ)), está caracterizada por ser la solución del siguiente problema de valor inicial matricial lineal con coeficientes variables, no autónomo,    Φ0(t) = Dxf(x(t))Φ(t), Φ(¯ t) = ¯ Φ, (3.6) donde Dxf(x(t)),(3.7) es la matriz jacobiana de fevaluada en x≡x(t, (t0, x0)). Hacemos notar que la función en (3.5) adquiere el valor 1en t=¯ t. De cara a estudiar la convergencia (o divergencia) de ambas órbitas realizamos, para tsuficientemente pequeño, la siguiente identificación kΦ(t, (¯ t, Id))¯vk ≡ eλi(t−¯ t)(3.8) que de algún modo define a λi≡λi(¯ t, ¯v). Esto motiva la siguiente definición Definición 3.1. Dado ¯ ty¯vdefinimos el exponente instantáneo de Lyapunov mediante λi(¯ t, ¯v) := 1 kΦ(t, (¯ t, Id))¯vk d dtkΦ(t, (¯ t, Id))¯vkt=¯ t .(3.9) Pero como kΦ(t, (¯ t, Id))¯vk= 1 entonces d dtkΦ(t, (¯ t, Id))¯vkt=¯ t=λi(¯ t, ¯v).(3.10) 3.1. MOTIVACIÓN E INTRODUCCIÓN 17 Nótese que además, kh(t)k=ph(t)th(t), d dtkh(t)k=1 2√hth2hth0, por lo que, λi(¯ t, ¯v) = ¯vtΦ(¯ t, (¯ t, Id))Φ0(¯ t, (¯ t, Id))¯v= ¯vtDxf(¯x)¯v= ¯vtDxf(x(¯ t))¯v. (3.11) Esta magnitud nos da una idea de cuanto divergen en un instante dos órbitas perturbadas por una cantidad infinitamente pequeña en la dirección de ¯v. A continuación introducimos una magnitud que indique la divergencia en media a lo largo de toda una trayectoria. Consideramos ahora una perturbación en el instante inicial dada por un vector v0. Esta perturbación inicial generará una perturbación en un tiempo futuro ¯ tdada por: ¯v=Φ(¯ t, (t0,Id))v0 kΦ(¯ t, (t0,Id))v0k= ¯v(v0).(3.12) El exponente de Lyapunov para esta dirección, en ¯ t, vendría dado por (ver definición 3.1 con ¯vdado por 3.12): λi(¯ t, ¯v) = 1 kΦ(¯ t, (t0,Id))v0k2vt 0(Φ(¯ t, (t0,Id)))tDxf(¯x)Φ(¯ t, (t0,Id))v0.(3.13) Este exponente depende en realidad de (¯ t, v0)con lo cual podríamos introducir la siguiente definición: ¯ λi(¯ t, v0) := λi(¯ t, ¯v(v0)).(3.14) Para estimar la sensibilidad de la órbita a pequeñas perturbaciones en la dirección v0 promediamos en el tiempo estos exponentes de Lyapunov, con lo cual tenemos: Definición 3.2. Dada una direción v0y un tiempo tse define el exponente característico promediado como λ(v0) := l´ım sup t→∞ 1 t−t0Zt t0 ¯ λi(s, v0)ds (3.15) Veamos a continuación que esta definición coincide con la definición habitual en la literatura [7]. Proposición 3.3. Tenemos que λ(v0) = l´ım sup t→∞ 1 t−t0 ln kΦ(t, (t0,Id))v0k.(3.16) 18 CAPÍTULO 3. CÁLCULO DE LOS EXPONENTES DE LYAPUNOV Demostración. Desarrollamos a continuación la expresión (3.15). Para ello introducimos la función g(t) := 1 2ln kΦ(t, (t0,Id))v0k2= ln kΦ(¯ t, (t0,Id))v0k.(3.17) Está claro que g(t0) = 1 2ln kv0k2=1 2ln(1) = 0, g0(t) = 1 2kΦ(t, (t0,Id))v0k22vt 0(Φ(t, (t0,Id)))tΦ0(t, (t0,Id))v0. (3.18) Empleando (3.13), (3.15) y (3.18) deducimos una expresión alternativa para el exponente característico promediado: λ(v0) = l´ım sup t→∞ 1 t−t0 ln kΦ(t, (t0,Id))v0k(3.19) que es la definición habitual de exponente característico de Lyapunov. Observación: Podríamos llegar a esta misma expresión empleando otros argumentos más intuitivos pero menos rigurosos. Partiendo de la solución del problema de valor inicial no lineal (3.1) construimos la solución del sistema linealizado con condición inicial v0:    z0(t) = Dxf(x(t))z(t), z(t0) = v0. (3.20) Nótese que: z(t)≡Φ(t)v0. Nos interesamos por la evolución de su norma kz(t)kbuscando un comportamiento exponencial: kz(t)k=e˜ λ(t−t0),(3.21) expresión que define ˜ λ(t)por ˜ λ(t, v0) = 1 t−t0 ln kz(t)k=1 t−t0 ln kΦ(t, (t0,Id))v0k(3.22) a comparar con (3.16). 3.2. Exponentes característicos de Lyapunov A continuación veremos que la función λ(·)(exponente promediado de Lyapunov) alcanza tan solo nvalores siendo nla dimensión de la EDO. Estos nvalores recibirán el nombre de exponentes característicos de Lyapunov. (Ver definición 3.7) 3.2. EXPONENTES CARACTERÍSTICOS DE LYAPUNOV 19 Para demostrar el resultado, para un λ∗∈Rfijado introducimos el conjunto V(λ∗) := {v∈Rn/λ(v)⩽λ∗}.(3.23) Para el cual tenemos los siguientes resultados: Proposición 3.4. Los conjuntos V(λ∗)definidos en (3.23) verifican: Si λ∗ 1⩽λ∗ 2⇒V(λ∗ 1)⊂V(λ∗ 2) V(λ∗)es un subespacio vectorial de Rn Demostración. El primer punto resulta evidente. Para demostrar el segundo punto probaremos que dados v1∈V(λ∗), v2∈V(λ∗)yµ∈R entonces λv ∈V(λ∗)yv1+v2∈V(λ∗). λ(µv1) = l´ım sup t→∞ 1 tln kΦ(t)(µv1)k= l´ım sup t→∞ 1 t{ln |µ|+ ln kΦ(t)(v1)k} = l´ım sup t→∞ 1 tln kΦ(t)(v1)k=λ(v1)⩽λ∗ (3.24) Por otra parte, λ(v1+v2) = l´ım sup t→∞ 1 tln kΦ(t)(v1+v2)k⩽l´ım sup t→∞ 1 tln (kΦ(t)(v1)k+kΦ(t)(v2)k) ⩽l´ım sup t→∞ 1 tln (2 m´ax{kΦ(t)(v1)k,kΦ(t)(v2)k}) = l´ım sup t→∞ 1 tln (m´ax{kΦ(t)(v1)k,kΦ(t)(v2)k}) = l´ım sup t→∞ 1 tm´ax{ln (kΦ(t)(v1)k),ln (kΦ(t)(v2)k)} ⩽m´ax l´ım sup t→∞ 1 tln (kΦ(t)(v1)k),l´ım sup t→∞ 1 tln (kΦ(t)(v2)k) = m´ax{λ(v1), λ(v2)}⩽λ∗ (3.25) Como queríamos demostrar. Observación: Nótese que la definición de exponente promediado no necesitaría a priori que el vector que lo define tenga norma 1, sin embargo el exponente instantáneo sí lo necesitaría. Para ser más precisos tenemos la siguiente proposición: Proposición 3.5. La aplicación λ(·)→dim(V(λ∗)) es creciente. Por otra parte existen ¯nvalores λ1< λ2< ... < λ¯ncon ¯n⩽ndonde es discontinua. Por último, para cada λi existe un vi∈Rntal que λ(vi) = λi. 20 CAPÍTULO 3. CÁLCULO DE LOS EXPONENTES DE LYAPUNOV Demostración. De lo visto anteriormente ya tendríamos probado las dos primeras afirmaciones, veamos la última. Para ε > 0suficientemente pequeño, los espacios vectoriales Vi−1=V(λi−ε)(V(λi+ε) = Vi(3.26) no coinciden. Por ello debe existir vi∈V(λi+ε)tal que vi/∈V(λi−ε)(3.27) y por ello λ(vi)⩽λi+εyλ(vi)> λi−ε. (3.28) Pero debido a (3.26) este vector no depende de ε(con tal de que sea suficientemente pequeño) por lo que de (3.28) se deduce: λ(vi)⩽λiyλ(vi)≥λi(3.29) es decir, λ(vi) = λi. Por último vemos que λ(·)alcanza solamente esos valores. Proposición 3.6. La aplicación λ(·)alcanza solamente los valores λi, i ∈ {1, ..., ¯n}, de la Proposición 3.5. Demostración. Supongamos que se alcanza algún otro valor ¯ λ, entonces debe existir ¯v∈Rn/λi<¯ λ=λ(¯v)con ¯ λ<λj∀j > i. (3.30) Sabemos que V(λ∗) = Vi∀λ∗∈(λi,¯ λ+ε)(3.31) al menos para εsuficientemente pequeño. En particular V(µ)≡V(¯ λ)≡Vi∀µ∈(λi,¯ λ). Pero, λ(¯v) = ¯ λ>µ⇒¯v∈V(¯ λ)y¯v /∈V(µ)(3.32) lo cual es una contradicción. Todo esto motiva a la siguiente definición. 3.2. EXPONENTES CARACTERÍSTICOS DE LYAPUNOV 21 Definición 3.7. Llamaremos exponentes característicos de Lyapunov a los números λ1, λ2, ..., λ¯n introducidos en la proposición anterior. Proposición 3.8. Supongamos f∈ C1y que la órbita no tiende a un punto de equilibrio y está acotada, entonces sean λ1, λ2, ..., λ¯nlos exponentes característicos de Lyapunov, se tiene que existe un λi, i ∈ {1, ..., ¯n}tal que, λi= 0. Demostración. La demostración detallada se puede ver en [6], aquí mostraremos un esbozo de la misma. Consideramos la siguiente EDO    x0(t) = f(x(t)), x(t0) = x0, (3.33) y el sistema lienalizado asociado    z0(t) = Dxf(x(t))z, z(t0) = v0. (3.34) Se define el exponente característico de Lyapunov como λ(v0) = l´ım sup t→∞ 1 t−t0 ln kzk.(3.35) Nuestro objetivo será construir una solución del linealizado z(·)tal que l´ım sup t→∞ 1 t−t0 ln kzk= 0. Tomamos z(t) := x0(t)entonces tenemos z0(t) = d dtf(x(t)) = Dxf(x(t))x0=Dxf(x(t))z En consecuencia, z(t)es solución del sistema linealizado. Sea la condición inicial v0:= x0(t0) = f(x(t0)) = f(x0)6= 0 entonces nuestro objetivo será probar que λ(v0)⩽0. Por hipótesis tenemos que fes continua y se anula en un número finito de puntos y que x(·)es acotada y no tiende a un punto de equilibrio. Entonces, kf(x(t))k⩽cte, ∀t⩾0 22 CAPÍTULO 3. CÁLCULO DE LOS EXPONENTES DE LYAPUNOV lo cual implica que |z|⩽cte ya que, |z|=|x0(t)|=|f(x(t))|⩽cte. Como consecuencia tenemos que λ(v0) = l´ım sup t→∞ 1 t−t0 ln kzk⩽l´ım t→∞ 1 t−t0 ln(cte) = 0. Entonces, λ(v0)⩽0.Veamos que es igual a cero. Sea ε > 0existirá un ¯ t⩾t0suficientemente grande tal que si t > ¯ ttendremos que 1 t−t0 ln kzk=1 t−t0 ln kx0k< λ +ε < 0. Tomando λ0:= λ+εy v(t) = v0e−|λ0|(t−t0) = v0eλ0(t−t0) tenemos que kx0k=kzk⩽kv0ke−|λ0|(t−t0),∀t⩾¯ t y esto nos conduce a que kx0(t)k=kf(x(t))k=kz(t)kt→∞ →0. Es decir, x(t)tiende a un punto de equilibrio, en contradicción a las hipótesis propuestas. Concluimos así que λ(v0)=0. Para finalizar esta sección veamos como calcular los exponentes característicos de Lyapunov. 3.2.1. Cálculo de los LCE Veremos que esta tarea no será trivial. Las dificultades que surgirán motivarán la sección posterior. Para calcular estos números parece claro que debemos evaluar λ(·)para algún vector v∈Rn, usando, por ejemplo, la fórmula (3.16). Esta operación puede resultar complicada debido a las dificultades numéricas. En efecto, las magnitudes kΦ(t, (t0,Id))v0kpodrían crecer o decrecer exponencialmente rápido, lo cual conduciría a un problema de overflow ounderflow rápidamente. Para evitar este contratiempo podemos ver lo siguiente. 3.2. EXPONENTES CARACTERÍSTICOS DE LYAPUNOV 23 Sea T > 0arbitrario, un tiempo característico, usando la fórmula (3.15) tenemos: Zt0+¯ kT t0 ¯ λi(s, v0)ds = ¯ k X k=1 Zt0+kT t0+(k−1)T ¯ λi(s, v0)ds = ¯ k X k=1 (g(t0+kT)−g(t0+ (k−1)T)) = ¯ k X k=1 ln kΦ(t0+kT)v0k kΦ(t0+ (k−1)T)v0k = ¯ k X k=1 ln     Φ(t0+kT)v0 1 kΦ(t0+ (k−1)T)v0k    = ¯ k X k=1 ln kΦ(t0+kT, (t0+ (k−1)T, Id))vk−1k. (3.36) Donde vkes el vector unitario definido por: vk:= Φ(t0+kT, (t0,Id))v0 kΦ(t0+kT, (t0,Id))v0k.(3.37) En consecuencia, podemos resolver el problema (3.1) para a continuación calcular zk−1(·) solución de:    z0 k−1(t) = Dxf(x(t))zk−1(t), zk−1(t0+ (k−1)T) = vk−1. (3.38) Tendremos finalmente: Zt0+¯ kT t0 ¯ λi(s, v0)ds = ¯ k X k=1 ln kzk−1(t0+kT)k(3.39) y entonces: λ(v0) = l´ım sup t→∞ 1 t−t0  t T X k=1 ln kzk−1(t0+kT)k+ ln kzt T(t)k (3.40) donde zes solución de (3.38). Dado que este razonamiento se puede realizar para cualquier T > 0, bastará tomarlo de tal modo que las soluciones zk−1(t) adquieran valores razonables en el intervalo [t0+ (k−1)T, t0+kT]. 24 CAPÍTULO 3. CÁLCULO DE LOS EXPONENTES DE LYAPUNOV Es decir, que las magnitudes kΦ(t, (t0,Id))v0kcrezcan o decrezcan moderadamente. Una segunda dificultad que se nos presenta en el cálculo de estos exponentes, λi, i ∈ {1, ..., ¯n}es que tendríamos que seleccionar un vi∈Vi≡V(λi)tal que vi/∈Vi−1=V(λi−1) y evaluar λ(vi). En ese caso λi−1< λ(vi)⩽λi⇒V(λi) = λi. Dado que estos espacios Vi=V(λi)no se conocen a priori esto muestra que si tomamos un v∈Rnarbitrario, digamos de forma aleatoria, tendremos que λ(v) = λn, siendo λn el mayor de todos. Esto se dará con probabilidad uno, ya que la dimensión de V¯n−1es estrictamente inferior a la dimensión de V¯n(V¯n−1tendría medida nula como subvariedad de V¯n). Veámoslo con la siguiente ilustración: Figura 3.1: Representación esquemática de los espacios V(λ)en un caso de dimensión 3. Si la perturbación estuviera en la recta tendríamos el menor exponente, en cambio si la perturbación estuviese en el plano pero no en la recta entonces obtendríamos el segundo exponente mientras que si la perturbación no está en el plano, y por tanto tampoco en la recta, entonces obtendremos el mayor exponente de todos. Este es el motivo por el que necesitamos introducir los exponentes característicos de Lyapunov de orden superior. 3.3. Exponentes de Lyapunov de orden superior De la definición vista en la sección anterior (3.7) se deduce entonces que el número de exponentes de Lyapunov no es superior a n, la dimensión del espacio, y que los V(λ∗)a 3.3. EXPONENTES DE LYAPUNOV DE ORDEN SUPERIOR 25 partir de unos valores aislados de λ∗cambian de dimensión, lo que hace que el cálculo sea muy difícil para calcular los exponentes característicos de Lyapunov que no sean el mayor. Entonces, para calcular los demás vamos a introducir los llamados exponentes de orden superior (p > 1). Se trata de reproducir el proceso explicado anteriormente con la diferencia de que en vez de observar como cambian las distancias entre las órbitas ahora veremos como cambian las áreas perturbadas al propagarlas mediante el flujo del linealizado (de modo más genérico, como cambian los p-volúmenes). Consideramos las condiciones iniciales infinitamente pequeñas en un politopo engendrado por pvectores, que propagaremos por el linealizado (ver figura 3.2). x y z · ·¯ P P0v0 w0x0=x(t0,(t0,x0)) ¯x=x(¯ t,(t0,x0)) ¯v=Φ(¯ t,(t0,Id))v0 ¯w=Φ(¯ t,(t0,Id))w0 Figura 3.2: Imagen por el linealizado de un paralelogramo para p= 2. Vamos a presentar un par de propiedades que cumplirán estos p-volúmenes. Propiedad 1: Si se amplian los vectores v0yw0(ver figura 3.2), es decir, si hacemos µv0yηw0con µ, η > 0entonces los correspondientes ¯vy¯wtambién se amplían, sin embargo, el cociente entre el área generada por ¯ PyP0no varía, se mantiene constante. Propiedad 2: Si se toman otros dos vectores ¯v0y¯w0unitarios que están en el mismo espacio vectorial que generan v0yw0entonces las áreas de los paralelogramos pueden 32 CAPÍTULO 4. APLICACIONES Dado que ∂x ∂x0(t, 0, x0)es solución de la ecuación variacional    Ψ0=Df(x)Ψ Ψ(0) = I (4.7) tenemos usando la identidad de Abel-Jacobi det(∂x ∂x0 (t, 0, x0)) = det(∂x ∂x0 (0,0, x0))eRt 0tr(Df(x(s)))ds =e−(a+b+1)t, ya que la jacobiana del campo que define el sistema tiene una traza constante, es decir, independiente de la solución de la ecuación. Finalmente obtenemos que Vx=V0e−(a+b+1)t, es decir, el volumen tiende exponencialmente a cero con el tiempo. En consecuencia, toda órbita que comience en esta esfera tiende a acercarse exponencialmente rápido con el tiempo al conjunto R∞=R0∩R1∩... ∩Rk∩... que tiene volumen nulo. Este resultado tiene especial interés pues nos permite relacionar el volumen de la región Rtcon los exponentes característicos de Lyapunov, pues la suma de los tres exponentes coincide con el valor −(a+b+ 1), que es la traza de la matriz jacobiana (4.2). Este resultado nos da una herramienta para comprobar si efectivamente los cálculos realizados son correctos, lamentablemente esto solo se produce para el sistema de Lorenz, cuya traza de la matriz jacobiana es constante. Para otros sistemas saber con exactitud si los cálculos son correctos se hace más complicado. Un caso de interés, debido a la riqueza de los resultados sería hacer depender el sistema de un único parámetro, es decir, si hacemos bparámetro, entonces pondremos aycen función de b. [3] Entonces, sea el sistema de Lorenz,          x0=a(y−x) y0=cx −xz −y z0=xy −bz (4.8) considerando bcomo parámetro tal que b∈[0.138,0.148] y a=b+1+p2(b+ 1)(b+ 2) 4.1. ATRACTOR DE LORENZ 33 c=aa+b+ 3 −(−a+b+ 1) Tenemos que la gráfica para los tres exponentes de Lyapunov y la suma de ambos es la siguiente, 0.138 0.139 0.14 0.141 0.142 0.143 0.144 0.145 0.146 0.147 0.148 -5 -4 -3 -2 -1 0 1 Figura 4.1: Expontenes de Lyapunov y su suma para Lorenz Si observamos la figura 4.2 en la que hemos hecho un zoom de los dos primeros exponentes podemos observar una gran riqueza de resultados. Vamos a analizarlos. 0.138 0.139 0.14 0.141 0.142 0.143 0.144 0.145 0.146 0.147 0.148 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 (a) Primer y segundo exponente. 0.1445 0.145 0.1455 0.146 0.1465 0.147 0.1475 0.148 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 (b) Zoom para los valores (0.143, 0.148). Figura 4.2: Zoom para los dos primeros exponentes. 34 CAPÍTULO 4. APLICACIONES Para valores de b∈[0.138,0.1441] no obtenemos caos, pues ningún exponente de Lyapunov es positivo. Por otra parte vemos claramente la presencia de un superatractor para el valor b= 0.1416 cuya figura del atractor se puede ver en la figura 4.3. Finalmente para el valor b= 0.1397 hay un ciclo de período dos que atrae a las órbitas. (Ver figura 4.4). Figura 4.3: Superatractor de Lorenz para b= 0.1416.Exponentes: λ1=−0.0004, λ2= −0.1427, λ3=−4.3513. Figura 4.4: Atractor de Lorenz (ciclo de orden 2) para b= 0.1397.Exponentes: λ1= −0.0004, λ2=−0.0031, λ3=−4.4843. Para valores de b∈[0.1441,0.148] vemos que el primer exponente de Lyapunov se vuelve positivo, y por lo tanto comenzamos a observar caos. Un ejemplo del atractor de Lorenz para un valor caótico se puede ver en la figura 4.5. 4.1. ATRACTOR DE LORENZ 35 Figura 4.5: Atractor de Lorenz para b= 0.147 Exponentes: λ1= 0.0371, λ2= −0.0000, λ3=−4.5503. Finalmente para terminar esta sección es necesario mostrar el caso más nombrado en la bibliografía [3] y en cualquier documento que consultemos sobre el sistema de Lorenz. La llamada mariposa (ver figura 4.6) para un caso de caos resultante de darle los siguientes valores a los parámetros del sistema: a= 10, b =8 3, c = 28 Figura 4.6: Atractor de Lorenz para a= 10, b =8 3yc= 28.Exponentes: λ1= 0.8987, λ2= 0.0021, λ3=−14.5675. Observación: De acuerdo con la proposición vista en el capítulo anterior (3.8) podemos ver en los anteriores ejemplos que si tenemos un órbita acotada entonces uno de los exponentes característicos de Lyapunov siempre es 0. 36 CAPÍTULO 4. APLICACIONES 4.2. Rössler El sistema de Rössler es un sistema tridimensional no lineal estudiado por Otto Rössler en 1976. Estas ecuaciones diferenciales definen un sistema dinámico que muestra situaciones caóticas asociadas con el atractor. El documento original de Rössler [4] afirma que el atractor de Rössler estaba destinado a comportarse de manera similar al atractor de Lorenz (ver sección 4.1), pero también a ser más fácil de analizar cualitativamente. Una órbita dentro del atractor sigue una espiral exterior cerca de un punto fijo inestable, una vez que el gráfico gira en espiral lo suficiente, un segundo punto fijo influye en el gráfico, causando un aumento y giro en el eje z. En el dominio del tiempo, se hace evidente que aunque cada variable está oscilando dentro de un rango fijo de valores, las oscilaciones son caóticas. Otto Rössler diseñó el atractor pero luego se descubrió que las ecuaciones originalmente teóricas eran útiles para modelar el equilibrio en reacciones químicas. El sistema de Rössler es el siguiente          x0=−y−z y0=x+ay z0=b+z(x−c) (4.9) cuya matriz jacobiana es la siguiente     0−1−1 1a0 z0x−c     (4.10) Observación: A diferencia del sistema de Lorenz vemos que la traza de la matriz jacobiana no es constante, con lo que el argumento de la traza para comprobar que el cálculo de los exponentes es correcto no es posible en este caso, con lo que su verificación se hace más complicada. Ahora bien, al igual que en la sección anterior vamos a considerar el siguiente caso de interés. Consideramos los parámetros a, b fijos y hacemos variar cdel siguiente modo: a= 0.1, b = 0.1, c ∈[2,18] cuya gráfica de los tres exponentes de Lyapunov y su suma se puede ver en la figura 4.7. Nuevamente si hacemos un zoom a los dos primeros exponentes tenemos la figura 4.8. Vamos a analizar los resultados vistos: 4.2. RÖSSLER 37 2 4 6 8 10 12 14 16 18 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 Figura 4.7: Exponentes de Lyapunov y su suma para Rössler Para valores de c∈[2,9) vemos que e primer exponente es prácticamente nulo, con lo cual no tenemos caos. Como ocurría en el sistema de Lorenz también tenemos un superatractor para un valor aproximado a c= 2.4. Por otra parte destacar que para valores de c= 4 tenemos un ciclo de período 1, para c= 6 un ciclo de período 2, para c= 8.5un ciclo de período 4 y para c= 8.7 un ciclo de período 8. (Ver figura 4.9). Para valores de c∈[9,18] vemos que el primer exponente se vuelve notablemente positivo y por tanto tenemos una zona caótica. ( Ver figura 4.10). Si observamos con atención la figura ampliada 4.8 vemos que para el valor c= 18 hay una zona bastante más caótica que para el valor c= 9 (Ver 4.11). Cabe resaltar que para valores de ccomprendidos entre [12,13] el primer exponente vuelve a valer prácticamente 0volviendo así a la existencia de ciclos de período 3 y 6 para valores de cigual a 12 y12.6respectivamente. (Ver figura 4.12). Para finalizar, en la múltiple bibliografía que se puede consultar, entre ellas [4], para el sistema de Rössler los valores habituales de los parámetros para observar caos son los siguientes: a= 0.2, b = 0.2, c = 5.7 38 CAPÍTULO 4. APLICACIONES 2 4 6 8 10 12 14 16 18 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 (a) Primer y segundo exponente. 9 10 11 12 13 14 15 16 17 18 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 (b) Zoom para los valores (9,18). Figura 4.8: Zoom para los dos primeros exponentes. Para los cuales tenemos la figura 4.13. 4.3. HyperChaos Un atractor hipercaótico se define usualmente como un sistema que presenta un comportamiento caótico con al menos dos exponentes característicos de Lyapunov positivos. Además de poseer un exponente nulo y otro negativo para así asegurar el límite de la solución. Por lo tanto es importante resaltar que la dimensión mínima para tener un sistema hipercaótico continuo es dimensión 4. El primer sistema tetradimensional fue propuesto en 1979 por Rössler [5],                x0=−y−z y0=x+ay +w z0=b+xz w0=−cz +dw (4.11) cuya matriz jacobiana es la siguiente       0−1−1 0 1a0 0 z0x0 0 0 −c d       (4.12) del mismo modo que para el sistema de Rössler vemos que la traza de la matriz jacobiana 4.3. HYPERCHAOS 39 (a) c= 4, período 1. (b) c= 6, período 2. (c) c= 8.5, período 4. (d) c= 8.7, período 8. Figura 4.9: Distintos tipos de ciclos para el sistema de Rössler no es constante, esto nos lleva a deducir la gran importancia y peculiaridad que el sistema de Lorenz aporta. Finalmente en la literatura [5] nos podemos encontrar con los siguientes valores para los parámetros del sistema (4.11) en el que obtenemos un comportamiento hipercaótico. Estos son a= 0.25, b = 3, c = 0.5, d = 0.05 para los cuales obtenemos los cuatro exponentes característicos de Lyapunov correspondientes λ1= 0.0930, λ2= 0.0211, λ3= 0.0004, λ4=−22.6686 Las distintas representaciones tridimensionales del atractor hipercaótico son las siguientes 40 CAPÍTULO 4. APLICACIONES Figura 4.10: Atractor de Rössler para un valor caótico c= 9. Exponentes: λ1= 0.0303, λ2= 0.0017, λ3=−8.8669. (ver figura 4.14). 4.3. HYPERCHAOS 41 Figura 4.11: Atractor de Rössler para un valor caótico c= 18. Exponentes: λ1= 0.0961, λ2= 0.0055, λ3=−17.9360. (a) c= 12, período 3. (b) c= 12.6, período 6. Figura 4.12: Ciclos para el sistema de Rössler 48 ÁLGEBRA EXTERIOR cada vector vjpuede escribirse como una combinación lineal de los vectores base ei. Usando la bilinealidad del producto exterior, esto puede expandirse a una combinación lineal de productos exteriores de esos vectores base y cualquier producto exterior en el que el mismo vector base aparezca más de una vez es cero. Además cualquier producto exterior en el que los vectores base no aparezcan en el orden adecuado se pueden reordenar. En consecuencia, dim( p ^(V)) = n p!.(13) Ahora vamos a introducir el producto escalar entre elementos de Vp(V) < v1∧... ∧vp, w1∧... ∧wp>p=det     < v1, w1>··· < v1, wp> . . ..... . . < vp, w1>··· < vp, wp>     .(14) Esto permite calcular la norma de un elemento de este espacio, kv1∧...∧vpk2 p=det     < v1, v1>··· < v1, vp> . . ..... . . < vp, v1>··· < vp, vp>     =det     kv1k2··· < v1, vp> . . ..... . . < vp, v1>··· kvpk2     . (15) finalmente, tenemos que la raíz cuadrada de esta norma coincide con el área del p-volumen correspondiente. Observación: Cabe mencionar que si los vectores vianteriores son ortogonales tendremos que kv1∧... ∧vpk2 p=det     kv1k2 ... kvpk2     =kv1k2·... ·kvpk2.(16) Códigos de programación .0.1. Cálculo de los exponentes %%%%%%%%%%%%%%%%%%%%%%%% % Exponentes .m % Programa que c a lc u la l o s CLE con ODE45. % Aurea Mene Hevia % Universidad de Santiago de Compostela %%%%%%%%%%%%%%%%%%%%%%%% % Condicion i n i c i a l . x0 = repmat ( [ 0 ; 1 ; 0 ] , 1 , 3 ) ; t0 = 0 . 0 ; z0 = [ 1 . 0 0.0 0 . 0 ; . . . 0.0 1.0 0 . 0 ; . . . 0.0 0.0 1 . 0 ] ; % Parametros del sistema de Lorenz . %a l f a [ 1 0 , 1 6 ] , beta [ 8 / 3 , 4 ] , ro (0 , i n f i n i t y ) % c = [10;... % 8 / 3 ; . . . % 28] alfamin = 10; alfamax = 10; 49 50 CÓDIGOS DE PROGRAMACIÓN betamin = 8/3; betamax = 8/3; romin = 28; romax = 28; cBounds = [ alfamin betamin romin ; . . . alfamax betamax romax ] ; nb_alfa = 1; nb_beta = 1; nb_ro = 1; a l f a = linspace ( alfamin , alfamax , nb_alfa ) ; beta =linspace ( betamin , betamax , nb_beta ) ; ro = linspace ( romin , romax , nb_ro) ; options = odeset ( ’ RelTol ’ ,1 . e−5) ; % Numero de i t e r a c i o n e s . DeltaT = 0 . 2 5 ; T = 2000.0; Npas = T/DeltaT ; m = size ( z0 , 1 ) ; % Numero de condic ion es i n i c i a l e s x=x0 ; t=t0 ; z=z0 ; SUM = zeros(m, 1 ) ; CLE = zeros( nb_alfa , nb_beta , nb_ro , 3 ) ; MLCE = zeros(m, 1 ) ; figure (1) ; clf hold on 51 figure (2) clf hold on for i_alfa = 1: nb_alfa i_al fa for i_beta = 1: nb_beta for i_ro = 1: nb_ro x=x0 ; t=t0 ; z=z0 ; c = [ a l f a ( i_al fa ) ; beta( i_beta ) ; ro ( i_ro ) ] ; f = @( t , xz ) LorenzfDxf ( t , xz , c ) ; SUM = zeros(m, 1 ) ; LE = [ Inf ;Inf ;Inf ] ; for k=1:Npas tspan = [ ( k−1)∗DeltaT+t0 , k∗DeltaT+t0 ] ; for i = 1:m xz = [ x ( : , i ) ; z ( : , i ) ] ; [ t , xz ] = ode45( f , tspan , xz , options ) ; i f k>7500 figure (1) plot ( t , xz ( : , 1 ) , ’b−’ , t , xz ( : , 2 ) , ’ r−’ , t , xz ( : , 3 ) , ’g−’ , ’ Linewidth ’ ,3) figure (2) plot3( xz ( : , 1 ) , xz ( : , 2 ) , xz ( : , 3 ) , ’b−’,’ 52 CÓDIGOS DE PROGRAMACIÓN Linewidth ’ ,3) drawnow end t = t (end) ; xz = xz (end , : ) . ’ ; z ( : , i ) = xz (m+1:end) ; x ( : , i ) = xz ( 1 :m) ; end % Calculos de l a s c on tr ib uci on es a l o s exponentes x = repmat ( xz ( 1 :m, 1 ) ,1 ,3) ; z = GramSchmidtOrth(z) ; r = sqrt(sum( z .^2 ,1) ) ; z = z ./ repmat ( r , size ( z , 1 ) ,1) ; SUM = SUM + log ( r ) . ’ ; CLE( i_alfa , i_beta , i_ro , : ) = SUM/( tspan (2)−t0 ) ; % Test de parada por s i ya alcanzamos la estabilidad t o l = 1 . e −6; LU = SUM/( tspan (2)−t0 ) ; i f norm((LU−LE) ./(1+norm(LU, Inf) ) , Inf)<t o l [ k , Npas ] disp ( ’ Se␣ha␣ cumplido ␣ e l ␣ t e s t ␣de␣ parada ’ ) break ; 53 else LE = LU; end end [ i_alfa , i_beta , i_ro ,CLE( i_alfa , i_beta , i_ro , 1 ) ,CLE( i_alfa , i_beta , i_ro , 2 ) ,CLE( i_alfa , i_beta , i_ro , 3 ) ] end end end % Representacion g r a f i c a de l o s exponentes y su suma figure (3) clf; hold on ; % Si vari a e l parametro a l f a plot ( a lfa ,CLE( : , 1 , 1 , 1 ) , ’b−o ’ , ’ Linewidth ’ ,4) plot ( a lfa ,CLE( : , 1 , 1 , 2 ) , ’ r−o ’ , ’ Linewidth ’ ,4) plot ( a lfa ,CLE( : , 1 , 1 , 3 ) , ’g−o ’ , ’ Linewidth ’ ,4) plot ( a lfa ,CLE( : , 1 , 1 , 1 )+CLE( : , 1 , 1 , 2 )+CLE( : , 1 , 1 , 3 ) , ’m−o ’ , ’ Linewidth ’ ,4) plot ( a lfa , a l f a ∗0 , ’k−−’ , ’ Linewidth ’ ,4) % Si va ria e l parametro beta plot (beta ,CLE( 1 , : , 1 , 1 ) , ’b−o ’ , ’ Linewidth ’ ,4) plot (beta ,CLE( 1 , : , 1 , 2 ) , ’ r−o ’ , ’ Linewidth ’ ,4) 54 CÓDIGOS DE PROGRAMACIÓN plot (beta ,CLE( 1 , : , 1 , 3 ) , ’g−o ’ , ’ Linewidth ’ ,4) plot (beta ,CLE( 1 , : , 1 , 1 )+CLE( 1 , : , 1 , 2 )+CLE( 1 , : , 1 , 3 ) , ’m−o ’ , ’ Linewidth ’ ,4) plot (beta ,beta∗0 , ’k−−’ , ’ Linewidth ’ ,4) % Si va ria e l parametro ro v1=zeros(4 ,1) ; v2=zeros(4 ,1) ; v3=zeros(4 ,1) ; vsum=zeros(4 ,1) ; for i = 1 : nb_ro v1 ( i ) = CLE(1 ,1 , i , 1 ) ; v2 ( i ) = CLE(1 ,1 , i , 2 ) ; v3 ( i ) = CLE(1 ,1 , i , 3 ) ; vsum( i ) = sum(CLE(1 ,1 , i , : ) ) ; end plot ( ro , v1 , ’b−o ’ , ’ Linewidth ’ ,4) plot ( ro , v2 , ’ r−o ’ , ’ Linewidth ’ ,4) plot ( ro , v3 , ’g−o ’ , ’ Linewidth ’ ,4) plot ( ro , vsum , ’m−o ’ , ’ Linewidth ’ ,4) plot ( ro , ro ∗0 , ’k−−’ , ’ Linewidth ’ ,4) 55 .0.2. Algoritmo Gram-Schmidt modificado function v = GramSchmidtOrth(v) %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % GramSchmidtOrth .m % Programa que r e a l i z a e l metodo GramSchmidtOrth % Aurea Mene Hevia % Universidad de Santiago de Compostela %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Ejemplo : % % consideramos 2 ve ct or es v1 =[3; 1] , v2 = [2 ;2 ] % Primero l o s colocamos en una matriz , ya que l o s v ec to re s se establecen en forma de columna % v = [3 2;1 2] % %A = GramSchmidtOrth(v) % %A = % 0.9487 −0.3162 % 0.3162 0.9487 % % prueba para asegurarse de que es co r r ec t o % % dot (A( : , 1 ) ,A( : , 2 ) ) % % ans = % 0 % k = size (v , 2 ) ; for i i = 1 : 1 : k for j j = i i +1:1:k v ( : , j j ) = v ( : , j j ) −proj (v ( : , i i ) ,v ( : , j j ) ) ; end end 56 CÓDIGOS DE PROGRAMACIÓN function w = proj (u , v) % Esta funcion proyecta e l v ector v en e l vector u w = (dot(v , u) / dot(u , u) ) ∗u ; end end Bibliografía [1] Perko Lawrence, Differential Equations and Dynamical Systems, 2nd ed. [2] Wanner G., Hairer E., Solving Ordinary Differential Equations V.I, 2nd ed. [3] Lorenz Edward, Deterministic Nonperiodic Flow, Volumen 20. Massachusetts Institute of Technology (1963). [4] Rössler O.E, An Equation For Continuous Chaos, Volumen 57A, No5. Germany (1976). [5] Rössler O.E, An Equation For Hyperchaos, Volumen 71A, No2,3. Germany (1979). [6] Haken H., At Least One Lyapunov Exponent Vanishes If The Trajectory Of An Attractor Does Not Contain A Fixed Point, Volumen 94A, No2. Germany (1982). [7] Barreira Luís, Lyapunov Exponents, Birkhäuser. [8] Benettin Giancarlo, Galgani Luigi, Giorgilli Antonio, Strelcyn Jean-Marie. Lyapunov Characteristic Exponents For Smooth Dynamical Systems And For Hamiltonian Systems; A Method For Computing All Of Them., (1980). 57