Continuación numérica de soluciones periódicas en sistemas tridimensionales lineales a trozos

Equation Chapter 1 Section 1 Trabajo Fin de Grado Ingeniería Aeroespacial Continuación Numérica de Soluciones Periódicas en Sistemas Tridimensionales Lineales a Trozos Autor: Antonio José Aguilera Albendín Tutor: Enrique Ponce Nuñez Dep. de Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 ii Trabajo Fin de Grado Ingeniería Aeroespacial Continuación Numérica de Soluciones Periódicas en Sistemas Tridimensionales Lineales a Trozos Autor: Antonio José Aguilera Albendín Tutor: Enrique Ponce Nuñez Catedrático Dep. de Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 iv Trabajo Fin de Grado: Continuación Numérica de Soluciones Periódicas en Sistemas Tridimensionales Lineales a Trozos Autor: Antonio José Aguilera Albendín Tutor: Enrique Ponce Nuñez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2017 vi El Secretario del Tribunal A mi familia A mis maestros viii Trabajo Fin de Grado Ingeniería Aeroespacial: Continuación Numérica de Soluciones Periódicas en Sistemas Tridimensionales Lineales a Trozos Antonio José Aguilera Albendín 10 de julio de 2017 8ÍNDICE DE FIGURAS Capítulo 1 Introducción Este proyecto surge con el interés de estudiar el comportamiento de un sistema dinámico no lineal y poder analizar su comportamiento en función de las diferentes variables que modelan dicho sistema. El proyecto sigue el trabajo realizado por Carlos Vallet Burguillos en su Proyecto de Fin de Grado “Análisis de la dinámica de un oscilador electrónico tipo Bonhoeffer-Van der Pol”, tutelado por Javier Ros Padilla. Para completar el análisis se partirá de una serie de valores fijos de los parámetros del oscilador y determinar mediante continuación numérica algunas de las curvas de bifurcación y ciclos límites del sistema. El método de continuación y las distintas simulaciones se llevarán a cabo mediante el software Mathematica, utilizando un código existente, realizado por los profesores E. García Medina y F. Fernández Sánchez. 9 10 CAPÍTULO 1. INTRODUCCIÓN Capítulo 2 Sistemas dinámicos 2.1. Introducción a los sistemas dinámicos Los sistemas que pueden ser modelados mediante ecuaciones algebraicas no lineales o ecuaciones diferenciales no lineales son aquellos que son denominados sistemas no lineales. Ejemplos de dichos sistemas están presentes en muchas disciplinas de la ingeniería y de la ciencia. Este serie de ecuaciones lineales dependen de ciertos parámetros que representan los posibles valores de las magnitudes que rigen la dinámica del sistema. Por tanto, la elección de los valores de estos parámetros determina el comportamiento del sistema y los posibles estados de equilibrio del mismo, así como su respuesta ante pequeñas perturbaciones que puedan producirse. Un sistema dinámico es aquel cuyo estado evoluciona con el tiempo t. Esta evolución está gobernada por una serie de reglas (no necesariamente ecuaciones) que especifica el estado del sistema tanto para valores continuos como discretos de t. Una evolución discreta en el tiempo es descrita normalmente mediante un sistema de ecuaciones algebraicas, mientras que una evolución continua en el tiempo es descrita mediante un sistema de ecuaciones diferenciales. El comportamiento asintótico de un sistema dinámico cuando t→ ∞ es llamado el estado estacionario del sistema. A menudo, este estado estacionario puede corresponder a un conjunto delimitado, el cual puede tener tanto una solución estática como dinámica. El comportamiento del sistema dinámico antes de alcanzar el estado estacionario es denominado estado transitorio, y consecuentemente, la correspondiente solución del sistema dinámico es lla- 11 12 CAPÍTULO 2. SISTEMAS DINÁMICOS mada solución transitoria. La solución de un sistema dinámico puede ser tanto constante en el tiempo como variante. Puntos fijos, soluciones de equilibrio y soluciones estacionarias son otros nombres para las soluciones constantes, mientras que las soluciones dinámicas es otro nombre para las soluciones variantes en el tiempo. 2.2. Ciclos límites En los sistemas no lineales existen soluciones periódicas que se caracterizan por tener estados invariantes en el tiempo. La solución periódica es una respuesta dinámica caracterizada por una frecuencia básica F. Matemáticamente, la solución x(t) es periódica para todo tiempo t, si existe algún T que satisfaga: x(t+T) = x(t) 0 < t < T. A una solución periódica se le llama ciclo límite si no existen otras soluciones periódicas cercanas a ella (es decir que se encuentra aislada). De forma similar a los puntos de equilibrio, los ciclos límites pueden ser estables o inestables. Un ciclo límite es asintóticamente estable si todas las órbitas, comenzando cerca de éste, alcanzan el ciclo límite conforme el tiempo incrementa. Un ejemplo claro es la figura 2.1 donde se ve que las órbitas convergen al ciclo límite estable. Figura 2.1: Ciclo límite estable. 2.3. ESTABILIDAD 13 Un ciclo límite es inestable si todas las órbitas comenzando cerca de éste, se alejan o divergen del mismo. Esto se muestra en la figura 2.2 donde se muestra como se alejan las órbitas. Figura 2.2: Ciclo límite inestable. 2.3. Estabilidad 2.3.1. Estabilidad estructural Nuestro sistema a estudiar está sujeto a pequeñas perturbaciones, lo cual hace necesario estudiar la estabilidad de la estructura orbital de dicho sistema dinámico. Por estructura orbital entendemos las diferentes órbitas del sistema considerado y su estabilidad para un conjunto de valores para los diferentes parámetros que definen dicho sistema. El Teorema de Hartman-Grobman establece que un campo vectorial es estructuralmente estable en la proximidad de puntos hiperbólicos fijos. Por otro lado, en la proximidad de un punto de bifurcación, la estructura orbital es inestable y el campo vectorial asociado es estructuralmente inestable. 14 CAPÍTULO 2. SISTEMAS DINÁMICOS 2.3.2. Estabilidad de bifurcaciones frente a perturbaciones Una bifurcación es estable frente a una perturbación si no hay diferencias cualitativas entre el diagrama de bifurcaciones en el sistema perturbado y en el sistema no perturbado 2.4. Bifurcaciones El concepto bifurcación es utilizado para indicar un cambio cualitativo en las características de un sistema. Estos cambios cualitativos están asociados al cambio en el número, tipo y estabilidad de los puntos de equilibrio u órbitas periódicas del sistema, bajo la variación de uno o más parámetros de los que depende el sistema. Las bifurcaciones se pueden considerar locales o globales. Se considera como locales aquellas en las que el cambio cualitativo ocurre en la proximidad de un punto fijo de una solución periódica del sistema. Cualquier otro cambio cualitativo es considerado como bifurcación global. En problemas de bifurcaciones es útil considerar un espacio formado por las variables de estado y los parámetros de control, llamado espacio estadocontrol. En este espacio, los puntos en los que se producen las bifurcaciones son llamados puntos de bifurcación. Muchas ramas de similares y/o diferentes soluciones se unen o emergen de puntos locales de bifurcación. Una bifurcación que requiere al menos m parámetros de control para que se produzca es llamada bifurcación de dimensión m. A menos que se indique lo contrario, debe entenderse que los parámetros de control son modificados muy lentamente para que sus valores instantáneos puedan ser considerados constantes. 2.4.1. Bifurcaciones locales Supongamos que uno o más parámetros de control son lentamente modificados, un punto fijo se convierte en no hiperbólico en un cierto punto del espacio estado-control. Entonces, si el cambio antes y después de este pun- to del espacio de estados es cualitativo, a este punto lo consideramos como punto de bifurcación. Un punto fijo puede perder estabilidad debido a una de las siguientes bifurcaciones: bifurcación silla-nodo, bifurcación pitchfork, bifurcación trans- 2.4. BIFURCACIONES 15 crítica y bifurcación Hopf. Las bifurcaciones silla-nodo, pitchfork y transcríticas son clasificadas como bifurcaciones estáticas, debido a que a ellas sólo convergen las soluciones asociadas a puntos de equilibrio, es decir, soluciones estáticas. En cambio, en la bifurcación Hopf convergen tanto las soluciones correspondientes a puntos fijos y soluciones periódicas, por lo que esta bifurcación es clasificada como bifurcación dinámica. Bifurcación Silla-Nodo Este tipo de bifurcación se produce debido a la aparición y desaparición de puntos de equilibrio conforme varía el parámetro de control. Por ejemplo, consideremos el siguiente sistema: ˙x =F(x;λ) = −λ−x2,(2.1) donde λes un escalar que hace de parámetro de control. Las soluciones de este sistema se muestran en la figura 2.3: Figura 2.3: Las gráficas corresponden a las soluciones del sistema para los casos λ<0, λ=0 y λ>0, respectivamente. En este caso tenemos un espacio bidimensional. Para λ<0 existen dos puntos de equilibrio en el sistema, uno estable y otro inestable. Cuando el parámetro λse anula los puntos estable e inestable chocan y se produce el punto de bifurcación Silla-Nodo, en el cual el Jacobiano del sistema se vuelve singular J(x∗) = 0, y no es posible obtener una solución de la ecuación - x2−λ= 0 . Por último, para el caso λ>0, no existen puntos de equilibrio. La respuesta del sistema es mostrada en la figura 2.4: en la que se ha re- 16 CAPÍTULO 2. SISTEMAS DINÁMICOS Figura 2.4: Diagrama bifurcación Silla-Nodo. presentado en trazo continuo los puntos estables y con línea discontinua los puntos inestables. Este tipo de bifurcaciones existe siempre que se cumplan las siguientes condiciones: f(x, λ)=0;∂f ∂xx∗ = 0;∂f ∂λ 6= 0;∂2f ∂x2x∗ = 0 La primera ecuación implica que la solución es un punto de equilibrio y la segunda que es un punto no hiperbólico. La tercera ecuación significa que coexiste en el punto de bifurcación (x∗, λ)una función continua diferenciable h(x)que intersecta a la línea x=x∗de manera transversal. Por último, la cuarta ecuación indica que el punto de equilibrio permanece de manera local a un lado de la línea λ=λ0. Estas dos últimas ecuaciones se les denomina condiciones de transversalidad. Bifurcación Pitchfork En este tipo de bifurcación, conforme se varía el parámetro de control un punto de equilibrio estable se convierte en inestable, y simultáneamente se originan dos puntos de equilibrios estables. De esta forma podemos decir que la bifurcación se asemeja a un “tridente”. Esta bifurcación puede ser de dos tipos: subcrítica y supercrítica. 2.4. BIFURCACIONES 17 La bifurcación supercrítica está relacionada con la creación de dos puntos de equilibrio estables a partir de de otro que se volvió inestable conforme varía el parámetro de control λ, tal y como se puede ver en la figura 2.5: Para el Figura 2.5: Creación de dos puntos fijos estables a partir de uno inestable. caso λ<0, sólo existe un punto de equilibrio estable. Cuando λse anula, este punto es críticamente estable debido a que la linealización desaparece, aún así sigue siendo estable. Para el caso que λ>0, el punto de equilibrio estable se convierte en inestable y se obtienen dos puntos estables simétricos. Un ejemplo de un sistema en el que se produzca este tipo de bifurcación es el siguiente: ˙x =F(x;λ) =λx−x3,(2.2) En la figura 2.6 podemos ver la respuesta del sistema en función del parámetro de control λ: en la que se ha representado en trazo continuo y grueso los puntos estables y en línea discontinua los puntos inestables. Se puede observar que para λ<0, el punto de equilibrio es inestable. Por el contrario, para λ>0, el punto de equilibrio se vuelve inestable y aparecen dos puntos estables de forma simétrica. Bifurcación Hopf Este tipo de bifurcación involucra a un punto de equilibrio y a un ciclo límite. La bifurcación se produce debido al cruce de dos autovalores complejos con el eje imaginario. Dicha condición la podemos ver en la figura 2.7: 24 CAPÍTULO 2. SISTEMAS DINÁMICOS Figura 2.15: Bifurcación Neimar-Sacker de Ciclos Límites. ambos ciclos) y subcrítica (se repelen ambos ciclos). Duplicación de Periodo o Flip Al contrario de la bifurcación tipo Silla-Nodo, donde sólo existen dos ciclos límites, estas bifurcaciones empiezan por tener un ciclo límite estable para un determinado valor inicial del parámetro de bifurcación dentro de un plano de fases. En el instante en el que se presenta una variación en el parámetro de control, las órbitas asociadas a las variables del sistema comienzan a dividirse en un par de órbitas, formándose una primera bifurcación de periodo doble. Pero si el parámetro sigue cambiando, en el sistema pueden ocurrir otras bifurcaciones de periodo doble. Por tanto, cada una de las órbitas anteriores se volverá a dividir en dos órbitas, originándose cuatro órbitas. Al seguir variándose el parámetro de bifurcación se pueden seguir formando más pares de órbitas en cascada dentro del sistema, conduciendo al sistema a un estado caótico. En la figura 2.17 aparece una bifurcación de duplicación de periodo típica. Para el valor crítico de bifurcación α=0 existe un ciclo límite estable, que corresponde a un punto fijo estable de la aplicación de Poincaré (nodo). Para α>0, el punto fijo pierde la estabilidad y se convierte en un punto de silla, correspondiendo a un ciclo límite inestable mientras aparece un ciclo límite estable de periodo doble. Para α<0 el ciclo límite desaparece. 2.4. BIFURCACIONES 25 Figura 2.16: Actractor extraño del Sistema de Lorenz. Figura 2.17: Bifurcación Flip de Ciclos Límites. 26 CAPÍTULO 2. SISTEMAS DINÁMICOS 2.5. Osciladores electrónicos 2.5.1. Osciladores y ciclos límites Los osciladores son circuitos electrónicos que pueden presentar oscilaciones de amplitud y periodos fijos sin excitación exterior. Estas oscilaciones se denominan ciclos límites u oscilaciones automantenidas. Un buen oscilador debe presentar al menos dos características para un correcto funcionamiento: 1) Una vez transcurrido el régimen transitorio inicial, los valores de amplitud y periodo de los ciclos límite no deben depender de las condiciones iniciales del sistema. 2) Las oscilaciones han de existir para unos rangos de los parámetros del sistema suficientemente amplios. Debido a que los sistemas lineales no presentan estas características, cualquier oscilador tiene que ser un sistema no lineal y ser modelado como tal para poder establecer analíticamente la amplitud y el periodo de las oscilaciones. La existencia de al menos un dispositivo no lineal en un circuito convierte a éste en no lineal. 2.5.2. Circuitos electrónicos no lineales Los circuitos electrónicos no lineales incluyen, además de algunos dispositivos lineales (resistencias, bobinas y condensadores), ciertos dispositivos no lineales usualmente basados en materiales semiconductores (diodos, transistores y amplificadores operaciones principalmente). Mientras que los dispositivos lineales no presentan ninguna dificultad y pueden modelarse mediante ecuaciones algebraicas o diferenciales lineales bien conocidas, para los dispositivos no lineales es necesario establecer unos modelos matemáticos adecuados que describan su funcionamiento de forma suficientemente aproximada. 2.5.3. Circuitos electrónicos lineales a trozos Una posibilidad muy utilizada es considerar que los dispositivos no lineales presentan varios modos de funcionamiento y establecer para cada uno de ellos un modelo lineal. Así se obtiene para cada dispositivo no lineal un modelo lineal a trozos. Cada modo de funcionamiento está asociado a un rango de las variables utilizadas para describir el comportamiento del circuito, o 2.5. OSCILADORES ELECTRÓNICOS 27 expresado de otra forma, a cada modo corresponde una región en un espacio de fase. De esta forma la resolución de un circuito no lineal puede hacerse resolviendo varios circuitos lineales, cada uno de los cuales tiene validez sobre periodos de tiempo que llamamos intervalos. Un intervalo finaliza y comienza el siguiente cuando las variables asociadas a los dispositivos no lineales pasan de un rango de valores a otro, o equivalentemente, cuando el punto representativo del estado del sistema en el espacio de fase pasa de una región a otra. En ese instante decimos que ha cambiado el modo de funcionamiento del circuito. 28 CAPÍTULO 2. SISTEMAS DINÁMICOS Capítulo 3 Análisis del sistema dinámico 3.1. Descripción del sistema El estudio de los ciclos límite de un sistema diferenciable es, después del análisis de soluciones en equilibrio, uno de los problemas más importantes en el análisis cualitativo de sistemas dinámicos. Diferentes herramientas provenientes de la teoría de bifurcación están disponibles para garantizar la bifurcación y la existencia de ciclos límite en sistemas lineales, sin embargo en este proyecto, trataremos sistemas a trozos lineales no diferenciables, debido a la aplicación a la vida real, en concreto, con osciladores electrónicos. Consideramos el sistema de la familia de sistemas lineales a trozos escritos en la forma de Luré: ˙ X=F(X) = ARX+bsat(x),(3.1) donde X= (x, y, z)T∈R3,y la función de saturación viene dada por: sat(x) =      1si x > 1, xsi |x|⩽1, −1si x < −1. La matriz ARtendrá la siguiente expresión: AR=   t−1 0 m0−1 d0 0    El vector bestá definido por: 29 30 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO b=   T−t M−m D−d  , de manera que, para |x| ≤ 1la dinámica está gobernada por el sistema: ˙ X=ACX=   T−1 0 M0−1 D0 0      x y z  ,(3.2) mientras que ARrepresenta la dinámica del sistema en las zonas externas y ACla dinámica del sistema en la zona central. El parámetro Trepresenta la traza, Mes la suma de los adjuntos de los elementos de la diagonal y D es el determinante de la matriz AC. Analógamente,trepresenta la traza, m la suma de los adjuntos de los elementos de la diagonal y del determinante de la matriz AR. Toda la información sobre la dinámica del sistema en la zona lineal está contenida en los autovalores λ1,λ2,λ3de ACque satisfacen el polinomio característico de la matriz. El polinomio característico de AC viene dado por: pAC(λ) = det(AC−λI) = λ3−Tλ2+Mλ −D. (3.3) Estamos interesados en la generación de oscilaciones automantenidas a partir de una situación estable. Cuando el sistema posee un autovalor real estable (negativo) y dos autovalores complejos conjugados que pasan de estables a inestables, es decir, la parte real de los autovalores complejos pasa de negativa a positiva, se producen oscilaciones periódicas cuya amplitud viene limitada por la saturación existente en el sistema. Equivalentemente, podemos decir que buscamos en que condiciones los dos autovalores complejos de la matriz ACcruzan el eje imaginario del plano complejo. Sabemos que el sistema se hará estable para: T < 0, D < 0, MT −D < 0.(3.4) Se tiene que AC=AR+beT 1, donde eT 1={1,0,0}, dando como resultado la familia de sistemas llamados en su forma generalizada de Liénard, el cuál, es análogo a la forma canónica desarrollada en teoría de control. Por tanto, bajo estas condiciones para cada sistema de la forma expuesta, haciendo un cambio de variables obtendríamos lo expuesto en (3.4). 3.2. ANÁLISIS DEL SISTEMA DINÁMICO 31 Considerando una estructura general de los autovalores para realizar el análisis dinámico, introducimos εcomo el parámetro de bifurcación tal que los tres autovalores de la matriz ACson λ(ε)para autovalores reales y σ(ε)±iω(ε) para autovalores con parte real y compleja. Los parámetros λ,σ,ωtienen las siguientes expresiones: λ(ε) = λ0+λ1ε+O(ε2), σ(ε) = σ1ε+O(ε2), ω(ε) = ω0+ω1ε+O(ε2), (3.5) donde se asumirá siempre ω0>0. Además, asumiremos que σ1yλ1no se anulan. Usando las expresiones anteriores resultan las siguientes expresiones para T,MyD: T(ε) = λ(ε)+2σ(ε), M(ε) = σ2(ε) + ω2(ε)+2λ(ε)σ(ε), D(ε) = λ(ε)[σ2(ε) + ω2(ε)], (3.6) donde tenemos que: T0=T(0) = λ0, M0=M(0) = ω2 0, D0=D(0) = λ0ω2 0. (3.7) Sabemos que ante las siguientes transformaciones no se verá afectado el sistema propuesto, conservándose invariante: (x, y, z, τ, t, m, d, ε)−→ (x, −y, z, −τ, −t, m, −d, −ε)(3.8) Esta propiedad simétrica es útil para simplificar el análisis de este tipo de sistemas. Después, sintetizaremos todas las afirmaciones fundamentales acerca del posible equilibrio de este tipo de sistema. 3.2. Análisis del sistema dinámico Asumimos en lo que sigue λ0= 0, es decir que al pasar εpor 0, tanto el autovalor real como la parte real del par complejo conjugado se anulan. Analizamos primero el número de equilibrios que puede presentar nuestro sistema. Si nos fijamos en la tercera ecuación del sistema, es decir, ˙z= dx + (D−d)sat(x), vemos que los equilibrios deben cumplir la condición 32 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO dx + (D−d)sat(x)=0. Si d=D6= 0, entonces la única solución es x= 0. Si D−d6= 0, podemos dividir por esta expresión y escribir: d d−Dx= sat(x)⇔1 1−D d x= sat(x), Si dyDtienen signos distintos (D/d < 0), la fracción 0<1 1−D d <1, signo de que hay tres soluciones, en concreto x= 0 yx=±(1 −D d), con |x|>1. Si dyDtienen el mismo signo (D/d > 0), 1−D d<1y puede ocurrir que: a) 0 <1−D d<1, b) 0 <1−D d= 0, c) 1 −D d<0. En los tres casos es fácil ver que sólo hay equilibrio si x= 0. Finalmente si D= 0 (que es lo que ocurre cuando ε= 0) queda la condición d[x−sat(x)] = 0, por tanto, todo valor de x∈[−1,1] serviría como solución. Para εpequeño en valor absoluto, se denota que el signo de Des igual al signo de λ1ε, por tanto, se desprenden las siguientes afirmaciones: (a) Si dλ1ε > 0el único punto de equilibrio es el origen. (b) Si dλ1ε < 0el equilibrio está en el origen y los dos puntos X+ Eq =1−D(ε) d, T(ε)−tD(ε) d, M(ε)−mD(ε) dT , X− Eq =−X+ Eq. (c) Si ε= 0 todos los puntos del siguiente segmento son equilibrios del sistema. {(x, y, z)T∈R3: (x, y, z)T=µ(1,0, ω2 0)T,|µ|⩽1} Si ademas d6= 0, el segmento contendría todos los puntos del equilibrio. A partir de la conclusión c) vemos que para ε= 0 se produce lo que se denomina bifurcación pitchfork degenerada, que se observa en la figura 3.1: 3.2. ANÁLISIS DEL SISTEMA DINÁMICO 33 Figura 3.1: El diagrama de bifurcación del equilibrio del sistema dibujado cuando d > 0. Sólo se ha representado la coordenada xdel equilibrio de puntos. Nótese que para dλ1ε < 0los puntos X± Eq son puntos convergentes del campo de vectores correspondientes a |x|>1pero estos puntos están fuera de sus correspondientes zonas. Estos puntos no constituyen un equilibrio real, aunque organizan la dinámica en las regiones externas. Estos puntos de equilibrio son comúnmente llamados puntos de equilibrio virtual. Para estudiar la existencia de órbitas periódicas en este tipo de sistemas, comenzaremos considerando la zona central del sistema con ε= 0. En este caso, una computación directa muestra que la solución de nuestro sistema a estudiar partiendo de un punto arbitrario tal como (x0, y0, z0)Tpara τ= 0 viene dado por la solución: ˙ X(τ) = e(ACτ)   x0 y0 z0   =1 ω2 0   ω2 0cosω0τ−ω0sinω0τ1−cosτ ω3 0sinω0τ ω2 0cosω0τ−ω0sinω0τ 0 0 ω2 0      x0 y0 z0   . (3.9) Si |x0|<1la correspondiente órbita es un arco de la siguiente elipse: 40 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO 3.3.2. Órbitas periódicas trizonales Asumimos la existencia de una órbita periódica existente intersectando los planos frontera x=−1yx= 1 en cuatro puntos, cuyos puntos tienen la siguiente expresión: X0=   1 y0 z0   ,X1=   1 y1 z1   ,X2=   −1 y2 z2   ,X3=   −1 y3 z3   . Como se muestra en la figura 3.4, se representan los puntos X0,X1,X2 yX3del ciclo límite trizonal y sus tiempos de vuelo en la parte central (τC), en la parte izquierda (τL) y en la parte derecha (τR). Figura 3.4: Ciclo límite trizonal en los planos frontera x=−1yx= 1 La simetría de la órbita periódica trizonal del sistema implica que X2=−X0y X3=−X1, 3.3. ANÁLISIS DE ÓRBITAS PERIÓDICAS MEDIANTE ECUACIONES DE CIERRE41 por tanto, el estudio de la trayectoria entre X0yX2es suficiente para determinar por completo la órbita. En la zona central tenemos que ˙x =ACx. Tomando como punto inicial X0tenemos como solución en la zona derecha x(τ) = XR+eARτ(X0−XR).(3.19) El punto X1de la trayectoria donde τ=τRsatisface que X1=XR+eARτR(X0−XR),(3.20) proporcionando la primera ecuación de cierre del ciclo límite trizonal. En la parte central tenemos que ˙x =ACx. Tomando como punto inicial X1 tenemos como solución x(τ) = eACτX1,(3.21) Para el punto X2donde τ=τCtenemos que X2=eACτCX1.(3.22) Sabemos que X2=−X0, por tanto, sustituyendo en 3.22 tenemos que X0=−eACτCX1.(3.23) La ecuación 3.23 proporciona la segunda ecuación de cierre para el ciclo límite trizonal, por consiguiente, tenemos el sistema no lineal (X1=XR+eARτR(X0−XR), X0=−eACτCX1, donde se tienen dos ecuaciones vectoriales de tres componentes teniendo como incógnitas τC,τR,y0,z0,y1,z1. Realizando una serie de operaciones reducimos el sistema a uno con 4 ecuaciones escalares y 4 incógnitas, las cuales son: y0,z0,τC,τR. El sistema reducido resultante a resolver queda de la forma:          eT 1e−ACτCX0+ 1 = 0 eT 1(XR+eARτR(X0−XR)) −1=0 eT 2(XR+eARτR(X0−XR+e−ACτCX0+)) = 0 eT 3(XR+eARτR(X0−XR+e−ACτCX0+)) = 0 (3.24) siendo los vectores eT 1={1,0,0},eT 2={0,1,0},eT 3={0,0,1}. 42 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO Figura 3.5: Interpretación gráfica de la continuación parámetrica natural. 3.4. Continuación Consideramos el cálculo de familias uniparamétricas de equilibrios de x0(t) = f(x(t), µ), es decir, soluciones de: f(x, µ) = 0. A un conjunto de estas soluciones se le denomina rama de soluciones. Con µ∈RyX=X(x, µ), la ecuación se puede reescribir como: f(X)=0, donde f:Rn+1 →Rn. Una solución X0= (x0, µ0)de f(X) = 0 se considera regular si fx(x0, µ0)tiene rango máximo, es decir, si rango(fx(x0, µ0)) = n. Se verifica que cerca de una solución regular no existe sino una única rama de soluciones. 3.4.1. Continuación paramétrica natural Se toma µcomo parámetro de continuación y suponemos que tenemos una solución x0de f(x, µ) = 0 en µ0, así como sus derivadas ˙x0respecto al parámetro µ, y que queremos calcular la solución x1en µ1=µ0+4µ(ver figura 3.5 para interpretación gráfica). 3.4. CONTINUACIÓN 43 Para hallar x1, resolvemos f(x1, µ1)=0para x1por el Método de Newton: fx(x(k) 1, µ1)4x(k) 1=−f(x(k) 1, µ1) x(k+1) 1=x(k) 1+4x(k) 1 k= 0,1,2, ... (3.25) Siendo x(0) 1=x0+4µ˙x0. Si fx(x1, µ1)es no singular y 4µsuficientemen- te pequeÃśo, la Teoría del Método de Newton garantiza que este proceso iterativo converge. Despue’s de la convergencia, el nuevo vector ˙x1se puede obtener resolviendo: fx(x1, µ1) ˙x1=−fµ(x1, µ1)(3.26) Esta ecuación resulta de derivar fx(x(µ), µ) = 0 con respecto a µen µ=µ1. 3.4.2. Método de pseudo longitud de arco de Keller Este método fue desarrollado por Keller (1977, 1987) y es utilizado como método de continuación en el software AUTO desarrollado por Doedel (1986). Permite la continuación de cualquier solución regular y es considerado geométricamente como el método de continuación más natural. Supongamos que tenemos una solución (x0, µ0)de f(x, µ)=0, así como del vector normalizado ( ˙x0,˙µ0)de la rama de soluciones en (x0, µ0). El método consiste en resolver las siguientes ecuaciones para x1yµ1: f(x, µ) = 0 (x1−x0) ˙x0+ (µ1−µ0) ˙µ0−4s= 0 (3.27) Resolviendo por Método de Newton: (f1 x)(k)(f1 µ)(k) ˙x0˙µ0! 4x(k) 1 4µ(k) 1!=− f(x(k) 1, µ(k) 1) (x(k) 1−x0) ˙x0−(µ(k) 1−µ0) ˙µ0−4s! (3.28) Determinamos el vector de dirección resolviendo: f1 xf1 µ ˙x0˙µ0! ˙x1 ˙µ1!= 0 1!(3.29) Posteriormente se normaliza dicho vector. La orientación se conserva si 4s es suficientemente pequeño. El nuevo vector de dirección se reescala para que sea unitario: ˙x12+ ˙µ12= 1. En la práctica, el paso 4s se actualiza durante el cálculo de la rama de soluciones. En el caso más simple, la elección del nuevo paso 4s se basa en la convergencia del Método de Newton. 44 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO Figura 3.6: Interpretación gráfica del Método de Keller. 3.4.3. Algoritmo utilizado en el Método de pseudo longitud de arco de Keller A partir de las series proporcionadas en los trabajos [2] y [3] obtenemos aproximaciones de x0yµ0tales que f(x0, µ0)=0. A partir de estos resultados, aplicamos el Método de Newton para hallar valores fiables de x0yµ0. Una vez que tenemos el punto X0= (x0, µ0), calculamos el nuevo punto X01 a partir de una pequeÃśa variación en la coordenada µ→µ+hy volvemos a aplicar el Método de Newton para hallar dicho nuevo punto X01. A partir de estos dos puntos, X0yX01, calculamos el vector V0mediante la siguiente expresión: V0=X01 −X0 h(3.30) y posteriormente se normaliza V0=V0 |V0|(3.31) 3.4. CONTINUACIÓN 45 Figura 3.7: Valores de x0yµ0obtenido por el Método de Newton. Este primer vector no es necesario para el método de continuación, pero sí para el cálculo del vector dirección necesario para el funcionamiento del algoritmo. Figura 3.8: Cálculo de V0. 46 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO El vector al que llamamos V0se corresponde en el apartado anterior con ˙ X0= ( ˙x0,˙µ0), que se añade como fila extra al Jacobiano para el cálculo del nuevo vector de dirección V= ( ˙x1,˙µ1)en las ecuaciones descritas en método de pseudo longitud de arco. f1 xf1 µ ˙x0˙µ0! ˙x1 ˙µ1!= 0 1!(3.32) Predicción En este paso definimos el nuevo punto cercano a la curva, al cual llamaremos “predicción”, de manera que si llamamos hal paso de integración se tiene Predicción =X0+h·V(3.33) También debemos calcular la matriz Jacobiana en dicho punto, siendo su última fila el vector ~ V, que corresponde a la primera matriz de la siguiente ecuación: (f1 x)(k)(f1 µ)(k) ˙x0˙µ0! 4x(k) 1 4µ(k) 1!=− f(x(k) 1, µ(k) 1) (x(k) 1−x0) ˙x0−(µ(k) 1−µ0) ˙µ0−4s! (3.34) Figura 3.9: Cálculo del vector “predicción”. 3.4. CONTINUACIÓN 47 Correción El objetivo de este paso es calcular la “corrección” necesaria para cada iteración del método. Para ello es necesario calcular previamente el incremento : 4x(k) 1 4µ(k) 1! Correccion=Prediccion-Incremento (3.35) A partir del paso de continuación 4s y de los valores x(k) 1yµ(k) 1hallamos el valor de (x(k) 1−x0) ˙x0−(µ(k) 1−µ0) ˙µ0−4s. La resolución del sistema de ecuaciones lineales (3.34) nos proporciona el incremento (4x(k) 1,4µ(k) 1)necesario para el posterior cálculo de la “corrección”. Esta operación corresponde a la resolución por el Método de Newton, acercándose más a la curva de bifurcación tras cada iteración. La “predicción” es actualizada en cada iteración, salvo en el punto 1 de la curva de bifurcación (ver figura 3.10). Las Figura 3.10: Cálculo de X1tras el proceso iterativo. predicciones ˙x1y˙µ1se determinan resolviendo el sistema (3.32), que nos proporciona el nuevo vector de dirección. A partir de este momento, comienza el proceso iterativo “corrección-predicción” recorriendo la curva de bifurcación. Como se puede observar en la figura (3.10), la predicción utiliza la dirección a la tangente a la curva, mientras que la corrección se lleva a cabo sobre el hiperplano ortogonal a dicha tangente. 48 CAPÍTULO 3. ANÁLISIS DEL SISTEMA DINÁMICO Capítulo 4 Oscilador electrónico BVP 4.1. Obtención del sistema modelado del oscilador electrónico BVP Para la aplicación de los resultados teóricos de los capítulos anteriores, consideramos el oscilador electrónico de Bonhoeffer-Van der Pol, que consta de dos condensadores (C1yC2de igual capacitancia C), una inductancia L, una resistencia ry una conductancia no lineal g, con una región activa, como se muestra en la figura 4.1: Figura 4.1: Esquema electrónico del oscilador BVP La dinámica del oscilador viene definida por las siguientes ecuaciones diferenciales: 49 56 CAPÍTULO 4. OSCILADOR ELECTRÓNICO BVP De las expresiónes t=d < 0ymt −d < 0obtenemos las condiciones que ha de cumplir el circuito: (α < γ (α−γ)[1 −αγ]<0(4.17) De la segunda expresión de (4.17) sabemos que (α−γ)es negativo debido a la primera condición, por tanto, el término [1 −αγ]debe ser positivo. Así obtenemos las condiciones que debe cumplir el circuito para que el equilibrio en las zonas estables sea estable: α < γ < 1 α(4.18) En la figura (4.4) se han representado ambas condiciones: 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 Figura 4.4: Condiciones que ha de cumplir el circuito para que el equilibrio en las zonas externas sea estable. El área de color azul del plano (α−γ) corresponde a la zona de estabilidad. 4.2.4. Zona correspondiente a tres equilibrios En este apartado se va a estudiar la zona del plano (α−γ)en el que coexisten tres equilibrios. Para ello se debe cumplir la condición det(AC)det(AR)< 4.2. ESTUDIOS DE EQUILIBRIOS (ESTABILIDAD) 57 0, por tanto dD < 0(4.19) Sustituyendo en (4.19) los valores de los invariantes del circuito, obtenemos la condición para que coexistan los tres equilibrios: (α−γ)(α+β−γ)<0(4.20) Para que esta inecuación se cumpla tenemos dos posibles opciones: 1. (α−γ)<0;(α+β−γ)>0 2. (α−γ)>0;(α+β−γ)<0. Este caso no es posible debido a que β > 0y por tanto, (α−γ)>0y(α+β−γ)<0son incompatibles. Finalmente obtenemos las condiciones que debe cumplir el circuito para que coexistan tres equilibrios: α < γ < α +β(4.21) En la figura (4.5) se han representado ambas condiciones: 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 Figura 4.5: Condiciones que ha de cumplir el circuito para que coexistan tres equilibrios. El área de color azul del plano (α−γ)corresponde a la zona de existencia. 58 CAPÍTULO 4. OSCILADOR ELECTRÓNICO BVP 4.3. Estudio de bifurcaciones En este apartado estudiaremos las diferentes bifurcaciones que se producen en el sistema en el plano (α−γ). Para ello se trabajará con los matrices ACyARy las diferentes propiedades que deben cumplir para cada tipo de bifurcación. 4.3.1. Bifurcación Foco-Centro-Ciclo-Límite La bifurcación Foco-Centro-Ciclo-Límite (FCCL) asociada al origen se origina cuando el par de autovalores complejos conjugados σ±iw cruzan el eje imaginario para el valor crítico de bifurcación, aparte de un autovalor real λ. Para que se produzca el cruce de los autovalores complejos por el eje imaginario, la matriz ACdebe cumplir la condición: T=D/M, con M > 0 (ver Teorema 1.1 de [5])(4.22) Las condiciones (4.22) en función de los parámetros del circuito da lugar a las siguientes expresiones: (γ=1 α+β γ < 2 α+β (4.23) La segunda condición de se cumple trivialmente cuando se cumple la primera. Además, el coeficiente de no degeneración δ=D(M−m) + M(d−Mt)(ver Teorema 1.1 de [5]), que en función de los parámetros del circuito queda δ=−βγ(α+β−γ)+(α−γ)(−1 + γ(α+β))(2γ(α+β)) (4.24) es distinto de cero sobre toda la curva de bifurcación γ=1 α+βexcepto en el punto (α= 1 −β, γ = 1). Dicho Teorema 1.1 estable que el Ciclo Límite existe cuando δ(T−D M)>0, condición que se verifica para los valores de γ por encima de la curva γ=1 α+β, para α6= 1 −β. En la gráfica (4.6) se han representado los punto en los que se produce el Ciclo Límite. 4.3. ESTUDIO DE BIFURCACIONES 59 Figura 4.6: Curva de bifurcación FCCL y ciclos límites existentes en sus cercanías. Aquellas zonas marcadas con “NCL” corresponde con aquellas en las que no se produce Ciclo Límite, mientras que aquellas marcadas con “CL” corresponden con las que sí se produce Ciclo Límite. 4.3.2. Bifurcación Hopf del infinito La bifurcación Hopf del infinito se produce cuando un par de autovalores complejos conjugados σ±iw , de la matriz ARde las zonas externas, cruzan el eje imaginario del plano complejo. Para que se produzca dicho cruce, la matriz ARdebe cumplir la condición: t=d/m, con m<0 (ver Teorema 5.10 de [8]) (4.25) Las condiciones (4.25) en función de los parámetros del circuito da lugar a las siguientes expresiones: (γ=1 α γ < 2 α (4.26) La segunda condición se cumple trivialmente cuando se cumple la primera. Además, el coeficiente de no degeneración δ=d(m−M) + m(D−mT)(ver Teorema 1 de [7]), que en función de los parámetros del circuito queda γ=βγ(α−γ)+(α+β−γ)(−1 + αγ)(2 −αγ)(4.27) 60 CAPÍTULO 4. OSCILADOR ELECTRÓNICO BVP es distinto de cero sobre toda la curva de bifurcación γ=1 αexcepto en el punto (α= 1, γ = 1). Dicho Teorema 1 estable que el Ciclo Límite existe cuando δ(t−d m)>0, condición que se verifica para los valores de γpor debajo de la curva γ=1 α, para α6= 1. En la gráfica (4.7) se han representado los punto en los que se produce el Ciclo Límite. Figura 4.7: Curva de bifurcación Hopf del infinito y ciclos límites existentes en sus cercanías. Aquellas zonas marcadas con “NCL” corresponde con aquellas en las que no se produce Ciclo Límite, mientras que aquellas marcadas con “CL” corresponden con las que sí se produce Ciclo Límite. 4.3. ESTUDIO DE BIFURCACIONES 61 4.3.3. Hopf Zero del Origen Para este caso estudiamos la matriz ACcorrespondiente a la zona central. Para poder realizar el estudio de esta bifurcación, tomaremos la siguiente estructura para los autovalores de AC(siguiendo la hipótesis de [4]): λ=− para el autovalor real, y ρ ±ωi para el par complejo imaginario, donde  será un parámetro de bifurcación. Así, los invariantes T,MyDde la matriz ACson de la forma:      T() = (2ρ−1) M() = ρ2(ρ−2) + ω2 D() = −(ρ2+ω2) Cuando el parámetro de bifurcación = 0, se produce el cruce por el eje imaginario del plano complejo de los autovalores. Esto implica que para el valor crítico de bifurcación los invariantes T,MyDsean igual a cero, dando lugar a (γ=α+β γ < 2 α+β (4.28) Combinando ambas ecuaciones de (4.28), obtenemos: γ=α+β < √2. Figura 4.8: Puntos correspondientes a la bifurcación Hopf Zero del Origen. 62 CAPÍTULO 4. OSCILADOR ELECTRÓNICO BVP 4.3.4. Hopf Zero del Infinito Para este caso estudiamos la matriz ARcorrespondiente a la zona externa. Por analogía al caso anterior (Hopf zero del origen), los invariantes t,myd de ARtienen que ser igual a cero, dando lugar a (γ=α γ < 2 α (4.29) Combinando ambas ecuaciones de (4.29), obtenemos: γ=α < √2. Figura 4.9: Puntos correspondientes a la bifurcación Hopf Zero del Infinito. 4.3.5. Triple zero del origen Este tipo de bifurcación se produce cuando los tres autovalores de la matriz ACcorrespodiente a la zona central son nulos para el valor crítico de bifurcación. Para ello, los invariantes T,MyDde ACdeben ser igual a cero, dando lugar a      γ=α+β γ=2 α+β γ=α+β (4.30) El anterior sistema (4.30) se trata en realidad de un sistema de 2 ecuaciones y 3 incógnitas, debido a que la primera y tercera ecuación son la misma. Por 4.3. ESTUDIO DE BIFURCACIONES 63 tanto, expresaremos αyγen función del parámetro β. Así, obtenemos: α=√2−β;γ=√2(4.31) Por lo tanto, el punto triple zero del origen se produce en las coordenadas (√2−β, √2) del plano (α−γ). 4.3.6. Triple zero de los equilibrios de la zona externa Para este caso estudiamos la matriz ARcorrespondiente a la zona externa. Por analogía al caso anterior (Triple zero del origen), los invariantes t,my dde ARdeben ser igual a cero, dando lugar a      γ=α γ=2 α γ=α (4.32) Igualando las expresiones de (4.32), obtenemos la solución α=√2,γ=√2. Finalmente representamos todas las curvas obtenidas anteriormente para así determinar el diagrama de bifurcaciones, el cual será de utilidad para el estudio de las curvas Silla-Nodo del sistema. 64 CAPÍTULO 4. OSCILADOR ELECTRÓNICO BVP Figura 4.10: Diagrama de bifurcaciones. 4.4. Análisis de bifurcaciones A partir de los resultados obtenidos en los apartados desarrollados anteriormente (4.2 y 4.3), analizaremos con más detalle las bifurcaciones que se producen en el oscilador electrónico BVP estudiado. Considerando el sistema expresado en la forma de Liénard: ˙x =   α−γ−1 0 2−αγ 0−1 α−γ0 0   x+   β −βγ β  sat(x), Asumiendo que los parámetros del circuito, α,βyγ, toman una serie determinada de valores, se obtienen los siguientes casos: 1. Con α > 0,β > 0y con γ0=1 α+βfijo: a) Para γ=γ0, se produce una bifurcación Foco-Centro-Ciclo-Límite (FCCL). 4.4. ANÁLISIS DE BIFURCACIONES 65 b) Para γ > γ0y suficientemente pequeño, aparece un ciclo límite. Si se cumple que α+β < 1, y por tanto, γ0>1, el ciclo límite es asintóticamente estable. 2. Con α > 0,β > 0,α+β6= 1 y0< γ0=α+β < √2fijo: a) Para γ > γ0, el origen es el único equilibrio del sistema. Además, si se cumple γγ0<1, éste es asintóticamente estable. b) Para γ= 0, el sistema experimenta una bifurcación propia de los sistemas diferenciales lineales a trozos, análoga a la bifurcación Hopf Zero. c) Para γ < γ0y suficientemente pequeño en valor absoluto, aparecen tres ciclos límite simultáneamente, uno de ellos trizonal y los otros dos bizonales, estando acompañados de dos puntos de equlibrio adicionales. 1) Si se cumple que γ0<1, el ciclo límite producido por la bifurcación trizonal es estable, mientras que el ciclo límite producido por la bifurcación bizonal es inestable. Los puntos de equilibrio que bifurcan serán estables siempre y cuando que γ0<1. 2) Si se cumple que 1< γ0<√2, el ciclo límite producido por la bifurcación trizonal es inestable, mientras que el ciclo límite producido por la bifurcación bizonal es estable. Los puntos de equilibrio que bifurcan serán estables siempre y cuando que γ0<1/α. En la figura (4.11) se han representado las zonas con los diferentes puntos de equilibrio del sistema y los distintos ciclos límites que se producen en el sistema. 72 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 γ 0.5 1.0 1.5 2.0 2.5 3.0 Figura 5.1: Evolución del tiempo de vuelo central τCpara α=0.9 en función de γ. El punto de Silla-Nodo (α, γ)=(0.9,0.972) es producido debido a la bifurcación Hopf Zero del origen. 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 γ 0.5 1.0 1.5 2.0 2.5 3.0 τr Figura 5.2: Evolución del tiempo de vuelo externo τRpara α=0.9 en función de γ. El punto de Silla-Nodo (α, γ)=(0.9,0.972) es producido debido a la bifurcación Hopf Zero del origen. En la figura (5.3) se han representado los ciclos límites obtenidos para diversos valores de los tiempos de vuelo τCyτR: 5.1. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITES TRIZONALES73 Figura 5.3: Ciclos Límites obtenidos para α=0.9 para diversos valores de τC yτR. Los ciclos límites atraviesan las tres zonas del sistema. A partir de los puntos de Silla-Nodo obtenidos producidos por la bifurcación Hopf Zero, añadimos la ecuación de Silla-Nodo a las ecuaciones de cierre del sistema y representamos la curva de continuación en el plano (α−γ). 2. Volvemos a repetir el mismo proceso que en el caso anterior, estudiando el ciclo límite correspondiente a (α, γ)=(0.9,1.028). En las figuras (5.5) y (5.6) se ha representado la evolución de los tiempos de vuelo correspondientes a α=0.9 en función de γ: 74 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP 0.85 0.90 0.95 1.00 α 0.975 0.980 0.985 0.990 0.995 1.000 Figura 5.4: Continuación numérica del ciclo límite trizonal producido por la bifurcación Hopf Zero a partir del punto Silla-Nodo (α, γ)=(0.9,0.972). 1.04 1.06 1.08 1.10 γ 0.5 1.0 1.5 2.0 2.5 3.0 3.5 τc Figura 5.5: Evolución del tiempo de vuelo central τCpara α=0.9 en función de γ. El punto de Silla-Nodo (α, γ)=(0.9,1.028) es producido debido a la bifurcación Hopf Zero del origen. En las figuras (5.7) se han representado los ciclos límites obtenidos para diversos valores de los tiempos de vuelo τCyτR: 5.1. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITES TRIZONALES75 1.04 1.06 1.08 1.10 γ 0.5 1.0 1.5 2.0 2.5 3.0 τr Figura 5.6: Evolución del tiempo de vuelo externo τRpara α=0.9 en función de γ. El punto de Silla-Nodo (α, γ)=(0.9,1.028) es producido debido a la bifurcación Hopf Zero del origen. Figura 5.7: Ciclos Límites obtenidos para α=0.9 para diversos valores de τC τR. Los ciclos límites atraviesan las tres zonas del sistema. A partir de los puntos de Silla-Nodo obtenidos producidos por la bifurcación Hopf Zero, añadimos la ecuación de Silla-Nodo a las ecuaciones de cierre del sistema y representamos la curva de continuación en el 76 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP plano (α−γ). 0.85 0.90 0.95 1.00 α 1.000 1.005 1.010 1.015 1.020 1.025 1.030 Figura 5.8: Continuación numérica del ciclo límite trizonal producido por la bifurcación Hopf Zero a partir del punto Silla-Nodo (α, γ)=(0.9,1.028). 3. El último ciclo límite estudiado es partiendo de un punto cercano a la degeneración triple zero. Dicho punto es (α, γ)=(√2−β, √2). En las figuras (5.9) y (5.10) se ha representado la evolución de los tiempos de vuelo para a α=√2−βen función de γ: 5.1. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITES TRIZONALES77 1.25 1.30 1.35 1.40γ 20 40 60 80 Figura 5.9: Evolución del tiempo de vuelo central τCpara α=√2−βen función de γ. Los puntos de Silla-Nodo presentes son producidos debido a la degeneración Triple Zero del Origen. 1.25 1.30 1.35 1.40γ 4.5 5.0 5.5 6.0 6.5 Figura 5.10: Evolución del tiempo de vuelo externo τRpara α=√2−βen función de γ. Los puntos de Silla-Nodo presentes son producidos debido a la degeneración Triple Zero del Origen. Las curvas obtenidas de la evolución de los tiempos de vuelo central τCy vuelo externo τRson típicas de una bifurcación homoclina tipo “Shilnikov” de un sistema tridimensional. Estas órbitas son trayectorias que parten de un punto de equilibrio X0Silla-Nodo con autovalores λ1(X0)>0> Reλ2,3(X0)>0y una bifurcación homoclina. Asumiendo 78 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP que se cumplen las condiciones: λ1(X0) + Reλ2,3(X0)>0λ2(X0)6=λ3(X0) Entonces, el sistema tiene un número infinito de ciclos límites en las proximidades del punto de equilibrio X0perteneciente a dicha órbita homoclina. Las órbitas nacen del punto de equilibrio y vuelven a dicho punto después de un tiempo infinito (Teorema 6.6 de [Kuznetsov]). Estas órbitas homoclinas están conectadas a la dinámica caótica de sistemas de diferentes áreas, como caos químico, fenómeno del ruido inducido y osciladores electroquímicos. En la figura (5.11) se han representado las órbitas producidas para el tiempo de vuelo central τCyτRdebido a la degeneración Triple Zero del Origen. Figura 5.11: Ciclos Límites obtenidos para α=√2−βpara diversos valores de τCτR. Los ciclos límites atraviesan las tres zonas del sistema. A partir de los puntos de Silla-Nodo obtenidos producidos por la bifurcación Hopf Zero, añadimos la ecuación de Silla-Nodo a las ecuaciones de cierre del sistema y representamos la curva de continuación en el plano (α−γ). 5.2. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITE BIZONALES79 1.25 1.30 1.35 1.40 α 1.390 1.395 1.400 1.405 1.410 Figura 5.12: Continuación numérica del ciclo límite trizonal producido por la degeneración Triple Zero a partir del punto Silla-Nodo (α, γ)=(√2− β,1.4037). 5.2. Continuación numérica de ciclos límite bizonales 5.2.1. Curvas Silla-Nodo correspondientes a ciclos límite bizonales Para el cálculo de estas curvas empleamos el algoritmo de detección de curvas Silla-Nodo explicado en el comienzo del capítulo. En este caso, al tratarse de ciclos l’mite bizonales, el apartado c)y < α +β−γdel algoritmo se tiene que cumplir. La zona estudiada corresponde al rombo del plano (α−γ)comprendido entre α∈[0.8,1] y γ∈[0.9,1.1]. En esta zona se continuará un ciclo límite de dos zonas correspondiente a la bifurcación Hopf Zero que se produce sobre γ=α+β. Manteniendo constante el valor de α, variamos el valor de γhasta encontrar un punto de bifurcación Silla-Nodo. Dicho punto es debido a un cambio del valor de los autovalores del producto M en el que uno de ellos cruza el círculo unidad por la parte positiva del eje real. Para el estudio de dichos ciclos límite, se analizarán tres casos distintos dentro de la zona en cuestión: 1. Para α=0.85 y γligeramente por debajo de γ=α+β. En las figuras 80 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP (5.13) y (5.14) se ha representado la evolución de los tiempos de vuelo correspondientes a α=0.85 en función de γ: 1.035 1.040 1.045 1.050 γ 4.0 4.5 5.0 5.5 6.0 6.5 Figura 5.13: Evolución del tiempo de vuelo central τCpara α= 0,85 en función de γ. El punto de Silla-Nodo (α, γ)=(0.85,1.033) es producido debido a la bifurcación Hopf Zero del origen. 1.035 1.040 1.045 1.050 γ 0.5 1.0 1.5 2.0 2.5 Figura 5.14: Evolución del tiempo de vuelo externo τRpara α= 0,85 en función de γ. El punto de Silla-Nodo (α, γ)=(0.85,1.033) es producido debido a la bifurcación Hopf Zero del origen. 5.2. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITE BIZONALES81 Como se ha comentado anteriormente, al tratarse de un ciclos límite bizonales, se tiene que cumplir la condición c)y<1.05-γ, condición que numéricamente se ha comprobado que deja de cumplirse cuando γ=1.05. Así, teniendo en cuenta esta restricción, se obtiene que el tiempo de vuelo central τCno llega a anularse, tal y como se puede observar en la figura (5.13). En la figura (5.15) se han representado los ciclos límite obtenidos para diversos valores de los tiempos de vuelo τCyτR: Figura 5.15: Ciclos Límites obtenidos para α=0.85 para diversos valores de τCyτR. Los ciclos límite atraviesan las dos zonas del sistema. 2. Para α=0.9 y γligeramente por encima de γ=1 α+β. En las figuras (5.16) y (5.17) se ha representado la evolución de los tiempos de vuelo correspondientes a α=0.9 en función de γ: 88 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP Figura 5.25: Ciclos Límites obtenidos para α=0.9 para diversos valores de τC yτR. Los ciclos límite atraviesan las dos zonas del sistema. A partir de los puntos de Silla-Nodo obtenidos producidos por la bifurcación Hopf Zero, añadimos la ecuación de Silla-Nodo a las ecuaciones de cierre del sistema y representamos la curva de continuación en el plano (α−γ). 0.85 0.90 0.95 1.00 α 0.96 0.97 0.98 0.99 1.00 γ Figura 5.26: Continuación numérica del ciclo límite bizonal producido por la bifurcación Hopf Zero. 5.2. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITE BIZONALES89 5.2.2. Curva Duplicación de Periodo correspondiente a ciclos límite bizonales Para estudiar este tipo de curva de bifurcación utilizamos el mismo procedimiento empleado en las curvas de bifurcación Silla-Nodo. En este caso, la zona a estudiar corresponde a la vertical de α=1.2 y con γligeramente por debajo de la curva γ=α+β. En esta zona se continuará un ciclo límite de dos zonas correspondiente a la bifurcación Hopf Zero que se produce sobre γ=α+β. Manteniendo constante el valor de α, variamos el valor de γhasta detectar un punto de bifurcación. Dicho punto es debido a un cambio del valor de los autovalores del producto M en el que uno de ellos cruza el círculo unidad por la parte negativa del eje real. En las figuras (5.27) y (5.28) se ha representado la evolución de los tiempos de vuelo correspondientes a α=1.2 en función de γ: 1.25 1.30 1.35 1.40 γ 5 10 15 20 25 30 τ Figura 5.27: Evolución del tiempo de vuelo central τCpara α=1.2 en función de γ. 90 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP 1.25 1.30 1.35 1.40 γ 2 4 6 8 Figura 5.28: Evolución del tiempo de vuelo externo τRpara α=1.2 en función de γ. En la figura (5.29) se han representado los ciclos límites obtenidos para diversos valores de los tiempos de vuelo τCyτR: Figura 5.29: Ciclos Límite obtenidos para α=1.2 para diversos valores de τC yτR. Los ciclos límites atraviesan las dos zonas del sistema. 5.2. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITE BIZONALES91 Añadiendo la ecuación de la curva Duplicación de Periodo a las ecuaciones del cierre del sistema y representamos la curva de continuación en el plano (α−γ). 1.15 1.20 1.25 1.30 1.35 1.40 α 1.30 1.35 1.40 1.45 Figura 5.30: Continuación numérica del ciclo límite bizonal producido por la curva Duplicación de Periodo. 5.2.3. Curva Bifurcación Toro correspondiente a ciclos límite bizonales Para este tipo de curva de bifurcación la zona a estudiar corresponde a la vertical de α=1.1 y con γligeramente por debajo de la curva γ=α+β. En esta zona se continuará un ciclo límite de dos zonas correspondiente a la bifurcación Hopf Zero que se produce sobre γ=α+β. Manteniendo constante el valor de α, variamos el valor de γhasta detectar un punto de bifurcación. Dicho punto es debido a un cambio del valor de los autovalores del producto M en el que dos de ellos cruzan el círculo unidad. En las figuras (5.31) y (5.32) se ha representado la evolución de los tiempos de vuelo correspondientes a α=1.1 en función de γ: 92 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP 1.15 1.20 1.25 1.30 γ 2 4 6 8 10 Figura 5.31: Evolución del tiempo de vuelo central τCpara α=1.1 en función de γ. 1.15 1.20 1.25 1.30 γ 1 2 3 4 5 6 7 Figura 5.32: Evolución del tiempo de vuelo externo τRpara α=1.1 en función de γ. En la figura (5.33) se han representado los ciclos límite obtenidos para diversos valores de los tiempos de vuelo τCyτR: 5.2. CONTINUACIÓN NUMÉRICA DE CICLOS LÍMITE BIZONALES93 Figura 5.33: Ciclos Límite obtenidos para α=1.1 para diversos valores de τC yτR. Los ciclos límite atraviesan las dos zonas del sistema. Añadiendo la ecuación de la curva Bifurcación Toro a las ecuaciones del cierre del sistema y representamos la curva de continuación en el plano (α−γ). 0.95 1.00 1.05 1.10 1.15 α 1.05 1.10 1.15 1.20 1.25 Figura 5.34: Continuación numérica del ciclo límite bizonal producido por la curva de Bifurcación Toro. 94 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP 5.3. Análisis de las curvas de continuación numérica obtenidas En este apartado se van a analizar las curvas de continuación numérica de ciclos límite trizonales y bizonales. En la figura (5.35) se encuentran representadas todas las curvas obtenidas mediante continuación numérica en el plano (α−γ). 0.6 0.8 1.0 1.2 1.4 α 0.6 0.8 1.0 1.2 1.4 Figura 5.35: Diagrama de bifurcaciones ampliado. Las curvas Silla-Nodo de color rojo, azul y verde corresponden a la continuación numérica de ciclos límites trizonales, y las curvas Silla-Nodo de color naranja y amarillo a ciclos límite bizonales. La curva de color morado corresponde a la curva de Duplicación de Periodo y la curva de color negro a la curva de Bifurcación Toro, ambas de ciclos límite bizonales. 5.3. ANÁLISIS DE LAS CURVAS DE CONTINUACIÓN NUMÉRICA OBTENIDAS95 En la figura (5.36) se ha indicado el número de cicles límites presentes en la zona del plano (α−γ)para α∈[0.80,1] y γ∈[0.90,1.10]: 2 CL 2 CL 1 CL 1 CL 0 CL 3 CL 3 CL 1 CL 1 CL 0 CL 0 CL 0 CL 0 CL 0 CL 0 CL 0.8 0.9 1.0 1.1 α0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 Figura 5.36: Diagrama de bifurcaciones ampliado. En la zona del plano que aparece la etiqueta “2CL” se produce tanto un ciclo límite bizonal como un ciclo límite trizonal. Estos ciclos límite se deben a la bifurcación Hopf Zero del origen o bifurcación FCCL y a la bifurcación Hopf del infinito o la bifurcación Hopf Zero del infinito. En cambio, la etiqueta “1CL” corresponde a un ciclo límite trizonal, el cual está producido por la bifurcación Hopf Zero del origen o bifurcación FCCL. 96 CAPÍTULO 5. CONTINUACIÓN NUMÉRICA EN OSCILADOR BVP En la figura (5.37) se ha indicado el número de cicles límites presentes en la zona del plano (α−γ)para α∈[1,1.50] y γ∈[1.20,1.55]: En la zona del 0 CL 0 CL 3 CL 2 CL 1.1 1.2 1.3 1.4 α1.1 1.2 1.3 1.4 1.5 Figura 5.37: Diagrama de bifurcaciones ampliado. plano que aparece la etiqueta “2CL” se produce tanto un ciclo límite bizonal como un ciclo límite trizonal. Estos ciclos límite se deben a la bifurcación Hopf Zero y a la degeneración Triple Zero. Finalmente, hay que aclarar que el análisis realizado durante este trabajo no comprende todas las posibles curvas de bifurcación posibles, por lo que es importante no trabajar fuera de las zonas no estudiadas del plano (α−γ) para poder evitar posibles errores. Bibliografía [1] Ali H. Nayfeh, Balakumar Balachandran. Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods. [2] Yuri A. Kuznetsov. Elements of Applied Bifurcation Theory,Second Edition. Springer. [3] Enrique Ponce, Javier Ros, Elísabet Vela. A Unified Approach to Piecewise Linear Hopf and Hopf-Pitchfork Bifurcations. Springer, Chapter 12. 2014. [4] Enrique Ponce, Javier Ros, Elísabet Vela. Unfolding the fold-Hopf bifurcation in piecewise linear continuous differential systems with symmetry. 2013. [5] Emilio Freire, Enrique Ponce, Javier Ros. The focus-center-limit cycle bifurcation in symmetric 3D piecewise linear systems. 2005. [6] Y. Nishiuchi, T. Ueta, H. Kawakami. Stable torus and its bifurcation phenomena in a simple three-dimensional autonomous circuit 2006. [7] V. Carmona, E. Freire, E. Ponce, J. Ros, F. Torres. Limit cycle bifurcation in 3D continuous piecewise linear systems with two zones. Application to Chua‘s circuit. 2005. [8] Jabier Ros Padilla. Estudio del comportamiento dinámico de sistemas autónomos tridimensionales lineales a trozos 2003. [9] Carlos Vallet Burguillos. Análisis de la dinámica de un oscilador electrónico tipo Bonhoeffer-Van der Pol 2014. [10] Luis Gómez de la Herrán Ramírez. Continuación numérica de ciclos límites en osciladores electrónicos. 2016. 97