Estimación de regiones de estabilidad en un modelo de
péndulo elástico

Stability zone estimation in an elastic pendulum
model

Nestor Rodrigo Condori U. & Diego Sanjinés Castedo ]]> Carrera de Física Universidad Mayor de San Andresc. 27 Cota-Cota, Casilla de Correos 8635
La Paz, Bolivia
Recibido: 28 de julio de 2019    aceptado: 15 de septiembre de 2019

Resumen

Se investiga el sistema dinamico del péndulo elástico con dos grados de libertad (radial y angular) a traves del formalismo de Lagrange que conduce a un sistema de dos ecuaciones diferenciales no-lineales acopladas. Este sistema se transforma en una ecuacion de Hill para la proyección horizontal de la coordenada radial, cuyo coeficiente periodico en el tiempo se modela por una función constante por tramos a fin de aplicar un criterio de estabilidad (traza de la matriz de evolucion); este criterio permite obtener numericamente un diagrama de estabilidad para los parámetros característicos del sistema sin invocar la aproximacion de pequeñas oscilaciones. Así, la simulacion de la solucion del sistema de ecuaciones diferenciales con los parámetros elegidos según dicho diagrama de estabilidad permite verificar la precision de su cálculo, encontrándose resultados bastante buenos. Para pequenas oscilaciones se recupera los resultados conocidos de la ecuación de Mathieu y su correspondiente diagrama de estabilidad (diagrama de Strutt). Finalmente se discute algunas perspectivas interesantes.

Codigo(s) PACS: 05.45.-a - 45.20.Jj - 43.20.Ks

Descriptores: Dinamica no-lineal - Mecánica lagrangiana y hamiltoniana - Conservación de la energía

Abstract

Subject headings: Nonlinear dynamics - Lagrangian and Hamiltonian mechanics - Energy conservation

1. INTRODUCCIÓN

El problema del pendulo elástico ha sido ampliamente estudiado en regiones de oscilaciones pequenas (ángulos y amplitudes despreciables) con lo cual el sistema no-lineal de ecuaciones diferenciales ordinarias acopladas se hace mas sencillo de tratar matemáticamente. Los primeros que aplicaron el modelo del pendulo elástico a un problema fueron Vitt & Gorelik (1933); ellos propusieron el modelo para describir la resonancia parametrica (esto es, el fenómeno de transferencia de energía entre dos modos de oscilacion) en la molécula del dioxido de carbono y hallar sus líneas espectrales. Simplificando el problema y usando la teorıa de perturbaciones Gorelik y Vitt encontraron soluciones periodicas aproximadas. En 1976 Olsson (Olsson 1976) en estudio el problema del pendulo elástico en un sistema cartesiano; realizo un desarrollo de Taylor del lagrangiano y desprecio términos superiores bajo la aproximación de oscilaciones pequenas llegando así a la ecuacion de Mathieu para la proyeccion de la longitud del resorte en el eje horizontal. Olsson describio también la resonancia parametrica y concluyó que no hay soluciones analíticas expresables mediante el principio de superposicion en funciones armonicas simples de los modos de oscilación. En 2002 Peter Lynch (Lynch 1996; Lynch 2002) trato el problema del pendulo elástico en dos y en tres dimensiones en coordenadas cartesianas, cilindricas y esfericas bajo la aproximación de pequeñas oscilaciones. Aunque Lynch deriva la ecuacion de Mathieu, no se preocupa tanto por las regiones de estabilidad, si no mas bien por hallar las expresiones analíticas de las soluciones en funcion al tiempo. En 2005 Sanjinés (Sanjinés 2005) estudio la relación entre las soluciones del péndulo elastico (de nuevo, bajo la aproximación de pequeñas oscilaciones) y su aplicacion a la teoría de bandas en un modelo unidimensional de solido cristalino; así, una vez que relaciona la ecuacion de Mathieu con la ecuación de Schrodinger, utiliza el criterio de la traza para hallar el diagrama de Strutt (que indica las regiones de estabilidad e inestabilidad) y relaciona las regiones estables con un estado físico de aislante electrico y las regiones inestables con un estado físico de conductor eléctrico.

