%      SyncVanDerPol.M
%   
%  Synchronisation of van-der-Pol oscillators
%  
% J. Lunze
% 
% 5.9.2018
% Version v. 5.9.2018
% for 2nd edition: 3.6.2021
%
echo off
clear
close all

%%  Behaviour of van-der-Pol oscillators (with u=0)
%
figure(1)
subplot(4,6,[1 2 7 8]);
epsilon=1.5;
Tend=20;
dt=0.01;
Time=[0:dt:Tend];
%
X0=[0.2 0];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-4 4 -4 4]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latextitle('$${\tt SyncVanDerPol:}$$ A single van-der-Pol oscillator');
%latexxlabel('$$x$$');
%latexylabel('$$\dot{x}\;$$');
%rotateY;
%latextext(-3.7, 4.5, ['$$\varepsilon=\; $$', num2str(epsilon)]);
hold off
%
subplot(4,6, [3 4 9 10]);
X0=[0.92 3.243];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-4 4 -4 4]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latexxlabel('$$x$$');
yleer;
hold off
%
subplot(4,6, [5 6 11 12]);
X0=[2 -3];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-4 4 -4 4]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latexxlabel('$$x$$');
yleer;
hold off


%%  Synchronisation of van-der-Pol oscillators
%   3 oscillators in a directed ring
%
figure(2)
Tend=80;
dt=0.002;
Time=[0:dt:Tend];
epsilon=1.5;
c=0.12;    % 125
%   initial states
X10=[-3 1];
X20=[2 3];
X30=[3 0];  % 0 3
X=[];
X(1,:)=[X10 X20 X30];   % row vector includes the states of all oscillators
%
for k1=2:length(Time)
    X(k1,1:2)=X(k1-1,1:2)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
    X(k1,3:4)=X(k1-1,3:4)+...
        [X(k1-1,4)*dt -epsilon*X(k1-1,4)*(X(k1-1,3)*X(k1-1,3)-1)*dt-X(k1-1,3)*dt];
    X(k1,5:6)=X(k1-1,5:6)+...
        [X(k1-1,6)*dt -epsilon*X(k1-1,6)*(X(k1-1,5)*X(k1-1,5)-1)*dt-X(k1-1,5)*dt];
    %  controller
    u1=-c*(X(k1-1,2)-X(k1-1,6))*dt;   % u depends upon xdot
    u2=-c*(X(k1-1,4)-X(k1-1,2))*dt;
    u3=-c*(X(k1-1,6)-X(k1-1,4))*dt;
    X(k1,:)=X(k1,:)+[0 u1 0 u2 0 u3];
end
plot(X(:,1), X(:,2), 'b', X(:,3), X(:,4), 'r', X(:,5), X(:,6), 'k');
hold on
plot(X(1,1), X(1,2), 'bo', X(1,3), X(1,4), 'ro', X(1,5), X(1,6), 'ko');
axis([-3.5 3.5 -3.5 3.5]);
axis('square');
yticks([-3 0 3]);
xticks([-3 0 3]);
latextitle('$${\tt SyncVanDerPol:}$$ Synchronisation of van-der-Pol oscillators');
latexxlabel('$$x$$');
latexylabel('$$\dot{x}$$');
rotateY;
latextext(-3, 3, ['$$\varepsilon=\; $$', num2str(epsilon)]);
latextext(-3, 2.5, ['$$c=\; $$', num2str(c)]);
hold off
%
%   Synchronisation error
figure(3)
subplot(2,1,1)
plot(Time, X(:,1)-X(:,3), 'b');
latexylabel('$$x_1-x_3$$');
rotateY;
xleer;
axis([0 Tend -5 4]);
yticks([-4 0 4]);
subplot(2,1,2)
plot(Time, X(:,1), 'b');
latexylabel('$$x_1\;$$');
rotateY;
latexxlabel('$$t$$');
axis([0 Tend -3 3]);
yticks([-3 0 3]);
xticks([0 20 40 60 80 100]);


