Estudio de las soluciones numéricas del

oscilador de doble pozo

David Casas García-Minguillán

Edi...cio C-2, Campus de Rabanales, Universidad de Córdoba, 14071 España

9 de febrero de 2005

Resumen

En este artículo he calculado las soluciones numéricas, con distintas

condiciones iniciales, de la ecuación de Du� ng que rige el movimiento

del oscilador de doble pozo. He utilizado el conocido método de Runge-

Kutta de orden cuatro, programado en Fortran 90. Finalmente, comento

de forma cualitativa el paso de oscilaciones regulares a caóticas en este

oscilador no lineal.

1.

Introducción

El oscilador de doble pozo es muy referenciado, como ejemplo de oscilador

no lineal, en la literatura cientí...ca dedicada al estudio del caos, atractores ex-

traños, soluciones aperiodicas, etc. Para determinadas soluciones se asemeja a

un sistema cáotico, sensible a pequeñas perturbaciones externas y por tanto

se comporta de forma impredecible, a pesar de estar de...nido por ecuaciones

deterministas, como veremos a continuación[1].

En este artículo he hecho una simulación por ordenador del comportamiento

del oscilador, sin embargo, aparte de los estudios computacionales, hay estu-

dios experimentales hechos sobre este oscilador los cuales con...rman el compor-

tamiento caótico de este sistema[2].

A continuación expondré la descripción matemática del problema, es decir,

el sistema de ecuaciones diferenciales con sus condiciones iniciales. Despues co-

mentaré el método de Runge-Kutta, el cual lo he programado usando Fortran

90. Explicaré los programas y su funcionamiento. Calcularé tres casos particu-

lares, y hablaré sobre los datos y grá...cas obtenidos, los cuales dan un indicio

del paso del orden al caos y rehaciendo los calculos para un rango de valores

de una variable, mostraré donde cambia bruscamente de un comportamiento a

otro. Este salto se visualiza mucho mejor en los videos que se encuentran en

el CD-ROM que acompaña a este artículo. Son animaciones hechas a partir de

muchas grá...cas individuales1 , generadas con aproximadamente medio millón

1 Realizadas a mano por el propio autor, lo que coloquialmente se denomina "trabajo de chinos".

de puntos distintos, en total. Finalmente, en el apéndice, he dejado escritos los

códigos de los programas. En el CD-ROM, aparte de los videos, está grabado

este mismo artículo, los dos programas en código fuente, los resultados de cada

apartado en sus carpetas correspondientes y dos archivos que corresponden a

las referencias 1 y 2 de este artículo.

2.

Descripción matemática del problema

El oscilador de doble pozo puede describirse matemáticamente mediante la
llamada ecuación de Duffing:

2. Dado que la frecuencia es el inverso del tiempo, la t que aparece dentro del coseno es adimensional: t(u:a:) = 1(T 1 ) t(T ) donde u.a. son unidades arbitrarias y T es la dimensión tiempo.

3.

3.1.

Métodos de Runge-Kutta

Introducción

El método de Runge-Kutta es uno de los existentes dentro de la clase de

métodos de un solo paso, los cuales se pueden describir genericamente mediante

las ecuaciones:

donde es una función que caracteriza el método.

Para estimar la rapidez de convergencia de un método u otro, se dice que

(4) es de orden p si para todo xi axib, y para todo h su...cientemente

pequeño, existen constantes C y p tales que si

La constante C depende, en general, de la solución y(x), de sus derivadas y de

la longitud del intervalo sobre el que se busca la solución, pero es independiente

de h[4].

Dentro de estos métodos están los de la serie de Taylor de la solución y(x).

Si y p+1) (x) es continua en [a; b] entonces el desarrollo viene dado por:

3.2.

Runge-Kutta

La idea general de este método es sustituir las derivadas de f presente en

(8), por combinaciones lineales de la propia función f (en diferentes puntos)

para obtener el orden requerido.

3

Por ejemplo, un algoritmo de Runge-Kutta de orden 2 es el llamado método

de Euler mejorado o de Heun:

La mayor limitación de las fórmulas de Runge-Kutta es la cantidad de tra-

bajo que requieren; este trabajo se mide en términos del número de veces que la

función f se evalua. Para fórmulas de alto orden, el trabajo crece excesivamente;

para los métodos de orden p con p = 1; 2; 3 y 4 el número de evaluaciones de f

es justamente p, pero para p = 5, es necesario hacer 6 evaluaciones, para p = 6

hay que hacer 7, para p = 7 ya son 9 las que se necesitan, para p = 8 son 11,

etc.

Como consecuencia, los métodos de orden menor a cinco con tamaño de paso

