Circuito de chua lineal a trozos: órbitas homoclinas y comportamiento caótico

GRADO EN INGENIER´ IA EL´ ECTRICA Escuela Polit´ecnica Superior Curso acad´emico 2022-2023 Trabajo fin de grado CIRCUITO DE CHUA LINEAL A TROZOS: ´ ORBITAS HOMOCLINAS Y COMPORTAMIENTO CA ´ OTICO Autor: Alexandro Moya Montes Tutor: Victoriano Carmona Centeno Departamento: Matem´atica Aplicada II A lo largo de mi vida, has sido mi mayor apoyo. Tus palabras de aliento y tu amor incondicional han sido fundamentales en mi camino acad´emico. Gracias por creer en m´ı, por sacrificarte, por ense˜narme el valor del esfuerzo y la perseverancia y por ense˜narme a quedarme siempre con el lado bueno de las cosas. Este logro no hubiera sido posible sin ti. Te dedico este trabajo como muestra de mi profundo agradecimiento y amor. Te quiero Mam´a. Agradecimientos A mi tutor, Victoriano, por su orientaci´on experta y su invaluable conocimiento. A mi padre, Salvador, por brindarme la oportunidad de cumplir mis metas. A mi hermano, Adri´an, por iluminarme siempre que se me oscurec´ıa el porvenir. A mi apoyo incondicional, Cristina, por creer en mi cuando yo mismo no lo hac´ıa. A mi compa˜nero de batallas, Camacho, por acabar juntos lo que empezamos juntos. Y a mi. Sevilla, 21 de Julio de 2023 Alexandro Moya Montes Resumen Este trabajo se centra en el estudio de un sistema din´amico espec´ıfico, el circuito de Chua lineal a trozos. En los diferentes cap´ıtulos, se aborda la descripci´on detallada del circuito, incluyendo su historia y las ecuaciones diferenciales que lo rigen. Se presentan las ecuaciones adimensionalizadas que simplifican el an´alisis y se exploran los sistemas din´amicos lineales a trozos en general, destacando sus propiedades y elementos b´asicos. Se introduce el principio de Shilnikov y se profundiza en la existencia de ´orbitas homocl´ınicas y su relaci´on con el caos. Finalmente se prueba mediante m´etodos num´ericos la existencia de una homoclina. Adem´as, se incluye un cap´ıtulo dedicado a la implementaci´on pr´actica de Matlab, donde se proporcionan los c´odigos utilizados para el estudio. El objetivo principal de este trabajo es profundizar en el an´alisis de los sistemas din´amicos lineales a trozos, buscando obtener resultados significativos y contribuir al conocimiento en este campo de estudio. iv ´ Indice general 1. Introducci´on 1 2. Circuito de Chua lineal a trozos 3 2.1. Sistemas ca´oticos y sus caracter´ısticas . . . . . . . . . . . . . . . . . . . 3 2.1.1. Caracter´ısticas del caos . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2. Sistemas ca´oticos . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2. Breve revisi´on hist´orica . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3. G´enesisdelcircuito............................. 6 2.4. Ecuaciones diferenciales dimensionadas . . . . . . . . . . . . . . . . . . 7 2.4.1. Leyes de Ohm y Kirchoff . . . . . . . . . . . . . . . . . . . . . . 7 2.4.2. Elementos almacenadores de energ´ıa . . . . . . . . . . . . . . . 8 2.4.3. DiododeChua ........................... 8 2.4.4. Ecuaciones del circuito . . . . . . . . . . . . . . . . . . . . . . . 10 2.5. Ecuaciones diferenciales adimensionadas . . . . . . . . . . . . . . . . . 11 3. Sistemas din´amicos lineales a trozos 16 3.1. Definici´on del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2. Propiedades del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1. Existencia y unicidad . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3. Puntosdeequilibrio............................. 20 3.4. Elementos din´amicos b´asicos y estabilidad de los puntos de equilibrio . 20 3.5. Semiaplicaciones de Poincar´e . . . . . . . . . . . . . . . . . . . . . . . . 21 4. Resoluci´on de los sistemas lineales asociados 23 4.1. Puntosdeequilibrio............................. 23 4.2. Variedades invariantes . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3. Soluci´on del sistema diferencial (lineal a trozos) . . . . . . . . . . . . . 28 5. ´ Orbitas homoclinas: Comportamiento ca´otico y Teorema de Shilnikov 31 5.1. Puntosdetangencia ............................ 32 v ´ INDICE GENERAL vi 6. Existencia de orbitas homoclinas 34 6.1. ContextoyObjetivos............................ 34 6.2. Resultados y Metodolog´ıa . . . . . . . . . . . . . . . . . . . . . . . . . 34 6.3. Limitaciones................................. 36 6.3.1. Limitaciones de la funci´on ODE45 . . . . . . . . . . . . . . . . . 36 6.3.2. Necesidad de utilizar el m´etodo ode78 . . . . . . . . . . . . . . . 37 6.3.3. Consideraciones sobre la fiabilidad de los resultados . . . . . . . 37 6.4. Discusi´on y Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . 39 7. C´odigos en Matlab 41 7.1. Coeficientes zona derecha . . . . . . . . . . . . . . . . . . . . . . . . . . 41 7.2. Matriz A zona derecha . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 7.3. Coeficientes zona central . . . . . . . . . . . . . . . . . . . . . . . . . . 42 7.4. Matriz A zona central . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 7.5. Calculo de punto de interseccion . . . . . . . . . . . . . . . . . . . . . . 43 7.6. Representacion gr´afica zona central . . . . . . . . . . . . . . . . . . . . 43 7.7. Representacion gr´afica zona derecha . . . . . . . . . . . . . . . . . . . . 43 7.8. Representacion orbita homoclina . . . . . . . . . . . . . . . . . . . . . 44 7.9. Solucion del sistema zona derecha . . . . . . . . . . . . . . . . . . . . . 45 7.10. Solucion del sistema zona izquierda . . . . . . . . . . . . . . . . . . . . 45 7.11. Sistema para ODE derecho . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.12. Sistema para ODE central . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.13.SistemaparaODE ............................. 46 7.14.ODE45.................................... 46 7.15. ODE45 con baja tolancia . . . . . . . . . . . . . . . . . . . . . . . . . . 47 7.16.ODE78.................................... 48 7.17. C´alculo de tiempo. Zona central 1 . . . . . . . . . . . . . . . . . . . . . 48 7.18. C´alculo de tiempo. Zona central 2 . . . . . . . . . . . . . . . . . . . . . 49 7.19. C´alculo de tiempo. Zona derecha . . . . . . . . . . . . . . . . . . . . . 49 7.20.Comprobaciones............................... 49 7.21. Puntos de tangencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 7.21.1. Punto de tangencia T . . . . . . . . . . . . . . . . . . . . . . . . 50 7.21.2. Punto de tangencia S . . . . . . . . . . . . . . . . . . . . . . . . 50 Bibliograf´ıa 52 ´ Indice de figuras 2.1. Comportamiento ca´otico del sistema de Lorenz. Mariposa de Lorenz. . 4 2.2. Diciembre de 1993. Leon Ong Chua sosteniendo el primer modelo de su circuito. ................................... 6 2.3. Esquema del circuito de Chua. . . . . . . . . . . . . . . . . . . . . . . . 7 2.4. Curva caracter´ıstica del diodo de Chua: intensidad IDfrente a tensi´on V1. 10 2.5. Nudo1. ................................... 10 2.6. Nudo2. ................................... 11 3.1. Representaci´on gr´afica de los 3 semiespacios RR,RC,RLy los planos de separaci´on ΣRy ΣC............................. 18 3.2. Ejemplo gr´afico de la primera semiaplicaci´on de Poincar´e definida para nuestrosistema3.6.............................. 22 3.3. Ejemplo gr´afico de la segunda semiaplicaci´on de Poincar´e definida para nuestro sistema 3.6. N´otese el cambio en los semiespacios RCyRR respecto a la primera figura. . . . . . . . . . . . . . . . . . . . . . . . . 22 5.1. Ejemplo gr´afico de una ´orbita homoclina. . . . . . . . . . . . . . . . . . 31 6.1. α= 14,43746643008159, β= 19. Dos ´orbitas homoclinas (sim´etricas) tipo Shilnikov. Cada color representa un cambio de semiespacio. Prestar especial atenci´on al origen de coordenadas donde vemos que ambas trayectorias se ”enrollan” lo que quieren a decir que tienden al infinito en dicho punto. Ver Figura 6.2 para apreciar una vista detallada . . . . 35 6.2. Imagen detallada del origen de coordenadas de la Figura 6.1 donde se aprecia claramente la tendencia de ambas ´orbitas (celeste y amarilla) a ”enrollarse” en el origen de coordenadas. . . . . . . . . . . . . . . . . 35 6.3. Resultado de utilizar la funcion ode45 con tolerancia por defecto. Es evidente que el resultado no es el esperado, pues se espera una ´orbita que retorne a su punto de origen. . . . . . . . . . . . . . . . . . . . . . 36 vii ´ INDICE DE FIGURAS viii 6.4. Resultado mejorado usando tolerancia de 1−15. Se puede apreciar claramente como la ´orbita no tiende al punto de equilibrio. . . . . . . . 37 6.5. Soluci´on obtenida usando ode78. A simple vista parece satisfactorio,sin embargo, en la Figura 6.5. podemos apreciar como las trayectorias no tienden a infinito en el punto de origen. . . . . . . . . . . . . . . . . . . 38 6.6. Detalle de la soluci´on obtenido usando ode78. Se puede observar claramente como no se tiende al punto de equilibrio, pues ambas orbitas (amarilla y morada) deber´ıa ”enrollarse” al origen de coordenadas. . . 38 Cap´ıtulo 1 Introducci´on El mundo actual est´a en constantes cambios, ya sean naturales o artificiales. Todo fen´omeno natural sufre cambios constantes. Algunos de estos cambios son de f´acil percepci´on, como por ejemplo, la migraci´on de las aves o el cambio en la duraci´on de los d´ıas. Sin embargo, otros son infinitamente m´as complicados, como el transporte de energ´ıa o la propagaci´on de una enfermedad, fen´omenos a la orden del d´ıa, que se sit´uan a la cabeza de la problem´atica mundial. Para ordenar y sistematizar todas estas formas de cambio con el fin de estudiar y predecir su comportamiento, ser´ıa ideal poder representarlos de una manera comprensible. Fue en este punto donde naci´o una rama de la matem´atica conocida como los sistemas din´amicos. Todo proceso en el que hay movimiento, entendido como variaci´on a lo largo del tiempo, puede ser considerado como un sistema din´amico. Esto deja en evidencia que el cambio y el movimiento son sin´onimos y no pueden existir el uno sin el otro. Para construir un sistema din´amico a partir de un fen´omeno en evoluci´on a lo largo del tiempo es necesario seguir un proceso cient´ıfico que consiste, principalmente en: Detectar las variables que influyen en el fen´omeno de inter´es Analizar su comportamiento relativo, es decir, determinar c´omo influyen unas sobre otras y c´omo se comportan a lo largo del tiempo. Como resultado de este proceso se obtendr´a una ecuaci´on diferencial (o un sistema de ellas). En este trabajo, se abordar´a el estudio de un sistema din´amico particular: el circuito de Chua lineal a trozos. En el Cap´ıtulo 2, se describir´a en detalle este circuito, incluyendo su historia y las ecuaciones diferenciales que lo gobiernan. Adem´as, se presentar´an las ecuaciones adimensionalizadas que permiten un an´alisis simplificado, 1 CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 8 donde nes el n´umero total de corrientes que entran o salen de un nodo e Ii representa la corriente en´esima que entra o sale del nodo: La ley de Kirchhoff de mallas establece que la suma algebraica de las ca´ıdas de tensi´on en un lazo cerrado de un circuito debe ser igual a cero: m X i=1 Vi= 0,(2.3) donde mes el n´umero total de caidas de tensi´on en una malla cerrada y Vi representa el voltaje en´esimo en la malla. A pesar de ser menos popular en aplicaciones en el mundo real, es especialmente ´util para analizar circuitos con varias fuentes de alimentaci´on y bucles de diferentes tama˜nos. 2.4.2. Elementos almacenadores de energ´ıa Ser´a necesario definir las ecuaciones que rigen el comportamiento de los elementos que conforman nuestro circuito. Para ello, tambi´en nos basaremos en los principios de conservaci´on de la carga (condensador) y energ´ıa (bobina). La ecuaci´on del condensador, muestra que este se opone al cambio en la tensi´on a trav´es de ´el y almacena energ´ıa en forma de campo el´ectrico: IC=CdVC dt ,(2.4) donde ICes la corriente que fluye a trav´es del condensador, Ces la capacidad del condensador y dVC dt es la tasa de cambio de la tensi´on a trav´es del condensador en el tiempo. La ecuaci´on de la Bobina, muestra que una bobina se opone al cambio en la corriente que fluye a trav´es de ella y almacena energ´ıa en forma de campo magn´etico: VL=LdI dt ,(2.5) donde VLes la ca´ıda de voltaje a trav´es de la bobina, Les la inductancia de la bobina ydI dt es la tasa de cambio de la corriente a trav´es de la bobina en el tiempo. 2.4.3. Diodo de Chua El diodo de Chua, es en realidad un conjunto de diodos y amplificadores operacionales que se comporta como un sistema lineal a trozos. Su funci´on es CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 9 retroalimentar el circuito y mantenerlo oscilando actuando como una resistencia negativa. La curva de tensi´on-intensidad puede tener diversas formas, pero el circuito de Chua original especifica una funci´on continua definida a trozos, donde en cada uno de los tramos la funci´on es lineal. La expresi´on de la funci´on tensi´on-intensidad depende del valor de la variable tensi´on V1, y se determina teniendo en cuenta que −V0yV0son los puntos de quiebre de la funci´on lineal a trozos del diodo de Chua. En la Figura 2.4 podemos observar como en el tramo central, cuando −V0≤V1≤ V0, la funci´on es una recta de pendiente m0<0 que pasa por el origen de coordenadas. Por tanto, la funci´on en este tramo es un segmento de la recta, ID=m0V1, cuyos puntos terminales son (−V0,−m0V0) y (V0, m0V0). Para los otros dos tramos restantes, la funci´on lineal es continua y las pendientes de las rectas coinciden, siendo su valor m1<0. Cuando el valor de la tensi´on es menor que −V0, la recta de pendiente m1debe pasar por el punto (−V0,−m0V0). Sustituyendo dicho punto en la expresi´on ID=m0V1+nse obtiene el t´ermino n=m1−m0. Por lo que, cuando V1≤ −V0, la expresi´on de la recta es ID=m1V1−(m0−m1)V0. Por otro lado, cuando el valor de la tensi´on es mayor o igual que V0, la recta de pendiente m1debe pasar por el punto (V0, m0V0). Por lo que, en este tramo, la funci´on tensi´on-intensidad viene dada por la expresi´on ID=m1V1+ (m0−m1)V0, para todo V1≥V0. Teniendo en cuenta la definici´on de la funci´on tensi´on-intensidad en cada uno de los tramos, la funci´on no lineal IDest´a definida como: ID(V1) =      m1V1−(m0−m1)V0si V1≤ −V0, m0V1si |V1| ≤ V0, m1V1+ (m0−m1)V0si V1≥V0. (2.6) donde m0ym1son las pendientes de las rectas en los diferentes tramos de linealidad, y−V0yV0son los puntos de quiebre de la funci´on lineal a trozos del diodo de Chua. CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 10 Figura 2.4: Curva caracter´ıstica del diodo de Chua: intensidad IDfrente a tensi´on V1. 2.4.4. Ecuaciones del circuito Una vez que conocemos el valor que toma la funci´on intensidad frente a tensi´on (2.6) en el diodo de Chua, el siguiente paso es obtener el sistema de ecuaciones din´amicas que describe al circuito. Para ello, vamos a aplicar todas las ecuaciones que hemos estado estudiando en apartados anteriores. Se ha denotado por I1eI2a las intensidades que pasan por los condensadores C1yC2, respectivamente. Por IRse denota a la intensidad que pasa por la resistencia, por ILa la intensidad que pasa por la bobina y por IDa la intensidad que hay en el diodo de Chua. En primer lugar, vamos a analizar el nudo 1. Aplicando la ley de corrientes de Kirchhoff (2.2), se deduce que la suma de las corrientes que pasan por el nudo 1 (ver Figura 2.5) es igual a cero I2−IR+ID= 0,(2.7) donde IR,I2eIDvienen dadas por las expresiones (2.1), (2.6) y (2.4), respectivamente. Las expresiones de las intensidades IReI2en (2.7) obtenemos la siguiente ecuaci´on Figura 2.5: Nudo 1. diferencial: C2 dV2 dt =1 R(V1−V2)−ID.(2.8) CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 11 Aplicando de nuevo la ley de corrientes de Kirchoff en el nudo 2 (ver Figura 2.6): IL−I1−IR= 0,(2.9) donde las expresiones de las intensidades IReI1(2.1) y (2.4), y la intensidad que pasa por la bobina, IL, satisface la relaci´on (2.5). Reemplazando en la ecuaci´on (2.9), Figura 2.6: Nudo 2. obtenemos: C1 dV1 dt =1 R(V2−V1) + IL.(2.10) Por ultimo, aplicando la ley de tensiones de Kirchhoff (2.3) en el bucle que contiene al condensador C1y la bobina, podemos determinar que la tensi´on total. Esto es: dIL dt L=−V1.(2.11) Por lo tanto, el sistema de ecuaciones diferenciales que describe el circuito de Chua est´a compuesto por las ecuaciones (2.8), (2.10) y (2.11)            C2 dV2 dt =1 R(V1−V2)−ID, C1 dV1 dt =1 R(V2−V1) + IL, dIL dt L=−V1. (2.12) 2.5. Ecuaciones diferenciales adimensionadas Es evidente que dada la complejidad del sistema (2.12), hay que tomar acci´on. Para ello procedemos a simplificarlo, de est´a manera nos ahorraremos tiempo en c´alculos posteriores y el manejo de unidades. Denotaremos a los potenciales de los condensadores C1yC2como V1=xeV2=y, respectivamente, y se considerar´a una tensi´on ficticia que va a ser V3=ILR=z. Una vez definido las tensiones, se procede a adimensionalizar las nuevas variables introducidas dividiendo por V0, valor del punto de quiebre de la funci´on no lineal del diodo de Chua. De este modo, se obtienen las siguientes expresiones: x=V1 V0 ,(2.13) CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 12 y=V2 V0 ,(2.14) z=V3 V0 =ILR V0 .(2.15) El sistema evoluciona a lo largo del tiempo, por lo que deberemos tenerlo en cuenta a la hora de adimensionalizar. Para ello, realizamos la siguiente transformaci´on τ=t RC2 ,(2.16) donde τes el tiempo caracter´ıstico del circuito y ten tiempo en segundos. N´otese que la unidad de medida de la resistencia es el ohmio y la del condensador es el faradio, midi´endose el producto de ambas en segundos. Una vez obtenido el tiempo caracter´ıstico del sistema, vamos a obtener las ecuaciones que rigen el sistema en funci´on de la variable temporal τ. Para ello, derivaremos las expresiones dadas respecto a τ. Por un lado, se calcula la derivada de la variable zrespecto a τutilizando la regla de la cadena: dz dτ =dz dIL dIL dt dt dτ .(2.17) La expresi´on de la derivada de zrespecto a la intensidad de la bobina es inmediata y viene dada por dz dIL =R V0 .(2.18) A partir de la expresi´on (2.11), se obtiene el valor de la derivada de la intensidad IL respecto a la variable temporal t dIL dt =−V2 L. Por otro lado, a partir de la expresi´on (2.16) podemos obtener la relaci´on t=RC2τ, cuya derivada respecto a la variable τes inmediata y viene dada por la expresi´on dt dτ =RC2. Sustituyendo las ambas expresiones en (2.17), obtenemos dz dτ =−V2 V0 ·R2C2 L. Teniendo en cuenta la definici´on dada en (2.14), denotando por β=R2C2 L, se transforma en la siguiente dz dτ =−β·y. (2.19) Calculando la derivada de la variable yrespecto a τ, utilizando la regla de la cadena dy dτ =dy dV2 dV2 dt dt dτ , CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 13 donde la derivada de la variable temporal trespecto de τviene dada por la expresi´on (2.16). Es inmediato comprobar que las expresiones de las otras dos derivadas vienen dadas por: dy dV2 =1 V0 , dV2 dt =1 RC2 (V1−V2) + IL C2 , Sustituyendo las expresiones de las derivadas dadas, se tiene, realizando algunas manipulaciones, la siguiente expresi´on dy dτ =V1 V0 −V2 V0 +ILR V0 . Ahora, teniendo en cuenta la notaci´on introducida, se obtiene que la expresi´on se transforma en la siguiente dy dτ =x−y+z. (2.20) Por ´ultimo, vamos a calcular la derivada de la variable xrespecto a τ. Para ello, usando la regla de la cadena se obtiene dx dτ =dx dV1 dV1 dt dt dτ .(2.21) donde la derivada de la variable temporal trespecto de τviene dada por la expresi´on (2.16), la derivada de xrespecto de la variable V1se obtiene directamente a trav´es de la definici´on de la funci´on dx dV1 =1 V0 . Y la derivada del voltaje V1respecto de tse calcula directamente, siendo su valor dV1 dt =1 RC1 ·(V2−V1) + ID C1 . Sustituyendo cada una de estas derivadas en la expresi´on (2.21), se deduce tras realizar algunas simplificaciones la siguiente relaci´on: dx dτ =C2 C1 ·V2 V0 −V1 V0 +R V0 ID. Denotando por α=C2 C1 y teniendo en cuenta la definici´on de las variables dadas, la expresi´on se transforma en la siguiente: dx dτ =α·(y−x+RID V0 ),(2.22) CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 14 donde la funci´on no lineal del diodo de Chua, ID, est´a dada por la expresi´on (2.6). Para obtener la expresi´on de la derivada de xrespecto a la variable temporal en cada una de las zonas de linealidad, introduciremos la siguiente funci´on: k(x) = x−RID V0 .(2.23) En la zona central comprendida entre los dos puntos de quiebre de la funci´on tensi´on-intensidad, la funci´on no lineal est´a definida por ID=m0V1. Por lo que, la expresi´on de la funci´on k(x) en esta zona ser´a: k(x) = x−R(m0V1) V0 . Teniendo en cuenta la expresi´on de la variable xdada, se tiene que: k(x) = x−R(m0V1) V0 = (1 −m0R)x=ax, donde hemos llamado a= 1 −m0R. En la zona derecha, cuando la tensi´on V1≥V0, la funci´on no lineal tensi´onintensidad est´a definida mediante ID=m1V1+ (m0−m1)V0. Siguiendo el mismo procedimiento, a partir de (2.24) se obtiene la expresi´on k(x) = x−R[m1V1+ (m0−m1)V0] V0 = (1 −m1R)x−R(m0−m1) = bx +a−b, siendo b= 1 −m1R. En el tramo izquierdo,cuando V1≤ −V0, la funci´on no lineal est´a definida mediante ID=m1V1−(m0−m1)V0. Sustituyendo dicha expresi´on en (2.24) se obtiene que k(x) = x−R[m1V1−(m0−m1)V0] V0 = (1 −m1R)x+ (m0−m1)R=bx −a+b. Teniendo en cuenta la definici´on de la funci´on k(x) en las tres zonas de linealidad k(x) =      bx −a+bsi x≤ −1, ax si |x| ≥ 1, bx +a−bsi x≥1. (2.24) podemos reescribir la expresi´on de la derivada de xrespecto a la variable temporal en funci´on de dichas zonas como dx dt =α(y−k(x)). CAP´ ITULO 2. CIRCUITO DE CHUA LINEAL A TROZOS 15 Finalmente, utilizando las ecuaciones (2.19), (2.20) y (2.5) y las variables adimensionadas, se obtiene un sistema de ecuaciones diferenciales no lineal de primer orden      ˙x=α(y−k(x)), ˙y=x−y+z, ˙z=−βy, (2.25) donde los par´ametros α=C2 C1 yβ=R2C2 Lson estrictamente positivos. Cap´ıtulo 3 Sistemas din´amicos lineales a trozos 3.1. Definici´on del sistema Para definir el tipo de sistema din´amico que consideramos en la memoria, en primer lugar debemos introducir una notaci´on general. En todo lo que sigue, ˙ xdenota la derivada de xrespecto de la variable temporal t,eiel i-´esimo vector de la base can´onica de Rn, y Mn(R) el conjunto de las matrices cuadradas de orden ncon elementos en R. Mediante ⟨·,·⟩ denotamos al producto escalar usual de Rny por ∥·∥a la norma eucl´ıdea, asociada a dicho producto escalar. En general, un sistema din´amico n-dimensional se puede representar por la ecuaci´on ˙ x=F(x), donde x= (x1, . . . , xn)T∈RnyF(x) = (f1(x), . . . , fn(x))T:Rn→Rn. Conocido esto, un sistema lineal puede formularse mediante ˙ x=Ax+b=F(x). donde A∈ Mnyb∈ Rn. En el caso de que un sistema din´amico est´e representado por una o varias ecuaciones diferenciales ordinarias aut´onomas o no forzadas, es decir que no depende de la variable temporal, de la forma ˙x=F(x), el sistema din´amico se dice que es aut´onomo. Por el contrario, si la ecuaci´on que modela al sistema es no aut´onoma o forzada ˙x=F(x, t), es decir depdende de su variable temporal, el sistema din´amico es no aut´onomo. Teniendo en cuenta lo anterior mencionado, podemos definir el sistema que ser´a objeto de estudio en este trabajo. Definici´on 3.1 Decimos que la ecuaci´on diferencial aut´onoma ˙ x=F(x),con x= (x, y, z)Tdefine un sistema din´amico continuo lineal a trozos tri-zonal en R3si existen cuatro vectores b1,b2,b3,v∈R3, con v=0, tres matrices A1, A2, A3∈ M3(R)y dos 16 CAP´ ITULO 3. SISTEMAS DIN ´ AMICOS LINEALES A TROZOS 17 escalares δ1, δ2∈R, con δ1< δ2, tales que ˙ x=     A0x+b0si ⟨v,x⟩+δ1<0, A1x+b1si −δ1≤ ⟨v,x⟩ ≤ δ2, A2x+b2si ⟨v,x⟩+δ2>0. (3.1) Los planos de ecuaci´on {⟨v,x⟩+δi= 0}(i= 1,2) se denominan planos de separaci´on. Los planos de separaci´on dividen el espacio en tres regiones, en cada uno de las cuales el sistema (3.1) es lineal. Adem´as, por ser el campo vectorial del sistema una funci´on continua se verifica A0x+b0=A1x+b1,∀x∈ {⟨v,x⟩ − δ1}, A1x+b1=A2x+b2,∀x∈ {⟨v,x⟩ − δ2}. Realizando un adecuado cambio de variable podemos transformar los planos de separaci´on en los planos de ecuaci´on {x= 1}y{x=−1}y, por tanto, el sistema (3.1) puede escribirse de la siguiente forma ˙ x=f(x) =      ALx+bLsi x < −1, ACx+bCsi −1≤x≤1, ARx+bRsi x > 1. (3.2) donde AL, AC, AR∈ M3(R) y bL,bC,bRson vectores de R3. Los planos x=−1 y x= 1, dividen al espacio en tres semiespacios, que denotamos por RL,RCyRR, denominamos semiespacio izquierdo, central y derecho respectivamente. RL=(x, y, z)T∈R3:x < −1 RC=(x, y, z)T∈R3:−1≤x≤1 RR=(x, y, z)T∈R3:x > 1 Los planos de separaci´on se denotar´an por ΣL=(x, y, z)T∈R3:x=−1y ΣR= (x, y, z)T∈R3:x= 1(ver Figura 3.1). Definici´on 3.2 Un sistema ˙ x=F(x),es sim´etrico respecto al origen si F(x) = −F(−x)para cualquier x∈Rn. Obviamente, el sistema ser´a sim´etrico respecto al origen si y solo si AR=AL,bR=−bL ybC= 0. CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 24 componente de los puntos de equilibrio tambi´en debe ser nula y el par´ametro tampoco αes nulo, para que exista un punto de equilibrio en el sistema, la funci´on k(x) (2.24), debe anularse en dicho punto. Como ya sabemos, esta funci´on es lineal a trozos, por lo que calcularemos sus valores en la zona central y en una de las zonas exteriores, pues recordemos que nuestro sistema es sim´etrico. En el semiespacio central, tenemos que k(x) = ax, donde a= 1 −m0R= 0. Como hemos dicho, para que exista un punto de equilibrio en este semiespacio, la primera componente de ´este debe anularse. Por tanto, teniendo en cuenta los razonamientos anteriores, se deduce que el punto de equilibrio en RCes el origen de coordenadas, al que denominaremos p0= (0,0,0). En el semiespacio derecho RR, la funci´on no lineal est´a definida por k(x) = bx+a−b, donde a= 1 −m0Ryb= 1 −m1R. Teniendo en cuenta que aybson distintos de 0 e imponiendo k(x) = 0 se obtiene x= 1 −a b. El punto de equilibrio en el semiespacio derecho queda definido como pR=     b−a b 0 a−b b      , donde se tiene que verificar que 1 −a b>1. Teniendo en cuenta la ya mencionada simetr´ıa 3.2 del sistema, en el semiespacio izquierdo RL, el punto de equilibrio queda definido como pL=     a−b b 0 b−a b      , donde se tiene que verificar que a b−1<−1. Para determinar la configuraci´on local de dichos puntos de equilibrio, tendremos que calcular los autovalores de las matrices de coeficientes asociadas al sistema (3.2). De nuevo, recordar que debido a que el sistema es sim´etrico, los autovalores asociados a la matriz de coeficientes en el semiespacio derecho e izquierdo coinciden. El polinomio caracter´ıstico de la matriz AR=  −αb α 0 1−1 1 0−β0  est´a dado por PAR(λ) = −λ3+p1λ2+p2λ+p3,(4.1) CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 25 donde p1=−1−αb,p2=α(1 −b)−β, y p3=−αβb. Notese que el t´ermino independiente del polinomio caracter´ıstico no se anula, lo que quiere decir que no pueden existir autovalores nulos. Este polinomio caracter´ıstico tiene grado impar y coeficientes reales, por lo que al menos tiene una ra´ız real. Las otras dos ra´ıces del polinomio caracter´ıstico las supondremos complejas conjugadas [11]. A continuaci´on vamos a determinar el valor de las ra´ıces del polinomio PARen el semiespacio derecho RR. En esta zona, la ra´ız real del polinomio caracter´ıstico la denotaremos por λR, y las ra´ıces complejas conjugadas de dicho polinomio por σR±ωRi. Para determinar dichas ra´ıces, se utilizar´an algunas propiedades que verifican los autovalores asociados a la matriz AR. En primer lugar, sabemos que la traza de la matriz AR,−1−αb, tiene que coincidir con la suma de los autovalores, λR+ 2σR. A partir de aqu´ı, se deduce que el valor de σRen funci´on de λRy el par´ametro b σR=1 2(−λR−αb −1).(4.2) Por otro lado, la parte imaginaria de los autovalores complejos se determina utilizando que el determinante de A, cuyo valor −αβb, ha de coincidir con el producto de los autovalores, λR(σ2 R+ω2 R). De esta relaci´on, se determina que el valor de ωRes el siguiente: ωR=rαβ λR b−σ2 R.(4.3) En el semiespacio central RC, la ra´ız real del polinomio caracter´ıstico la denotaremos por λC, y por σC±ωCia las ra´ıces complejas conjugadas de dicho polinomio. Determinaremos dichas ra´ıces utilizando las propiedades que tienen que verificar los autovalores asociados a la matriz AC. AC=  −αa α 0 1−1 1 0−β0  El polinomio caracter´ıstico de la matriz es PAC(λ) = −λ3+p1λ2+p2λ+p3,(4.4) donde p1=−1−αa,p2=α(1 −a)−β, y p3=−αβa. Notese que el t´ermino independiente del polinomio caracter´ıstico no se anula, lo que quiere decir que no pueden existir autovalores nulos. CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 26 En primer lugar, sabemos que la traza de la matriz, −1−αa, tiene que coincidir con la suma de los autovalores, λC+ 2σC. A partir de aqu´ı, se deduce que el valor de σCen funci´on de λCy el par´ametro b σC=−1 2(λC+αb + 1).(4.5) Por otro lado, la parte imaginaria de los autovalores complejos se determina utilizando que el determinante, cuyo valor −αβa, ha de coincidir con el producto de los autovalores, λC(σ2 C+ω2 C). De esta relaci´on, se determina que el valor de ωCes el siguiente: ωC=r−αβ λC a−σ2 C.(4.6) 4.2. Variedades invariantes De cara a obtener una soluci´on, es necesario obtener las variedades invariantes asociadas a los puntos de equilibrio del sistema (3.2). En concreto nos centraremos en la zona central, pues es la ´unica que nos resulta de inter´es. Ver el trabajo de M. Camacho [1], donde es calculan las variedades de las zonas externas. Recordemos que ACposee un par de autovalores complejos conjugados, σC±ωCi, que determinan la variedad bidimensional, y un autovalor real λC, que determina la variedad unidimensional. Denotaremos ambas variedades por PF(p0) y RC(p0), respectivamente. A partir de este punto tomaremos, como se hace en Medrano [8], a=−8 7yb=−5 7. La variedad RC(p0) es localmente una variedad lineal unidimensional la regi´on central RC. M´as concretamente, en este semiespacio dicha variedad es una recta que pasa por el punto de equilibrio p0y tiene como vector director un autovector asociado al autovalor real λC. El primer paso para determinar las ecuaciones de esta variedad ser´a calcular un autovector de ACasociado al autovalor λC. Para ello, resolveremos el sistema de ecuaciones (AC−λCI)v = 0, siendo Ila matriz identidad,    α 7−λCα0 1−1−λC−β 0−β−λC     v1 v2 v3 = 0.(4.7) Se obtiene que el subespacio propio asociado al autovalor λCes  V(λC) =        αλC βα 7−λC −λC β 1        CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 27 Recordando el car´acter linealmente dependiente de los autovectores, podemos expresarlo como  V(λC) =      αλC α 7−λC −λC β      (4.8) Se determinan las ecuaciones de RC(p0) localmente en RCsabiendo que es una variedad lineal que pasa por el origen de coordenadas que est´a generada por el autovector  V(λC). RC(p0) = p0+µ V(λC) =  0+     αλC α 7−λC µ −λCµ βµ      donde − α 7−λC αλC ≤µ≤ α 7−λC αλC (4.9) Observemos que en el caso de que el valor de µsea igual a cero obtendremos el punto de equilibrio p0, y si el valor de µcoincide con el l´ımite superior de la desigualdad, la recta interseca al plano de separaci´on ΣRen el punto IR=         1 α 7−λC α βα 7−λC αλC         . La variedad PF(p0), est´a localmente contenida en el semiplano que pasa por el origen de coordenadas y que est´a generado por los autovalores complejos asociados a la matriz AC,σC±ωCi. Teorema 4.1 Dada una matriz A∈ M3(R), con autovalores conocidos λ∈Ry σ±ωi, el subespacio propio de ATasociado a λes ortogonal al plano generado por los autovectores de Aasociados a los autovalores complejos, y por lo tanto un vector normal al plano. A la luz del Teorema 4.1, el primer paso para obtener la ecuaci´on general de dicho semiplano focal PFes determinar el subespacio propio de AT Casociado al autovalor real λC, el cual est´a generado por los autovectores que verifican el sistema (AT C−λCI)nλC=   α 7−λC1 0 α−1−λC−β 0 1 −λC     v1 v2 v3 = 0.(4.10) Resolviendo el sistema de ecuaciones se obtiene f´acilmente que CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 28 nλC=     −λC α 7−λc λC 1      es un autovector de AT Cque genera el subespacio propio asociado a λC. Este autovector es ortogonal a PFy por tanto, su vector director. La ecuaci´on del semiplano focal central se determina imponiendo que pase por el punto p0y tenga por vector normal anλC. De ese modo, se obtiene que su ecuaci´on es PF=    −λC α 7−λc x+λCy+z= 0   donde −1≤x≤1 (4.11) Es obvio que, la intersecci´on del semiplano focal central PFcon el plano de separaci´on ΣRes una recta dada por RI1=     1 r λC +1 α 7−λC r      , donde r∈R. De forma similar, podr´ıamos obtener todas las variedades invariantes para el resto de zonas. V´ease el trabajo de M. Camacho [1]. 4.3. Soluci´on del sistema diferencial (lineal a trozos) Esta secci´on esta dedicada a describir detalladamente la soluci´on anal´ıtica del sistema diferencial en la zona central RR. En el trabajo realizado por M. Camacho[1], se encuentra la soluci´on anal´ıtica para la zona central. Vamos a determinar la soluci´on del problema de valores iniciales asociado al sistema 3.2 en el semiespacio derecho RRcon una condici´on inicial arbitraria, X(0) = p= (xp, yp, zp).Esto es (X′=ARX(t) + bR X(0) = pdonde X(t) = X(t)SGH +X(t)SP NH En primer lugar para c´alcular X(t)SGH , el sistema a resolver ser´ıa    −2 7α α 0 1−1 1 0−β0     x y z =  0 0 0 (4.12) CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 29 Recordemos que este sistema posee un autovalor real denominado λR, y dos autovalores complejos conjugados, los cuales son σR±iωR. A partir de estos autovalores, obtendremos el autovector asociado al autovalor λR vR=      αλR α1 + −2 7α+λR λR −β       , y un autovector complejo asociado al par de autovalores complejos conjugados, que denominaremos γR=uR±qRi uR±qRi=  (1 + σR)σR−ω2 R σR −β +i  ωR(1 + 2σR) ωR 0 . Una vez conocidos los vectores vR,uRyqR, la expresi´on de la soluci´on particular no homog´enea del sistema ser´a X(t) = c1X1(t) + c2X2(t) + c3X3(t).(4.13) donde c1, c2yc3son coeficientes constates y X1(t) = eλRtvR, X2(t) = eσRt(cos(ωRt)uR−sen(ωRt)qR), X3(t) = eσRt(cos(ωRt)qR+sen(ωRt)uR). Obteniendo X=c1X1(0) + c2X2(0) + c3X3(0) . Las soluciones particulares X1, X2yX3en el instante t= 0 ser´an X1(0) = vR, X2(0) = cos(0)uR−sin(0)qR=uR=uR, X3(0) = cos(0)qR+ sin(0)uR=qR=qR. Por otro lado, X(t)SP NH es la soluci´on del sistema AX +B= 0   −α(1 + b)α0 1−1 1 0−β0   x y z +  α(b−a) 0 0 = 0,(4.14) Cuya soluci´on, no es ni mas ni menos que el punto de equilibrio en la zona derecha pR=3 2,0,−3 2T . CAP´ ITULO 4. RESOLUCI ´ ON DE LOS SISTEMAS LINEALES ASOCIADOS 30 Entonces,imponiendo nuestra condici´on inicial, X(0) = p, la soluci´on del PVI vendr´a dada por la expresi´on c1vC+c2uC+c3qC+pR=p. Despejando pRy desarollando     αλR α(1 + b) + λR (1 + σR)σR−ω2 RωR(1 + 2σR) λRσRωR −β−β0      c1 c2 c3 =     xp−3 2 yp zp+3 2      .(4.15) De forma an´aloga, podr´ıamos proceder con los sistemas diferenciales, lineales en las zonas restantes [1]. Finalmente , es notable que el c´alculo de los coeficientes posee un grado de dificultad elevado, incluso podr´ıamos concluir que dicho sistema es imposible de resolver de forma convencional. Por lo que se ha recurrido, a un m´etodo alternativo: desarrollar una serie de scripts en MatLab (ver Cap´ıtulo 7) con la finalidad de abordar el sistema usando metodolog´ıa num´erica. Cap´ıtulo 5 ´ Orbitas homoclinas: Comportamiento ca´otico y Teorema de Shilnikov Anteriormente hemos mencionado en varias ocasiones que nuestro trabajo esta enfocado enfocado en la identificaci´on de orbitas homocl´ınicas de tipo Shilnikov, en este cap´ıtulo definiremos que son las orbitas homoclinas, exhibir la estrecha relaci´on que poseen con el caos y como probar su existencia. Teorema 5.1 (Teorema de Shilnikov [11]) 1 Supongamos que el sistema diferencial ˙ x=f(x),x∈R, siendo funa funci´on continua, con derivadas parciales de 1ºy 2ºorden continuas, posee una orbita homoclina Hasociada a un punto de equilibrio cuya matriz de linealizaci´on posee un autovalor real λy un par de autovalores complejos conjugados σ±ωi. Si se satisface la condici´on |λ|>|σ|,(5.1) entonces el sistema posee un r´egimen ca´otico en un entorno de la homoclina. La condici´on (5.1) se denomina condici´on de Shilnikov. Figura 5.1: Ejemplo gr´afico de una ´orbita homoclina. 1C.Tresser prob´o que este teorema aplica en sistemas lineales a trozos continuos[15] [14]. 31 CAP´ ITULO 5. ´ ORBITAS HOMOCLINAS Y TEOREMA DE SHILNIKOV 32 En la secci´on 6.2, hemos obtenidos para un valor fijado de β= 19 que α= 14,43746643008159. Utilizando la funci´on autmatAC (ver 7.4) para estos valores y un comando b´asico de MatLab, podremos verificar f´acilmente c´omo se cumple la condici´on 5.1 en nuestro sistema central, donde se encuentra el punto de equilibrio. >>abs ( lambdaC ) ans = 3.5181 >>abs( real ( autovalorcomplejoC )) ans = 1.2278 5.1. Puntos de tangencia En el plano de separaci´on ΣR, existe una recta donde las ´orbitas del sistema diferencial son tangentes al plano ΣR. Tomemos como Tal punto de intersecci´on de esa recta de tangencia con el plano focal de la zona central PF. Entendemos una homoclina directa como aquella que corta al plano de separaci´on RRy excatamente en dos puntos IRy ΠR R(IR) . Recordando las semiaplicaciones defendidas en el Cap´ıtulo 3 (3.6 y 3.7), tomemos S= ΠC R(T)−1. Existir´a una homoclina directa, si y solo si ΠR R(IR)∈ST. Para corroborar esto, partiendo de AC=   α 7α0 1−1−β 0−β0   ,(5.2) y teniendo en cuenta que nuestro punto T est´a contenido en el plano de separaci´on RR. Es f´acil concluir que la recta de tangencia RTviene dada por y=−x 7. Por ´ultimo, es obvio que la intersecci´on de RTcon PFser´a un punto, al que denominaremos T. T=      1 −1 7 λC α 7−λC +λC 7       . CAP´ ITULO 5. ´ ORBITAS HOMOCLINAS Y TEOREMA DE SHILNIKOV 33 Como hemos mencionado previamente, en la secci´on 6.2, hemos obtenido para un valor fijado de β= 19 un valor de α= 14,43746643008159. Haciendo uso de la funci´on puntodetangencia 7.21.1 para estos valores, obtenemos T=  1 0,1429 1,1082 . Por otro lado, utilizando la funci´on solucionzonaCT 7.21.2 obtenemos el punto S. S=  1 −1,2917 2,1297  . Conocidos estos puntos es f´acilmente demostrable que ΠR R(IR)∈ST. Visto que ambos requisitos, cuando β= 19 y α= 14,43746643008159 quedan satisfechos, podemos concluir que existe una homoclina tipo Shilnikov que induce al sistema a un comportamiento ca´otico para estos valores. CAP´ ITULO 6. EXISTENCIA DE ORBITAS HOMOCLINAS 40 obtener resultados confiables. Es importante resaltar que este trabajo no solo nos ha permitido profundizar en el campo de estudio elegido, sino que tambi´en hemos adquirido habilidades y conocimientos t´ecnicos en diversas ´areas, como el an´alisis num´erico, la programaci´on y la interpretaci´on de resultados. Adem´as de la redacci´on de texto y el uso de LaTex. Estas habilidades son valiosas y nos han preparado para futuras investigaciones y desaf´ıos acad´emicos. En resumen, este Trabajo de Fin de Grado ha sido una experiencia enriquecedora y gratificante, donde hemos aplicado y ampliado nuestros conocimientos en el campo de la ciencia y la tecnolog´ıa. Nos sentimos orgullosos de haber alcanzado los objetivos planteados y esperamos que este trabajo contribuya de manera significativa al avance y la comprensi´on en esta ´area de estudio. Cap´ıtulo 7 C´odigos en Matlab 7.1. Coeficientes zona derecha function sol = coefc1c2c3terminado (alpha ,beta ) %Calculo de los coef . c1 , c2 y c3 [~ ,vR ,u,w]= autmatrizAR (alpha ,beta ); %obtenemos autovalores zona R %lambdaC = autmatrizAC (alpha ,beta ); Obtenemos autovalores zona C (en desuso ya que se opto por calcular lambdaC dentro de " Interseccion ") M=[ vR u w]; % matriz con autovalores en forma de columna [I1 ]= Interseccion (alpha ,beta ); %calculamos la recta interseccion "I1" b=I1 -[3/2 0 -3/2] ';% [3/2 0 -3/2] punto de equilibrio sol=M\b; %solucionamos el sistema end 7.2. Matriz A zona derecha function [ lambdaR ,vR ,u ,w, autovalorcomplejoR ]= autmatrizAR ( alpha ,beta) %AUTMATRIZAR Autovalor real de matriz derecha AR en Chua %Construimos la matriz AR =[- alpha *(2/7) ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; %Calculo de los autovalores de la matriz de la parte R [V,D]= eig (AR); %V contiene columnas que son autovectores . El primero es el real y los %demas complejos %D contiene autovalores en la diagonal , fuera de esta son 0 lambdaR =D (1 ,1); % Guardo autovalor real autovalorcomplejoR =D (2 ,2) ; % Guardo autovalor complejo vR=V(: ,1) ; % guardo el primer autovector ( real) u= real(V (: ,2)); % guardo el segundo autovector ( complejo )[ parte real ] w= imag(V (: ,2)); % guardo el segundo autovector ( complejo )[ 41 CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 42 parte imaginaria ] %Forma alternativa %aut= sort(aut ,'ComparisonMethod ','real '); %lambdaR = aut (1) ; %vR= null (AR - lambdaR * eye (3) ); end 7.3. Coeficientes zona central function sol = coefc1c2c3central (alpha , beta ) %Calculo de los coef . c1 , c2 y c3 en la zona central [~,vC ,uC ,wC ]= autmatrizAC (alpha , beta ); M=[ vC uC wC]; %[p]= fsolve (@ sisparahomo ,[1.1915 , alpha ]); [p]= fsolve (@ sisparahomo ,[1.1871 , alpha ]); Estrella = solucionzonaR (p(1 ,1) ,p (1 ,2) ,beta ) % t12= tiempovueloenzonaR (12.165457244 ,19 ,1.4) ; % Estrella = solucionzonaR (t12 ,12.165457244219103 ,19) ; sol =M\ Estrella ; end 7.4. Matriz A zona central function [lambdaC ,vC ,uc ,wc , autovalorcomplejoC ]= autmatrizAC ( alpha ,beta) %AUTMATRIZAR Autovalor real de matriz central AC en Chua %Construimos la matriz AC =[ alpha /7 ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; %Calculo de los autovalores de la matriz de la parte C [V,D]= eig (AC); %V contiene columnas que son autovectores . El primero es el real y los %demas complejos %D contiene autovalores en la diagonal , fuera de esta son 0 %%almaceno los valores obtenidos lambdaC =D (1 ,1); autovalorcomplejoC =D (2 ,2) ; vC=V(: ,1) ; uc= real (V (: ,2)); wc= imag (V (: ,2)); end CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 43 7.5. Calculo de punto de interseccion %function [ I1 ]= Interseccion ( alpha ,beta , lambdaC ) function [I1 ]= Interseccion ( alpha , beta ) lambdaC = autmatrizAC ( alpha ,beta ); %calculo autovalor real zona C %Calculo del punto de interseccion de vC con el plano x=1 I1 =[1; -(( alpha /7) -lambdaC )/ alpha ;beta *(( alpha /7) -lambdaC )/( alpha * lambdaC )]; 7.6. Representacion gr´afica zona central function dibujosolzonaC = dibsolucionzonaCmejorado ( intervalo , alpha ,beta) % Dibuja la solucion en la zona C del sistema de Chua l= length ( intervalo ); % Obtenemos la longitud del intervalo % Prealocar espacio para las orbitas orbita1 = zeros ( size ( intervalo )); orbita2 = zeros ( size ( intervalo )); orbita3 = zeros ( size ( intervalo )); for i=1:l %Recorremos cada valor en el intervalo orbita = solucionzonaC ( intervalo (i),alpha , beta ); orbita1 (i)= orbita (1) ; orbita2 (i)= orbita (2) ; orbita3 (i)= orbita (3) ; end 7.7. Representacion gr´afica zona derecha function dibujosolzonaR = dibsolucionzonaRmejorado ( intervalo , alpha ,beta) % Dibuja la solucion en la zona R del sistema de Chua % Dado un intervalo y valores para los parametros alpha y beta , se obtiene % la solucion en la zona R del sistema de Chua para cada valor del % intervalo y se almacena en las variables orbita1 , orbita2 y orbita3. % Luego se dibuja la solucion en 3D utilizando la funcion plot3 . l= length ( intervalo ); % Obtenemos la longitud del intervalo % Prealocar espacio para las orbitas CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 44 orbita1 = zeros ( size ( intervalo )); orbita2 = zeros ( size ( intervalo )); orbita3 = zeros ( size ( intervalo )); for i=1:l % Recorremos cada valor en el intervalo orbita = solucionzonaR ( intervalo (i),alpha , beta ); % Obtenemos la solucion en la zona R orbita1 (i)= orbita (1) ; % Almacenamos la posicion x de la solucion orbita2 (i)= orbita (2) ; % Almacenamos la posicion y de la solucion orbita3 (i)= orbita (3) ; % Almacenamos la posicion z de la solucion end dibujosolzonaR =[ orbita1 ; orbita2 ; orbita3 ]; % Guardamos la solucion en una matriz plot3(orbita1,orbita2 ,orbita3 ,'LineWidth',2) % Graficamos la solucion en 3D 7.8. Representacion orbita homoclina function dibujoalfa ( alpha ,beta ,tini ,op) % lambdaC = autmatrizAC (alpha , beta ); [I1 ]= Interseccion (alpha ,beta ); if op == 1 plot3 ([0 , I1 (1) ] ,[0, I1 (2) ],[0 ,I1 (3)],'LineWidth',2); hold on; dibsolucionzonaRmejorado (0:0.01: tini +0.01 , alpha , beta ); dibsolucionzonaCmejorado (0:0.01:3.3 , alpha , beta ); hold off elseif op == -1 plot3 ([0 ,- I1 (1) ],[0,- I1 (2) ],[0,- I1 (3) ], 'LineWidth',2); hold on; dibsolucionzonaRsime (0:0.01: tini +0.01 , alpha ,beta ); dibsolucionzonaCsime (0:0.01:3.3 , alpha ,beta ); hold off else plot3 ([0 , I1 (1) ] ,[0, I1 (2) ],[0 ,I1 (3)],'LineWidth',2); hold on; dibsolucionzonaRmejorado (0:0.01: tini +0.01 , alpha , beta ); dibsolucionzonaCmejorado (0:0.01:3.3 , alpha , beta ); plot3 ([0 ,- I1 (1) ],[0,- I1 (2) ],[0,- I1 (3) ], 'LineWidth',2); dibsolucionzonaRsime (0:0.01: tini +0.01 , alpha ,beta ); dibsolucionzonaCsime (0:0.01:3.3 , alpha ,beta ); hold off end end CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 45 7.9. Solucion del sistema zona derecha function orbita = solucionzonaR (t,alpha , beta ) %Solucion del sistema en la zona R con c.i. I1. I1 es la recta %interseccion del plano con el limite de la zona R [lambdaR ,vR ,u ,w, autovalorcomplejoR ]= autmatrizAR ( alpha ,beta ) ;%obtengo autovalores % " Formulas " soluciones EDO : X1= exp ( lambdaR *t)* vR; sigmaR =real ( autovalorcomplejoR ); omegaR =imag ( autovalorcomplejoR ); X2= exp ( sigmaR *t) .*( cos ( omegaR *t)*usin ( omegaR *t)*w); X3= exp ( sigmaR *t) .*( sin ( omegaR *t)*u+ cos ( omegaR *t)*w); %Calculo coeficientes para satisfacer c.i: sol = coefc1c2c3terminado ( alpha ,beta ); c1= sol (1); c2= sol (2); c3= sol (3); %Obtengo solucion ( Solucion Homo + Solucion Nohomo ) orbita = c1 *X1+ c2 *X2+ c3 *X3 +[3/2;0; -3/2]; 7.10. Solucion del sistema zona izquierda function orbita = solucionzonaC (t,alpha , beta ) %Solucion del sistema en la zona C con c.i. Estrella . %Consultar " solucionzonaR .m para una mejor explicacion " [lambdaC ,vC ,uC ,wC , autovalorcomplejoC ]= autmatrizAC (alpha , beta ); X1= exp ( lambdaC *t)* vC; sigmaC =real ( autovalorcomplejoC ); omegaC =imag ( autovalorcomplejoC ); X2= exp ( sigmaC *t) .*( cos ( omegaC *t)*uC - sin ( omegaC *t)*wC ); X3= exp ( sigmaC *t) .*( sin ( omegaC *t)* uC +cos ( omegaC *t)*wC ); sol = coefc1c2c3central (alpha , beta); c1= sol (1); c2= sol (2); c3= sol (3); orbita = c1 *X1+ c2 *X2+ c3 *X3; CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 46 7.11. Sistema para ODE derecho function Y= sisdifzonaC (~,X, alpha , beta ) %Termino derecho del sistema diferencia en la zona R. Se usa para ODE45 . %Matriz A en la zona central AC =[ alpha /7 ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; Y=AC*X; %Y=AR*X +[3* alpha /7;0;0]; 7.12. Sistema para ODE central function Y= sisdifzonaR (~,X, alpha , beta ) %Termino derecho del sistema diferencia en la zona R. Se usa para ODE45 AR =[- alpha *(2/7) ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; Y=AR*X +[3* alpha /7, 0, 0] '; 7.13. Sistema para ODE function Y= sisdifzonasRyC (~,X, alpha ,beta ) %Termino derecho del sistema diferencial en la zona R if X(1) >=1 AR =[- alpha *(2/7) ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; Y=AR*X +[3* alpha /7, 0, 0] '; else AC =[ alpha /7 ,alpha ,0;1 , -1 ,1;0 , - beta ,0]; Y=AC*X; end end 7.14. ODE45 function [t ,solR ,solC , solC2 ]= Solv2 ( alpha ,beta , tiempo ) % %Zona Central ini % [~, solC ]= ode45 (@(t,X) sisdifzonaC (t,X, alpha , beta ) ,[0, tiempo ] ,0); % orbita1 = solC (: ,1) ; % orbita2 = solC (: ,2) ; % orbita3 = solC (: ,3) ; CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 47 % plot3(orbita1,orbita2 ,orbita3,'LineWidth ',2) % hold off [I1 ]= Interseccion (alpha ,beta ); plot3 ([0 , I1 (1) ] ,[0, I1 (2) ],[0 ,I1 (3)],'LineWidth',2); hold on %Zona derecha [~ , solR ]= ode45 (@(t ,X) sisdifzonaR (t ,X,alpha , beta ) ,[0, tiempo ],I1); orbita1R = solR (: ,1) ; orbita2R = solR (: ,2) ; orbita3R = solR (: ,3) ; plot3 (orbita1R , orbita2R , orbita3R ,'LineWidth',2) hold on %solR (end ,:) %Zona Central 2. tiempovueloenzonaC (14.43746643008159,19,1.1871) = 0.1151 %tc2= tiempovueloenzonaC (14.43746643008159 ,19 ,1.1871) ; Estrella =[1.0000 , -0.2734 , -1.4550] '; [~, solC2 ]= ode45 (@(t,X) sisdifzonaC (t,X,alpha , beta ) ,[0 ,2.8] , solR (end ,:) '); orbita1C = solC2 (: ,1); orbita2C = solC2 (: ,2); orbita3C = solC2 (: ,3); plot3 (orbita1C , orbita2C , orbita3C ,'LineWidth',2) hold off end 7.15. ODE45 con baja tolancia function [t ,solR ,solC , solC2 ]= Solv3 ( alpha ,beta , tiempo ) % %Zona Central ini % [~, solC ]= ode45 (@(t,X) sisdifzonaC (t,X, alpha , beta ) ,[0, tiempo ] ,0); % orbita1 = solC (: ,1) ; % orbita2 = solC (: ,2) ; % orbita3 = solC (: ,3) ; % plot3(orbita1,orbita2 ,orbita3,'LineWidth ',2) % hold off [I1 ]= Interseccion (alpha ,beta ); plot3 ([0 , I1 (1) ] ,[0, I1 (2) ],[0 ,I1 (3)],'LineWidth',2); pause hold on plot3 ([0 ,- I1 (1) ],[0,- I1 (2) ],[0,- I1 (3) ], 'LineWidth',2); CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 48 %Zona derecha y central options = odeset ( 'RelTol',1e-14 , 'AbsTol',1e -15) ; [~ , solR ]= ode78 (@(t ,X) sisdifzonasRyC (t,X, alpha , beta ) ,[0 ,4.63] , I1 , options ); orbita1R = solR (: ,1) ; orbita2R = solR (: ,2) ; orbita3R = solR (: ,3) ; plot3 (orbita1R , orbita2R , orbita3R ,'LineWidth',2) plot3 (- orbita1R ,- orbita2R ,-orbita3R ,'LineWidth',2) hold off %solR (end ,:) % %Zona Central 2. tiempovueloenzonaC (14.43746643008159,19,1.1871) = 0.1151 % %tc2= tiempovueloenzonaC (14.43746643008159 ,19 ,1.1871) ; % Estrella =[1.0000 , -0.2734 , -1.4550] '; % [~, solC2 ]= ode45 (@(t,X) sisdifzonaC (t,X,alpha , beta ) ,[0 ,2.8] , solR (end ,:) '); % orbita1C = solC2 (: ,1); % orbita2C = solC2 (: ,2); % orbita3C = solC2 (: ,3); % plot3 ( orbita1C ,orbita2C , orbita3C ,'LineWidth ',2) % hold off % end 7.16. ODE78 function [t ,sol ]= Solv (alpha , beta , tiempo ) Estrella =[1.0000,-0.17104641,-0.74666314]'; [t,sol ]= ode45 (@(t,X) sisdifzonasRyC (t,X,alpha , beta ) ,[0, tiempo ], Estrella ); orbita1 = sol (: ,1) ; orbita2 = sol (: ,2) ; orbita3 = sol (: ,3) ; plot3(orbita1,orbita2 ,orbita3 ,'LineWidth',2) end 7.17. C´alculo de tiempo. Zona central 1 function t1= tiempovueloenzonaC ( alpha ,beta , tiempoini ) %I1= Interseccion ( alpha ,beta ); Estrella =[1.0000,-0.17104641,-0.74666314]'; function sol11fin = primeracompsolzonaC (alpha ,beta , tiempo ) [~ , sol1 ]= ode45 (@(t ,X) sisdifzonaC (t ,X,alpha , beta ) ,[0, tiempo ], Estrella ); sol11 =sol1 (: ,1) -1; CAP´ ITULO 7. C ´ ODIGOS EN MATLAB 49 sol11fin = sol11 (end ); end t1= fzero (@( tiempo ) primeracompsolzonaC ( alpha ,beta , tiempo ) ,[0.0001 , tiempoini ]); end 7.18. C´alculo de tiempo. Zona central 2 function t1= tiempovueloenzonaC1 ( alpha , beta , tiempoini ) function sol11fin = primeracompsolzonaC (alpha ,beta , tiempo ) [~ , sol1 ]= ode45 (@(t ,X) sisdifzonaC (t ,X,alpha , beta ) ,[0, tiempo ],0); sol11 =sol1 (: ,1) -1; sol11fin = sol11 (end ); end t1= fzero (@( tiempo ) primeracompsolzonaC ( alpha ,beta , tiempo ) ,[0, tiempoini ]); end 7.19. C´alculo de tiempo. Zona derecha function t1= tiempovueloenzonaR ( alpha ,beta , tiempoini ) %I1= Interseccion ( alpha ,beta ); I1 =[0.999996373164299,-0.107934733707622,-0.439167178205539]'; function sol11fin = primeracompsolzonaR (alpha ,beta , tiempo ) [t, sol1 ]= ode45 (@(t,X) sisdifzonaR (t,X, alpha ,beta ) ,[0, tiempo ],I1); sol11 =sol1 (: ,1) -1; sol11fin = sol11 (end ); end t1= fzero (@( tiempo ) primeracompsolzonaR ( alpha ,beta , tiempo ) ,[0.0001 , tiempoini ]); end 7.20. Comprobaciones function [ myvR1 ]= myvR (alpha ,beta , lambdaR ,z) %Valor esperado de VR segun mis calculos myvR1 =[( alpha* lambdaR *z) /( beta *( - alpha *(1 -5/7) - lambdaR ));(- lambdaR *z)/ beta ;z]; end function x1= prisolucionzonaR (t,alpha , beta)