Contents

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

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

s=3; b=1; r=20;
Tfin=1000; h=0.02; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);
tmin=900; tmax=1000; nmin=find(t>=tmin,1,'first'); 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')

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 =
       25001

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)
        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
L_m=length(m)
% implemento la interpolación parabólico en todos los cálculos, al haber
% comprobado que da mejor resultado.
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 =
       25001
L_m =
   319

Doblamiento de periodo

close all
s=3; b=1; r=99.0;
Tfin=100; h=0.01; t1=0:h:Tfin;
[t1,solnum1] = rungekut('laser',t1,[1,1,1]);

tmin=80; tmax=Tfin; nmin=find(t1>=tmin,1,'first'); 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')

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)
%Se ve la generación del 2º armónico al haberse duplicado el periodo

m=[];
for k=10: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
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}')
% Se observan puntos sueltos, lo cual es lógico al no estar en una zona
% caótica o autopulsada, sino que hay una oscilación de amplitud y
% frecuencias bien definidas.

L_FT1 =
        2001
L_m =
    28

close all

Published with MATLAB® 7.11