La determinacion de las regiones de estabilidad e inestabilidad en el pendulo elástico tiene diversas aplicaciones, por ejemplo, en dinamica atmosférica (Lynch 2002), construccion de barcos, aeroplanos, sistemas roboticos, pues se requiere un estrecho rango de variacion de parámetros característicos para el correcto funcionamiento de estos dispositivos. En este trabajo u-saremos la teorıa de Floquet (Sanjinés 1990; Magnus & Winkler 1966) para calcular los parametros para los cuales hay estabilidad en las soluciones de una ecuacion diferencial con coeficientes periodicos; los resultados permitiran elaborar el diagrama de estabilidad correspondiente sin invocar la aproximacion de pequeñas oscilaciones.

La organizacion de este artículo es la siguiente. En la Seccion II se describe la dinámica del péndulo elástico en coordenadas polares sin invocar la aproximacion de pequenas oscilaciones; a partir del sistema acoplado de ecuaciones diferenciales no-lineales se deduce una ecuacion diferencial de segundo orden en variable compleja, lo que es posible gracias a un cambio de variable muy oportuno que no es trivial. En la Seccion III se supone que la coordenada radial r(t) varía de manera armonica, por lo que la ecuación de movimiento del pendulo se convierte en una ecuación de Hill (esto es, una ecuacion de segundo orden con un coeficiente periodico); en el límite de pequenas oscilaciones se define las constantes apropiadas para realizar luego la comparacion con una ecuación de Mathieu. En la Seccion IV se aplica la teoría de Floquet para segmentar el coeficiente periodico de la ecuación de Hill y elaborar el diagrama de estabilidad, que en el lımite de pequenas oscilaciones corresponde al conocido diagrama de Strutt, lo que es consistente con el hecho de que en ese lımite se obtiene la ecuación de Mathieu a partir de la ecuacion de Hill. En la Sección V se realiza las simulaciones numericas de las ecuaciones de movimiento e-xactas del sistema (i.e., sin aproximar el coeficiente de la ecuacion de Hill por una función constante por tramos); se obtiene un acuerdo bastante bueno entre las predicciones analıticas realizadas para el modelo aproximado y las simulaciones numericas para el sistema exacto. Finalmente, en los Apendices se trata los temas complementarios sobre: las aproximaciones para la funcion de la traza y el diagrama de estabilidad (Apéndice A), la ecuación radial (Apéndice B) y el caos en el péndulo elástico (Apéndice C).

2. EL PÉNDULO ELÁSTICO

En general el pendulo elástico es un problema tridimensional cuyo movimiento con tres grados de libertad es complicado. Se ha visto que ocurren fenomenos tales como la flexion del resorte o la alternancia de regımenes de orden y caos, y aunque hay simulaciones numericas de estos fenómenos, aún no se ha conseguido una solucion analítica completa a mas de 300 años desde la deduccion de sus ecuaciones de movimiento (Fitch 2009; Van Der Weele & De Kleine 1996). Como aproximaciones adicionales a las ya referidas, se supondra en este trabajo que el pendulo oscila en un plano (dos grados de libertad) y que la coordenada radial depende del tiempo de manera armónica.

A continuacion deduciremos las ecuaciones de movimiento del sistema (Fig. 1) para las coordenadas radial r(t) y angular Φ(t). La energía cinetica es

y la energía potencial es

así que la energía total es

donde r₀ es la longitud de relajamiento del resorte, m es la masa puntual de la carga y k es la constante elástica.

La función lagrangiana del sistema es L = T - U:

así que de las ecuaciones de Euler-Lagrange se obtiene:

Este es un sistema no-lineal de ecuaciones diferenciales acopladas que no tiene una solucion analítica cerrada, pues posee solo una constante de movimiento que es la energía total; sin otra constante de movimiento no es posible resolver analíticamente el sistema (Van Der Weele & De Kleine 1996). Multiplicando (6) por cos Φ, (7) por sen Φ, y sumando ambas ecuaciones, se obtiene:

De manera similar multiplicamos (6) por -sen Φ y (7) por cos Φ, de donde obtenemos:

Aunque parece que las nuevas ecuaciones (8) y (9) no son mas sencillas que (6) y (7) podemos sin embargo trabajar con ellas. Para ello realizamos el cambio de variable ξ= re^iΦ,con ello el sistema (8) y (9) se reduce a:

Esta es una ecuación diferencial de segundo orden en variable compleja, cuya parte real es (9) y cuya parte imaginaria es (8). Para dilucidar las propiedades dinamicas del sistema se puede graficar r(Φ) o bien proponer una función r(t) adecuada y resolver (10) para Φ(t). La ecuacion para ε se construye tomando el complejo conjugado de (10):

