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

: Santos Bravo Yuste. email:  santos@unex.es

El fractal de las condiciones iniciales en el oscilador de Duffing

A continuación vamos a discutir la importancia de las condiciones iniciales en el comportamiento final del oscilador de Duffing.
El oscilador de Duffing que vamos a estudiar viene dado por los parámetros d=0.29,  γ = 0.1 y  ω = 1.4 .

Para las condiciones iniciales  x0=-1, v0=2, la solución final oscila con  periodo 2π/ω  alrededor del mínimo x=1.
Empezamos resolviendo numéricamente la ecuación diferencial con las condiciones de contorno anteriores.  La solución la llamamos sol4.

Representamos de forma paramétrica la solución en el plano (x,v) para tiempos muy grandes, en particular, para 700≤ t≤800 y comproamos que la solución final es oscilatoria alrededor del pozo con mínimo en x=1

Sin embargo, para  x0=-1 y  v0=-1, la solución final es también  oscilatoria  con  periodo 2π/ω pero ahora oscilando alrededor del otro mínimo situado en x=-1.
Empezamos resolviendo numéricamente la ecuación diferencial con las condiciones de contorno anteriores.  La solución la llamamos sol4.

Representamos de forma paramétrica la solución en el plano (x,v) para tiempos muy grandes, en particular, para 700≤ t≤800 y comprobamos que la solución final es oscilatoria alrededor del pozo con mínimo en x=-1

Hemos pues comprobado que el movimiento final de la partícula depende de sus condiciones iniciales (x0,v0).
¿Es posible predecir en qué pozo va a terminar oscilando la partícula? Veremos que esta cuestión no es nada sencilla.  Para intentar responder a esta pregunta empezamos definiendo una función llamada pozo que es igual a +1 si la partícula,  tras el tiempo tmax=numper*T=numper*2*Pi/ω,  oscila en el pozo de la derecha, y -1 si oscila en el pozo de la izquierda. Esta circunstancia  se determina simplemente comprobando si  x(tmax) es mayor o menor que cero.
Escogemos numper=200, de modo que tmax=Floor[2*Pi/ω*200]=897

La función pozo es

He aquí  un par de ejemplos que están de acuerdo con los resultados obtenidos más arriba. De paso mostramos el tiempo que tarda el ordenador en hacer estos dos cálculos.  

Generamos una matriz  en la que el elemento (m,n) representa el valor de la función pozo para
(x0=x0min+m*increx0, v0=v0min+m*increv0). A estas matrices las llamamos basinN, con N=5, 15, 30, 40, 50,...
A continuación calculamos esta matriz para una discretización muy gruesa: increx0=increv0=0.5;

Dibujamos la matriz anterior asignando el color rojo a la celdilla (m,n) si la trayectoria que partió del punto (x0=x0min+m*increx0, v0=v0min+m*increv0) temina oscilando alrededor del pozo centrado en x=1; si termina en el pozo x=-1, pintamos la casilla de verde.  Dibujamos la matriz tranpuesta para que la coordenada x esté sobre las abcisas y v sobre la ordenada.

Ahora usamos una discretización más fina:  increx0=increv0=0.25;.  Esto implica más tiempo de cálculo.  Este tiempo puede estimarse multiplicando el número de veces que hay que calcular la función pozo por el tiempo que tarda en evaluarse esta función (este último tiempo se cálculo más arriba). Probablemente podrás  tomarte un café antes de que salga el resultado.

Ahora usamos una discretización aún  más fina: increx0=increv0=0.1.  Esto implica aún más tiempo de cálculo.

Si no tienes paciencia para esperar el resultado final, la gráfica que resultaría es la siguiente

ListDensityPlot[basin30//Transpose, Mesh -> True,
    ColorFunctionScaling->False,
    ColorFunction->(Hue[(#+2)/3]&)];

ListDensityPlot[basin30//Transpose, Mesh -> False,
    ColorFunctionScaling->False,
    ColorFunction->(Hue[(#+2)/3]&)];

La discretización  increx0=increv0=0.01 es extraordinariamente fina y el cómputo se demora durante decenas de horas.   No es por tanto sensato ejecutar  

y esperar el resultado. .Lo mejor es ver directamente el resultado que  obtuve hace unos días (el cálculo de essta matriz tardó unas 45 horas) y que puedes descubrir haciendo doble click en la celda siguiente

ListDensityPlot[basin50//Transpose, Mesh -> False,
    ColorFunctionScaling->False,
    ColorFunction->(Hue[(#+2)/3]&)];

Esta figura es fractal y, por tanto, autosemejante (en sentido estadístico). Esto puede verse en la siguiente figura que no es más que la ampliación de la zona del espacio de fases definida por -2≤x0≤-1, -2≤v0≤-1 (cada cajita tiene un ancho  increx0=increv0=0.0025).   La matriz de puntos  se generó mediante las instrucciones

x0min=-2;v0min=-2; increx0=increv0=0.0025;
x0max=v0max=-1;
basin60=Table[pozo[xx,vv],{xx,x0min,x0max,increx0},{vv,v0min,v0max,increv0}];

Este cálculo requiere mucho tiempo de computación (varias horas, como poco) por lo que es recomendable que veas directamente el resultado haciendo doble click en la celda siguiente

ListDensityPlot[basin60//Transpose, Mesh -> False,
    ColorFunctionScaling->False,
    ColorFunction->(Hue[(#+2)/3]&)];

Podemos seguir con este proceso de ampliaciones sucesivas y siempre encontramos una figura muy parecida a la primera.  Ahora vamos a ampliar la región del plano de fases  -1.75≤x0≤-1.5, -1.75≤v0≤-1.5 (cada cajita tiene un ancho  increx0=increv0=0.000625).   La matriz de puntos  se generó mediante las instrucciones

x0min=-1.75; v0min=-1.75; increx0=increv0=0.000625;
x0max=v0max=-1.5;
basin70=Table[pozo[xx,vv],{xx,x0min,x0max,increx0},{vv,v0min,v0max,increv0}];

Por supuesto, este cálculo también requiere mucho tiempo de computación (en principio debe tardar tanto como el empleado para hallar basin60; ¿por qué?) por lo que es recomendable que veas directamente el resultado haciendo doble click en la celda siguiente

ListDensityPlot[basin70//Transpose, Mesh -> False,
    ColorFunctionScaling->False,
    ColorFunction->(Hue[(#+2)/3]&)];