function [X, outcome] =AndrzejakChaosMain(N, b, alphaD, Deltaalpha, epsilon, couplingonset, X0, WindowLength, NumberWindows, dt) % % This source code allows you to run the simulations presented in our paper % [1] Andrzejak, Ralph Gregor, Giulia Ruzzene, and Irene Malvestio. % Generalized synchronization between chimera states." Chaos; 27 (5): 053114. (2017). % % In general, we follow the notation of [1]. Only the auxiliary response is called Theta instead of Psi' % % Input % N: Number of oscillators in each network (50 in [1]) % b: broadness of coupling kernel (18 in [1]) % alphaD: phase lag parameter of the driver network (1.46 in [1]) % Deltaalpha: difference between the phase lag parameter of the driver and response network. alphaR = alphaD+Deltaalpha (ranged from 0 to 0.08 in [1]) % epsilon: coupling strength (ranged from 0 to 1 in [1]) % couplingonset: number of sample points at which the coupling is turned on from zero to epsilon (100,0000 in [1]) % X0: initial conditions in a vector of size (3*N,1). (random intial conditions uniformly distributed between 0 and 2*pi in [1], using rand(150,1)*2*pi) % % NumberWindows: Total number of windows -of WindowLength samples each- used for simulation (250 in [1]) % WindowLength: Number of samples in each window. The last window is given as output X. (20,000 used in [1]) % dt sampling time used for the numerical integration (0.01 in [1]) % Taken these three parameters, we have 250 windows * 20,000 samples/window = 5,000,000 samples. Corresponding to 50,000 dimensionless time units % % Output % X: phases of all three network in one joint matrix of size: (3N,WindowLength); % In X: 1:N phases of driver; N+1:2*N phases of response; 2*N+1:3*N phases of auxiliary response. % outcome: string explaining what happened for this realization % The three networks (driver Phi, response Psi, and auxiliary response Theta) % are generated using one set of differential equations. We here set the % indices to address them: Phi = 1:N; Psi = N+1:2*N; Theta = 2*N+1:3*N; % Writing the phase lag parameters into a matrix of size (3*N,3*N). % This matrix representation is used for computational efficiency. alphaR = alphaD+Deltaalpha; alphaV = zeros(3*N); alphaV(Phi,Phi) = alphaD; alphaV(Psi,Psi) = alphaR; alphaV(Theta,Theta) = alphaR; % Defining the coupling kernel matrix G Rectangularwindow = zeros(N,1); Rectangularwindow(1:b+1)=1; Rectangularwindow(N-b+1:N)=1; G = zeros(3*N,3*N); for i = 1:round(N) G(i,1:N)=circshift(Rectangularwindow,i-1); G(i+N,Psi)=G(i,Phi); G(i+2*N,Theta)=G(i,Phi); end G = G/2/b; % Allocating the memory space for the phases of all three networks. X = zeros(WindowLength,N*3); % Setting the last state of the network to the initial conditions. x = X0; % Specifying the function which contains the differential equation of the % entire network (driver & response & auxiliary response). Included below. deriv = @AndrzejakChaosNetworkDifferentialEquation; % default entry for outcome outcome = 'No chimera collapsed. No synchronization occurred'; TotalSamples = WindowLength*NumberWindows; for j =1:TotalSamples % one-time activation of coupling if j==couplingonset; for i = 1:round(N) G(i+N,i)=epsilon; G(i+2*N,i)=epsilon; end end % Runge-Kutta integration of order four k1 = dt*feval(deriv,x,G,alphaV); k2 = dt*feval(deriv,x+k1/2,G,alphaV); k3 = dt*feval(deriv,x+k2/2,G,alphaV); k4 = dt*feval(deriv,x+k3,G,alphaV); x = x + k1/6+k2/3+k3/3+k4/6; % wrapping the phases x = mod(x,2*pi); % Writing the phases to the output X using the modulus of Windowlength. % In this way, always the last WindowLength samples are in X. X(mod(j-1,WindowLength)+1,:)=x; % What follows all serves to check for possible chimera collapses in one of the % networks or onsets of synchronization. This is done only every 100th % sample to reduce the computational load. % Order parameter of driver network to check if its chimera has collapsed Z = 1/N*abs(sum(exp(1i*x(Phi)))); if Z >0.9999 outcome='Driver network collapsed to fully synchronized state.'; break end % Order parameter of response network to check if its chimera has collapsed Z = 1/N*abs(sum(exp(1i*x(Psi)))); if Z >0.9999 outcome='Response network collapsed to fully synchronized state.'; break end % Order parameter of auxiliary response network to check if its chimera has collapsed Z = 1/N*abs(sum(exp(1i*x(Theta)))); if Z >0.9999 outcome='Auxiliary response network collapsed to fully synchronized state.'; break end deltaphipsi=mean(abs(sin((x(Phi)-x(Psi))/2))); if deltaphipsi == 0 outcome='Identical synchronization between driver and response.'; break end deltaphitheta=mean(abs(sin((x(Phi)-x(Theta))/2))); if deltaphitheta == 0 outcome='Identical synchronization between driver and auxiliary response.'; break end deltapsitheta=mean(abs(sin((x(Psi)-x(Theta))/2))); if deltapsitheta == 0 outcome='Identical synchronization between response and auxiliary response => Generalized synchronization between driver and response.'; break end end % Make sure that the last point is always at the right end of the window. % See above how the data is written circularly into X. X = circshift(X,-(mod(j,WindowLength))); disp(outcome); function dy = AndrzejakChaosNetworkDifferentialEquation(x,G,alpha) diffm=bsxfun(@minus,x,x'); dy = -sum(sin(diffm+alpha).*G,2);