pequeño son preferibles a los métodos de orden superior usando un tamaño de

paso más grande.

Los métodos más utilizados son los de orden p = 4 entre ellos la fórmula

clásica de Runge-Kutta de orden 4, esta se escribe en la forma

�

Este método es muy utilizado por su rapidez de convergencia y por su sen-

cillez de fórmula. De hecho, este método es capaz de conseguir precisiones altas

sin tener que tomar h tan pequeño como para hacer excesiva la tarea computa-

cional o que los cálculos de redondeo planteen serias di...cultades[5].

Cuando tenemos ecuaciones de orden superior a 1, se pueden reescribir en

forma de sistemas de ecuaciones de primer orden, si tenemos derivadas de y

hasta el orden m, el sistema sería:

Escrito en notación vectorial, el algoritmo de Runge-Kutta de…nido en (10),
con Y = (y1; :::; ym), F = (f1; :::; fm),  = (1; :::; m), se escribiría[3]:

Este será el método, que programaré para resolver numéricamente la ecuación del oscilador de doble pozo.

4.

4.1.

Resultados

Programa RK4Oscilador

Este sencillo programa genera dos ficheros, posición.txt y velocidad.txt con

los resultados de los ídem. Los datos que hay que introducir son el tiempo de ini-

cio, el de finalización, el numero de divisiones del intervalo (con estos tres datos

hemos establecido la longitud del paso, h), las condiciones iniciales de posición

y velocidad, el coe...ciente de amortiguamiento y la amplitud de la fuerza. El

comando TIMEF() sirve para evaluar el tiempo que tarda en ejecutarse la parte

del programa situada entre las dos llamadas de esta orden. Además, al acabar,

da en pantalla el ultimo valor del intervalo de la variable independiente 3 , la

posición, la velocidad, el valor del paso, h y el tiempo que ha tardado la parte

ejecutora del programa. A continuación expondré los tres casos que he estudia-

do con este programa, explicando que valor he tomado para todas las variables,

comentando las gra...cas del espacio de fases, de la posición y de la velocidad.

4.1.1.

Primer caso

Los valores de las variables son los de la siguiente tabla:

He elegido este primer caso para comprobar que no había cometido ningún
error de bulto4 al escribir el algoritmo. Para ello elegí un software comercial que
resolviese una EDO de forma numérica y contrasté los resultados con los de mi
programa. Usé el Scientific WorkPlace 5.0, el mismo que me ha servido para
redactar este artículo, y esta es la tabla que compara uno y otro:

3 Es la variable temporal del sistema físico. No hay que confundirla con el tiempo que tarda en ejecutarse el programa, tal y como se ha comentado antes.
4 En una prueba preliminar con otra función puse en el algoritmo h=2 donde tendría que haber puesto h. El resultado no fue catastro…co, pero convergía a la solución verdadera mucho más lentamente.

desconozco que método numérico usa el SWP, hay que comentar, como se
verá a continuación, que con estas condiciones el sistema oscila de forma estable
alrededor del punto de equilibrio 1, de forma que las pequeñas diferencias de
método y precisión no afectan significativamente al resultado. Posiblemente, en
una región caótica no deban coincidir, aunque no lo he comprobado.
Las grá…cas y la pantalla de ejecución son las siguientes:

4.2.

Programa RK4ODP

Quería visualizar como ocurría el cambio del orden al caos con los datos de los

casos dos y tres, variando únicamente la fuerza. Por eso modi...qué el programa,

aunque el algoritmo de Runge-Kutta es el mismo, lo único es que se evalua

8

más veces. Introduzo un intervalo de variación de la fuerza, y el programa hace

los mismos cálculos que el anterior para cada valor de fuerza. Luego graba los

resultados en matrices de tres dimensiones (el anterior los grababa en matrices

de dos) y escribe los resultados en los siguientes ...cheros: fuerza.txt, posición.txt,

velocidad.txt y fases.txt. Tanto este programa como el anterior graba todos los

números en un único vector columna, ya que así se pueden grabar grandes

cantidades de números sin problema. Luego el programa con el que construyo

las grá...cas, el Microcal Origin 6.0, al importar el ...chero ASCII tiene la opción

de interpretar los caracteres no numericos como inicio de una nueva columna,

de forma que los ordena en formato de hoja de datos.

Ejecuté el programa con un valor de inicio de la fuerza de 0;20 y un valor

...nal de 0;45 en 50 intervalos, con lo cual variaba en pasos de 0;005. En las

siguientes cuatro grá...cas vemos como salta del orden al caos, en las posiciones,

entre F = 0;26 y F = 0;265:

5.

Conclusiones

He obtenido las soluciones numéricas de la ecuación de Du� ng aplicadad al

