Práctica de Métodos de la Física Matemática.
Área de Física Teórica. Dpto. Física. UEx

[Graphics:Images/Duffing.nb_gr_1.gif]: Santos Bravo Yuste. email:  santos@unex.es

La ecuación de Duffing

Cuaderno Mathematica basado en el cuaderno duffing.nb de  Peter Young (htttp://bartok.ucsc.edu/peter/115)

Introducción

Sistemas con comportamiento caótico son más abundantes de lo que habitualmente se piensa.  Este comportamiento puede surgir en sistemas descritos mediante (i)  ecuaciones iterativas o en diferencias (el iterador logístico[Graphics:Images/Duffing.nb_gr_2.gif] es un ejemplo bien conocido), (ii) ecuaciones diferenciales ordinarias (ejemplos más famosos: ecuaciones de Lorentz y ecuación de Duffing) y (iii) ecuaciones en derivadas parciales (ejemplo: flujo turbulento).
En este cuaderno vamos a discutir brevemente la ecuación de  Duffing

La ecuación de  Duffing  es un oscilador no lineal que describe el  movimiento de una partícula clasica dentro del potencial
                [Graphics:Images/Duffing.nb_gr_3.gif]
con A>0 y B>0.  Este potencial se llama de doble pozo porque tiene dos mínimos simétricamente situados respecto una barrera de potencial centrada en x=0.  Vamos a escoger A=B=1 de modo que los mínimos  de valor V=-1/4 están situados en  x = ± 1 (esto es equivalente a decir que escogemos las unidades de longitud y energía de modo que el minimo está situado en   x = ± 1 y su profundiad es -1/4).

Dibujemos este potencial:

[Graphics:Images/Duffing.nb_gr_4.gif]

[Graphics:Images/Duffing.nb_gr_5.gif]

La fuerza correspondiente a este potencial es F(x) = -dV/dx = x - [Graphics:Images/Duffing.nb_gr_6.gif]  La trayectoria x(t) que sigue la partícula sometida a esta fuerza se obtiene resolviendo la ecuación diferencial de segundo orden  F = [Graphics:Images/Duffing.nb_gr_7.gif] . Esta ecuación la escribiremos como un sistema de dos ecuaciones de primer orden:   v = dx/dt,         mdv/dt = F.   En lo que sigue tomaremos la masa de la particula igual a la unidad de modo que  la ecuación a resolver es

[Graphics:Images/Duffing.nb_gr_8.gif]

Empezaremos usando la función NSolve de  Mathematica para resolver la ecuación cuando la partícula está situada sobre el mínimo de la derecha en  x=1 y tiene una velocidad inicial igual a 1. Esta velocidad es suficiente para hacer que la partícula supere la barrera de potencial  en x=0, pase al pozo de la izquierda , vuelva a superar la barrera de potencial, retorne al seno de la derecha y así continue oscilando indefinidamente.  Esto siempre se producirá  si la velocidad inicial de la partícula es  mayor que  [Graphics:Images/Duffing.nb_gr_9.gif], como es fácil deducir  si se iguala la energía de la barrera de potencial a la energía cinética incial de la partícula que parte de uno de los mínimos.

[Graphics:Images/Duffing.nb_gr_10.gif]


Representamos la posición de la partícula en función del tiempo, x(t), para t entre 0 y 10.

[Graphics:Images/Duffing.nb_gr_11.gif]

[Graphics:Images/Duffing.nb_gr_12.gif]



Otro modo muy conveniente de representar la solución es usando el plano de fases donde la abcisa es x y la ordenada es v. Nuevamente representamos la solución (trayectoria en el plano de fases) con t entre 0 y 10.

[Graphics:Images/Duffing.nb_gr_13.gif]

[Graphics:Images/Duffing.nb_gr_14.gif]

Se ve que la trayectoria es cerrada. Esto debe verificarse dado que la energía se conserva, es decir,   [Graphics:Images/Duffing.nb_gr_15.gif], de modo que para un x dado la velocidad está, salvo por el signo, también dada. Como la posición y la velocidad en un instante dado es todo lo que se necesita para determinar la posición y velocidad en un instante posterior, la conservación de la energía exige que el movimiento sea periódico. Puedes experimentar cambiando el tiempo hasta que obtengas justamente aquel en el que la trayectoria se cierrra por primera vez. Este tiempo no es nada más que el periodo de oscilación. Puedes comprobar que dejando correr la solución durante más tiempo la trayectoria en el plano de fases  no cambia porque el movimeiento simplemente se repite.

Fricción y amortiguamiento; órbitas periódicas

Debido a que la energía se conserva no puede esperarse encontrar movimiento caótico en sistemas con un sólo grado de libertad (en nuestro ejemplo, la posición es el único grado de libertad del oscilador unidimesional).  Por este motivo  vamos a complicar el movimiento de la partícula que se mueve sometida al potencial de Duffing   introduciendo  amortiguamiento (rozamiento)  y una fuerza motriz externa periódica.  (Si sólo introdujeramos amortiguamiento, el comportamiento no sería muy divertido: al final la particula acabaría siempre quieta sobre uno de los dos pozos de potencial).  La ecuación del movimiento es ahora

[Graphics:Images/Duffing.nb_gr_16.gif]

donde γ es el coeficiente de amortiguamiento y d es la magnitud (amplitud) de la fuerza motriz cuya  frecuencia de oscilación es   ω.   Veremos que a medida que la fuerza externa va aumentando en su magnitud, el movimiento de la partícula se irá haciendo más complejo hasta convertirse en  caótico.  Fijaremos γ = 0. y  ω = 1.4 en todo lo que sigue, y empezaremos con la amplitud de la fuerza externa d = 0.1 (veremos que esto hace que estemos fuera del regimen caótico).

[Graphics:Images/Duffing.nb_gr_17.gif]

Ahora consideraremos la situación en la que la partícula está inicialmente en reposo (velocidad v=0)  en x=0, es decir,  en la cima de la barrera de potencial.  Integraremos la ecuación del movimiento hasta  t = 200.  Para alcanzar este tiempo integrando con una precisión aceptable hay que asignar a la opción MaxSteps de NDSolve un valor más grande que el que Mathematica usa por defecto (Mathematica se quejará si no se hace así; puedes hacer la prueba).

[Graphics:Images/Duffing.nb_gr_18.gif]

Mostramos a continuación el aspecto que toma la trayectoria en el espacio de fases:

[Graphics:Images/Duffing.nb_gr_19.gif]

[Graphics:Images/Duffing.nb_gr_20.gif]

Su aspecto es complicado,  pero eso se debe a que se ha representado el "transitorio", es decir, la solución en el intervalo temporal en el cual la solución  se está acercando a su comportamiento final, el cual es mucho más simple. Es fácil ver que esto es efectivamene lo que sucede sin más que representar la solución a partir de un tiempo suficientemente grande. Por ejemplo, si representamos la trayectoria en el plano de fases con t entre 150 y 200  se obtiene una figura cuyo aspecto es mucho más sencillo que el de la figura anterior:

[Graphics:Images/Duffing.nb_gr_21.gif]

[Graphics:Images/Duffing.nb_gr_22.gif]

Lo que vemos ahora es simplemente  una órbita periódica, una  trayectoria cerrada.  Su periodo es justamente igual al periodo T de la fuerza motriz:  T=2 π / ω = 4.48799.... Esto se puede comprobar integrando para diferentes intervalos temporales. La figura siguiente muestra el resultado de la integración durante un intervalo de duración 4. Se ve , como era de esperar, que  la trayectoria no se ha completado, no se ha cerrado.

[Graphics:Images/Duffing.nb_gr_23.gif]

[Graphics:Images/Duffing.nb_gr_24.gif]

Transición al caos

A continuación vamos a indagar cómo cambian  las trayectorias en el plano de fases  cuando cambia  la amplitud de la fuerza motriz. Es conveniente definir una función que llamaremos "solutcion[d, tmax]" , la cual integra las ecuaciones del movimiento del oscilador de Duffing que estamos considerando desde t=0 hasta t=tmax para un valor dado d de la amplitud de la fuerza externa motriz.  (Usamos la asignación demorada ":=" porque queremos que sólo se evalúe  cuandomás tarde  se llame con ciertos valores particulares de  d y tmax.)

[Graphics:Images/Duffing.nb_gr_25.gif]
[Graphics:Images/Duffing.nb_gr_26.gif]

Empezamos eligiendo los valores  d = 0.32 y  tmax = 800.

[Graphics:Images/Duffing.nb_gr_27.gif]

Para dibujar la solución en el plano de fases es conveniente definir la función "graph[tmin, tmax]" que dibuja la solución entre  tmin y tmax:

[Graphics:Images/Duffing.nb_gr_28.gif]

Hasta t = 200 el movimiento de la particula tiene un aspecto muy complicado:

[Graphics:Images/Duffing.nb_gr_29.gif]

[Graphics:Images/Duffing.nb_gr_30.gif]

Fíjate que la partícula se ha movido dentro de los dos senos del pozo doble. Sin embargo, de nuevo, esta  complejidad se debe a que se dibuja el transitorio: la complejidad desaparece si trazamos la trayectoria sólo para tiempos más grandes.  Por ejemplo, la solución para t entre 750 y  800,  da lugar a

[Graphics:Images/Duffing.nb_gr_31.gif]

[Graphics:Images/Duffing.nb_gr_32.gif]

Esta gráfica muestra que la partícula, pasado el transitorio, se limita a oscilar dentro del pozo de la izquierda (alrededor del mínimo situado en  x = -1) y da dos vueltas completas (va, viene,va y viene) antes de volver al punto desde el cual partió.  Ahora el periodo se ha doblado para alcanzar el valor de 2T= 4π/ω  , tal como debe comprobarse mediante prueba y error.
Incrementemos ahora la magnitud de la fuerza externa hasta d= 0.339.

[Graphics:Images/Duffing.nb_gr_33.gif]
[Graphics:Images/Duffing.nb_gr_34.gif]

[Graphics:Images/Duffing.nb_gr_35.gif]


Ahora la órbita periódica, antes de volver al punto de partida y repetirse,  gira cuatro veces en torno al mínimo en x=1.
Parece natural pensar que el periodo se ha duplicado de nuevo hasta alcanzar el valor de 4T= 8π/ω . Esto es efectivamente lo que sucede  tal como uno puede comprobar mediante prueba y error. Lo que estamos viendo no es más que el comienzo de una secuencia de duplicación de periodo que es una forma característica mediante la cual los sistemas alcanzan el régimen caótico. Veamos qué pasa cuando la amplitud de la fuerza motriz se incrementa hasta d=0.35.

[Graphics:Images/Duffing.nb_gr_36.gif]
[Graphics:Images/Duffing.nb_gr_37.gif]

[Graphics:Images/Duffing.nb_gr_38.gif]


No se observa que la trayectoria tienda a ninguna estructura periódica (observa que la trayectoria es abierta).  ¿Podría suceder que estemos aún observando el regimen  transitorio? Dibujemos la trayectoria para tiempos mucho más grandes.

[Graphics:Images/Duffing.nb_gr_39.gif]
[Graphics:Images/Duffing.nb_gr_40.gif]

[Graphics:Images/Duffing.nb_gr_41.gif]

No se ve ningún indicio de que las trayectorias tiendan hacia una órbita periódica.  Podríamos aumentar más y más los tiempos de integración sin alcanzar nunca un régimen periódico.  Un estudio detallado de los valores en los cuales se producen las duplicaciones de los periodos mostraría que la secuencia [Graphics:Images/Duffing.nb_gr_42.gif]
tiende a la constante de Feigenbaum  δ=4.6692.... cuando n->∞, siendo [Graphics:Images/Duffing.nb_gr_43.gif] el valor de d para el que  se produce la n-sima duplicación de periodo. Esta constante es "universal", es decir, aparece en  los escenarios en los que se alcanza el régimen caóticos a través de bifurcaciones (duplicaciones) de periodo.

Secciones de Poincaré

Un modo útil de analizar el paso del comportamiento periódico al comportamiento caótico es mediante un diagrama conocido como sección de Poincaré. En vez de dibujar la trayectoria  completa [es decir, los puntos solución (x(t),v(t)) para todo instante] tal como hemos hecho hasta ahora,  en la sección de Poincaré se representan los puntos (x(t),v(t)) de la trayectoria para valores discretos del tiempo separados por el periodo T=2π/ω de la fuerza externa motriz. Es decir, se representan los puntos   (x(t),v(t)) para  t =T, 2T, 3T,...  Si consideramos los diagramas de fases que hemos obtenido anteriormente, es claro que para una órbita periódica de periodo 2π/ω  (por ejemplo, la que obtuvimos con d = 0.1) la sección de Poincaré consiste en un único punto; cuando el periodo se duplica (por ejemplo, para  d = 0.32) en la sección de Poincaré aparecen dos puntos;  cuando se vuelve a duplicar  el periodo (para d=0.338, por ejemplo) aparecen cuatro puntos, etc.

A continuación definimos una función llamada "poincare" la cual dibuja la sección de Poincaré correspondiente a los coeficientes  d, γ, y ω. Esta función asume que los primeros  "ndrop" periodos corresponden al transitorio inicial de modo que no son dibujados y sólo se representan los puntos (x(t),v(t)) correspondientes a los restantes "nplot" periodos, es decir, se representan los puntos (x(t),v(t)) con t=(ndrop+1)*T,  (ndrop+2)*T ,... nplot*T.

[Graphics:Images/Duffing.nb_gr_44.gif]

Si usamos los coeficientes γ = 0.1,  ω = 1.4, d = 0.339  correspondientes a la órbita de periodo cuatro obtenemos justamente cuatro puntos, como era de esperar.

[Graphics:Images/Duffing.nb_gr_45.gif]

[Graphics:Images/Duffing.nb_gr_46.gif]

Sin embargo, en el régimen caótico, la sección de Poincaré toma un aspecto mucho más imponente, sugerente, sorprendente...

[Graphics:Images/Duffing.nb_gr_47.gif]

[Graphics:Images/Duffing.nb_gr_48.gif]


Este diagrama tan extraño se llama atractor extraño. Es el conjunto de puntos hacia el cual tienden las órbitas caóticas, sean cuales sean las condiciones iniciales.  La estructura es complicada pero no amorfa, no aleatoria. Si uno amplia y observa con detalle cualquier zona del atractor, se encuentra una estructura similar a la anterior. Esta estructura es por tanto autosemejante: es un  fractal.


Ejercicios:
(1)  Dibuja las gráficas de las posiciones x frente al tiempo para los  distintos casos  discutidos anteriormente y compara con las trayectorias del diagrama de fase.  
(2*) Usa distintos valores de
[Graphics:Images/Duffing.nb_gr_49.gif] obtén los valores de la amplitud [Graphics:Images/Duffing.nb_gr_50.gif] en  los que se producen duplicaciones de periodo. (Ejercicio difícil)

*=Opcional


Converted by Mathematica      November 12, 2002