%%  Synchronisation of van-der-Pol oscillators
%   3 completely coupled oscillators
%
figure(4)
Tend=100;
dt=0.001;
Time=[0:dt:Tend];
epsilon=0.4;
c=0.5;
%   initial states
X10=[3 -3]/10;
X20=[-3 -1]/10;
X30=[0 3]/10;
X=[];
X(1,:)=[X10 X20 X30];   % row vector includes the states of all oscillators
%
for k1=2:length(Time)
    X(k1,1:2)=X(k1-1,1:2)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
    X(k1,3:4)=X(k1-1,3:4)+...
        [X(k1-1,4)*dt -epsilon*X(k1-1,4)*(X(k1-1,3)*X(k1-1,3)-1)*dt-X(k1-1,3)*dt];
    X(k1,5:6)=X(k1-1,5:6)+...
        [X(k1-1,6)*dt -epsilon*X(k1-1,6)*(X(k1-1,5)*X(k1-1,5)-1)*dt-X(k1-1,5)*dt];
    %  controller
    u1=-c*(X(k1-1,1)-X(k1-1,5)+X(k1-1,1)-X(k1-1,3))/2*dt;
    u2=-c*(X(k1-1,3)-X(k1-1,1)+X(k1-1,3)-X(k1-1,5))/2*dt;
    u3=-c*(X(k1-1,5)-X(k1-1,3)+X(k1-1,5)-X(k1-1,1))/2*dt;
    X(k1,:)=X(k1,:)+[0 u1 0 u2 0 u3];
end
plot(X(:,1), X(:,2), 'b', X(:,3), X(:,4), 'r', X(:,5), X(:,6), 'k');
hold on
plot(X(1,1), X(1,2), 'bo', X(1,3), X(1,4), 'ro', X(1,5), X(1,6), 'ko');
%axis([-4 4 -4 4]);
axis('square');
yticks([-3 0 3]);
xticks([-3 0 3]);
latextitle('$${\tt SyncVanDerPol:}$$ Synchronisation of van-der-Pol oscillators');
latexxlabel('$$x$$');
latexylabel('$$\dot{x}$$');
rotateY;
hold off
%
%   Synchronisation error
figure(5)
plot(Time, X(:,1)-X(:,3), 'b');
hold off

%%  Synchronisation of linear oscillators
%   3 oscillators in a directed ring
%
figure(6)
Tend=70;
dt=0.001;
Time=[0:dt:Tend];
epsilon=0;
c=0.1;    % 125
%   initial states
X10=[3 -3];
X20=[-3 -1];
X30=[0 3];  % 0 3
X=[];
X(1,:)=[X10 X20 X30];   % row vector includes the states of all oscillators
%
for k1=2:length(Time)
    X(k1,1:2)=X(k1-1,1:2)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*dt-X(k1-1,1)*dt];
    X(k1,3:4)=X(k1-1,3:4)+...
        [X(k1-1,4)*dt -epsilon*X(k1-1,4)*dt-X(k1-1,3)*dt];
    X(k1,5:6)=X(k1-1,5:6)+...
        [X(k1-1,6)*dt -epsilon*X(k1-1,6)*dt-X(k1-1,5)*dt];
    %  controller
    u1=-c*(X(k1-1,2)-X(k1-1,4))*dt;   % u depends upon xdot!
    u2=-c*(X(k1-1,4)-X(k1-1,2))*dt;
    u3=-c*(X(k1-1,6)-X(k1-1,4))*dt;
    X(k1,:)=X(k1,:)+[0 u1 0 u2 0 u3];