oscilador de doble pozo usando el método de Runge-Kutta de orden 4. Apli-

candolo a distintos casos, con diferentes valores de las variables y condiciones

10

iniciales, he mostrado como las soluciones cambian de la periodicidad a la ape-

riodicidad e impredictibilidad. Cuando el sistema se encuentra en el régimen

caótico, es muy sensible a las condiciones iniciales, de forma que hacía falta

dar pasos muy ...nos para precisar donde ocurria. Eso es lo que he hecho, sin

embargo hay que matizar que lo he hecho para una sección muy concreta del

espacio de fases total, donde están representadas todas las variables y condi-

ciones iniciales del problema, como son t, x, v, x (t0 ), v (t0 ), b, F , h, t0 , t1 e

incluso ! que aunque la hemos obviado, también interviene en la ecuación orig-

inal. Intuitivamente se puede adivinar el comportamiento del sistema, ya que

si sabemos que parte de un punto estable, con poca velocidad inicial, con un

coe...ciente amortiguador distinto de cero y una fuerza no muy grande, oscilará

en torno al punto de equilibrio. Podemos por tanto, intentar adivinar que zonas

pueden ser de interés para estudiarlas, ya que la simulación total, variando todo

lo que puede variar, puede generar una cantidad astronómica de números, que

para ser tratados necesitarían un gran tiempo de computación incluso en un

superordenador.

Sin necesidad de ser tan ambicioso este artículo muestra de forma didáctica

y aproximada el paso del orden al caos en un oscilador no lineal.

6.

6.1.

Apéndice

Código del programa RK4Oscilador

PROGRAM RK4Oscilador

USE PORTLIB

!Este es un programa para probar el metodo

!de Runge-Kutta de orden 4 en el caso del

!Oscilador de Doble Pozo

IMPLICIT NONE

REAL::b,h,f,t,x,v,tiempo,t0,t1

INTEGER::i,j,n

REAL,ALLOCATABLE,DIMENSION(:,:)::posres,velres

REAL,DIMENSION(4)::k,l

PRINT *, "Introduce los extremos del intervalo"

PRINT *, "t0 y t1:"

READ *,t0

READ *,t1

PRINT *, "Introduce las condiciones iniciales x(0) y v(0):"

READ *,x

READ *,v

PRINT *, "Introduce las divisiones del intervalo, N"

READ *,n

PRINT *, "Introduce el coe...ciente de amortiguamiento, b"

READ *,b

PRINT *, "Introduce la amplitud de la fuerza, F"

11

READ *,f

ALLOCATE(posres(n+1,2),velres(n+1,2))

h=(t1-t0)/n

t=t0

posres(1,1)=t

posres(1,2)=x

velres(1,1)=t

velres(1,2)=v

tiempo=TIMEF()

DO i=1,n

k(1)=v

l(1)=-b*v-x*x*x+x+f*cos(t)

DO j=1,2

k(j+1)=v+l(j)*h/2

l(j+1)=-b*(v+l(j)*h/2)-(x+k(j)*h/2)*(x+k(j)*h/2)*(x+k(j)*h/2)+(x+k(j)*h/2)+f*cos(t+h/2)

END DO

k(4)=v+l(3)*h

l(4)=-b*(v+l(3)*h)-(x+k(3)*h)*(x+k(3)*h)*(x+k(3)*h)+(x+k(3)*h)+f*cos(t+h)

x=x+(h/6)*(k(1)+2*k(2)+2*k(3)+k(4))

v=v+(h/6)*(l(1)+2*l(2)+2*l(3)+l(4))

t=t0+h*i

posres(i+1,1)=t

posres(i+1,2)=x

velres(i+1,1)=t

velres(i+1,2)=v

END DO

tiempo=TIMEF()

PRINT *,t,"t"

PRINT *,x,"x"

PRINT *,tiempo,"segundos"

PRINT *,h,"h"

OPEN(8,...le="posicion.txt")

DO j=1,2

DO i=1,n+1

WRITE (8,*) posres(i,j)

END DO

WRITE (8,*) 'a'

END DO

CLOSE(8)

OPEN(9,...le="velocidad.txt")

DO j=1,2

DO i=1,n+1

WRITE (9,*) velres(i,j)

END DO

WRITE (9,*) 'a'

END DO

12

CLOSE(9)

END PROGRAM RK4Oscilador

6.2.

Código del programa RK4ODP

PROGRAM RK4ODP

USE PORTLIB

!Este es un programa para probar el metodo

!de Runge-Kutta de orden 4 en el caso del

!Oscilador de Doble Pozo

IMPLICIT NONE

REAL::b,h,f,f0,f1,p,t,x,v,tiempo,t0,t1