y multiplicando por ξ* para obtener:

Analogamente, multiplicamos (11) por —ξ:

y se suma (12) y (13):

Esta ecuacion para ξ y su complejo conjugado equivale a (6). De manera similar, restando (12) de (13) se obtiene el equivalente de (7). Vemos pues que el pendulo elástico se describe efectivamente por (10) ya que de aquí se deducen (6) y (7).

3. ECUACIÓN DE HILL

En la ecuacion (10) es posible desarrollar (en principio) la funcion periódica r(t) en una serie de Fourier, por lo que (10) tendra un coeficiente periódico. Ésta es la ecuacion de Hill (Magnus & Winkler 1966). No obstante no seguiremos ese procedimiento en este trabajo sino que se propondra un modelo simple para r(t):

donde es la amplitud de oscilacion (que no se considera pequena) y ω, como se vera, es la frecuencia de oscilación de la coordenada radial. Reemplazando (15) en (10) se obtiene:

Esta ecuacion se reduce a una ecuación de Mathieu para el caso de pequenas oscilaciones. Tal reducción nos permitira relacionar las variables de los ejes del diagrama de Strutt con las correspondientes del diagrama de estabilidad para nuestra ecuacion de Hill (sección V).

A continuacion proyectamos ξ(t) en coordenadas cartesianas:

derivando dos veces obtenemos:

de donde se obtiene un sistema de ecuaciones diferenciales acopladadas al sustituir (17) y (18) en (16):

Si hacemos a 0 en la ecuacion de Hill (20) y desarrollamos el denominador en serie para alrededor de cero se obtiene:

Esta es la ecuación de Mathieu que posee como coeficiente periodico una función armónica en el numerador. Comparando con las halladas en Olsson (1976); San-jines (2005) se ve que poseen la misma estructura y las mismas constantes, con la siguiente identificacion: δ = g/p y λ = aω²r₀ /p² = kr₀ /mp² (δ y λ son constantes); así, identificamos ω² = k/m y se prueba que nuestra ecuacion de Hill se reduce a una ecuación de Mathieu. De manera equivalente, este resultado se obtiene tambien multiplicando el segundo término de (20) en el numerador y denominador por
p + cos ωt:

Si a tiende cero se desprecia el termino cuadrático y se obtiene nuevamente (21). Ahora que hemos encontrado la constante ω, volvamos a la ecuacion de Hill (20). Factorizando p en el denominador del coeficiente, obtenemos dicha ecuacion donde se introducen las constantes α y β:

donde

En adelante identificaremos por Q(τ) al coeficiente periodico de la ecuación de Hill (23).

4. APROXIMACION DE LA FUNCIÓN Q(τ)

Consideremos a continuacion la función Q(τ):

con n un numero entero y

En la Fig. 2 (a) y 2 (b) se verifica un comportamiento periodico (cosenoidal) cuando se dan las condiciones expresadas arriba; si dichas condiciones no se cumplen entonces Q(τ) es singular (Fig. 2 (c)) o posee una forma periodica complicada. A continuación segmentemos Q(τ) modelandola por la siguiente función constante por tramos

La condicion (27) parecería restringir el valor de a de tal forma que no sería posible determinar las regiones de estabilidad sin la aproximacion de oscilaciones pequenas pues posee una cota superior. Veamos que este no es el caso. A partir de (15) se tiene que, para los valores extremos de cos ωt, r = r₀+ mg/k ± , con r₀ + mg/k (condicion (27)), así que:

5. FUNCION DE LA TRAZA Y DIAGRAMA DE ESTABILIDAD

Para elaborar el diagrama de estabilidad seguimos el criterio de la traza que se deduce de la teoría de Floquet (Sanjines 2005, 1990; Magnus & Winkler 1966). Expresamos (23) en forma matricial:

El vector X definido así en (32) contiene a la funcion x(τ) y a su derivada. La solucion formal para X(τ) se puede expresar en terminos de la condición inicial X(τ₀) y de la matriz de evolución T(τ, δτ₀):

Combinando (33) y (32) se obtiene la ecuacion diferencial:

La matriz de evolucion (matriz de Floquet) es unimo-dular y consiste de un mapeo biparametrico o aplicación lineal que posee las propiedades de grupo abeliano; esta matriz cumple la propiedad de composicion:

