synchro.m

%Do the simulation and compute a few things from  the traces

Ten = 24;%entrainement period
Tin = 12;%intrinsic period

m = @(th,t) 2*pi/Tin + 0.2*sin(th) * (0.5*cos( 2*pi/Ten*t) + 0.5)^10;

Nt = 1e5;
dt = 1e-2;
T = linspace(0,Nt*dt,Nt)';

th = zeros(Nt,1);
th(1) = 2*pi*rand;

phi = 2*pi/Ten*T;

for t=1:Nt-1
        
    th(t+1) = th(t) + dt*m(th(t),dt*t);        
end

%compute effective period : w = lim theta / t
weff = mean(th(end-1000:end)./T(end-1000:end));
Teff = 2*pi/weff;

%is the phase difference growing like T ?
phaseDiffVsT = abs(th-phi)./T;
if(  mean(phaseDiffVsT(end-100:end)) > 1e-2)
   fprintf('\nNo 1:1 sychronization!\n') 
else
   fprintf('\n 1:1 Sychronization!\n') 
end

phaseDiff = mod( mean(th(end-1000:end)-phi(end-1000:end)), 2*pi);

fprintf('effective period:\t %2.2f\n',Teff)
fprintf('rotation number:\t  %2.2f\n',Teff/Ten)
fprintf('mean phase difference:\t %2.2f\n',phaseDiff)

clf
subplot(2,1,1)
plot(T,[th phi])
subplot(2,1,2)
semilogy(T,(abs(th-phi)./T))


%% Compute the map F

% F gives the value of theta_n+1 given theta_n. 
% n correspond to the time t_n : phi(t_n) = 2*n*pi

dt = 5e-3;
Nt = round( Ten/dt );

N = 200;
initialConds = linspace(0,2*pi,N);

F = zeros(N,1);
th = zeros(Nt,1);

for k=1:N

    th(1) = initialConds(k);

    for t=1:Nt-1

        th(t+1) = th(t) + dt*(m(th(t),dt*t));

    end
    
    F(k) = th(end);
    
end

clf
plot(initialConds,F)
xlabel('theta n')
ylabel('theta n+1')

%% Iterate the map F with different initial values

Nit = 50;

Fs = zeros(Nit,1);

clf; hold on

for trial=1:20

    Fs(1) = 2*pi*rand;

    for k=2:Nit

        Fs(k) = interp1( initialConds,F, Fs(k-1) );
        Fs(k) = mod(Fs(k),2*pi);

    end

    plot(Fs)

end


%%