Órbitas Canard en sistemas lineales a trozos y aplicaciones

Trabajo de Fin de Grado: Órbitas Canard en Sistemas Lineales a Trozos y Aplicaciones Alejandro Bandera Moreno Doble Grado en Física y Matemáticas Tutelado por: Soledad Fernández García Macarena Gómez Mármol Junio de 2019 DPTO. ECUACIONES DIFERENCIALES Y ANÁLISIS NUMÉRICO - US Resumen El estudio de los sistemas dinámicos continúa teniendo una gran importancia en el mundo actual, ya que estos sistemas modelan múltiples fenómenos que evolucionan con el tiempo. Dichos fenómenos provienen de campos tan diversos como la Medicina, la Física o la Economía. Una clase de estos sistemas especialemente interesante es la de los sistemas con múltiples escalas de tiempo, en los que aparecen fenómenos sumamente complejos. En este trabajo se realiza un estudio sobre órbitas canard, presentes en sistemas dinámicos con varias escalas de tiempo, desde el punto de vista de la modelización lineal a trozos. Veremos que los sistemas lineales a trozos son capaces de reproducir fenómenos canard clásicos, reduciendo la complejidad de los modelos originales. Además, comprobaremos los resultados teóricos presentados en un modelo aerodinámico de pérdida de estabilidad de los aviones realizando un giro en la tierra, que sabemos que exhibe el fenómeno canard. Abstract Nowadays, the study of dynamical systems continues having a major importance, because of the great quantity of phenomena that evolves with time that they model. These phenomena come from fields so diverse as Medicine, Physics or Economy. A class of this systems specially interesting are systems with multiple time scales, in which extremely complex phenomena appear. In this work we will consider a study of canard orbits, appearing in dynamical systems with multiple time scales, from the point of view of piecewise linear modeling. We will see that piecewise linear models are able to reproduce classic canard phenomena, reducing the complexity of original models. In addition, we will check the theoretical results of an aerodynamic model of loss of lateral stability of an aircraft turning on the ground, which we know that exhibits the canard phenomenon. Índice general Introducción 3 1. Resultados Preliminares 5 1.1. Sistemas de Ecuaciones Diferenciales Ordinarias . . . . . . . . . . . . 5 1.1.1. Existencia y Unicidad de Solución . . . . . . . . . . . . . . . . 5 1.1.2. Órbitas y Estabilidad . . . . . . . . . . . . . . . . . . . . . . . 9 1.2. Teorema de la Función Implícita . . . . . . . . . . . . . . . . . . . . . 15 1.3. Bifurcaciones ............................... 16 1.3.1. Bifurcaciones en Dimensión Uno . . . . . . . . . . . . . . . . . 16 1.3.2. Bifurcación de Poincaré-Andronov-Hopf . . . . . . . . . . . . 19 2. Sistemas diferenciables Lento-Rápido 25 2.1. Sistemas Dinámicos con Múltiples Escalas de Tiempo . . . . . . . . . 25 2.1.1. La Variedad Crítica y el Flujo Lento . . . . . . . . . . . . . . 26 2.1.2. Explosión Canard . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2. La Ecuación de Van der Pol . . . . . . . . . . . . . . . . . . . . . . . 31 2.2.1. La Variedad Crítica y el Flujo Lento . . . . . . . . . . . . . . 32 2.2.2. Explosión Canard . . . . . . . . . . . . . . . . . . . . . . . . . 34 3. Aproximación Lineal a Trozos 37 3.1. Existencia y Unicidad de Solución . . . . . . . . . . . . . . . . . . . . 37 3.2. PrimerosModelos............................. 40 3.2.1. Sistema Lineal a Trozos con una Esquina . . . . . . . . . . . . 40 3.2.2. Sistema Lineal a Trozos con dos Esquinas . . . . . . . . . . . 41 3.3. Sistemas Lineales a Trozos con Tres Zonas . . . . . . . . . . . . . . . 42 3.3.1. Primeras Propiedades . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.2. Resultados principales . . . . . . . . . . . . . . . . . . . . . . 46 3.4. Similitudes y Diferencias con el Caso Diferenciable . . . . . . . . . . . 57 3.4.1. Primeros Modelos . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.2. Sistemas Lineales a Trozos con Tres Zonas . . . . . . . . . . . 58 3.5. Conclusiones................................ 60 4. Estudio de un ejemplo concreto: Ciclos Canard en la Dinámica de Aeronaves en Tierra 63 4.1. Presentación del Caso . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.1.1. Introducción............................ 64 4.1.2. Modelo de la Aeronave . . . . . . . . . . . . . . . . . . . . . . 64 4.2. Resultados del Caso Diferenciable . . . . . . . . . . . . . . . . . . . . 66 4.2.1. Investigación Uniparamétrica . . . . . . . . . . . . . . . . . . 66 4.2.2. Identificación de Ciclos Canard . . . . . . . . . . . . . . . . . 67 4.2.3. Modelo Plano de la Explosión Canard . . . . . . . . . . . . . 69 4.3. Construcción del Modelo Lineal a Trozos . . . . . . . . . . . . . . . . 70 4.3.1. Estudio de la Variedad Crítica del Sistema Diferenciable . . . 71 4.3.2. Modelo Lineal a Trozos con Tres Zonas . . . . . . . . . . . . . 72 4.3.3. Modelo Lineal a Trozos con Cuatro Zonas . . . . . . . . . . . 75 4.4. Simulación y Validación . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.4.1. Caso Reversible con Pendiente uno . . . . . . . . . . . . . . . 79 4.4.2. Caso Reversible con Pendiente distinta de uno . . . . . . . . . 84 Bibliografía 95 2 Introducción El estudio de los sistemas dinámicos continúa teniendo una gran importancia en el mundo actual. Esto se debe, entre otros factores, a que modelan múltiples fenómenos que evolucionan con el tiempo. Dichos fenómenos provienen de campos muy diversos, como pueden ser la Neurociencia, la Física o la Economía. Dentro de los sistemas dinámicos, el estudio de los sistemas con múltiples escalas de tiempo goza de gran interés. En los sistemas con múltiples escalas de tiempo o sistemas lento-rápido, aparecen fenómenos sumamente complejos que requieren del conocimiento de métodos de la Teoría Geométrica de las Perturbaciones Singulares para su estudio. En los últimos años se ha comprobado que los sistemas lineales a trozos son útiles para reproducir el comportamiento de los sistemas lento-rápido de una manera más sencilla. No obstante, la teoría clásica de ecuaciones diferenciales no se puede aplicar a estos sistemas, por lo que es necesaria la construcción de una nueva teoría que nos permita estudiar los sistemas lineales a trozos de forma sencilla. Uno de los fenómenos más interesantes que aparecen en algunos sistemas con múltiples escalas de tiempo es el fenómeno canard, objeto en el que centraremos nuestro estudio durante este trabajo. La característica principal asociada al fenómeno canard es la explosión canard, que se manifiesta como una transición muy rápida, con variaciones pequeñas de los parámetros, desde pequeños ciclos provenientes de una bifurcación de Hopf a oscilaciones de relajación, ciclos límite de gran amplitud formados por segmentos rápidos y lentos alternantes. Esta transición ocurre a través de una familia de ciclos límite cuya propiedad característica es que contienen segmentos canard, segmentos de la órbita que siguen la variedad lenta repulsiva. Este tipo de fenómenos se encuentra, por ejemplo, en modelos de la actividad neuronal, como el de FitzHugh-Nagumo [11]. Para el análisis de los modelos que vamos a estudiar en este trabajo es necesario recordar algunos resultados fundamentales de la teoría de ecuaciones diferenciales, como son la existencia y unicidad de solución de un sistema diferencial [22], y algunos resultados sobre órbitas y su estabilidad [21], [23]; recordaremos, a su vez, el Teorema de la Función Implícita [20] y enunciaremos algunas herramientas de la teoría de bifurcaciones en una y dos dimensiones, centrándonos en la Bifurcación de Poincaré- 3 Andronov-Hopf [12] y [19]. Todo esto se situará en el Capítulo 1. En el Capítulo 2, vamos a presentar el fenómeno de la explosión canard en los sistemas diferenciables lento-rápido, tomando como ejemplo el sistema de Van der Pol, analizando las condiciones necesarias para que el sistema exhiba este tipo de comportamiento que vamos a estudiar [4], [10]. El Capítulo 3 está dedicado íntegramente, siguiendo con el sistema de Van der Pol como ejemplo, a presentar una aproximacion lineal a trozos de los sistemas originales. En primer lugar, presentaremos un resultado sobre existencia y unicidad de solución en este tipo de sistemas [9], y continuaremos estudiando los primeros modelos que se propusieron para modelar el comportamiento de la explosión canard [2] y [13]. Veremos que el modelo estudiado en [2] reproduce de forma fiel ciertas caracterísiticas de las órbitas canard, pero aún no cuenta con todas las propiedades del fenómeno clásico. Así, será necesario añadir una zona más de linealidad para poder reproducir el fenómeno en su totalidad. Esto es lo que se realiza en [9], como se presenta al final de este capítulo. Realizaremos también una aplicación a un modelo físico en el Capítulo 4. Presentaremos un modelo sobre la pérdida de estabilidad lateral de un avión cuando realiza un giro mientras está posado en tierra [17], en el que se da el fenómeno de la explosión canard, para posteriormente construir un modelo lineal a trozos asociado y estudiar si efectivamente se puede reproducir la dinámica observada en el modelo original. Para terminar el capítulo nos centraremos en visualizar numéricamente ("cazar") una de las órbitas canard que están presentes en nuestro modelo lineal a trozos. 4 Capítulo 1 Resultados Preliminares En este capítulo veremos algunos resultados importantes que son necesarios para el desarrollo de nuestro trabajo. Algunos de ellos se han presentado en asignaturas del Grado en Matemáticas (como son EDO [22], AED [21] y DFVV [20]), mientras que otros resultados no se han tratado durante el Doble Grado en Física y Matemáticas (Teoría de Bifurcaciones [12] y [19]). 1.1. Sistemas de Ecuaciones Diferenciales Ordinarias Para lograr un correcto entendimiento de los resultados que presentamos en este trabajo, es necesario recordar algunos conceptos de cursos anteriores (como son el de función globalmente y localmente lipchitziana,solución e intervalo maximal, entre otros). Para el desarrollo de esta sección se han seguido principalmente las referencias [21] y [22]. 1.1.1. Existencia y Unicidad de Solución Funciones Lipchitzianas Consideremos fijado un conjunto no vacío, Ω⊂Rn+1, cuyos puntos denotamos por (t, x)con t∈Rex∈Rn. Definición 1 Se dice que f: Ω →Rnes una función globalmente lipschitziana respecto de la variable xen Ω, si existe una constante L > 0(dependiente de f), tal que: |f(t, x1)−f(t, x2)| ≤ L|x1−x2|,∀(t, x1),(t, x2)∈Ω. 5 Estabilidad Evidentemente, en la estimación del estado inicial de un sistema pueden producirse errores, por lo que resulta fundamenteal conocer cómo afectará este pequeño error a la verdadera solución del problema de Cauchy (PC). En esta sección empezaremos dando nociones sobre la estabilidad de los puntos de equilibrio para acabar proporcionando un método para conocer el tipo topológico de los mismos. Para la elaboración de esta sección se ha seguido la referencia [12]. Definición 14 Un punto de equilibrio x0del sistema (1.1) se dice estable si dado ε > 0, existe un δ > 0tal que, para cualquier xcon kx−x0k< δ, la solución φ(t, x) del sistema (1.1) que pasa por xen t= 0 satisface la desigualdad kφ(t, x)−x0k< ε, para todo t≥0. El punto de equilibrio x0se dirá inestable si no es estable. Definición 15 Un punto de equilibrio x0se dice que es asintóticamente estable si es estable y, además, si existe r > 0tal que kφ(t, x)−x0k → 0cuando t→+∞ para todo xsatisfaciendo kx−x0k< r. En el siguiente teorema, estudiaremos la estabilidad de los puntos de equilibrios de los sistemas lineales, siendo A la matriz de coeficientes del sistema. Teorema 5 Si todos los autovalores de la matriz de coeficientes Aen el sistema lineal ˙x=Ax tienen parte real negativa, entonces el punto de equilibrio x0= 0 es asintóticamente estable. Además, existen constantes positivas Kyαtales que keAtxk ≤ Ke−αtkxk,∀t≥0, x ∈Rn. Si uno de los autovalores de la matriz de coeficientes A tiene parte real positiva, entonces el punto de equilibrio x0= 0 es inestable.  El tipo de estabilidad de un punto de equilibrio es una propiedad local. Por ello, es de esperar que bajo ciertas condiciones el tipo de estabilidad de x0se pueda aproximar el campo vectorial fpor sus derivadas, que se trata de un campo vectorial lineal. Suponemos que fde (1.1) es C1(Ω) y sea Df(x) =            ∂f1 ∂x1 (x)∂f1 ∂x2 (x)··· ∂f1 ∂xn (x) ∂f2 ∂x1 (x)∂f2 ∂x2 (x)··· ∂f2 ∂xn (x) . . .. . ..... . . ∂fn ∂x1 (x)∂fn ∂x2 (x)··· ∂fn ∂xn (x)            (1.4) la matriz Jacobiana de fen el punto x. 12 Definición 16 Si x0es un punto de equilibrio de ˙x=f(x), entonces la ecuación diferencial lineal ˙x=Df(x0)x se denomina ecuación lineal variacional o la linealización del campo vectorial fen el punto de equilibrio x0. De igual forma que para los sistemas lineales, se tienen dos resultados sobre la estabilidad de la solución de la ecuación diferencial (1.1) en función del signo de la parte real de los autovalores de Df(x0): Teorema 6 Sea f∈C1(Ω). Si todos los autovalores de la matriz Jacobiana Df(x0)tienen parte real negativa, entonces el punto de equilibrio x0de la ecuación diferencial (1.1) es asintóticamente estable.  Teorema 7 Sea f∈C1(Ω). Si al menos uno de los autovalores de la matriz Jacobiana Df(x0)tienen parte real positiva, entonces el punto de equilibrio x0de la ecuación diferencial (1.1) es inestable.  En la siguiente definición y el resultado posterior vamos a destacar aquellos puntos de equilibrio en los que se espera reproducir la dinámica no lineal mediante la linealización. Definición 17 Un punto de equilibrio x0del sistema (1.1) se dice hiperbólico si todos los autovalores de la matriz Jacobiana Df(x0)tienen parte real no nula. En estos puntos de equilibrio la linealización nos da toda la información, como se recoge en el siguiente teorema: Teorema 8 (Grobman-Hartman) Si x0es un punto de equilibrio hiperbólico de (1.1), entonces existe un entorno de x0en el cual fes topológicamente equivalente al campo vectorial ˙x=Df(x0)x. Restringiéndonos a dimensión dos, vamos a presentar una clasificación del tipo topológico de los equilibrios en función de su matriz de coeficientes, esta clasificación servirá tanto para sistemas lineales como sistemas linealizados en el caso en el que el equilibrio sea hiperbólico. Además, también nos dará una expresión de los autovalores de la matriz de coeficientes. Observación 18 Se considera el sistema lineal ˙x=Ax con x∈R2, y sean δ= det(A)yτ=traza(A). Entonces el punto de equilibrio x0= 0 será: Un punto de silla si δ < 0. 13 Un nodo estable si δ > 0,τ2−4δ≥0yτ < 0. Un nodo inestable si δ > 0,τ2−4δ≥0yτ > 0. Un foco estable si δ > 0,τ2−4δ < 0yτ < 0. Un foco inestable si δ > 0,τ2−4δ < 0yτ > 0. Un centro si δ > 0,τ= 0. Además, los valores propios de Ase pueden escribir de la forma: λ=τ±√τ2−4δ 2.  De esta forma podremos saber, calculando el determinante y la traza de la matriz de coeficientes, cuál es el tipo topológico de un equilibrio y calcular los autovalores de la matriz de coeficientes. Para terminar esta sección, realizaremos un estudio de las órbitas cercanas a una órbita periódica y utilizaremos la aplicación de Poincaré para determinar su estabilidad. Sea φ(t, p)una solución periódica con periodo mínimo Tde la ecuación diferencial ˙x=f(x)y denotemos la órbita periódica correspondiente por Γ. Escogemos ahora un vector v∈R2tal que, vy el vector tangente f(p)de Γen psean linealmente independientes. De esta forma, podemos definir un segmento: Lε={x∈R|x=p+av, 0≤ |a| ≤ ε}, que se denomina sección transversal de la órbita periódica Γen el punto p. A continuación, se define una aplicación en un subconjunto de Lεinducida por el flujo. Se elige un εtan pequeño que Lεintersecta a la curva Γen un solo punto p, y que todas las órbitas que crucen Lεlo hagan en la misma dirección. Como φ(T, p) = py sus soluciones dependen de forma continua del valor inicial, hay un δ > 0tal que, si x0∈Lδ, entonces hay un primer instante T(x0)>0en el que φ(T(x0), x0)∈Lε. Definición 19 La aplicación de Poincaré en un entorno de una órbita periódiva Γ se define como: Π : Lδ−→ Lε x07−→ φ(T(x0), x0) 14 Los puntos de la sección transversal Lεtienen un orden natural: dos puntos x0=p+a0vyx1=p+a1vcumplen que x0≥x1en Lδ, si y solo si, a0≥a1. Por tanto, una aplicación de Poincaré Πse dice que es monótona si x0≥x1en Lδ implica que Π(x0)≥Π(x1). Teorema 9 La aplicación de Poincaré tiene las siguientes propiedades: 1. La aplicación de Poincaré Πcerca de la órbita periódica Γes una aplicación monótona de C1. 2. La órbita γ(x0)de un punto x0∈Lδes una órbita periódica si y solo si es punto fijo de la aplicación de Poincaré, es decir, Π(x0) = x0. 3. La órbita periódica Γ, con p∈Γ, es asintóticamente estable si Π0(p)<1, e inestable si Π0(p)>1. 1.2. Teorema de la Función Implícita En algunos casos necesitamos saber si una ecuación o conjunto de ecuaciones de varias variables permite definir a una de ellas o varias de ellas como función de las demás. Para el desarrollo de esta sección se ha seguido la referencia [20]. Una función y(x)está dada de forma implícita cuando está definida de la forma F(x, y)=0en lugar de la forma habitual (explícita). Dada la ecuación F(x, y)=0 (lo que se conoce como función implícita), bajo ciertas exigencias sobre la derivada de Fpodríamos, al menos localmente, despejar y=f(x). El siguiente Teorema nos da las condiciones suficientes para poder despejar una de las variables: Teorema 10 (de la Función Implícita) Sea F: Ω ⊂Rm+n→Rntal que F∈C1(Ω), y sea (a, b)∈Rm×Rnun punto tal que F(a, b)=0. Escribimos entonces DF(a, b)=[DxF(a, b), DyF(a, b)] y supongamos que DyF(a, b)es invertible. Entonces existe U⊂Rm+nyW⊂Rmcon (a, b)∈Uya∈Wtales que para cada x∈W, existe un único ytal que (x, y)∈Uy F(x, y)=0 esto define una función G:W→Rn, con G∈C1(W)y que satisface: F(x, G(x)) = 0,∀x∈W. Además, DG(x) = −[DyF(x, G(x))]−1DxF(x, G(x)) ∀x∈Wy evidentemente G(a) = b. 15 1.3. Bifurcaciones La teoría de bifurcaciones es un campo matemático centrado en el estudio de los cambios en la estructura cualitativa o topológica del comportamiento de un conjunto de ecuaciones diferenciales. Esta teoría tiene una importancia práctica muy importante en ingeniería y física. La teoría de bifurcaciones estudia el comportamiento de familias de soluciones de ecuaciones diferenciales. Generalmente en referencia a sistemas dinámicos, una bifurcación se da cuando una pequeña variación en los valores de los parámetros de un sistema (parámetros de bifurcación) causa un brusco cambio cualitativo o topológico en su comportamiento. A los valores de los parámetros en los que ocurren se les llama puntos de bifurcación. Las bifurcaciones pueden producirse tanto en sistemas continuos como en sistemas discretos. Para realizar esta sección se han seguido principalmente las referencias [12] y [19]. Antes de empezar, demos una definición sobre la forma más sencilla que adquiere una familia de sistemas. Definición 20 Se llamará forma normal de una familia de ecuaciones diferenciales (o de sistemas de ecuaciones diferenciales) con el mismo comportamiento cualitativo a la ecuación diferencial (o sistema de ecuaciones diferenciales) más simple representativa de la familia de ecuaciones diferenciales (o de sistemas de ecuaciones diferenciales). 1.3.1. Bifurcaciones en Dimensión Uno La dinámica de los sistemas unidimensionales autónomos dependientes de un parámetro es sencilla. Sin embargo, son interesantes debido a la dependencia de los parámetros, ya que la estructura cualitativa del flujo puede cambiar a medida que varían los parámetros. En particular, los puntos críticos pueden ser creados o destruidos, o pueden cambiar su estabilidad. Empecemos por una de las bifurcaciones más conocidas. Bifurcación Silla-Nodo Una bifuración silla-nodo es un mecanismo básico por el que los puntos de equilibrio se crean o se destruyen. A medida que el parámetro varía, dos puntos críticos se van acercando, colisionan y se aniquilan mutuamente. La forma normal de una bifurcación silla-nodo está dada por la ecuación diferencial de primer orden: ˙x=r+x2,(1.5) donde res un parámetro, que puede ser positivo, negativo, o cero (ver Figura 1.3). 16 Figura 1.3: Plano de fases del sistema (1.5) para distintos valores de r. La curva ˙x=r+x2 se muestra en azul y los puntos críticos en rojo. Cuando res negativo, hay dos puntos críticos, uno estable y otro inestable. A medida que rse aproxima a 0 por la izquierda la parábola se desplaza hacia arriba y los puntos críticos se acercan. Cuando r= 0, los puntos críticos se funden en un punto crítico semiestable en x= 0. Este tipo de punto crítico es extremadamente delicado, se desvanece en cuanto r > 0, y ahora ya no hay ningún punto crítico. En este ejemplo, decimos que la bifurcación ocurre en r= 0, ya que los campos vectoriales para r < 0yr > 0son cualitativamente diferentes. Bifurcación Transcrítica Hay ciertas situaciones científicas donde debe existir un punto de equilibrio para todos los valores del parámetro y nunca pueden ser destruidos. Sin embargo, el punto de equilibrio puede cambiar su estabilidad si el parámetro varía. La bifurcación transcrítica es el mecanismo estándar de este cambio en la estabilidad. Figura 1.4: Plano de fases del sistema (1.6) para distintos valores de r. La curva ˙x= rx −x2se muestra en azul y los puntos críticos en rojo. 17 La forma normal de la bifurcación transcrítica es: ˙x=rx −x2,(1.6) que se trata de la ecuación logística, pero ahora no se imponen restricciones sobre el signo de xyr. La Figura 1.4 muestra el campo vectorial a medida que varía el parámetro r. Notemos que hay un punto crítico en x= 0 para todos los valores de r. Para r < 0, hay un punto crítico inestable en x=ry un punto crítico estable en x= 0. Conforme va creciendo r, el punto crítico inestable se acerca al origen y se fusiona con el cuando r= 0. Finalmente, cuando r > 0, el origen se ha vuelto inestable y x=res ahora estable. Bifurcación Tridente Ahora pasamos a un tercer tipo de bifuración, la llamada bifurcación tridente. Este tipo de bifurcación es bastante común en problemas físicos que presentan simetría. En estos casos, los puntos de equilibrio tienden a aparecer y desaparecer por pares simétricos. Hay dos tipos de bifurcación tridente. Al tipo más simple se le llama supercrítico y al otro tipo se le denomina subcrítico. Discutiremos primero el caso supercrítico. Bifurcación Tridente Supercrítica La forma normal de una bifurcación tridente supercrítica es: ˙x=rx −x3.(1.7) Notemos que esta ecuación es invariante bajo el cambio de variables x→ −x. La invariancia es la expresión matemática de la simetría izquierda-derecha mencionada anteriormente. Figura 1.5: Plano de fases del sistema (1.7) para distintos valores de r. La curva ˙x= rx −x3se muestra en azul y los puntos críticos en rojo. 18 La Figura 1.5 muestra el campo vectorial para distintos valores de r. Cuando r < 0, el origen es el único punto de equilibrio, y es estable. Cuando r= 0, el origen sigue siendo estable, pero de forma más débil ya que el término lineal ha desaparecido. Finalmente, cuando r > 0, el origen se convierte en inestable y dos nuevos puntos críticos aparecen a ambos lados del origen, simetricamente localizados en x=±√r. Bifurcación Tridente Subcrítica En el caso supercrítico, el término cúbico es estabilizador. Si en cambio el término cúbico fuese desestabilizador, como en ˙x=rx +x3,(1.8) tendríamos una bifurcación tridente subcrítica. En la Figura 1.6 se muestra el plano de fase para este sistema: Figura 1.6: Plano de fases del sistema (1.8) para distintos valores de r. La curva ˙x= rx +x3se muestra en azul y los puntos críticos en rojo. Ahora los puntos críticos no nulos, x=±√−rson inestables y existen sólo ‘debajo’ de la bifurcación (r < 0), que motiva el término subcrítico. Más importante es que el origen es estable para r < 0e inestable para r > 0, como en el caso supercrítico, pero ahora la inestabilidad para r > 0no es opuesta al término cúbico. De hecho, este efecto lleva a explosión en tiempo finito si tenemos una condición inicial x06= 0. 1.3.2. Bifurcación de Poincaré-Andronov-Hopf Pasamos ahora a estudiar sistemas en dimensión 2. Suponemos un sistema bidimensional que tiene un punto crítico estable. ¿Cuáles son las formas posibles en las que puede perder la estabilidad a medida que varía un parámetro µ? Los autovalores del Jacobiano son la clave para encontrar la respuesta. Si el punto crítico es estable, 19 los autovalores λ1,λ2tienen que pertenecer al semiplano izquierdo, Re(λ)<0. Como los autovalores deben satisfacer una ecuación cuadrática con coeficientes reales, hay dos posibilidades: 1. Ambos autovalores son reales y negativos. 2. Ambos autovalores son complejos conjugados. Para desestabilizar el punto crítico, necesitamos que uno (o ambos) autovalores crucen al semiplano derecho, Re(λ)>0, conforme µvaría. En el apartado anterior hemos explorado los casos en los que un autovalor real pasa por λ= 0. Ahora consideramos el otro escenario posible, en los que dos autovalores complejos conjugados cruzan simultáneamente el eje imaginario hacia el semiplano derecho. Se sigue del Teorema de la Función Implícita (10) que bajo pequeñas variaciones del campo vectorial, el punto crítico persiste y no se crean nuevos puntos críticos. Sin embargo, si el tipo de estabilidad del punto crítico cambia cuando está sujeto a variaciones, entonces este cambio está normalmente acompañado de la aparición o desaparición de pequeñas órbitas periódicas que circunscriben al punto crítico. En particular, se tiene el siguiente resultado: Teorema 11 (de Poincaré-Andronov-Hopf) Sea ˙x=A(λ)x+F(λ, x) un campo vectorial plano Ck, con k≥3, que depende del parámetro escalar λy tal que F(λ, 0) = 0 yDxF(λ, 0) = 0 para todo |λ|lo suficientemente pequeño. Asumimos que la parte lineal A(λ)en el origen tiene los autovalores α(λ)±iβ(λ)con α(0) = 0 yβ(0) 6= 0. Además, suponemos que los autovalores cruzan el eje imaginario con velocidad no nula, esto es, dα dλ(0) 6= 0. Entonces, en un entorno Udel origen en R2y dado cualquier λ0>0 hay un ˆ λcon ˆ λ < λ0tal que la ecuación diferencial ˙x=A(ˆ λ)x+F(ˆ λ, x) tiene una órbita periódica no trivial en U. Dependiendo de si la órbita periódica no trivial que nace tras la bifurcación es estable o inestable, la bifurcación de Poincaré-Andronov-Hopf (a partir de ahora, bifurcación de Hopf por brevedad) será supercrítica o subcrítica respectivamente. 20 Bifurcación de Hopf Supercrítica Una bifuración de Hopf supercrítica ocurre cuando una espiral estable se transforma en una espiral inestable rodeada de un ciclo límite pequeño cuasi-elíptico. La bifurcación de Hopf puede ocurrir en espacios de fases para cualquier dimensión n≥2, pero en este trabajo nos restringiremos a dos dimensiones. Un ejemplo simple de una bifurcación de Hopf supercrítica está dado por el siguiente sistema: ˙r=µr −r3, ˙ θ=ω+br2.(1.9) Hay tres parámetros: µcontrola la estabilidad del punto crítico en el origen, ω da la frecuencia de las oscilaciones infinitesimales, y bdetermina la dependencia de la frecuencia con la amplitud para oscilaciones de gran amplitud. En la Figura 1.7 se muestra el plano de fases para µantes y después de la bifurcación. Cuando µ < 0el origen r= 0 es una espiral estable cuyo sentido de rotación depende del signo de ω. Para µ= 0 el origen sigue siendo una espiral estable, aunque muy débil. Finalmente, para µ > 0hay una espiral inestable en el origen y un ciclo límte circular estable en r=√µ. Figura 1.7: Plano de fases (en cartesianas) del sistema (1.9) con b= 1 yω= 1 para distintos valores de µ. Los puntos iniciales son (en coordenadas polares) (1/2,0) para la curva azul, (1, π/2) para la curva roja y (3/2, π)para la curva verde. El único punto crítico del sistema (0,0) se presenta en negro. Veamos qué le ha ocurrido a los autovalores del jacobiano del sistema en este proceso. Cambiando a coordenadas cilíndricas y evaluando en el origen: µ−ω ω µ →λ=µ±iω. Como esperábamos, los autovalores cruzan el eje imaginario desde la izquierda a la derecha a medida que µcrece desde valores negativos a positivos. 21 Teorema 12 (de Fenichel) Suponemos que M=M0es una subvariedad compacta normalmente hiperbólica (posiblemente con frontera) de la variedad crítica Sdel sistema (2.1) y que f,g∈Crcon r < +∞. Entonces para ε > 0lo suficientemente pequeño se tiene que: 1. Existe una variedad localmente invariante Mεdifeomorfa a M0. 2. Mεtiene una distancia Hausdorf de O(ε)desde M0. 3. El flujo de Mεconverge al flujo lento con ε→0. 4. Mεes Cr-diferenciable. 5. Mεes normalmente hiperbólica y tiene las mismas propiedades de estabilidad respecto a las variables rápidas de M0. 6. Usualmente Mεno es única. En las regiones que se mantienen a una distancia fija de la frontera de Mε, todas las variedades que satisfacen las propiedades 1-5 se encuentran a una distancia de Hausdorff de orden O(e−K/ε)las unas de las otras para algún K>0 con K=O(1). La variedad normalmente hiperbólica M0tiene asociadas variedades localmente estables e inestables. Ws loc(M0) = [ p∈M0 Ws loc(p)Wu loc(M0) = [ p∈M0 Wu loc(p)(2.5) donde Ws loc(p)son las variables locales estables de pcomo equilibrio hiperbólico de las ecuaciones de capa (2.2) mientras que Wu loc(p)son las inestables. Estas variedades también persisten para ε > 0, el resultado se matiene si intercambiamos MεyM0por Ws loc(Mε)yWs loc(M0). Llamaremos a Mεla variedad de Fenichel, subclase de las variedades lentas, variedades invariantes en las que el campo vectorial tiene una velocidad que tiende a 0 en la escala temporal rápida con ε→0. Geométricamente, la variedad estable Ws loc(Mε)de una variedad de Fenichel Mεconsiste en los puntos cuyas trayectorias se aproximan a Mεcon el tiempo creciente; análogamente, Wu loc(Mε)consiste en los puntos cuyas trayectorias se aproximan a Mεcon el tiempo decreciente. 2.1.2. Explosión Canard La dinámica de los sistemas lento-rápido en el entorno de los puntos de la variedad crítica en los que se pierde la normalidad hiperbólica puede ser sorprendentemente complicada y nada parecido a lo que conocemos con una única escala temporal. 28 Dentro del estudio de los sistemas con múltiples escalas temporales nos centraremos en el fenómeno canard. Su característica asociada, la explosión canard se manifiesta como una transición muy rápida, tras una pequeña variación de los parámetros que determinan el sistema, desde pequeños ciclos provinientes de una bifurcación de Hopf a oscilaciones de relajación, que son ciclos límite de gran amplitud formados por segmentos rápidos y lentos alternantes. Esta transición ocurre a través de una familia de ciclos cuya propiedad característica es que contienen segmentos canard, segmentos de la órbita que siguen a la variedad lenta repulsiva durante un tiempo O(1). El crecimiento rápido de la amplitud cerca de la bifurcación de Hopf es una explosión canard, más generalmente, uno se refiere a una trayectoria como órbita canard si sigue la variedad repulsiva para un tiempo O(1) en la escala temporal lenta. Una órbita canard se llama canard maximal si conecta las variedades lentas atractiva y repulsiva. Esta definición depende de la selección específica de las variedades lenta y rápida. Este fenómeno se estudió por primera vez en el contexto del sistema de Van der Pol con un forzado constante [4] (que estudiaremos más adelante en la sección siguiente). La principal condición geométrica es que la variedad crítica (2.4) tenga un fold (x0, y0)para el λ=λ∗tal que g(x0, y0, λ∗,0) = 0. Sin pérdida de generalidad podemos asumir que el punto fold se encuentra en el origen y que λ∗es el vector de parámetros nulo. Para ε= 0, un punto fold de la variedad crítica la divide localmente en dos partes, una rama atractiva, Sa, y otra rama repulsiva Sr, con respecto al flujo del sistema rápido (2.2). Lejos del límite singular y para valores de εlo suficientemente pequeños y positivos, estas dos variedades invariantes del sistema rápido persisten como variedades lentas localmente invariantes (Sε aySε r). Las variedades de Fenichel están bien definidas hasta los puntos fold donde la hiperbolicidad normal de la variedad crítica falla. Éstas se pueden extender por el flujo y se comportan de forma diferente cerca de un fold genérico y cerca de un punto canard (un fold que satisface g(0,0,0,0) y las dos condiciones para la nodegeneración gx(0,0,0,0) 6= 0 ygλ(0,0,0,0) 6= 0). Cerca de un fold genérico. Una variedad lenta de Fenichel atractiva puede seguir muy de cerca la rama atractiva Sa, pasar por un entorno del fold y continuar siguiendo aproximadamente las dinámicas rápidas, pudiendo aparecer oscilaciones de relajación, ciclos límite de gran amplitud formados por segmentos rápidos y lentos alternantes. Cerca de un punto canard. Una variedad lenta de Fenichel atractiva puede seguir muy de cerca la rama atractiva Sa, pasar por un entorno del fold y continuar siguiendo muy de cerca la rama repulsiva Sr. Esto sugiere la existencia 29 de solución en el sistema original (para 0< ε 1) conteniendo segmentos canard. Pasemos a enunciar un teorema clásico que resume cuando se va a producir una explosión canard dentro de un sistema plano lento-rápido, y cuya demostración se puede consultar en [15]: Teorema 13 (Explosión Canard en R2)Suponemos un sistema plano lento-rápido que tiene un punto tipo fold genérico p∗= (xp, yp)∈S, esto es, f(p∗, λ, 0) = 0,∂ ∂xf(p∗, λ, 0) = 0,∂2 ∂x2f(p∗, λ, 0) 6= 0,∂ ∂yf(p∗, λ, 0) 6= 0. Asumimos que la variedad crítica es localmente atractiva para x < xpy repulsiva para x>xpy que existe una singularidad fold para λ= 0 en p∗, es decir, g(p∗,0,0) = 0,∂ ∂xg(p∗,0,0) 6= 0,∂ ∂λg(p∗,0,0) 6= 0. Entonces un bifurcación de Hopf aparece en: λH=H1ε+O(ε3/2) y la Explosión Canard ocurre para: λC= (H1+K1)ε+O(ε3/2) Donde los coeficientes H1yK1se pueden calcular explícitamente.  En el límite singular tenemos λH=λC. Para cualquier ε > 0lo suficientemente pequeño el sistema lineal en el punto de la bifurcación de Hopf tiene un par de autovalores singulares. σ(λ;ε) = α(λ;ε) + ıβ(λ;ε) con α(λH;ε)=0,∂ ∂λα(λH;ε)6= 0 y l´ımε→0β(λH;ε) = ∞en la escala temporal lenta τy l´ımε→0β(λH;ε)=0en la escala temporal rápida t. Hay que resaltar que las variedades lentas típicamente no son únicas, sino que están exponencialmente cerca unas de otras. Esto permite el estudio de las posiciones relativas de las variedades lentas atractivas y repulsivas, y mostrar la existencia de canards. La presencia de términos exponencialmente pequeños en la ε-expansión de las variedades lentas implican que su posición respectiva puede cambiar con una variación exponencialmente pequeña del parámetro. 30 Cuando Sε aestá más cerca de la variedad crítica que Sε ren la región repulsiva, se fuerza la existencia de canards sin cabeza. Cuando Sε restá mas cerca de la variedad crítica que Sε aen la región repulsiva, pueden existir canards con cabeza. La transición de un tipo a otro de canard ocurre cuando Sε aestá conectada con Sε r, esto sucede a lo largo de una curva en el plano de los parámetros (ε, λ) y la solución canard asociada se llama canard maximal. Los ciclos canard se desarrollan a lo largo de una rama nacida en una bifurcación de Hopf y la explosión canard ocurre a una distancia O(ε)del punto de Hopf (muy cerca de la bifurcación), antes de la explosión, los ciclos tienen las características de un ciclo de Hopf típico. Pasemos ahora a un sencillo caso de aplicación de estos resultados que nos servirá de base para el desarrollo de los resultados posteriores en los Sistemas Lineales a Trozos. 2.2. La Ecuación de Van der Pol El contenido de esta sección, incluidas las gráficas y el tratamiento matemático, ha sido obtenido de [4]y[10]. Uno de los sistemas más simples en el que se manifiestan estos conceptos, siendo además el primer sistema donde se estudió el fenómeno canard, es la ecuación de Van der Pol con una constante λ∈Rque indica que está forzado: ε˙x=y−1 3x3+x, ˙y=λ−x. (2.6) Este sistema lento-rápido tiene sólo una variable rápida y otra lenta, pero exhibe dinámicas complejas. Como    f(x, y) = 1 ε(y−1 3x3+x), g(x, y) = λ−x, son C∞, también son localmente lipschitzianas, f, g ∈Liploc(x, R), esto implica que f, g ∈C(R;R2)∩Liploc(x, R), por lo que el Teorema 4 (de Existencia y Unicidad de Solución) se aplica para todo (x0, y0)∈R2. Sin embargo, no podemos asegurar que la función sea globalmente Lipschitz, por tanto, podrían producirse fenómenos de explosión en tiempo finito. No obstante, dicho análisis queda fuera de los objetivos del presente trabajo, en el que nos centraremos en los fenómenos canard asociados a la naturaleza lento-rápida de este sistema. 31 2.2.1. La Variedad Crítica y el Flujo Lento Si fijamos ε= 0 en (2.6), obtenemos el sistema reducido con la ecuación algebraica que define la variedad crítica como la cúbica: S={(x, y)∈R2|y=c(x) := 1 3x3−x}.(2.7) Figura 2.1: Variedad Cúbica Crítica con un punto fold genérico (a la izquierda) y un punto canard (a la derecha). La variedad central (discontinua) es repulsiva, mientras que las dos exteriores son atractivas. Figura modificada de [10]. Ses normalmente hiperbólica lejos de los extremos locales de la cúbica (p±= (±1,∓2/3)) donde Stiene un fold respecto a la variable rápida x. Por tanto, estos puntos descomponen naturalmente a Sen tres ramas: S=Sa,−∪{p−}∪Sr∪{p+}∪Sa,+,(2.8) donde Sa,±son las ramas atractivas y Sres la rama repulsiva. Notemos que no dependen de λ. La dinámica de un punto que no esté en Sse encuentra totalmente determinadas por la dirección de la variable rápida x. Observamos que espacialmente la rama de Ssituada en medio es repulsiva y las dos ramas no acotadas son atractivas (como se muestra en la Figura 2.1). Para obtener el flujo lento en Sen la ecuación de Van der Pol, no es realmente necesario resolver y=c(x)para xen una de las ramas. Es más conveniente escribir el flujo lento en términos de la variable rápida x. 32 Diferenciamos f(x, y, λ, 0) = y−c(x) = 0 respecto a τy obtenemos: ˙y= ˙xx2−˙x= ˙x(x2−1),(2.9) que junto con la ecuación (2.6) para ˙ynos permite obtener una expresión para el flujo lento: ˙x=λ−x x2−1.(2.10) El flujo lento depende de λ, porque la dirección del flujo está parcialmente determinada por la posición del equilibrio x=λen S. El flujo lento está bien definido en las tres ramas salvo en x=±1. Podemos desingularizar el flujo lento cerca de x=±1reescalando el tiempo con el factor x2−1. Esto nos proporciona la ecuación del flujo lento desingularizado,˙x=λ−x. Este reescalado invierte la dirección del flujo en Sr, así que debemos tener cuidado cuando pasemos del espacio de fases de este sistema al del flujo lento original. Figura 2.2: Planos de fase de la ecuación de Van der Pol (2.6) para λ= 0 (a)y para λ= 1 (b). Se muestra la variedad crítica S(curva gris sólida) y la nulclinas en y (línea discontinua); las flechas dobles indican la dirección del flujo rápido y las fechas simples la del flujo lento. En el panel (a)se muestra una candidata a oscilación de relajación (en azul) rodeando un equilibrio inestable. El panel (b)es el momento de una relajación de Hopf singular con una sigularidad fold en el equilibrio local p+. Figura modificada de [4]. En primer lugar, vamos a centrarnos en el caso λ= 0 (representativo para |λ|<1 y lejos de λ≈1). En este caso, el origen es el único equilibrio, y se trata de una fuente para este valor de λ. La curva cerrada (Figura 2.2 (a)) es una órbita singular compuesta por dos trayectorias rápidas que empiezan en los fold (p±) concatenadas con segmentos de S. Estas concatenaciones continuas de las trayectorias de las ecuaciones de capa y el flujo lento son llamadas candidatas. La órbita singular sigue al flujo lento en Shasta el fold, luego hace una transición al segmento de la trayectoria rápida que fluye hasta la otra rama atractiva de S, y el mismo mecanismo devuelve la órbita singular a su rama inicial de S. 33 La órbita singular se perturba para ε > 0a una órbita periódica de la ecuación de Van der Pol que se encuentra a una distancia O(ε2/3)de la candidata. Van der Pol introdujo el término oscilación de relajación para describir las órbitas periódicas que alternan entre las fases de movimiento lento y de movimiento rápido. Pasemos ahora a ver que ocurre cuando |λ| ≈ 1, donde tiene lugar la explosión canard. 2.2.2. Explosión Canard El análisis de estabilidad lineal muestra que el único punto de equilibrio (x, y) = (λ, 1 3λ3−λ)es una fuente para |λ|<1pero un sumidero para |λ|>1, con una bifurcación de Hopf para λ=±1. Como resultado nos encontramos con órbitas periódicas emergiendo del punto de equilibrio, distinguimos dos casos genéricos en la bifurcación de Hopf: Supercrítica: Si las órbitas son estables. Subcrítica: Si las órbitas son inestables. El tipo de bifurcación de Hopf está determinado por el signo del coeficiente de Lyapunov (6= 0 en el caso genérico). Cerca de la bifurcación, la amplitud de las órbitas periódicas es comparable con la raíz cuadrada de la distancia entre el parámetro y el punto de bifurcación. Figura 2.3: Continuación numerica de las órbitas de la ecuación de Van der Pol (2.6) para ε= 0.05. El panel (a)muestra una selección de órbitas periódicas: la órbita verde oscura es un pequeño ciclo límite típico cerca de la bifurcación de Hopf en λ=λH, mientras que las órbitas en negro ocurren en un intervalo pequeño de variación del parámetro en λ≈0.993491. Los paneles (b)y(c)son esquemas de los diagramas de bifurcación correspondientes a las bifurcaciones de Hopf singulares supercrítica y subcrítitca respectivamente; aquí, A denota la amplitud del ciclo límite. Figura modificada de [4]. 34 Para la ecuación de Van der Pol en λH=±1ocurre una bifurcación de Hopf supercrítica y las órbitas periódicas existen en el intervalo del parámetro |λ|<1. El caso subcrítico lo podríamos observar en sistemas más complejos, como el de FitzHugh-Nagumo [11], el cual es una extensión del sistema de Van der Pol en el que añadimos otro parámetro que permite que varíe la pendiente de la y-nulclina. El análisis de cómo cambian las dinámicas estables observadas en la ecuación de Van der Pol con λde un foco estable a oscilaciones de relajación cuando ε > 0es pequeño fue un desarrollo importante en la teoría de los sistemas lento-rápido. Como se aprecia en la Figura 2.3, cerca de la bifuración de Hopf en λH= 1 la órbita periódica es pequeña. Sin embargo, a medida que λdecrece la amplitud de la órbita crece rápidamente hasta que alcanza la forma de una oscilación de relajación. Hay un intervalo pequeño de λen el que la órbita crece como una raíz cuadrada, como se espera cerca de una bifurcación de Hopf. Sin embargo, la amplitud después crece extremadamente rápido hasta que se estabiliza en oscilaciones de relajación. 35 36 Capítulo 3 Aproximación Lineal a Trozos Actualmente, en la comunidad científica hay un gran interés por el análisis de sistemas diferenciables a trozos, ya que modelan una amplia variedad de fenómenos, como, por ejemplo, osciladores de impacto o convertidores AC/DC [5]. Un objeto de estudio interensante dentro de este campo son los Sistemas Lineales a Trozos (PWL, por sus siglas en inglés). Su estudio viene motivado por el modelado de circuitos electrónicos, a principios del siglo XX [1]. Desde entonces, se ha visto que este tipo de sistemas reproduce todos los aspectos de la dinámica no lineal y el hecho de que se puedan obtener soluciones explícitas en cada zona de linealidad hace posible describir explícitamente algunos elementos básicos de la dinámica y la geometría del sistema en cuestión. Además, se ha observado que los sistemas PWL pueden mostrar nuevos comportamientos, imposibles de obtener bajo hipótesis de diferenciabilidad. Por ejemplo, el hecho de que la unión de dos sistemas lineales estables puede producir dinámica inestable. A pesar de estas ventajas, en estos sistemas no es posible obtener una solución general ni aplicar la teoría clásica de sistemas dinámicos diferenciables, por lo que se necesita una nueva teoría específica para los PWL. 3.1. Existencia y Unicidad de Solución Antes de ponernos a trabajar en este contexto, comprobaremos que los sistemas que vamos a estudiar tienen solución, que además esta es única y está definida en R. Así, una vez que obtengamos una propuesta de solución que cumpla con las ecuaciones, sabremos que se trata de la solución que estamos buscando. Para el desarrollo de esta sección se ha seguido principalmente la referencia [9]. 37 zona central es: AC=b1 −ε0. Por tanto, δC=det(AC)=+ε > 0, y τC=Traza(AC) = b, de los que obtenemos tres casos: Si b= 0,pCes un centro. Si b < 0,pCes un foco estable. Si b > 0,pCes un foco inestable. La orientación es la de las agujas del reloj en los tres casos. Aparte del punto de equilibrio pC, hay otros dos puntos que no son equilibrios reales del sistema, pero que influyen en la dinámica de la misma forma en la que lo haría un equilibrio. A estos equilibrios se les suele llamar equilibrios virtuales en el contexto de los sistemas PWL. Estos son pL= (a, −a−ε)ypR= (a, a −ε), puntos de equilibrio del sistema en las zonas izquierda y derecha, respectivamente, pero ambos localizados en la zona central. Es fácil obtener el tipo topológico de estos equilibrios, las matrices en las zonas respectivas son: AL=+1 1 −ε0, AR=−1 1 −ε0. de las que se obtiene que δL=det(AL) = ε=δR, y τL=traza(AL)=+1=−τR, por tanto el tipo topológico de estos puntos de equilibrio virtuales es: pLes un nodo inestable. pRes un nodo estable. Estudiemos los autovalores y el signo de los mismos en cada una de las zonas, así como la forma que adoptan los espacios invariantes asociados a cada autovalor. Los autovalores de la matriz de coeficientes de la zona izquierda están dados por: λL s=1−√1−4ε 2;λL f=1 + √1−4ε 2.(3.8) Es fácil ver que λL f> λL s>0, y que l´ım ε→0λL s= 0+; l´ım ε→0λL f= 1−. Las rectas µL s≡y=−ε(x−a) λL s−a−ε;µL f≡y=−ε(x−a) λL f−a−ε, (3.9) 44 son los espacios invariantes de los autovalores lento y rápido, respectivamente, para el punto de equilibrio virtual pL. Notemos que son variedades invariantes sólo para el flujo en la zona izquierda. La pendiente de µL ftiende a cero con ε→0y la pendiente de µL stiende a −1con ε→0. Consideramos ahora los autovalores de la matriz de coeficientes de la zona derecha: λR s=−1 + √1−4ε 2;λR f=−1−√1−4ε 2.(3.10) Es fácil ver que λR f< λR s<0, y que l´ım ε→0λR s= 0−; l´ım ε→0λR f=−1+, Las líneas rectas µR s≡y=−ε(x−a) λR s +a−ε;µR f≡y=−ε(x−a) λR f +a−ε, (3.11) son los espacios o variedades invariantes correspondientes al autovalor lento y rápido, respectivamente, para el punto de equilibrio virtual pR. Notemos que son variedades invariantes sólo para el flujo en la zona derecha. La pendiente de µR ftiende a cero con ε→0y la pendiente de µR stiende a 1con ε→0. Se observa fácilmente que se dan, además, las siguientes igualdades: λL s=−λR s yλL f=−λR f. Figura 3.4: Zonas de linealidad, nulclinas y equilibrios (reales y virtuales) del sistema (3.7) y las variedades lentas correspondientea a los nodos virtuales de las zonas exteriores. Figura modificada de [10]. 45 En nuestro caso, a diferencia del contexto clásico, sí podemos calcular las variedades canónicas de Fenichel de manera explícita. El espacio invariante de la zona derecha µR ses fuertemente atractivo y el espacio invariante de la zona izquierda µL s es fuertemente repulsivo. A medida que ε→0se aproximan a la variedad crítica del sistema (y=|x|). Se sigue entonces que (ver Figura 3.4): µR sproporciona una elección canónica de una variedad lenta de Fenichel atractiva, Sε a. µL sproporciona una elección canónica de una variedad lenta de Fenichel atractiva, Sε r. Además, la zona central permite al flujo pasar a través de ella, desde la variedad atractiva Sε ahasta la variedad repulsiva Sε r. Así, la clase de sistemas (3.7) posee en principio todos los ingredientes para obtener ciclos canard. 3.3.2. Resultados principales Nuestro objetivo principal es encontrar las condiciones bajo las cuales la familia de sistemas (3.7) posee una solución canard maximal. Para ello, empezamos analizando un sistema no genérico de la familia con a=b= 0. Figura 3.5: Plano de fase para el sistema (3.7) con a=b= 0 yε= 0,1. Además representamos las variables lentas µL syµR sen rojo y en línea gris gruesa las lineas de separación. Figura modificada de [10]. 46 Teorema 15 El sistema (3.7) con a=b= 0 y0< ε 1posee un continuo de órbitas periódicas acotado inferiormente por las variedades lentas µL syµR sen sus zonas correspondientes, y por la conexión de µR sa µL sen la zona central. Por debajo del continuo, las órbitas cruzan de la zona derecha a la central y después a la izquierda, en la que escapan al infinito.  La Figura 3.5 nos muestra el plano de fase del sistema (3.7) bajo estas hipótesis. Notemos que la situación no genérica descrita en el Teorema 15 es análoga a la del sistema Hamiltoniano encontrado en el análisis cerca del punto canard en el caso diferenciable [14]. En este caso encontramos una órbita especial que acota al continuo de órbitas periódicas. Esta órbita está formada por tres partes: Las variedades lentas invariantes µL syµR sen sus zonas correspondientes. Una conexión entre ellas en la zona central, como se aprecia en la Figura 3.6: Figura 3.6: Representación de un posible canard maximal en el sistema (3.7). Figura modificada de [10]. Como hemos descrito anteriormente, en nuestro sistema (3.7) las variedades lentas invariantes µL syµR sjuegan el papel de las variedades lentas invariantes repulsiva 47 y atractiva Sε rySε a, respectivamente. Entonces podemos decir que esta órbita especial es la órbita canard maximal y que las órbitas periódicas del continuo son las órbitas canard periódicas. La pregunta que naturalmente surge es si perturbando esta situación no genérica es posible encontrar valores apropiados para los parámetros de tal forma que la órbita canard maximal persista. En el siguiente teorema enunciamos la existencia de una curva en el espacio de parámetros a= ˜a(b, ε)tal que el canard maximal se mantiene tras la perturbación. Teorema 16 Existe una función a= ˜a(b, ε), analítica como función de (b, √ε), definida en un conjunto abierto U⊂R2conteniendo al origen y tal que para (b, √ε)∈U∩{ε > 0}, el sistema (3.7) posee una órbita conectando las variedades lentas de las zonas exteriores. Además, tenemos que el desarrollo en serie de Taylor con respecto a b cerca de b= 0 de la función aestá dada por: ˜a(b, ε) = ε((τ∗/2) −1)b/2 + O(b3) =π 4√ε+O(ε)b+O(b3),(3.12) donde τ∗=π−arcsin 2√ε √ε.(3.13)  Demostración: A partir de ahora denotaremos pR S(µ), donde µ= (a, b, ε), al punto de intersección entre la variedad lenta derecha, µR S, y la recta de separación x=ε/(1+b), y a pL S(µ)al punto de intersección entre la variedad lenta izquierda, µL S, y la línea de separación x=−ε/(1 −b). La existencia de una solución canard maximal se reduce a la existencia de una órbita conectando los puntos pR S(µ)ypL S(µ)(ver Figura 3.6). El conjunto de condiciones que caracterizan a esta conexión está dado por la existencia de: τ > 0,0< ε 1,|b|<2√εya∈(−ε/(1 −b), ε(b+ 1)),tales que, para pR S(µ) = (ε(b+ 1), a −ε−ε(ε/(1 + b)−a)/λR s)T:= (xR 0, yR 0)T,(3.14) las siguientes condiciones se satisfacen: xC(τ;µ, pR S) = −ε 1−b,(3.15) yC(τ;µ, pR S) = −ε λL S (xC(τ;µ, pR S)−a)−a−ε, (3.16) 48 xC(s;µ, pR S)∈−ε 1−b,ε 1 + b∀s∈(0, τ).(3.17) En una primera etapa, analizamos la existencia de soluciones de las ecuaciones (3.15)-(3.16), que son las ecuaciones de cierre asociadas a la conexión. Definimos las siguientes funciones:      F(τ, a, b, ε) = xC(τ;µ, pR S) + ε 1−b, G(τ, a, b, ε) = yC(τ;µ, pR S) + ε λL S (xC(τ;µ, pR S)−a) + a+ε. (3.18) Entonces, la solución de las ecuaciones de cierre es equivalente a la solución (τ, a, b, ε) del sistema F(τ, a, b, ε)=0, G(τ, a, b, ε)=0.(3.19) Además, la desigualdad (3.17) es equivalente a F(s, a, b, ε)∈0,ε 1−b+ε 1 + b∀s∈(0, τ).(3.20) Integrando el sistema lineal de la zona central con la condición inicial, pR S, obtenemos la siguiente expresión explícita de la solución:                  xC(τ;µ, pR S) = exp(bτ/2) ab +a−ε (b+ 1)√4ε−b2 ·(b−2λR S) sen(βτ)−2βcos(βτ)+a, yC(τ;µ, pR S) = exp(bτ/2) ab +a−ε (b+ 1)√4ε−b2 ·(b(λR S)+2ε) sen(βτ)+2β(b−λR S) cos(βτ)−ab, (3.21) donde β(b, ε) = √4ε−b2/2.(3.22) Tras reescalar los parámetros como sigue: b= 2ˆ b√ε;τ= ˆτ/ε;a= ˆaε, (3.23) encontraremos una rama de soluciones del sistema (3.19). Sean φ(ˆτ, ˆa,ˆ b, ε) = exp(ˆ bˆτ)(ˆa−1+2√εˆaˆ b) (1 + 2√εˆ b)p1−ˆ b2, θ(ˆ b, ˆτ) = ˆτp1−ˆ b2. (3.24) 49 Se sigue de la definición de F(3.18) y de la fórmula para xC(3.21) que F=εˆ F, donde ˆ Festá dado por ˆ F(ˆτ, ˆa,ˆ b, ε) = φ−p1−ˆ b2cos(θ)+(ˆ b+√ε) sen(θ) +1 1−2√εˆ b+O(ε).(3.25) Para Gesta estimación es más complicada. Primero, notemos que yC=ε3/2φ((2ˆ b+ε)q1−ˆ b2cos(θ)+(−2ˆ b2+√εˆ b) sen(θ)) −2ˆaˆ b+O(ε). En segundo lugar, reescribimos G como sigue: G=yC+ε λL S F−ε λL S ε 1−pεˆ b +εˆa!+εˆa+ε. Por las relaciones (3.8) es claro que ε/λL S= 1 + O(ε), por tanto, −ε λL S ε 1−√εˆ b+εˆa!+εˆa+ε=−ε 1−2√εˆ b+ε+O(ε2) = −2ε3/2+O(ε2). Se sigue entonces que G=ε3/2ˆ G+ (ε/λL S)ˆ F, donde ˆ Festá dado por (3.25) y ˆ G(ˆτ, ˆa,ˆ b, ε) = φ(2ˆ b+√ε)p1−ˆ b2cos(θ) + (1 −2ˆ b2+εˆ b) sen(θ) −2ˆaˆ b−2ˆ b+O(ε). Por tanto, resolver F=G= 0 para ε > 0es equivalente a resolver ˆ F=ˆ G= 0. Para ε= 0, tenemos: ˆ F=exp(ˆ bˆτ)(ˆa−1) p1−ˆ b2−q1−ˆ b2cos(θ) + ˆ bsen(θ)+ ˆa+ 1, ˆ G=exp(ˆ bˆτ)(ˆa−1) p1−ˆ b2(2ˆ b)q1−ˆ b2cos(θ) + (1 −2ˆ b2) sen(θ)−2ˆaˆ b−2ˆ b. (3.26) Para ˆa=ˆ b= 0, la ecuación ˆ F= 0 nos da la condición cos(ˆτ) = −1, de donde ˆτ=π. Además la ecuación ˆ G= 0 se satisface con ˆa=ˆ b= 0 yˆτ=π. Para aplicar el Teorema de la Función Implícita es necesario probar que det(J(π, 0,0,0)) 6= 0, donde J(ˆτ, ˆa,ˆ b, ε) =     ∂ˆ F ∂ˆτ(ˆτ, ˆa,ˆ b, ε)∂ˆ F ∂ˆa(ˆτ, ˆa,ˆ b, ε) ∂ˆ G ∂ˆτ(ˆτ, ˆa,ˆ b, ε)∂ˆ G ∂ˆa(ˆτ, ˆa,ˆ b, ε)     .(3.27) 50 Es fácil ver que ∂ˆ F ∂ˆτ(π, 0,0,0) = 0,∂ˆ F ∂ˆa(π, 0,0,0) = 2, ∂ˆ G ∂ˆτ(π, 0,0,0) = 1,∂ˆ G ∂ˆa(π, 0,0,0) = 0, y entonces det(J(π, 0,0,0)) = −26= 0. Por tanto, del Teorema de la Función Implícita concluimos que existe un conjunto U∈R2abierto que contiene a (ˆ b, ε) = (0,0) y dos funciones analíticas ¯τ(ˆ b, ε)y¯a(ˆ b, ε) definidas en U, tales que ¯τ(0,0) = π,¯a(0,0) = 0 y cuando ε > 0: (˜τ(b, ε),˜a(b, ε), b, ε) es una solución del sistema (3.19), donde ˜τ= ¯τ(b/(2√ε), ε)/√ε, ˜a= ¯a(b/(2√ε), ε)/ε. Para asegurar la existencia de la conexión, además, debemos probar que las soluciones de las ecuaciones de cierre (3.15) y (3.16) satisfacen la desigualdad (3.17), o equivalentemente, la condición (3.20). Necesitamos probar F(s, ˜a, b, ε)∈0,ε 1−b+ε 1 + b∀s∈(0,˜τ), Lo haremos por contradicción. Asumimos que existe τ1<˜τtal que F(τ1,˜a, b, ε)> ε/(1 −b) + ε/(1 + b), o equivalentemente, xC(τ1,˜η, pR S(˜η)) >−ε/(1 −b), donde ˜η= (˜a, b, ε). Entonces, para permitir al flujo volver a la recta x=−ε/(1 −b), es necesario que ˜τ≥2π/√ε≥2τ∗, donde τ∗=π−sen−1(2√ε) √ε, y corresponde a la solución para a=b= 0. Pero ˜τ(0, ε) = τ∗, lo que contradice ˜τ≥2τ∗. Un razonamiento análogo nos permite probar que F(s, ˜a, b, ε)>0∀s∈(0,˜τ). Por consiguiente, la prueba de la existencia de una subfamilia de sistemas PWL (3.7) que tienen una solución canard maximal está completa. Pasemos a demostrar la expresión en serie de Taylor para ˜a: Procederemos por diferenciación implícita. Sabemos que (˜τ(b, ε),˜a(b, ε), b, ε)es una solución del sistema (3.19), esto es, F(˜τ(b, ε),˜a(b, ε), b, ε)=0, G(˜τ(b, ε),˜a(b, ε), b, ε)=0,(3.28) 51 y además, se satisface que ˜τ(0, ε) = τ∗y˜a(0, ε) = 0. Denotamos p∗=pR S(0,0, ε) = (ε, λR Sε). Derivando en el sistema (3.28) con respecto a ben b= 0, es fácil ver que ∂˜τ ∂b (0, ε) y∂˜a ∂b (0, ε)deben satisfacer el siguiente sistema de ecuaciones:    ∂F ∂τ (τ∗,0,0, ε)∂F ∂a (τ∗,0,0, ε) ∂G ∂τ (τ∗,0,0, ε)∂G ∂a (τ∗,0,0, ε)     ∂˜τ ∂b (0, ε) ∂˜a ∂b (0, ε)  =   ∂F ∂b (τ∗,0,0, ε) ∂G ∂b (τ∗,0,0, ε)  .(3.29) Empecemos calculando ∂G ∂τ (τ, a, b, ε). Por la definición de la función G(ver 3.18), ∂G ∂τ (τ, a, b, ε) = ∂yC ∂τ (τ;η, pR S) + ε λL S∂xC ∂τ (τ;η, pR S).(3.30) Del sistema (3.7), en la zona central ˙x=y+by,˙y=ε(a−x), y entonces, ∂G ∂τ (τ, a, b, ε) = ε(a−xC(τ;η, pR S)) + ε λL SyC(τ;η, pR S) + bxC(τ;η, pR S).(3.31) Ahora, se sigue de xC(τ∗; 0, p∗) = −ε,yC(τ∗; 0, p∗) = λR Sε, y de la expresión (3.31) que ∂G ∂τ (τ∗,0,0, ε)=0. El cálculo de ∂F ∂τ (τ, a, b, ε)es similar al de ∂G ∂τ (τ, a, b, ε), ∂F ∂τ (τ, a, b, ε) = ∂xC ∂τ (τ;η, pR S) = yC(τ;η, pR S) + bxC(τ;η, pR S),(3.32) y entonces ∂F ∂τ (τ∗,0,0, ε) = yC(τ; 0, p∗) = λR Sε. (3.33) Ahora, obtenemos ∂G ∂a (τ, a, b, ε). De la definición (3.18), ∂G ∂a (τ, a, b, ε) = ∂yC ∂a (τ;η, pR S) + ε λL S∂xC ∂a (τ;η, pR S)−1+ 1.(3.34) 52 De la expresión explícita de la solución de en la zona central dada en (3.21), no es difícil llegar a ∂xC ∂a (τ∗; 0, p∗)=2, ∂yC ∂a (τ∗; 0, p∗)=2λR 2, (3.35) de donde obtenemos ∂G ∂a (τ∗,0,0, ε)=2. Ahora, para finalizar el cálculo de la matriz de coeficientes del sistema (3.29), sólo tenemos que encontrar ∂F ∂a (τ∗,0,0, ε). Tomando la derivada de la función F respecto a a, ∂F ∂a (τ, a, b, ε) = ∂xC ∂a (τ;η, pR S).(3.36) A la luz de la primera expresión de (3.35), la evaluación en (τ∗,0,0, ε), nos lleva a ∂F ∂a (τ∗,0,0, ε)=2.(3.37) Calculamos el segundo miembreo del sistema (3.29). De la definición de la función F, (ver 3.18) ∂F ∂b (τ, a, b, ε) = ∂xC ∂b (τ;η, pR S) + ε (1 −b)2,(3.38) y entonces, ∂F ∂b (τ∗,0,0, ε) = ∂xC ∂b (τ; 0, p∗) + ε. (3.39) De la expresión explícita de la solución en la zona central dada en (3.21), obtenemos ∂xC ∂b (τ∗; 0, p∗) = −(ετ∗)2, ∂yC ∂b (τ∗; 0, p∗) = ε(1 −τ∗/2), (3.40) y de la primera expresión, conseguimos ∂F ∂b (τ∗,0,0, ε) = ε(1 −τ∗/2).(3.41) Buscando ∂G ∂b (τ∗,0,0, ε), de la definición de la función G(ver 3.18), tenemos ∂G ∂b (τ, a, b, ε) = ∂yC ∂b (τ;η, pR S) + ε λL S ∂xC ∂b (τ;η, pR S).(3.42) 53 Figura 3.10: A la izquierda, diagrama de bifurcación para ˆ b > 0. Consideramos ε > 0 fijado y a > 0en el sector situado más a la derecha. Decreciendo a, toma lugar una bifurcación de Hopf, dando lugar a un pequeño ciclo límite estable. El ciclo límite crece a medida que adecrce. Cuando aalcanza la región sombreada en gris oscuro, el ciclo límite se convierte en un ciclo canard. A lo largo de la línea situada más a la izquierda la familia de ciclos canard terminan en una conexión de canard maximal. A la derecha, rama explosiva de ciclos límite obtenida mediante simulación directa cuando variamos el parámetro acon ε= 0,1yb= 0,009944 fijos. Figura modificada de [10]. La desaparición del ciclo está probablemente relacionada con una bifuración global en el infinito que podría analizarse compactificando el plano de fases. En el corolario anterior se aclara que podemos reproducir con la clase de Sistemas Lineales a Trozos (3.7) diagramas de bifurcación similares a los obtenidos en el caso diferenciable. 3.5. Conclusiones Como hemos visto, los sistemas PWL son capaces de reproducir dinamica de tipo canard, con la ventaja de que podemos elegir representantes canónicos de las variedades de Fenichel, Sε aySε r, debido a la posibilidad de calcular analíticamente las variedades lentas asociadas a los equilibrios en nuestro modelo PWL. A su vez, nuestros sistemas PWL también pueden repoducir otras dinámicas asociadas a fenómenos canard más complejos que se presentan en el caso diferenciable. Como son: 1. Si incluimos un parámetro, κ, que se encargue de romper la simetría de las zonas exteriores, podemos obtener la aparición de más de un ciclo límite ca- 60 nard, como además puede producirse en el caso [17] que estudiaremos en la siguiente sección. 2. Si ampliamos el número de dimensiones en el que estamos trabajando podemos crear un modelo PWL que reproduzca las llamadas Oscilaciones de Tipo Mixto (o MMO por sus siglas en inglés), que son un tipo de oscilaciones que presenta (al menos) oscilaciones de dos periodos y amplitudes distintas, cuya combinación produce la aparición de dinámicas complejas [3]. Destacamos que los sistemas con varias escalas de tiempo son los utilizados para estudiar este tipo de comportamientos y que las neuronas muestran este tipo de comportamiento mixto cuando se comunican entre ellas [6], [7]. 61 62 Capítulo 4 Estudio de un ejemplo concreto: Ciclos Canard en la Dinámica de Aeronaves en Tierra En este capítulo estudiaremos una versión PWL de un modelo de pérdida de estabilidad general de aeronaves posadas en el suelo presentado en [17]. En primer lugar, realizaremos una presentación del modelo original para establecer las propiedades fundamentales que debe reproducir nuestro modelo PWL. Figura 4.1: Representación esquemática de los ejes de referencia de una aeronave, junto a su centro de gravedad. 63 4.1. Presentación del Caso 4.1.1. Introducción La pérdida de estabilidad lateral de un avión de pasajeros cuando realiza un giro mientras está posado en el suelo bajo un empuje constante de sus motores está asociada con la saturación de las fuerzas de frenado. Esta propiedad física del sistema se puede relacionar directamente con una explosión canard. El modelo industrial de la aeronave no contiene explícitamente una separación de escalas de tiempo, aunque localmente, y cerca de la transición, la velocidad de la aeronave actúa como una variable lenta. Debido a la complejidad de las ecuaciones que se consideran en el modelo, no es práctico un estudio analítico directo. Sin embargo, una simulación numérica consigue identificar propiedades del sistema características de la separación de escalas de tiempo y oscilaciones de tipo canard. Este caso difiere de los ejemplos de explosión canard anteriores en que la variedad crítica (nulclina rápida) tiene un único pliegue y una asíntota en la dirección rápida, en lugar de ser una cúbica con dos pliegues. Para resaltar las principales diferencias cualitativas entre las explosiones canard clásicas de la literatura y la de la aeronave, en [18] se presenta un modelo bidimensional tratable analíticamente, mientras que en este capítulo, desarrollaremos una versión PWL de dicho modelo. 4.1.2. Modelo de la Aeronave El modelo matemático aquí presentado fue específicamente desarrollado en [18]. Un esquema de la disposición de los elementos y las fuerzas se puede consultar en esa referencia. Consideramos un modelo matemático completamente parametrizado de una aeronave para pasajeros de tamaño medio. La aeronave tiene una configuración de triciclo con dos trenes de aterrizaje principales situados en las alas y un tren de aterrizaje situado en la proa que se utiliza para realizar los cambios de dirección. La aeronave se representa con seis grados de libertad y los efectos no lineales importantes se capturan en las componentes del modelo para los neumáticos y la aerodinámica. Los tres primeros grados de libertad ya proporcionan una descripción razonable de la dinámica, estos son, la velocidad longitudinal vx, la velocidad lateral vy, y la velocidad de rotación a lo largo del eje vertical ωz. Como queremos modelar de forma precisa la asimetría de carga entre los trenes de aterrizaje, es necesario incluir otros tres grados de libertad, la velocidad vertical vz, la velocidad de alabeo (a lo largo del eje longitudinal) ωx, y la velocidad de cabeceo ωy(a lo largo del eje lateral). Los ejes pueden consultarse en la Figura 4.1. 64 Las ecuaciones de movimiento de este sistema se derivan del equilibrio de fuerzas y momentos en cada grado de libertad, y son las siguientes: m(˙vx+vzωy−vyωz) = FxT −FxR −FxL −FxN cos(δ)−FyN sen(δ)−FxA, m(˙vy+vxωz−vzωx) = FyR +FyL +FyN cos(δ)−FxN sen(δ) + FyA, m(˙vz+vyωx−vxωy) = FzW −FzR −FzL −FzN −FzA, Ixx ˙ωx−(Iyy −Izz)ωyωz=lyLFzL −lyRFzR −lzLFyL −lzRFzR −lzN FyN cos(δ) + lzN FxN sen(δ) + lzAFyA +MxA, Iyy ˙ωy−(Izz −Ixx)ωzωx=lxN FzN −lzN FxN cos(δ)−lzN FyN sen(δ) −lxRFzR −lzRFxR −lxLFzL −lzLFxL +lzT FxT +lxAFzA +lzAFxA +MyA, Izz ˙ωz−(Ixx −Iyy)ωxωy=lyRFxR −lyLFxL −lxRFyR −lxLFyL +lxN FyN cos(δ)−lxN FxN sen(δ) + lxAFyA +MzA. (4.1) En este capítulo identificaremos un caso específico, estableciendo los valores de algunos de los parámetros. Se fija la masa de la aeronave en m= 45420 Kg (caso ligero) y consideramos los momentos de inercia principales (Ixx, Iyy, Izz). El ángulo de dirección se aplica al tren de proa y se denota por δ. El peso de la aeronave es FzW =mg que actúa en el centro de masa que se sitúa a una distancia del 14 % del eje aerodinámico principal medido desde la proa y sobre ese mismo eje. Los elementos de fuerza ortogonales en cada uno de los neumáticos se denotan por F∗i con i=L, R, N, la fuerza aerodinámica y los elementos del momento actúan en el centro aerodinámico de la aeronave y se denotan por F∗AyM∗A, respectivamente. La fuerza de empuje se asume que actúa paralelamente al eje de la aeronave y se denota por FxT . Las distancias relativas al centro de gravedad de la aeronave se denotan por l∗. El modelo matemático anterior (4.1) ha sido validado en la industria aeronáutica (más detalles en [18]). Como se puede apreciar en las ecuaciones (4.1), las fuerzas laterales FxL, FyL, FxR yFyR, dependen de forma no lineal de la carga vertical en el neúmatico de proa, FxN , FyN y su ángulo de dirección del desplazamiento, δ. Para ángulos pequeños la relación entre el álgulo de desplazamiento y la fuerza lateral es lineal, pero si aumentamos el ángulo, la fuerza satura antes de decrecer. Esta característica no lineal juega un papel muy importante en los resultados que expondremos a continuación. 65 4.2. Resultados del Caso Diferenciable 4.2.1. Investigación Uniparamétrica Las soluciones circulares de radio constante son estados estacionarios del modelo que podemos seguir bajo variaciones paramétricas. Como el ángulo de desplazamien- to o la potencia del empuje, el resto de los parámetros se fijan en valores realistas. Figura 4.2: La Figura (a) muestra una rama (gris) de estados estacionarios (continua si es estable y discontinua si es inestable) representados por la velocidad longitudinal vx como función del ángulo de desplazamiento δ. Los cambios de estabilidad ocurren en los puntos de Hopf H1yH2, la velocidades máxima de vxa lo largo de una rama de soluciones periódicas se dibujan en negro y la mínima en verde. La Figura (b) muestra la proyección en el plano de fase (vx, vy)de la órbita periódica con δr= 6.5o. La Figura (c) muestra la trayectoria del centro de gravedad en el plano (X, Y ), con marcas a tiempos equiespaciados. Figura modificada de [17]. Desde el punto de vista de las aplicaciones, el objetivo es encontrar los límites de operación segura para una aeronave girando sobre el suelo, que corresponde a encontrar la región paramétrica donde la solución es estable. Fuera de los límites de esta región, el comportamiento inestable se traduce en que la aeronave no puede seguir la solución de radio constante. Los límites de estabilidad corresponden a bifurcaciones de Hopf que producen nuevas órbitas periódicas estables, que son soluciones de radio no constante. En la Figura 4.2(a) se muestra un análisis de una bifuración uniparamétrica del 66 ángulo de desplazamiento con la velocidad longitudinal utilizada como medida de la solución. El empuje de los motores se ha fijado en un 14 % del máximo valor posible, para el que existe dinámica inestable [18]. Los cambios de estabilidad ocurren en las bifurcaciones de Hopf H1para δH1≈2.94oyH2para δH2≈10.35o. Una rama de soluciones periódicas que conectan las dos bifurcaciones de Hopf se representa por dos curvas, una negra y otra verde, que muestran los valores máximos y mínimos de vx, respectivamente, la curva maximal sigue de cerca las soluciones estacionarias inestables. Cerca de ambos notamos una caída brusca del valor mínimo de vxjunto a la órbita periódica que alcanza un valor más bajo que cambia de forma gradual con δ. La Figura 4.2(b) muestra una proyección de una órbita periódica en el plano (vx, vy)con δr= 6,5oy un periodo de 81s. Se pasa mucho más tiempo cerca de vy= 0, lo que indica que vxes una variable lenta en esta parte del espacio de fase. Esta parte lenta está seguida por una transición rápida de vuelta a pequeños valores de vx, esto es, una oscilación de relajación con un segmento lento y otro rápido que corresponde con un giro inestable. En la Figura 4.2(c) se ilustra esto último. En ella se muestra la trayectoria del centro de gravedad de la aeronave sobre el plano de tierra (X, Y )(los marcadores estan equiespaciados temporalmente). La condición inicial de la órbita periódica se escoge donde la velocidad longitudinal es máxima. La aeronave intenta seguir un círculo de giro de radio constante pero empieza a cerrarse, resultando en una pérdida de la estabilidad lateral y un giro, que además, va acompañado con un descenso brusco de vx. La aeronave entonces se para, después de lo cual se acelera bajo el empluje constante de los motores para de nuevo intentar seguir un círuclo de radio constante y el proceso se repite. 4.2.2. Identificación de Ciclos Canard Nos concentramos en la naturaleza exacta de las caídas abruptas de vxcerca de los dos puntos de Hopf. Cerca de H2la rama de soluciones periódicas inicialmente crece con una progresión similar a una raíz cuadrada. Sin embargo, a medida que decrece δhay una caída brusca y repentina del valor mínimo de vxpara δc2≈9.83o. Cerca de H1encontramos una caída incluso más brusca de vxpara δc1≈2.92o. Este cambio brusco de la diferencia de amplitud para un rango pequeño del parámetro cerca de ambas bifurcaciones de Hopf es característico de una explosión canard. Los puntos δc1yδc2se determinan por ser los puntos de las ramas de órbitas periódicas con el mayor incremento de amplitud. Veamos que, efectivamente, son explosiones canard. Debido a la complejidad de las ecuaciones de movimiento no es posible identificar la separación de escalas temporales en términos de los parámetros del sistema. Sin embargo, la encontramos en la Figura 4.2(b) de dinámica lento-rápida en las soluciones periódicas de tipo relaja- 67 ción para δr, en la que vxactúa como la variable lenta y vycomo la variable rápida del sistema. Como la órbita periódica es altamente atractiva, las variables vz, ωx, ωy se mantienen prácticamente constantes, la velocidad rotacional ωzsí varía significativamente durante el movimiento, pero sigue de forma pasiva la dinámica de vy. Por tanto, podemos tomar esta última como representativa de las variables rápidas del sistema. Presentamos entonces la dinámica en el plano (vx, vy)y progresamos computando las variedades críticas cerca de H1yH2. Figura 4.3: Familia de órbitas canard en el plano (vx, vy)de la explosión canard para δc1 cerca de H1en (a) y para δc2cerca de H2en (b). La curva gruesa con la flecha es el canard maximal. También se representa la variedad crítica con su rama estable de forma continua y la inestable discontinua. Figura modificada de [17]. Añadimos la ligadura ˙vx= 0, así conseguimos congelar la dinámica de la variable lenta vxy trantándola como parámetro y computando los estados estacionarios del sistema reducido, podemos obtener las variedades críticas asociadas con las explosiones canard cerca de H1yH2. Cabe destacar que cerca del punto de Hopf y durante la explosión canard respectiva no se observa dependencia de la variedad crítica con el ángulo de desplazamiento. La Figura 4.3 nos muestra las variedades críticas para δc1yδc2. También se muestran las respectivas familias de órbitas periódicas cercanas a los puntos de Hopf. Ambas variedades críticas tienen una rama superior atractiva y una inferior repulsiva. La Figura 4.3 exhibe las principales características asociadas a los ciclos canard y la explosión canard. La primera de estas características es que muy cerca del equilibrio del sistema completo que pasa por una bifurcación de Hopf en H1yH2, hay puntos fold de una variedad crítica en F1yF2, respectivamente. La distancia entre estos puntos en el espacio de fase es del orden de la separación de escalas de tiempo, ε. La segunda de las características es que las órbitas periódicas crecen rápidamente en tamaño si variamos δalejándola de los puntos de Hopf, y siguen no sólo la parte estable de la variedad crítica sino también un segmento de su rama inestable, lo que claramente las identifica como órbitas periódicas canard. 68 En la Figura 4.3 identificamos un canard maximal que sigue la parte inestable durante el mayor tiempo. Las órbitas periódicas más pequeñas que el canard maximal suben desde la rama inestable hasta la estable, las llamadas canards sin cabeza, en contraste con las órbitas periódicas más grandes que el canard maximal, que inicialmente se mueven hacia abajo desde la rama inestable y vuelven a la rama estable después de pasar por valores negativos de vy, las llamadas canards con cabeza. Las órbitas canard existen para rangos exponencialmente pequeños del parámetro δ. Cabe destacar que la órbita más grande en la Figura 4.3(b) es realmente una oscilación de relajación, ya que no sigue la parte inestable de la variedad crítica (comparando con la Figura 4.2(b) vemos que tienen la misma forma). 4.2.3. Modelo Plano de la Explosión Canard Para capturar las características cualitativas clave de la explosión canard descrita anteriormente, en el trabajo [17], se presenta el modelo del campo vectorial plano (˙x=dx dτ ): (ε˙x=−y−(x−a)ex b, ˙y=x−α, (4.2) donde 0< ε 1es el parámetro de separación de las escalas temporales, y αes el parámetro de la bifurcación que determina el punto de equilibrio. Los parámetros a ybdeterminan la forma de la variedad crítica que no depende de α. Esta variedad tiene un fold en x=a−by la recta y= 0 es una asíntota horizontal donde xtiende a−∞. Hay un punto de Hopf para α=αH=a−b, exactamente en el fold de la variedad crítica. Figura 4.4: Diagrama de bifurcación uniparamétrica de xen función de α, en él se ven el equilibrio (en rojo, discontinuo) y el mínimo en x de la órbita periódica de la bifurcación (en azul, continuo). 69 Figura 4.8: Variedad crítica del modelo (4.15) para s= 40 en azul junto a su flujo y a la variedad crítica original en verde. También se marcan con un punto rojo los puntos de corte con los ejes. En la Figura 4.8 podemos apreciar una representación de esta nueva variedad crítica junto a su flujo. Aunque en las Figuras 4.7 y 4.8 parece que el flujo es el mismo, al hacer una ampliación de un entorno de la nueva zona se aprecia claramente que el flujo no es el mismo (ver Figura 4.9). Figura 4.9: Ampliación de la zona del máximo de las Figuras 4.7 y 4.8 Puntos de Equilibrio Añadir una nueva zona en el sistema lineal (4.15) tiene un efecto inmediato en los puntos de equilibrio con respecto a los de (4.9), ya que pasan de ser uno real y dos virtuales (en (4.9)) a ser uno real y tres virtuales (en (4.15)), lo que a su vez 76 conlleva a un cambio en la dinámica global del sistema. Los nuevos valores posibles para los puntos de equilibrio son: pLL ≡(α+b−a, 0), pL≡(α+b−a, β s(α+b−a+s)), pC≡(α+b−a, γ(α+b−a) + β−ε), pR≡(α+b−a, β(a−α) b). (4.16) De estos equilibrios, uno será el equilibrio real del sistema y los otros tres serán equilibrios virtuales, que, como hemos indicado en el Capítulo 3, influyen en la dinámica de la misma forma en la que lo haría un equilibrio. Los equilibrios virtuales se sitúan en la zona del equilibrio real. El equilibrio real será: pLL si α+b−a < −s, pLsi −s<α+b−a < ε γ−β s , pCsi ε γ−β s < α +b−a < ε γ+β b , pRsi α+b−a > ε γ+β b . Pasamos a obtener el tipo topológico de los equilibrios, calculando la matriz de coeficientes correspondiente a cada zona: ALL =0−1 ε0, AL=β s−1 ε0, AC=γ−1 ε0, AR=−β b−1 ε0, de las que obtenemos que δLL =δL=δC=δR=ε > 0,τLL = 0,τC=γ,τL=β/s yτR=−β/b. Como estamos suponiendo que β2−4εs2>0,β2−4εb2>0y γ2−4ε < 0, se tiene que τ2−4δ > 0en todas las zonas, menos en la central, por lo que el tipo topológico de los equilibrios es: pLL es un centro, pLes un nodo inestable, pRes un nodo estable, mientras que el tipo topológico de pCdepende del valor del parámetro γ, como sabemos que γ2−4ε < 0, se tiene la siguiente clasificación: 77 Si γ < 0,pCes un foco estable. Si γ > 0,pCes un foco inestable. Si γ= 0,pCes un centro. La orientación es la contraria a las agujas del reloj en los tres casos. En cuanto a los autovalores de la matriz de coeficientes y sus variedades invariantes correspondientes son: Para la zona izquierda: los autovalores (4.11) y las variedades invariantes (4.12). Para la zona derecha: los autovalores (4.13) y las variedades invariantes (4.14). Es decir, las mismas que para el modelo (4.9). De esta forma, podemos utilizar para el modelo (4.15) los resultados de las zonas izquierda y derecha obtenidos en el modelo (4.9). Una vez obtenidos los puntos notables para la dinámica del sistema y conocidos sus tipos topológicos, así como las variedades invariantes de las zonas izquierda y derecha, vamos a realizar simulaciones numéricas con el fin de comprobar que nuestros modelos reproducen la dinámica de la explosión canard del modelo de partida (4.2) y recoger las similitudes y diferencias que se presentan. 4.4. Simulación y Validación En esta sección vamos a realizar algunos tests numéricos que nos permitan validar nuestros modelos lineales a trozos (4.9) y (4.15). Para la resolución de los sistemas utilizaremos un método de Runge-Kutta de cuarto orden con paso h= 0.01, ya que las características asociadas a los modelos (definición a trozos, sistema lento-rápido) necesitan de una precisión elevada. 1. En primer lugar, simularemos un caso de tipo académico, correspondiente a un sistema reversible localmente (en un entorno de x= 0) con pendiente uno en valor absoluto. Estudiaremos el modelo para tres (4.9) y cuatro zonas (4.15) con el objetivo de reproducir la dinámica predicha por la teoría, similar a los casos de los modelos lineales a trozos (3.4) y (3.7) de la ecuación de Van der Pol, esto es, la superexplosión y la explosión canard, respectivamente. 2. En segundo lugar, con el propósito de recrear situaciones más realistas vamos a no considerar la restricción sobre la pendiente. De esta forma tenemos dos posibilidades, el caso nodo-nodo (cualitativamente análogo al caso anterior) y 78 el caso foco-foco (no representativo del fenómeno canard). De nuevo, estudiaremos ambos modelos con el objetivo de reproducir la dinámica predicha por la teoría para cada caso. 4.4.1. Caso Reversible con Pendiente uno En esta sección, vamos a realizar tests numéricos con el objetivo de reproducir los resultados de la teoría para el caso académico para ambos modelos, (4.9) y (4.15). Vamos a imponer que el sistema sea reversible localmente en un entorno de x= 0 y que la pendiente sea igual a uno en valor absoluto. Para que el sistema sea reversible localmente establecemos los siguientes valores para los parámetros: s=b= 10.(4.17) Ahora, para ajustar el valor absoluto de la pendiente igual a uno, fijamos: β= 10.(4.18) Elegimos el valor de εen 0.1para facilitar la visualización de los resultados posteriores y asegurarnos que nos situamos en el caso nodo-nodo (β2−4εs2>0y β2−4εb2>0) y para el modelo de las cuatro zonas (4.15), además se impone que γ= 0.01 (γ2−4ε < 0), para que la órbita que nazca tras la bifurcación de tipo Hopf de este modelo sea estable. De esta forma, la variedad crítica del modelo (4.9) y su flujo correspondiente tienen la forma que se aprecia en la Figura 4.10 y la variedad crítica del modelo (4.15) y su flujo correspondiente tienen la forma que se aprecia en la Figura 4.11. Figura 4.10: Varidad crítica (en línea contínua) del modelo (4.9) con las condiciones (4.17) y (4.18) junto a su flujo (flechas). 79 Tras una primera inspección, las variedades críticas parecen cualitativa y cuantitativamente iguales. Sin embargo, al realizar una ampliación de la zona del máximo veremos que realmente difieren, como se ha comentado en la sección anterior. Figura 4.11: Varidad crítica (en línea contínua) del modelo (4.15) con las condiciones (4.17) y (4.18) junto a su flujo (flechas). Ahora, situaremos el punto inicial en (b, 0.1), así el punto inicial está a la derecha del punto de equilibrio y se puede realizar un diagrama de bifurcación. Con este objetivo, vamos a realizar tests numéricos, realizando un barrido del parámetro α desde −12 hasta −15. De esta forma, vamos a evaluar las distintas etapas de nuestro sistema ya que analizamos los valores de αpara los que el punto de equilibrio cambia de zona, i.e. αH=−13 para el modelo de las tres zonas (4.9) y αL=−13+ε/(γ−β s) yαR=−13 + ε/(γ+β b)para el modelo de las cuatro zonas (4.15). Modelo de Tres Zonas Empezamos a trabajar con el modelo (4.9). Según la teoría, nuestro modelo debe exhibir una superexplosión canard cuando el equilibrio cambia de zona (para α=−13). Utilizamos los valores de los parámetros establecidos en el apartado anterior: s=b=β= 10 para la reversibilidad local y la restricción sobre la pendiente, y ε= 0.1para mejorar la visualización de los resultados y asegurar que nos encontramos en el caso nodo-nodo. Realizando el barrido mencionado anteriormente, obtenemos el diagrama de bifurcación que se representa en la Figura 4.12. En ella se observa claramente que 80 la transición a una oscilación de relajación es instantánea, lo que se conoce como super-explosión canard. Además, se trata del mismo comportamiento que se observaba para el caso nodo del modelo (3.5) y es una propiedad claramente diferente del caso diferenciable, donde la transición es abrupta, pero no instantánea. Figura 4.12: Diagrama de bifurcación uniparamétrica de x en función de α, en él se ven el equilibrio (en negro, discontinuo) y el mínimo en x de la órbita periódica de la bifurcación(en azul, continuo). Cabe destacar que la oscilación de relajación del sistema (4.9) no corresponde con una oscilación de relajación clásica. Para el sistema de Van der Pol, la oscilación de relajación tiene cuatro partes bien diferenciadas, dos lentas y dos rápidas, mientras que en nuestro caso tenemos una zona lenta y otra rápida. Con el objetivo de ilustrar esta comparativa realizamos un test numérico en las mismas condiciones que en las simulaciones anteriores para α1=−14, para estar en la etapa de oscilación de relajación, y representamos la variable xfrente al tiempo. Comparamos con la obtenida de forma teórica para la ecuación de Van der Pol en la Figura 4.13, y apreciamos claramente que ambas órbitas de relajación son cualitativamente diferentes, como hemos mencionado anteriormente. La forma de la "nueva" oscilación de relajación (concretamente su parte rápida) es la causante de la pérdida de estabilidad del avión y veremos que es similar para todos los casos que vamos a estudiar a continuación. 81 Figura 4.13: Comparativa de la oscilación de relajación del modelo PWL aerodinámico (izquierda) con la oscilación de relajación teórica para el sistema de Van der Pol (derecha). Concluimos que nuestro modelo PWL de tres zonas (4.9) reproduce una dinámica similar a la del sistema de Van der Pol (3.5), en particular, el fenómeno de la superexplosión, como estábamos esperando. Pasamos a trabajar con el modelo PWL de cuatro zonas (4.15). Modelo de Cuatro Zonas En esta sección vamos a trabajar con el modelo (4.15), según la teoría, este modelo ya no exhibe una superexplosión canard, pasando a reproducir el fenómeno de la explosión canard que estamos buscando. Elegimos los valores de los parámetros s=b=β= 10 para la restricción sobre la pendiente, ε= 0.1para mejorar la visualización de los resultados y asegurar que estamos en el caso nodo-nodo y γ= 0.01 para que la órbita que nazca tras la bifurcación de tipo Hopf sea estable. Realizamos el mismo barrido que en el apartado anterior para obtener el diagrama de bifurcación en este caso (ver Figura 4.14). A primera vista, los diagramas de bifurcación de las Figuras 4.12 y 4.14 parecen cualitativa y cuantitativamente idénticos. Sin embargo, al realizar un zoom a una distancia εdel valor del parámetro αpara el que se da la superexplosión o la explosión canard, respectivamente, se aprecian claramente las diferencias entre ambos diagramas (ver Figura 4.15). En el caso de las cuatro zonas, la caída es abrupta pero no instantánea como lo es en el caso de las tres zonas. Antes de la caida brusca hay un periodo en el que la curva del mínimo se aleja de la curva del equilibrio pero se mantiene cerca, esto se debe que la órbita periódica nace en una bifurcación tipo Hopf supercrítica. El comportamiento que tenemos es similar al de la Figura 3.10: tras la bifurcación, el ciclo crece de forma lineal, hasta que, a orden εde la bifurcación, se da la explosión canard. Esta es una propiedad del sistema diferenciable (ver Figura 4.4) y se trata justo de lo que estábamos buscando, conseguir reproducir diagramas de bifurcación similares a los 82 obtenidos en el caso diferenciable con nuestro modelo PWL. Podemos decir que la simulación numérica reproduce los resultados teóricos esperados. Figura 4.14: Diagrama de bifurcación uniparamétrica de x en función de α, en él se ven el equilibrio inestable (en negro, discontinuo) y el mínimo en x de la órbita periódica de la bifurcación(en azul, continuo; antes de la bifurcación, a la derecha respresenta el equilibrio estable). Figura 4.15: Ampliación de la zona de la explosión de los diagramas de las Figuras 4.12 y 4.14. En la siguiente sección vamos a aumentar la dificultad eliminando la restricción sobre la pendiente, lo que nos permite acercarnos un poco más a la situación real. 83 4.4.2. Caso Reversible con Pendiente distinta de uno En esta sección vamos a realizar tests numéricos con el objetivo de reproducir los resultados predichos por la teoría cuando nos encontramos en los casos foco-foco y nodo-nodo. Además para este último vamos a intentar visualizar numéricamente ("cazar") una órbita canard. Para la realización de estos tests numéricos vamos a seguir eligiendo los valores de los parámetros s=b= 10,ε= 0.1yγ= 0.01 por los motivos descritos en la sección anterior. Para variar el valor de la pendiente vamos a modificar el valor del máximo de la nulclina rápida, β, como sigue: 1. En primer lugar, fijamos β= 12 >10. De esta forma, tendremos un valor absoluto de las pendientes mayor que 1 y seguiremos en el caso nodo-nodo (β2−4εs2>0yβ2−4εb > 0). Esperamos unos resultados cualtitativamente similares a los del caso con pendiente uno, en el que también estábamos en el caso nodo-nodo. 2. En segundo lugar, fijamos β= 2.73. De esta forma, tendremos un valor absoluto de las pendientes menor que 1 y saldremos del caso nodo-nodo para entrar en el caso foco-foco (β2−4εs2<0yβ2−4εb < 0). A diferencia del apartado anterior, en lugar de situar el punto incial en un punto fijo y realizar un barrido para el parámetro α, realizaremos una representación del plano de fases para valores del parámetro αrepresentativos con órbitas que pasen por unos puntos determinados ((0, β/2),(0,0) y(0,−β)). De esta forma, tendremos una imagen más completa de la dinámica de los modelos (4.9) y (4.15), lo que complementa a los resultados obtenidos para el caso con pendiente igual a uno. Empezaremos estudiando el modelo de tres zonas (4.9), del que discutiremos ambos casos (nodo-nodo y foco-foco), para luego estudiar únicamente el caso nodonodo del modelo de cuatro zonas (4.15). El caso foco-foco del modelo de cuatro zonas (4.15) no se estudia ya que no es representativo para la explosión canard que queremos estudiar y las cualidades principales de la dinámica son semejantes al caso foco-foco para el modelo de tres zonas (4.9). Vamos a realizar los tests numéricos para valores representativos de cada etapa de la explosión (superexplosión en el caso de las tres zonas): para la etapa previa a la (super)explosión tomaremos α=−12, y para la etapa posterior, en la que existen las oscilaciones de relajación, tomaremos α=−14. Además, analizaremos el valor del parámetro α=−13 en el modelo de las tres zonas y algunos valores situados a menos de una distancia εdel valor del parámetro α=−13 para realizar una comparativa de la dinámica de ambos modelos. 84 Modelo de Tres Zonas Empezamos estudiando el modelo de las tres zonas (4.9) en el caso nodo-nodo con β= 12. Esperamos volver a obtener la superexplosión canard. Figura 4.16: Plano de fases del sistema (4.9) para α=−12 en el caso nodo-nodo. Para el caso previo a la superexplosión (con un valor representativo α=−12), el equilibrio real se sitúa en la zona lineal derecha. Este equilibrio es estable y es un atractor de todo el flujo. Esperamos que las trayectorias tengan al principio una etapa de acercamiento rápido a la variedad lineal lenta para posteriormente seguirla relativamente cerca hasta llegar a un entorno del punto de equilibrio. Todo esto se aprecia efectivamente en la Figura 4.16. Figura 4.17: Plano de fases del sistema (4.9) para α=−13 en el caso nodo-nodo. Para el valor del parámetro para el que se da la superexplosión (α=−13), esperamos tener un continuo de órbitas homoclinas como efectivamente se aprecia en la Figura 4.17. Debido a la existencia de una zona lineal extra (respecto al sistema (3.4)) en la parte izquierda del sistema, el contínuo está acotado por una homoclina que atrae el resto del flujo del sistema, i.e. es un atractor del sistema. Esta órbita será la candidata para convertirse en la oscilación de relajación en la siguiente etapa. 85 3. Repetimos el proceso hasta que la órbita siga la variedad repulsiva de forma que sea visible en nuestra representación. De esta forma estamos reduciendo, con cada iteración, el intervalo del parámetro para el que se da una órbita canard. En las Figuras 4.27 y 4.28 vemos la inicialización del proceso, correspondientes al Caso Pre y Caso Post. Figura 4.27: Órbita de la solución del sistema (4.15) (en negro, más grueso). Valor inicial (0, β −2ε)yαa=−12.995. Se muestra la variedad crítica en azul (menos grueso). Figura 4.28: Órbita de la solución del sistema (4.15) (en negro, más grueso). Valor inicial (0, β −2ε)yαb=−13.000. Se muestra la variedad crítica en azul (menos grueso). Calculamos el parámetro αm, punto medio de αayαb, resolvemos el sistema y pintamos la órbita obtenida. Figura 4.29: Órbita de la solución del sistema (4.15). Valor inicial (0, β −2ε)y αm=−12.9975. Figura 4.30: Órbita de la solución del sistema (4.15). Valor inicial (0, β −2ε)y αm=−12.99875. En la Figura 4.29 podemos observar que todavía nos encontramos en el Caso Pre, por tanto, actualizamos αacon αmy repetimos el proceso. En este caso, vemos 92 en la Figura 4.30 que nos encontramos en el Caso Post, entonces ahora actualizamos αbcon αm. Repitiendo este proceso de forma sucesiva, llegamos después de 8 iteraciones aαm=−12.998642578125. Para este valor, en la Figura 4.31 podemos apreciar claramente como la órbita representada sigue la variedad de la izquierda, que es repulsiva, durante un tiempo, que es la característica principal de las órbitas canard. Es decir, hemos conseguido visualizar una órbita canard como predecía la teoría desarrollada. Figura 4.31: Órbita de la solución del sistema (4.15). Valor inicial (0, β −2ε)yαm= −12.998642578125. Caso con Pendiente de los lados Izquierdo y Derecho distintas Terminamos este capítulo con un comentario sobre el caso en el que la pendiente de la variedad lineal izquierda, kL, y la pendiente de la variedad lineal derecha, kR, son distintas en valor absoluto (no reversible). En este caso tendríamos un diagrama de bifurcación más complejo con la posible aparición de más de una órbita canard, [15]. Ese escenario requiere de un estudio que queda fuera de los objetivos del presente trabajo, centrado en la explosión canard de una única órbita periódica. 93 94 Bibliografía [1] A. Andronov, A. Vitt, S. Khaikin, Theory of Oscillators, 1966, Pergamon Press. [2] M. Desroches, E. Freire, S.J. Hogan, E. Ponce y P. Thota, Canards in PWL Systems: Explosions and Super-explosions, 2013, Proceedings of the Royal S. [3] M. Desroches, A. Guillamon, E. Ponce, R. Prohens, S. Rodrigues, A. E. Teruel Canards, Folded Nodes, and Mixed-Mode Oscillations in Piecewise-Linear Slow- Fast Systems, 2016, Society for Industrial and Applied Mathematics. [4] M. Desroches et al., Mixed-Mode Oscillations with Multiple Time Scales, 2012, Society for Industrial and Applied Mathematics. [5] M. di Bernardo, CJ. Budd, AR. Champneys, et al. Piecewise-Smooth Dynamical Systems, 2008, Springer-Verlag London. [6] E. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, 1967, The MIT Press. [7] E. Izhikevich, Neural Excitability, Spiking, and Bursting, 1999, International Journal of Bifurcation and Chaos, 10(6):1171-1266. [8] N. Fenichel, Geometric Singular Perturbation Theory for Ordinary Differential Equations, 1979, J. Differential Equations 31:53-98. [9] S. Fernández-García, Bifurcations of Periodic Orbits and Invariant Sets in Piecewise Linear Dynamical Systems, 2012, Tesis Doctoral, Universidad de Sevilla. [10] S. Fernández-García, M. Desroches, M. Krupa y A.E. Teruel, Canard Solutions in Planar Piecewise Linear Systems with Three Zones, 2016, Dynamical Systems 31(2):173-197. [11] FitzHugh R. Impulses and Physiological States in Theoretical Models of Nerve Membrane, 1961, Biophysical J. 1:445–466. 95 [12] J. Hale y H. Kolçak, Dynamics and Bifurcations, 1991, Springer-Verlag New York. [13] M. Komuro y T. Saito, “Lost solution” in a PWL system, 1991, IEICE TRANS. [14] M. Krupa, P. Szmolyan, Extending Geometric Singular Perturbation Theory to Nonhyperbolic Points–Folds and Canard Points in Two Dimensions, 2001, SIAM J. Math Anal. 33(2):286-314. [15] M. Krupa, P. Szmolyan, Relaxation Oscillation and Canard Explosion, 2001, J. Differential Equations 174(2):312-368. [16] J. Llibre, E. Núñez y A.E. Teruel, Limit Cycles for Planar PWL Differential Systems, 2002, Qual. Theory Dyn. Syst. [17] J. Rankin, M. Desroches, B. Krauskopf y M. Lowenberg, Canard Cycles in Aircraft Ground Dynamics, 2011, Nonlinear Dynamics 66:681-688. [18] J. Rankin, B. Krauskopf, M. Lowenberg, E. Coetzee, Operational Parameter Study of Aircraft Dynamics on the Ground, 2010, J. Computational Nonlinear Dynamics 5, 021007. [19] S. Strogatz, Nonlinear Dynamics and Chaos With Applications to Physics, Biology, Chemistry, and Engineering, 1994, Perseus Books. [20] Departamento de Análisis Matemático, Apuntes de la Asignatura DFVV, 2015/16, Universidad de Sevilla. [21] Departamento de Ecuaciones Diferenciales y Análisis Numérico, Apuntes de la Asignatura AED, 2014/2015, Universidad de Sevilla. [22] Departamento de Ecuaciones Diferenciales y Análisis Numérico, Apuntes de la Asignatura EDO, 2013/14, Universidad de Sevilla. [23] Departamento de Ingeniería de Sistemas y Automática, Notas del Curso "Análisis de Sistemas No Lineales", 2007, Universidad de Sevilla. 96