Esta propiedad será útil para calcular la matriz de evolucion resultante de una región segmentada donde en cada intervalo (segmento) la integral (35) se pueda calcular de manera exacta (Magnus & Winkler 1966; Sanjines 2002). Luego, se aplica la condición de estabilidad usual en (33) con T(τ,τ₀) ya calculada: X₁(τ) y X₂(τ), con condiciones iniciales X₁(0) y X₂(0), son soluciones estables de (32) si para todo ε > 0 existe un δ > 0 de tal forma que si | X₁(0) - X₂(0) |< δ, entonces | X₁(τ) — X₂(τ) |< ε (τ > 0). Así, las regiones de estabilidad para los parametros α y β dados en (24) se determinan si la funcion de la traza S(α, β) ≡ Tr(T(L,0)) (-2, 2), donde L es el periodo de la matriz de evolucion; este es el criterio de la traza (Sanjinés 2005, 2002). En estas regiones la norma de X(τ) se mantiene acotada. Vale la pena mencionar que el criterio de la traza utilizado en este trabajo lleva al mismo resultado que el obtenido por el algoritmo de Meissner (Magnus & Winkler 1966). Ahora podemos aplicar el criterio de la traza a la función Q(τ) dada en (28). El resultado para S = Tr(T(τ₂ , τ₁)T( τ₁ , τ₀)) es:

con T(τ_i , τ_j) calculada segun (35) para los valores constantes de T(z) en (32) que corresponden a Q(τ) dado en (28). Si el valor de a es muy pequeno la ecuación (37) corresponde a las regiones de estabilidad para oscilaciones pequenas, lo que equivale a α ≈ 0 en el diagrama de estabilidad (apendice A). En la Fig. 3 se muestra el diagrama de estabilidad obtenido con la condición | S | ≤ 2.

Ahora volvamos al sistema (6)-(7) y trabajemos con las siguientes condiciones iniciales:

donde se ha utilizado (24) para relacionar β con la condicion inicial radial; r₀ se mide en metros; el angulo inicial es 1^o = 0.017 rad. (e iguales resultados se obtiene para Φ(0) < 1.53rad.) El regimen de inestabilidad corresponde al fenomeno de "resonancia paramétrica" definido como el regimen en el que la transferencia de energía es optima entre los modos radial y angular (Ols-son 1976); esta condicion ocurre cuando la frecuencia radial es el doble de la frecuencia angular, lo que corresponde a β = (π/2)² y, de manera equivalente, a que el cociente de las longitudes (en reposo) del resorte cargado y sin cargar sea r₀ /p = 4/3 (Sanjines 2005). En efecto, se puede verificar en la Figs. 4 (a) y 4 (b) la inestabilidad del sistema (resonacia parametrica). En las Figs. 4 (c) y 4 (d) se muestra casos del regimen de estabilidad donde los modos de oscilacion no se alternan, es decir, la energıa se mantiene en cada modo. Un caso extremo de esta situacion ocurre cuando Φ(0) = 0: la energía del sistema de mantiene en el modo radial y el modo angular no se excita pues nunca recibe energıa del modo radial, por lo que Φ(t) = 0 para t > 0. La simulacion numerica de este caso es trivial por lo que no se incluye en los casos de la Fig. 4; de todas maneras se realizo a fin de verificar la precision numérica del algoritmo para resolver el sistema (6) - (7).

Los graficos anteriores, que se obtuvieron resolviendo numericamente el sistema de ecuaciones (6)-(7), corresponden a los parametros α y β extraıdos del diagrama de estabilidad para valores arbitrarios de las coordenadas radial y angular.

Recordemos, sin embargo, que este diagrama de estabilidad es aproximado pues se elaboro con base en una función Q(τ) constante por tramos (seccion IV), pero para efectos practicos resulta bastante bueno. Debemos indicar que la condicion de resonancia paramétrica (β = (π/2)²) no es unica, pues otros valores de β corresponden al regimen de inestabilidad, sin embargo podemos afirmar cualitativamente que β = (π/2)² hace que el sistema trasfiera mas eficientemente la energía entre los modos radial y angular. Cuando β = (π/2)² no se puede dar la elongacion máxima al péndulo pues la masa choca con el punto de suspension; así, el valor máximo de elongacion inicial es r₀ + con mg/k . Para β = (2π/5) no ocurre tal choque y por lo tanto r₀ + mg/k puede tomar este valor maximo, lo que es consistente con el diagrama de estabilidad (ver Figs. 3, 4 (e), 4 (f), 4 (g) y 4 (h)). El diagrama de estabilidad aproximado que se obtuvo en este trabajo es bastante preciso en las regiones estables β > (π/2)² pero menos preciso en β < (π/2)². Notese que para β > α las areas de estabilidad son mas grandes que para β < α, lo que significa que los resortes con longitudes de relajacion pequeñas y valores grandes de la constante elastica tienden a ser inestables cuando son sometidos a elongaciones (o compresiones) grandes.