INTEGER::i,j,g,n,m

REAL,ALLOCATABLE,DIMENSION(:,:,:)::posres,velres,fases

REAL,ALLOCATABLE,DIMENSION(:)::fuerza

REAL,DIMENSION(4)::k,l

PRINT *, "Introduce los extremos del intervalo"

PRINT *, "t0 y t1:"

READ *,t0

READ *,t1

PRINT *, "Introduce las condiciones iniciales x(0) y v(0):"

READ *,x

READ *,v

PRINT *, "Introduce las divisiones del intervalo temporal, n"

READ *,n

PRINT *, "Introduce el coe...ciente de amortiguamiento, b"

READ *,b

PRINT *, "Introduce el intervalo de la fuerza, F:"

PRINT *, .a mplitud inicial y ...nal, f0 y f1:"

READ *,f0

READ *,f1

PRINT *, "Introduce las divisiones del intervalo virial, m"

READ *,m

ALLOCATE(posres(n+1,2,m+1),velres(n+1,2,m+1),fases(n+1,2,m+1),fuerza(m+1))

h=(t1-t0)/n

t=t0

p=(f1-f0)/m

DO g=1,m+1

posres(1,1,g)=t

posres(1,2,g)=x

velres(1,1,g)=t

velres(1,2,g)=v

fases(1,1,g)=x

fases(1,2,g)=v

END DO

tiempo=TIMEF()

13

DO g=1,m+1

f=f0+p*(g-1)

DO i=1,n

k(1)=v

l(1)=-b*v-x*x*x+x+f*cos(t)

DO j=1,2

k(j+1)=v+l(j)*h/2

l(j+1)=-b*(v+l(j)*h/2)-(x+k(j)*h/2)*(x+k(j)*h/2)*(x+k(j)*h/2)+(x+k(j)*h/2)+f*cos(t+h/2)

END DO

k(4)=v+l(3)*h

l(4)=-b*(v+l(3)*h)-(x+k(3)*h)*(x+k(3)*h)*(x+k(3)*h)+(x+k(3)*h)+f*cos(t+h)

x=x+(h/6)*(k(1)+2*k(2)+2*k(3)+k(4))

v=v+(h/6)*(l(1)+2*l(2)+2*l(3)+l(4))

t=t0+h*i

posres(i+1,1,g)=t

posres(i+1,2,g)=x

velres(i+1,1,g)=t

velres(i+1,2,g)=v

fases(i+1,1,g)=x

fases(i+1,2,g)=v

END DO

fuerza(g)=f

END DO

tiempo=TIMEF()

PRINT *,t,"t"

PRINT *,x,"x"

PRINT *,f,"fuerza"

PRINT *,tiempo,"segundos"

OPEN(8,...le="posicion.txt")

DO g=1,m+1

DO j=1,2

DO i=1,n+1

WRITE (8,*) posres(i,j,g)

END DO

WRITE (8,*) 'a'

END DO

WRITE (8,*) 'pppp'

END DO

CLOSE(8)

OPEN(9,...le="velocidad.txt")

DO g=1,m+1

DO j=1,2

DO i=1,n+1

WRITE (9,*) velres(i,j,g)

END DO

WRITE (9,*) 'a'

14

WRITE (9,*) 'vvvv'

END DO

END DO

CLOSE(9)

OPEN(10,...le="fases.txt")

DO g=1,m+1

DO j=1,2

DO i=1,n+1

WRITE (10,*) fases(i,j,g)

END DO

WRITE (10,*) 'a'

WRITE (10,*) '��'

END DO

END DO

CLOSE(10)

OPEN(11,...le="fuerza.txt")

DO g=1,m+1

WRITE (11,*) fuerza(g)

END DO

CLOSE(11)

END PROGRAM RK4ODP

Agradecimientos

Agradezco a Mercedes Marín Beltrán su ayuda y consejos en la elaboración

de este trabajo.

Referencias

[1] V Conferencia Nacional de Ciencias de la Computación (noviem-

bre 98, Potosí), Sistemas no lineales. Conceptos, algoritmos y aplicaciones,

José Manuel Gutiérrez

[2] Chaos, Solitons & Fractals, Vol. 8, No 4, pp. 681-697, Routes to

Chaos in a Magneto-Elastic Beam, R.C.A. Thompson, T. Mullin

[3] Tema 2, 26-27, Apuntes de Ampliación de Métodos Numéricos

(2005), M. Marín

[4] Tema 2, 4-6, Apuntes de Ampliación de Métodos Numéricos (2005),

M. Marín

[5] Tema 2, 11-14, Apuntes de Ampliación de Métodos Numéricos

(2005), M. Marín

15