Contents
Práctica sobre el MODELO DE LORENZ-HAKEN (2)
Programa traducido de Mathematica a Matlab 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]; % Umbral de la bifurcación de Hopf rHB=-s*(3+b+s)/(1+b-s);
Ecuaciones y resolución del modelo de Lorenz-Haken
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;%variables globales Tfin=1000;%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,[0.001,0.001,1]); tmin=900;%mínimo del eje de abcisas del plot tmax=1000;%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='Hard-mode excitation 1D'; figure('Name',titul),plot(t(nmin:nmax), solnum(nmin:nmax,1).^2); axis tight title(titul) xlabel('t') ylabel('F^2') titul='Hard-mode excitation 2D'; figure('Name',titul),plot(solnum(nmin:nmax,1), -solnum(nmin:nmax,3)); axis tight title(titul) xlabel('F') ylabel('-d') titul='Hard-mode excitation 3D'; figure('Name',titul),plot3(solnum(nmin:nmax,1), solnum(nmin:nmax,2),-solnum(nmin:nmax,3)); axis tight title(titul) xlabel('F') ylabel('P') zlabel('-d')
En esta ocasión, se fija un tmin>0 para saltarnos el transitorio hasta que se cae en el atractor caótico. De hecho, es fácil comprobar que si representamos de 0 a 100 en la escala de tiempos, se ven oscilaciones que van aumentando progresivamente (figura 1D) hasta llegar al régimen puramente caótico. Por tanto, el transitorio sería el intervalo entre 0 y 60 aproximadamente hasta ver el carácter autopulsado. Se observa un carácter similar autopulsado al del plot de Mathematica, pero con valores distintos dado que estamos en la solución caótica y hay dependencia grande con el método numérico utilizado. De hecho, si se cambia el paso h en Runge Kutta, el aspecto varía también considerablemente (análogo a la sensibilidad de las condiciones iniciales). En resumen, se observa claramente la estabilidad (la solución NO es estacionaria con un valor fijo, pero sí es estable u oscilante acotadamente) de la solución en torno al atractor caótico (atractor de Lorenz). Aparte, cabe apuntar que se ha elegido un paso de Runge-Kutta de 0.01 para que el cálculo no se haga muy largo (al menos en mi ordenador, que es poco potente). Se puede aumentar la precisión bajando h, pero se necesitaría más tiempo de cálculo. Por ejemplo, si quieres llegar a t=2000, conviene aumentar h hasta 0.04, sin perder excesiva resolución o precisión en la curva.
Comentarios físicos:
Cabe señalar que el atractor no "llena" todo el cubo del espacio de soluciones, sino que la región es aproximadamente plana, topología característica del atractor de Lorenz. De manera similar a los fractales, se podría definir una dimensión fraccionaria que lo describa. Aparte, las soluciones forman una especie de disco en cuyo centro debe estar la solución estacionaria, pero que está topológicamente aislada, por lo que no es alcanzable si el sistema ha caído en el atractor de Lorenz.
Transformada de Fourier
tmin=500; tmax=Tfin; nmin=find(t>=tmin,1,'first'); nmax=find(t<=tmax,1,'last'); FT=fft(solnum(nmin:nmax,1).^2); L_FT=length(FT) close all titul='Transformada de Fourier'; figure('Name',titul),semilogy(abs(FT(1:1000)).^2); axis tight title(titul)
L_FT =
50001
Se observa que el plot tiene los mismos picos que en Mathematica, como debería ser. Algunos valores se van ligeramente por la diferencia del método numérico, pero son valores sueltos y poco frecuentes respecto al comportamiento general. Cabe señalar que se escogen sólo las primeras frecuencias (de 1 a 1000), que es donde está lo interesante en la transformada. El hecho de que el espectro esté "lleno" de muchos picos es indicativo del comportamiento caótico del sistema.
Mapas de máximos
Verificar que has ejecutado previamente la orden solnum.
lista=solnum(nmin:nmax,1).^2; L=length(lista) m=[]; for k=2:L-1 if lista(k-1)<lista(k) && lista(k+1)<lista(k) m=[m,lista(k)]; end end L_m=length(m) close all titul='Mapa de máximos'; figure('Name',titul),plot(m(1:end-1),m(2:end),'.') xlabel('F^2max_i'), ylabel('F^2max_{i+1}') title(titul)
L =
50001
L_m =
318
Se observa que el plot es muy similar al obtenido en Mathematica. La posición de los picos no es exactamente la misma debido a las diferencias entre los distintos metódos numéricos (NDSolve frente a Runge-Kutta).
Comentarios físicos:
La forma de este mapa de máximos, con un pico pronunciado, es la signatura de que estamos estudiando el modelo de Lorenz. En la gráfica de F^2(t) se observa que los consecutivos máximos de oscilación parecen tener una cierta relación entre ellos (véase el comportamiento autopulsado). El diagrama pone de manifiesto esta relación entre máximos consecutivos, con una curva que es característica del atractor de Lorenz.
Doblamiento de periodo
close all s=3; b=1; r=140.0; Tfin=100;%tiempo hasta el que calculas h=0.01; t1=0:h:Tfin; %Define el paso h que utiliza Runge-Kutta [t1,solnum1] = rungekut('laser',t1,[1,1,1]); tmin=80;%mínimo del eje de abcisas del plot tmax=Tfin;%máximo del eje de abcisas del plot nmin=find(t1>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear nmax=find(t1<=tmax,1,'last'); titul='Doblamiento de periodo 1D'; figure('Name',titul),plot(t1(nmin:nmax), solnum1(nmin:nmax,1).^2); axis tight title(titul) xlabel('t') ylabel('F^2') titul='Doblamiento de periodo 2D'; figure('Name',titul),plot(solnum1(nmin:nmax,1), solnum1(nmin:nmax,2)); axis tight title(titul) xlabel('F') ylabel('P')
Se aprecia que la forma es muy similar al plot de Mathematica, al estar en una solución estable con periodo único (solución oscilante con misma amplitud y frecuencia). De nuevo, elegimos un valor de plot con tmin>0 para recortar el transitorio del cálculo.
tmin=80; tmax=Tfin; nmin=find(t1>=tmin,1,'first'); nmax=find(t1<=tmax,1,'last'); lista=solnum1(nmin:nmax,1).^2; L=length(lista); FT1=fft(solnum1(nmin:nmax,1).^2); L_FT1=length(FT1) close all titul='Transformada de Fourier'; figure('Name',titul),semilogy(abs(FT1(1:200)).^2); axis tight title(titul) m=[]; for k=60:L-1 if lista(k-1)<lista(k) && lista(k+1)<lista(k) m=[m,lista(k)]; end end L_m=length(m) titul='Mapa de máximos'; figure('Name',titul),plot(m(1:end-1),m(2:end),'.') title(titul); xlabel('F^2max_i'), ylabel('F^2max_{i+1}'); titul='Mapa de máximos (reescalado)'; figure('Name',titul),plot(m(1:end-1),m(2:end),'.') title(titul); xlabel('F^2max_i'), ylabel('F^2max_{i+1}') axis([0,500,0,500]);
L_FT1 =
2001
L_m =
34
Se comprueba que hay coincidencia entre estos plots y el de Mathematica. A continuación, se hace un nuevo cálculo de los máximos con ajuste parabólico, para mejorar la dispersión que se obtiene en el mapa de máximos no escalado. Aparte, se verifica que según se va duplicando el periodo, van apareciendo en la transformada de Fourier pico más pequeños, lo que indica que la función original es periódica (pero no armónica). A partir de multiplicar por 8 ó 16 el periodo inicial, es difícil distinguir cada vuelta en el plot 2D. La cola que aparece al final de la transformada de Fourier es ruido debido ya al método numérico, por lo que es conveniente cortarlo.
Comentarios físicos:
Bajando de r=140 hacia r=100, etc. iremos duplicando el periodo hasta llegar al comportamiento caótico del atractor (r_1=140, r_2=100, r_4=96.8, r_8=96.1, r_16=96.04). Este proceso se denomina ruta al caos por duplicación de periodo o de Feigenbaum. En la transformada de Fourier van apareciendo armónicos a distancias equiespaciadas (señal de movimiento uniforme), como múltiplos del fundamental. Esto es una prueba de que el periodo se duplica, cuadruplica, etc. También se puede razonar a partir de la intensidad del pico. Cuando bajas de r=91, pasas a la zona caótica, el periodo es infinito y se juntan todos los picos dando una transformada ruidosa. Si fuese igualmente probable (suceso azaroso), saldría una zona casi plana. Pero ése no es nuestro caso, sino que tenemos un comportamiento caótico, que muestra un cierto pico fundamental, correspondiente a lo que tarda en cambiar entre uno y otro atractor (las dos alas), y además modulado por las infinitas frecuencias. En r=78 se observa un nuevo atractor, pero con distinta topología. Hay 1 general, y otro "interior" en cada solución.
Cálculo con interpolación parabólica
lista=solnum1(nmin:nmax,1).^2; L=length(lista); m=[]; for k=60:L-1 if lista(k-1)<lista(k) && lista(k+1)<lista(k) parab=polyfit([-1;0;1], lista(k-1:k+1),2); x=-parab(2)/2/parab(1); y=polyval(parab,x); m=[m,y]; end end close all titul='Mapa de máximos con interpolación (reescalado)'; figure('Name',titul),plot(m(1:end-1),m(2:end),'.') title(titul); xlabel('F^2max_i'), ylabel('F^2max_{i+1}'); axis([395,397.2,396,397.2]);
Se ve, escogiendo la misma escala que en el mapa de máximos no escalado anterior, que la dispersión de los puntos es mucho menor, al haber utilizado un método más sofisticado de cálculo de los máximos (interpolación parabólica). Por tanto, el hecho de que haya dispersión de los puntos no se debe a que haya distintas frecuencias (desdoblamiento de periodo) sino al propio error numérico acumulado en el procedimiento de cálculo de máximos, que puede ser más o menos preciso.
De cara a una mejora del programa, sería apropiado definir ficheros específicos con las funciones de búsqueda de máximos, realización de transformada de Fourier, etc, para llamarlos desde este programa general, y así evitar tener que estar copiar y pegando líneas de código. Además, si quisiésemos que el mapa de máximos siempre se hiciese con interpolación parabólica, bastaría con cambiar dicho fichero aparte y no cada orden si se pega específicamente en cada apartado del programa.
close all