7. PENDULO ELÁSTICO AMORTIGUADO

Consideremos las ecuaciones de Euler-Lagrange para un sistema disipativo:

donde funcion de disipación de Rayleigh que corresponde a una fuerza disipativa proporcional a la velocidad (Goldstein 1980). Para el caso del pendulo elástico en coordenadas polares se obtiene el siguiente sistema de ecuaciones:

Para coeficientes de friccion iguales «i = n₂ = n las dos ecuaciones anteriores se expresan de manera compacta (con ξ = re^iΦ):

donde se han realizado operaciones analogas a las de la seccion II con un coeficiente de disipación constante. La anterior ecuacion es una ecuación de Hill con r(t) periodica según (15). En coordenadas cartesianas se obtiene:

Así, el movimiento en el eje x obedece la ecuacion:

que con el cambio de variable x(t) = e ^-Kt/2x(t) conduce a:

donde se definio el coeficiente de disipasión = kπ/ωα. Esta ecuacion posee la misma estructura que en el caso del movimiento del pendulo sin disipación, por lo tanto se puede repetir el criterio de la traza de la teoría de Floquet y obtener la traza de matriz de evolucion y luego el correspondiente diagrama de estabilidad. Así, (46) se reescribe (con ωt = πτ) como :

donde α y β tiene el mismo significado que en el caso del pendulo elástico sin fricción. Luego, la traza de la matriz de evolucion resulta ser:

donde:

8. CONCLUSIONES Y PERSPECTIVAS

En este trabajo se elaboro un procedimiento aproximado para hallar las regiones de estabilidad de un modelo de pendulo elástico con dos grados de libertad; no se invoco la aproximación de pequeñas oscilaciones, así que los resultados de este trabajo son una generalizacion de aquellos que se presentaron en Sanjines (2005). La formula (10) y (42) son novedosas y comprenden las ecuaciones dinamicas del péndulo elástico para el caso general de oscilaciones de amplitud arbitraria. Se logro aplicar el criterio de la traza de la matriz de evolucion a la ecuacion de Hill (deducida a partir de (10)) por medio de la segmentacion del coeficiente periódico Q(τ) en dos partes y de aquí se obtuvo el correspondiente diagrama de estabilidad (Fig. 3) que es una generalizacion del conocido diagrama de Strutt. En el límite de pequenas oscilaciones se recupera todos los resultados conocidos. Aunque la segmentacion elegida para Q(τ) ciertamente es burda, el diagrama de estabilidad (Fig. 3) es razonablemente bueno, al menos para los puntos senalados en los que el comportamiento dinamico del péndulo es el que se espera segun (10).

Se observa que las regiones de estabilidad se encogen a medida que aumenta la amplitud de las oscilaciones radiales hasta llegar a desaparecer completamente (en el regimen de resonancia paramétrica), lo que también se observa en el diagrama de Strutt, aunque este ultimo no permite una interpretacion válida para la estabilidad de oscilaciones de gran amplitud. En el diagrama de estabilidad obtenido en este trabajo se observa que no existe estabilidad para el valor β = (π/2)², lo que corresponde al "pendulo elástico 4/3" reportado en Sanjinés (2005), lograndose así una interpretacion física consistente del fenomeno de resonancia paramétrica para este sistema. Considerando que la elongacion máxima posible (en el modo radial) corresponde a la longitud de equilibrio del resorte cargado, se observa la predominancia de regiones de inestabilidad, salvo los casos extremos en donde la masa del sistema es muy pequena o la constante del resorte es muy grande.

