Contents

Práctica sobre el MODELO DE LORENZ-HAKEN (1)

Programa traducido de Mathematica a Matlab y comentado por Fernando Hueso González. Óptica cuántica (Seminarios de teoría semiclásica y simulaciones) - 4º de Grado de Física, 31 de mayo de 2011.

global s b r
syms d s b r
format compact

% Solución estacionaria (analítica)
d = 1;
P = sqrt(r-1);
F = sqrt(r-1);
solestac = [d,P,F];

% Estabilidad lineal (hasta bifurcación de Hopf rHB)
rHB=-s*(3+b+s)/(1+b-s);

COMPROBACIÓN DE:

a) Estabilidad lineal (condiciones iniciales próximas a la solución estacionaria)

Evaluación para valores concretos

RHB=subs(rHB,{s,b},[3,1])
Solestac=subs(solestac,r,22)

RHB =
    21
Solestac =
    1.0000    4.5826    4.5826

Comentarios físicos:

La solución estacionaria no es la única solución "estable" en la que el láser puede funcionar. Para r grande, mayor que la bifurcación de Hopf, tendremos soluciones caóticas en general, con posibles islas de estabilidad y soluciones oscilantes de frecuencia fija (Ver doblamiento de periodo el segundo archivo). También podremos abandonar la solución estacionaria si hacemos un encendido brusco (Hard-mode excitation), con condiciones iniciales lejanas a las de la solución estacionaria, de forma que ya no la "encuentra (se puede hacer la analogía mecánica de que intentas encestar una pelota en un pequeño pozo de potencial justo en la cúspide de una montaña). Por contra, si ponemos condiciones iniciales cercanas a la misma, ese pequeño ruido (cambio adiabático) no es suficiente para sacarlo de ella y cae a la línea de estabilidad.

Ecuaciones diferenciales del láser definidas en el archivo laser.m

%F'=s*(P-F)
%P'=-P+d*F
%d'=b*(r-d-F*P)
s=3; b=1; r=20; %no lo defino como vector "par" sino como variable global para ahorrar notación
Tfin=100;%tiempo hasta el que calculas

h=0.01; t=0:h:Tfin; %Define el paso h que utiliza Runge-Kutta
[t,solnum] = rungekut('laser',t,[4.58258,4.58258,1]);

solnum es una matriz con tres columnas ordenadas (F,P,d). Las filas marcan distintos valores del tiempo. Se fijan los valores iniciales [F0,P0,d0]=[4.58258,4.58258,1] por ejemplo. Se utiliza el método de Runge-Kutta porque el de ODE45 no se comporta adecuadamente, es decir, no converge a la solución estable sino que es infinitamente oscilante para r=20, s=3, b=1 debido a error numérico.

tmin=0;%mínimo del eje de abcisas del plot
tmax=Tfin;%máximo del eje de abcisas del plot

nmin=find(t>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear
nmax=find(t<=tmax,1,'last');
close all
titul='Estabilidad lineal';
figure('Name',titul),plot(t(nmin:nmax), solnum(nmin:nmax,1).^2);
%el 1 representa la primera columna, es decir, la F
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

Como se observa, el transitorio hasta llegar a la solución final estacionaria es muy largo. Si quisiésemos centrarnos en la estacionaria, deberíamos poner Tfin grande (500 aprox.) y establecer tmin=400. De esta manera podemos cortar el transitorio de los cálculos. Se definen en todos los plots estas variables para poder prescindir de los transitorios en caso de que no nos interesasen.

b) Hard-mode excitation -> biestabilidad generalizada

s=3; b=1; r=20;
Tfin=100;
h=0.01; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);%orden F0,P0,d0
close all
titul='Hard-mode excitation'
figure('Name',titul),plot(t, solnum(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

titul =
Hard-mode excitation

Comentarios físicos:

El Hard-mode excitation se alcanza dando unos valores iniciales alejados de la solución estacionaria. De esta manera, el sistema no es capaz de "encontrar" dicha solución estacionaria y cae al atractor caótico y permanece en ella de manera estable (ya no puede volver a la estacionaria). A diferencia de a), donde s,b y r son compartidos, el valor de F0 y P0 se cambia sustancialmente para sacar al sistema de la solución estacionaria. Para F0=P0=4.5826 (caso a), que son los valores de la solución estacionaria, el sistema es capaz de alcanzar la solución estacionaria pese al alto bombeo, ya cercano a rHB (bifurcación de Hopf), si se deja pasar suficientemente tiempo (t=300 ó 400). Por contra, en el caso b), cambiamos notablemente dichos valores iniciales, y provocamos el hard mode excitation, es decir, el salto a la solución caótica pese a estar por debajo de rHB, dado que la perturbación inicial es grande y se va amplificando cada oscilación en lugar de irse atenuando hacia la estable. Esta es la principal diferencia entre el cambio adiabático y el hard mode excitation. En otras palabras, si damos un "golpe" muy fuerte al láser (valores iniciales lejanos a la estacionaria) o subimos de golpe el bombeo manteniendo valores iniciales de otra solución, podemos salirnos de la solución estacionaria y caer a la caótica. De ahí la importancia del cambio adiabático (subir poco a poco, con condiciones iniciales siempre cercanas a la estacionaria) para mantenerse en la estacionaria hasta rHB.

