%% 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)
%% 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')
%% 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 plote