Una perspectiva interesante de este trabajo es refi-nar la dependencia temporal de r(t) dada por (15), en vista de los resultados de las Figs. 4(a) y 4(b), donde se aprecia una envolvente en las oscilaciones de r(t) con un periodo T = 2π/ω' que corresponde a la transicion entre los regímenes radial y angular. Se podria tratar por ejemplo, r(t) = p — α cos(ωt) (cos ω')t. Este refinamiento, así como una segmentacion más fina de Q(τ), seguramente mejoraran la precisión del diagrama de estabilidad (Fig. 3). Notemos, sin embargo, que a lo largo de la "lengua de estabilidad" (en la que se encuentra el punto (α,β) = (5,0.76)) la amplitud de las oscilaciones de r(t) se mantiene aproximadamente constante pues no se tranfiere energía al modo angular; esto es lo que se observa en las Figs. 4(c) y 4(d). Por lo tanto el refinamiento en el modelo de r(t) dado por (15) no afectara a la envolvente de las oscilaciones pero probablemente modificara su frecuencia, incluyendo el eventual caso de caos (Apendice C).

Por otra parte, de manera inesperada se observo un comportamiento interesante (Figs. 4(f), 4(h)): en el regimen de resonancia paramétrica (inestabilidad) correspondiente a β= (π/2)², para un angulo inicial del orden de 10~²⁰ rad (!) el pendulo tarda un tiempo apre-ciable de varios segundos (que es mayor mientras menor es el angulo inicial) en cambiar del modo de oscilación radial al modo angular; cuando el angulo inicial es cero el sistema aparentemente no cambia al modo angular, al menos dentro del intervalo temporal considerado. El que este fenomeno sea un resultado solamente de la simulacion computacional o pueda corresponder a un fenomeno físico bona fide queda como un tema abierto para una investigacion futura.

Conflicto de intereses

Los autores declaran que no hay conflicto de intereses con respecto a la publicacion de éste documento.

A. APROXIMACIONES PARA LA FUNCION DE LA TRAZA Y EL DIAGRAMA DE ESTABILIDAD

Hagamos el cambio de variable πτ = ωt en la ecuacion de Mathieu (21):

con q(τ) = (β — v cos πτ), β = gπ²/pω² y v = π²r₀ω²/p²ω². Si aproximamos q(τ) por una funcion constante por tramos:

se obtiene la siguiente funcion de la traza:

cuya grafica se muestra en la Fig. A1 y fue tomada de Sanjines (2005).

A continuacion veamos el caso particular de la función de la traza (37) para la ecuacion de Mathieu (caso de pequenas oscilaciones). Notemos que el eje de las ordenadas es el mismo en ambos diagramas pero no así el eje de las abscisas; manipulemos entonces el parametro a de (37) para lograr la coincidencia de las escalas de ese eje:

con v = r₀π²/p², donde se despreciaron potencias de a mayores a dos (según α y β dados en (24)). Se procede de manera similar para obtener el siguiente factor:

A continuacion trabajemos el término:

donde en el numerador se desprecio potencias de la amplitud mayores a dos y en el denominador se utilizo los resultados (A6) y (A7). La funcion de la traza que resulta es:

que es el mismo resultado que se obtuvo en Sanjines (2005) (Fig. A1).

El metodo de aproximar una función continua Q(τ) por una funcion constante por tramos (segmentación de la funcion) es más preciso mientras más funciones constantes por tramos se tomen, lo que mejora a su vez la precision del diagrama de estabilidad asociado a la ecuacion de Hill cuyo coeficiente periódico es Q(τ) (seccion IV). Para una secuencia infinita de funciones continuas por tramos el diagrama de estabilidad se aproxima mejor al diagrama verdadero, pero aun no lo es debido a la limitacion impuesta por la aproximación para r(t) en (15). Sin embargo, el calculo de la función de la traza a medida que aumenta la segmentacion de Q(τ) se hace demasiado complicado. En algunos casos es posible calcular directamente la funcion de la traza integrando Q(τ). Para ilustrar tal metodo obtendremos el diagrama de Strutt integrando de manera directa la matriz de transferencia a partir de la ecuacion de Mathieu (21):

que con el cambio de variable ωt= 2n adopta su forma estandar:

donde α' = 2r₀/p² y β' = 4g/pω². De manera equivalente, la solucion de (A11) se expresa en forma matricial como se hizo en (31)-(35):

con:

y con = ( β' - 2α' cos2n) ; así, el primer paso para encontrar la matriz de evolucion es resolver la integral:

donde p es el periodo, y M(α', β') es una funcion de dos variables. Para el caso de la ecuacion de Mathieu el resultado de la integracion anula la variable α', de tal manera que no es posible calcular la funcion de la traza con argumentos α' y β':