end
plot(X(:,1), X(:,2), 'b', X(:,3), X(:,4), 'r', X(:,5), X(:,6), 'k');
hold on
plot(X(1,1), X(1,2), 'bo', X(1,3), X(1,4), 'ro', X(1,5), X(1,6), 'ko');
%axis([-4 4 -4 4]);
axis('square');
yticks([-3 0 3]);
xticks([-3 0 3]);
latextitle('$${\tt SyncVanDerPol:}$$ Synchronisation of linear oscillators');
latexxlabel('$$x$$');
latexylabel('$$\dot{x}$$');
rotateY;
hold off
%
%   Synchronisation error
figure(7)
subplot(2,1,1)
plot(Time, X(:,3)-X(:,5), 'b');
hold off
subplot(2,1,2)
plot(Time, X(:,5), 'b');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Synchronisation of van-der-Pol oscillators with circular limit cycle
%  Behaviour of van-der-Pol oscillators (with u=0)
%
figure(8)
subplot(4,6,[1 2 7 8]);
epsilon=0.04;
Tend=30;
dt=0.1;
Time=[0:dt:Tend];
%
X0=[0.2 0];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-5 5 -5 5]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latextitle('$${\tt SyncVanDerPol:}$$ A single van-der-Pol oscillator');
%latexxlabel('$$x$$');
%latexylabel('$$\dot{x}\;$$');
%rotateY;
%latextext(-3.7, 4.5, ['$$\varepsilon=\; $$', num2str(epsilon)]);
hold off
%
subplot(4,6, [3 4 9 10]);
X0=[0.92 3.243];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-5 5 -5 5]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latexxlabel('$$x$$');
yleer;
hold off
%
subplot(4,6, [5 6 11 12]);
X0=[2 -3];   % row vector
X=[];
X(1,:)=X0;
for k1=2:length(Time)
    X(k1,:)=X(k1-1,:)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(1,1), X(1,2), 'bo');
axis([-5 5 -5 5]);
axis('square');
%yticks([-3 0 3]);
%xticks([-3 0 3]);
xleer;
yleer;
%latexxlabel('$$x$$');
yleer;
hold off



%%   3 oscillators in a directed ring
%
figure(9)
Tend=60;
dt=0.002;
Time=[0:dt:Tend];
c=0.28;    % 125
%   initial states
X10=[-3 1];
X20=[2 3];
X30=[3 0];  % 0 3
X=[];
X(1,:)=[X10 X20 X30];   % row vector includes the states of all oscillators
%
for k1=2:length(Time)
    X(k1,1:2)=X(k1-1,1:2)+...
        [X(k1-1,2)*dt -epsilon*X(k1-1,2)*(X(k1-1,1)*X(k1-1,1)-1)*dt-X(k1-1,1)*dt];
    X(k1,3:4)=X(k1-1,3:4)+...
        [X(k1-1,4)*dt -epsilon*X(k1-1,4)*(X(k1-1,3)*X(k1-1,3)-1)*dt-X(k1-1,3)*dt];
    X(k1,5:6)=X(k1-1,5:6)+...
        [X(k1-1,6)*dt -epsilon*X(k1-1,6)*(X(k1-1,5)*X(k1-1,5)-1)*dt-X(k1-1,5)*dt];
    %  controller
    u1=-c*(X(k1-1,2)-X(k1-1,6))*dt;   % u depends upon xdot
    u2=-c*(X(k1-1,4)-X(k1-1,2))*dt;
    u3=-c*(X(k1-1,6)-X(k1-1,4))*dt;
    X(k1,:)=X(k1,:)+[0 u1 0 u2 0 u3];
end
plot(X(:,1), X(:,2), 'b', X(:,3), X(:,4), 'r', X(:,5), X(:,6), 'k');
hold on
plot(X(1,1), X(1,2), 'bo', X(1,3), X(1,4), 'ro', X(1,5), X(1,6), 'ko');
axis([-3.5 3.5 -3.5 3.5]);
axis('square');
yticks([-3 0 3]);
xticks([-3 0 3]);
latextitle('$${\tt SyncVanDerPol:}$$ Synchronisation of van-der-Pol oscillators');
latexxlabel('$$x$$');
latexylabel('$$\dot{x}$$');
rotateY;
latextext(-3, 3, ['$$\varepsilon=\; $$', num2str(epsilon)]);
latextext(-3, 2.5, ['$$c=\; $$', num2str(c)]);
hold off
%
%   Synchronisation error
figure(10)
subplot(2,1,1)
plot(Time, X(:,1)-X(:,3), 'b');
latexylabel('$$x_1-x_3$$');
rotateY;
xleer;
axis([0 Tend -5 4]);
yticks([-4 0 4]);
subplot(2,1,2)
plot(Time, X(:,1), 'b');
latexylabel('$$x_1\;$$');
rotateY;
latexxlabel('$$t$$');
axis([0 Tend -3 3]);
yticks([-3 0 3]);
xticks([0 20 40 60 80 100]);



%%   Figures
epsfigc('SyncVanDerPol');