Cabe resaltar la diferencia entre estacionario y estable. La solución estacionaria indica que el valor de la variable es fijo en el tiempo, no cambia (derivada nula). La solución estacionaria es única y está calculada analíticamente (ver solestac). El resto de soluciones dinámicas del sistema son no estacionarias. Pero si el láser se mantiene en un subconjunto de soluciones (no estacionarias), oscilando por ejemplo alrededor de un valor medio, entonces decimos que la solución es estable. Es el caso de la forma autopulsada estable de la curva correspondiente al hard-mode excitation.

Representación de atractores

s=3; b=1; r=30;
Tfin=100;
h=0.01; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);%orden F0,P0,d0

close all
titul='Representación de atractores 2D';
figure('Name',titul),plot(solnum(:,1), -solnum(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('-d')

titul='Representación de atractores 3D';
figure('Name',titul),plot3(solnum(:,1), solnum(:,2),-solnum(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')

Comentarios físicos:

Aparecen dos "alas de mariposa" o cuencas de atracción, es decir, las soluciones dinámicas del sistema dan vueltas alrededor de una zona del espacio de soluciones. Esto es característico del atractor de Lorenz. El hecho de que aparezcan dos puntos de atracción se debe a que en los casos 2D y 3D se representa F y no F^2 como en el caso 1D. Si en el 1D representásemos F, se vería que saltaría a valores negativos y positivos aleatoriamente. Cabe señalar que la notación 1D, 2D y 3D se refiere a: 1D-> F(t); 2D->Representación paramétrica F(t), -d(t); 3D -> Paramétrica 3D F(t), P(t), -d(t).

Comprobación Sensibilidad a las Condiciones Iniciales

s=3; b=1; r=23;
Tfin=50;
h=0.01; t2=0:h:Tfin;
[t2,solnum2] = rungekut('laser',t2,[1,1,1]);%orden F0,P0,d0
% Evaluar los valores de F, P y D para el último t calculado:
F0=solnum2(end,1)
P0=solnum2(end,2)
d0=solnum2(end,3)

Tfin=20;
h=0.01; t3=0:h:Tfin; t4=0:h:Tfin;
[t3,solnum3] = rungekut('laser',t3,[F0,P0,d0]);
[t4,solnum4] = rungekut('laser',t4,[F0+0.01,P0+0.01,d0]);

close all
titul='Sensibilidad a las condiciones iniciales (s3)';
figure('Name',titul),plot(t3, solnum3(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

titul='Sensibilidad a las condiciones iniciales (s4)';
figure('Name',titul),plot(t4, solnum4(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

titul='Sensibilidad a las condiciones iniciales (3vs4)';
figure('Name',titul),plot(t3, solnum3(:,1).^2, t4, solnum4(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')
legend('solnum3','solnum4')

tmin=0;%mínimo del eje de abcisas del plot
tmax=2;%máximo del eje de abcisas del plot

nmin=find(t3>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear
nmax=find(t3<=tmax,1,'last');

titul='Sensibilidad a las condiciones iniciales (s3^2-s4^2)';
figure('Name',titul),plot(t3(nmin:nmax), solnum3(nmin:nmax,1).^2 - solnum4(nmin:nmax,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F_3^2-F_4^2')

F0 =
    5.7405
P0 =
    1.1309
d0 =
   -3.8313

Cabe señalar que estos últimos plots (y los de aquí en adelante) son ligeramente distinto al de Mathematica. La razón es que utilizamos un valor de F0, P0 y d0 calculado mediante Runge-Kutta que es notablemente distinto al obtenido mediante NDSolve. Esta diferencia se debe a que estamos en una zona caótica (solución no es estacionaria), con lo que cualquier mínima diferencia entre uno y otro método numérico genera variaciones grandes en la curva y por tanto en el valor de F, P y d en Tfin. Si se representa solnum2 en Matlab y Mathematica, se aprecia que su aspecto es similar pero desplazado en el tiempo, y a partir de cierto valor de t, el orden de picos es difícil seguirlo. De hecho, esta observación constituye una prueba más de la sensibilidad a las condiciones iniciales.

close all
titul='Sensibilidad a las condiciones iniciales 3D (s3)';
figure('Name',titul),plot3(solnum3(:,1), solnum3(:,2),-solnum3(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')

titul='Sensibilidad a las condiciones iniciales 3D (s4)';
figure('Name',titul),plot3(solnum4(:,1), solnum4(:,2),-solnum4(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')

titul='Sensibilidad a las condiciones iniciales 3D (3vs4)';
figure('Name',titul),plot3(solnum3(:,1), solnum3(:,2),-solnum3(:,3), solnum4(:,1), solnum4(:,2),-solnum4(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')
legend('solnum3','solnum4')

close all

Comentarios físicos:

Se aprecia que el sistema tiene un comportamiento caótico debido a la gran sensibilidad en las condiciones iniciales. Pese a poner una perturbación muy pequeña entre dos soluciones, ambas acaban divergiendo si se deja pasar el tiempo necesario (más tiempo cuanto menor sea la perturbación). Esto es característico de los sistemas caóticos clásicos (caos determinista), donde la incertidumbre en las condiciones iniciales te impiden predecir el devenir del sistema. Cualquier pequeña variación o ruido va a provocar que el cabo de un tiempo largo dé una solución radicalmente diferente a la obtenida sin dicha variación.

Published with MATLAB® 7.11