Esto ocurre por la simetrıa impar del resultado de integrar en nuestro caso particular. Una forma de evitar eso es aprovechar las propiedades de simetría de la siguiente forma:

de donde se obtiene:

Vemos pues que se recupero el parámetro α'. Ahora podemos hallar la funcion de la traza:

A continuacion apliquemos el anterior procedimiento para calcular el diagrama de estabilidad para el pendulo elastico sin la restricción de oscilaciones pequeñas. Partimos de (20):

y realizando el cambio de variables 2n/ = ωt:

donde α" = 2/ p y β' = 4g/ω²p; requerimos la integral:

donde p = π es el periodo de la funcion. Suponemos α < 2 que es la restriccion para que la longitud del péndulo elastico no sea negativa, lo que equivale a que la máxima elongacion del modo radial sea igual a la longitud de equilibrio del pendulo cargado; así, el denominador de la integral anterior se desarrolla en serie:

con lo que se obtendría -en principio- el diagrama de estabilidad. En la aproximacion a primer orden en la serie de potencias se genera la ecuacion de Mathieu y los ordenes superiores corrigen el diagrama de estabilidad. Sin embargo, para oscilaciones arbitrarias los ordenes superiores de la serie no se pueden despreciar. Se debe pues resolver analíticamente tal integral. Si esto no es factible entonces el metodo de segmentar la función continua es la mejor forma para calcular de manera aproximada el diagrama de estabilidad. En la Fig. A3 se compara el diagrama de Strutt exacto con el obtenido en Sanjines (2005).

B. ECUACION RADIAL

En este apendice veremos el tratamiento analítico para algunos casos particulares de ecuaciones radiales. La ecuacion (6) se puede re-escribir para obtener una aproximacion de Φ:

dividiendo entre r, aplicando la derivada del logaritmo e integrando se obtiene:

Si se conserva solo el primer término del desarrollo en serie de la exponencial se obtiene:

Esta aproximacion es razonable cuando la anterior integral tiene un valor pequeno, lo que puede ocurrir si Φ 0 es muy pequeno con respecto al producto en el denominador rΦ. A partir de (7) consideremos la ecuacion:

donde el caso de = 0 y Φ 0 corresponde a la propuesta para r(t) en (15) que se usa en este trabajo. Sustituyendo (B3) en (B4) se tiene que:

ya que cosΦ 1 para angulos pequeños, se obtiene una ecuacion para r:

Integrando:

( es la constante de integracion) y volviendo a integrar:

que es una integral abeliana (que hasta donde conocemos no es soluble). Si se recupera la propuesta para r(t) en este trabajo.

El paso de (B2) a (B3) es quiza poco natural, pues sólo se considera el primer termino del desarrollo del exponencial, no obstante podemos encontrar una ecuacion radial con la unica condición de Φ 0. Para tal objeto reordenemos (7):

reemplazando esta ecuacion en la ecuación de la energía (3) se obtiene:

Aproximando la anterior ecuacion para pequeños angulos y estableciendo ε = 2E/m - /m, (B10) queda como :

Esta ecuacion es general pues no se ha supuesto ángulos pequenos. Si se supone que Φ esta dada por (B3) se obtiene una ecuacion para r que no es soluble:

De las ecuaciones presentadas en esta seccion sólo (B8) es soluble cuando . Es probable que la solucion radial deba expresarse por una generalización de la integral elıptica tal como una integral abeliana (Barker 2010), sin embargo para el caso general de oscilaciones arbitrarias (amplitudes y angulos arbitrar-ios)una solucion analítica parece no existir por la mon-odromía del sistema Hamiltoniano en consideracion (Fitch 2009).

C. CAOS EN EL PENDULO ELÁSTICO

El estudio del caos en el pendulo elástico ha sido ampliamente tratado (ver, por ejemplo, Van Der Weele & De Kleine (1996); Núñez-Yépez et al. (1990); Carretero-Gonzalez et al. (1994); Anurag et al. (2019)) en vista de la cuestion natural de saber si en el régimen de resonancia parametrica (inestabilidad) determinado por β = (π/2)² 2.4674 (línea segmentada de la Fig. 3) alguna perturbacion caótica inicial se puede amplificar de la misma forma en que en dicho regimen el pequeño angulo inicial Φ(0) = 0.017 rad se amplifica hasta alcanzar valores significativos (Fig. 4). Dichos estudios reportan el fenomeno de transición orden-caos-orden como funcion de parámetros físicos fácilmente controlables, como por ejemplo, la energıa total del sistema. En un contexto mas general, el estudio del caos en el pendulo elástico está expuesto detalladamente en el trabajo de Lynch (2002), donde se invoca el teorema de Kolmogorov-Arnold-Moser (KAM) y su relacion con las condiciones iniciales lineal y no-lineal, aunque dicho estudio esta restringido a pequeñas oscilaciones. Recordemos que en nuestro trabajo estudiamos las regiones de estabilidad sin invocar la aproximacion de pequeñas oscilaciones.

Las cantidades α y β definidas en (24) son adecuadas para caracterizar la transicion orden-caos-orden pues α comprende la elongacion (o compresión) inicial = r₀ -r(0) + mg/k que se le da al pendulo en el modo radial, mientras que β comprende a los parametros propios del resorte (masa, constante elastica, longitud del resorte relajado). Así, combinando (3) y (24), la energía total E se puede expresar por

Podemos ver que E(α , β) es una funcion monótona creciente de α, por lo que se espera la transicion orden-caos-orden a medida que aumenta el valor de α sobre la línea segmentada de la Fig. 3, de manera consistente con lo reportado por lo autores citados en este apendice, en particular, por Núñez-Yépez et al. (1990). Aunque nues-tro proposito original no es investigar el caos, las simu-laciones que hicimos de r(t) y Φ(t) para diferentes valores de a en el regimen de resonancia paramétrica sugieren la transicion referida entre el orden y el caos (estas simulaciones no se muestran en nuestro trabajo); para efectos de caracterizar dichos estados de orden y caos solo acudimos a una inspección visual de los graficos r(t) y Φ(t) de la Fig. 4. Por otra parte, los casos (a), (b) y (g) de la Fig. 4 corresponden a un resorte con β = 2.47 que es un valor apenas mayor que el del regimen de resonancia paramétrica. En este caso podemos inferir que ocurre una transicion orden-caos para algun valor de α entre los valores correspondientes a los casos (a) y (b); no sabemos si ocurre una transicion caos-orden para algun valor de a mayor al que corresponde al caso (g). Otra observacion interesante de las simulaciones de la Fig. 4 es la de estados aparentemente caoticos (al menos para Φ(t) en los casos (c) y (d) que corresponden al regimen de estabilidad, lo que sugiere que el estado de caos no solamente estaría asociado al regimen de resonancia paramétrica. Estos aspectos, y seguramente otros, son posibles temas de futuras inves- tigaciones.

REFERENCIAS

Anurag, B. M., Bhattacherjee, J. K. & Chakraborty, S. (2019), Physica D (emprensa).        [ Links ]

Barker, A. L. (2010), Elliptic Functions, An Elementary Text-Book for Students of Mathematics, (Cornell University Library).        [ Links ]

Carretero-Gonzalez, R., Núñez-Yépez, H. N. & Salas-Brito, A. L. (1994), Eur. J. Phys. 15, 139.        [ Links ]

Fitch, N. J. (2009), Phys. Rev. Lett. 103, 034301.        [ Links ]

Goldstein, H. (1980), Classical Mechanics, (2nd ed., Addison-Wesley).        [ Links ]

Lynch, P. (1996), Journal of Non-Linear Mechanics 37, 345.        [ Links ]

Lynch, P. (2002), The Swinging Spring, A model for Atmospheric Balance, (Cambridge University Press: Large Scale Atmospheric-Ocean Dynamics, Vol II).        [ Links ]

Magnus, W. & Winkler, S. (1966), Hill's Equation, (Intercience Tracts in Pure and Applied Mathematics, Vol. 20, Wiley).         [ Links ]

Núñez-Yépez, H. N., Salas-Brito, A. L., Vargas, C. A. & Vicente, L. (1990), Phys. Lett. A 145, 101.         [ Links ]

Olsson, M. G. (1976), Am. J. Phys. 45, 1211.         [ Links ]

Sanjines, D. (1990), Rev. Bol. Fís. 36 (S1), 181.         [ Links ]

Sanjines, D. (2002), Rev. Bol. Fís. 8, 1.         [ Links ]

Sanjines, D. (2005), Rev. Bol. Fís. 11, 31.        [ Links ]

Van Der Weele, J. P. & De Kleine, E. (1996), Physica A 228, 245.         [ Links ]

Vitt, A. & Gorelik, G. (1971), Z. Tekh. Fiz. 3, 294.        [ Links ]