Le modele d'oscillateur le plus simple est l'oscillateur de phase $$\frac{d\phi}{dt} = \omega$$ La variabe $\phi$ représente la phase d'une trajectoire sur un cercle de rayon 1 en coordonnées polaires. En coordonnées cartésiennes, la trajectoire est $$(\cos(\omega t+\phi_{init}),\sin(\omega t+\phi_{init}))$$. L'oscillateur de phase a une periode $T = 2\pi/\omega$. Le paramètre $\omega$ est la fréquence naturelle ou intrinsèque de l'oscillateur. On peut rajouter un terme de forçage $I(t)$:
$$\frac{d\phi}{dt} = \omega + I(t)$$Si le terme $I(t)$ est suffisament petit, la solution reste oscillante ($\phi$ repasse toujours par zéro).
On s'intéresse maintenant à une population de $N$ oscillateurs de phase couplés entre eux. Chaque oscillateur a sa propre frequence naturelle $\omega_i$, pour $i=1,...,N$. Chaque oscillateur interagit avec tous les oscillateurs sinusoïdalement, pour donner
$$\frac{d\phi_i}{dt} = \omega_i + \frac{K}{N}\sum_{j=1}^{N} \sin(\phi_j-\phi_i)$$Le paramètre $K$ est la force de couplage et $N$ le nombre d'oscillateurs. Ce modèle est le modèle de Kuramoto.
On peut coupler chaque oscillateur avec seulement une partie des autres oscillateurs. En ce cas, la force de couplage devient une matrice $K$ de taille $N \times N$, avec $K_{i,j}$ la force de couplage de l'oscillateur $j$ sur l'oscillateur $i$ $$\frac{d\phi_i}{dt} = \omega_i + \sum_{j=1}^{N} K_{i,j} \sin(\phi_j-\phi_i).$$ Quand $K_{i,j} = K/N$, on retrouve le modèle de Kuramoto.
Le code kuramoto.m permet de simuler un système d'oscillateurs de phase et de tracer les solutions.
Le code est divisé en trois parties: initialisation, resolution, analyse. L'entête est standard.
function sol = kuramoto % KURAMOTO Systeme d'oscillateurs de phase couples % sol = kuramoto resoud un systeme de N d'oscillateurs de phase couples
% intervalle d'integration
t0 = 0;
tfinal = 500;
% nombre d'oscillateurs
N = 400;
% parametres des equations
% frequences intrinseques aleatoires
mean_freq = 1.0;
std_freq = 0.05;
omega = mean_freq + std_freq*randn(N,1);
% conditions initiales aleatoires uniforme sur le cercle
phi0 = 2*pi*rand(N,1);
% on definit la matrice de couplage
C = computeCouplingMatrix('local');
[t,phi]=solveKuramoto;
mean_field = mean(exp(1i*phi),1);
order_parameter = abs(mean(exp(1i*phi),1));
angle_parameter = angle(mean_field);
plotSolution('local')
% on met la solution dans une structure de sortie
sol.t=t;
sol.phi=phi;
sol.C = C;
sol.K = couplingStrength(t);
sol.mean_field = mean_field;
sol.order_parameter = order_parameter;
sol.angle_parameter = angle_parameter;
% FONCTIONS IMBRIQUEES----------------------------------------------------
solveKuramoto Résolution du système avec la méthode d'Euler. Quand la phase d'un des oscillateurs traverse $2\pi$, sa phase est remise à zéro. La fonction phaseEq est le membre de droite du système d'équations et la fonction coupling calcule le terme de couplage.
function [t,x] = solveKuramoto
% SOLVEKURAMOTO resoud le systeme avec Euler
% parametres de simulation
dt = 0.5;
nsteps = (tfinal-t0)/dt+1;
t = t0:dt:tfinal;
x = zeros(N,nsteps);
x(:,1) = phi0;
for ii = 1:nsteps-1
% integre avec Euler
x(:,ii+1) = x(:,ii) + dt * phaseEq(t(ii),x(:,ii));
% reinitialisation
ifire = (x(:,ii+1)>=2*pi);
x(ifire,ii+1) = x(ifire,ii+1) - 2*pi;
end
function dphidt = phaseEq(t,phi)
% PHASEEQ equations du systeme
dphidt = omega + couplingStrength(t)*coupling(phi);
end
function c = coupling(phi)
% COUPLING terme de couplage
phase_diff = bsxfun(@minus,phi',phi);
c = sum(C.*sin(phase_diff),2);
end
end
couplingStrength est la force de couplage $K$. Elle varie dans le temps, d'abord croissante jusqu'à atteindre la valeur maxK et rediminue ensuite.
function y = couplingStrength(t)
% COUPLINGSTRENGTH force de couplage variable
maxK = 0.4;
midT = (tfinal-t0)/2;
y = maxK/midT*t.*(t<=midT)+(2*maxK-maxK/midT*t).*(t>midT);
end
computeCouplingMatrix permet de générer une matrice de couplage. Pour $N$ oscillateurs, la matrice $C$ est une matrice $N \times N$. L'oscillateur $j$ affecte l'oscillateur $i$ avec une force $C_{i,j}$. Par défaut, la somme des lignes est normalisée à 1 si au moins un coefficient de la ligne est non-nul. Ce code n'est pas optimal pour générer un couplage de type local, mais permet de facilement construire d'autres types de couplages.
function C = computeCouplingMatrix(type)
% COMPUTECOUPLINGMATRIX calcule la matrice de couplage
% La matrice de couplage est une matrice NxN. L'entree
% M(i,j)=1 si l'oscillateur i recoit une connexion de l'oscillateur j. La
% matrice M est ensuite normalisee de telle sorte que la somme de chaque
% ligne soit 1, si il y a au moins un j non-nul, ou 0 si toute la ligne est
% nulle.
switch lower(type)
case 'local'
% Pour un couplage local, M est une matrice NxN
% On dispose les N oscillateurs sur une lattice (prendre N carre)
% On connecte deux oscillateurs si ils sont voisins (4 voisins, conditions
% de bord periodiques)
% xi position en x des oscillateurs;
% yi position en y des oscillateurs;
if sqrt(N)~=round(sqrt(N))
error('N doit etre carre');
end
sN = sqrt(N);
xi = mod((0:(N-1)),sN);
yi = reshape(repmat(0:(sN-1),sN,1),N,1);
distX = abs(bsxfun(@minus,xi',xi));
distX(distX==(sN-1))=1;
distY = abs(bsxfun(@minus,yi',yi));
distY(distY==(sN-1))=1;
dist = distX+distY;
M = dist<=1 & dist>0;
sumM = sum(M,2);
C = bsxfun(@rdivide,M,sumM);
C(isnan(C))=0;
C = sparse(C);
otherwise
% Pour un couplage global (par defaut), M reste
% un scalaire: M(i,j) = 1 pour tout i et j.
M=1;
C=M/N;
end
end
plotSolution trace la solution en temps réel. Pour le couplage local, trois panneaux sont affichés. À gauche, chaque oscillateur est un carré disposé sur une grille $N \times N$. La couleur de l'oscillateur représente sa phase: $\phi = 0$ en noir et $\phi = 2 \pi$ en blanc. Au milieu, la moyenne des sinus des oscillateurs en bleu et le paramètre d'ordre en rouge en fonction du temps. À droite, la position de chaque oscillateur sur le cercle unité en bleu et le centre de masse des oscillateurs en rouge. Les trois panneaux évoluent avec le temps.
function plotSolution(type)
% PLOTSOLUTION Trace la solution phi
figure(1); clf;
n = size(phi,2);
colormap gray;
switch lower(type)
case 'local'
sN = sqrt(N);
for i=1:n
subplot(131)
m = reshape(sin(phi(:,i)),sN,sN);
image((m+1)/2*64)
subplot(132)
hold off
plot(t(1:i),mean(sin(phi(:,1:i)),1))
hold on
plot(t(1:i),order_parameter(1:i),'r')
subplot(133)
theta = linspace(0,2*pi,100);
circle = exp(1i*theta);
hold off
plot(circle,'k')
hold on
plot(exp(1i*phi(:,i)),'o','MarkerFaceColor','b')
plot(mean_field(i),'*r');
pause(0.01)
end
otherwise
theta = linspace(0,2*pi,100);
circle = exp(1i*theta);
for i = 1:n
hold off
plot(circle,'k')
hold on
plot(exp(1i*phi(:,i)),'o','MarkerFaceColor','b')
hold on
plot(mean_field(i),'*r');
pause(0.0001)
end
end
end
% FIN FONCTIONS IMBRIQUEES---------------------------------------------------- end
Lancez une simulation depuis la Command Window.
sol = kuramoto;
sol =
t: [1x1001 double]
phi: [400x1001 double]
C: [400x400 double]
K: [1x1001 double]
mean_field: [1x1001 double]
order_parameter: [1x1001 double]
angle_parameter: [1x1001 double]
Effectuez des simulations avec les paramètres par défaut de kuramoto. Interprétez ce que vous voyez. Changez certain des paramètres: type de couplage, taille du système (pas trop gros !), etc.
Simulez le modèle en gardant la fréquence moyenne à 1.0.
% frequences intrinseques aleatoires mean_freq = 1.0; std_freq = 0.05; omega = mean_freq + std_freq*randn(N,1);
Maintenant fixez la fréquence moyenne à 0.0: mean_freq = 0.0;, et simulez. Qu'observez-vous? Gardez la fréquence moyenne à 0.0 pour les simulations suivantes.
Modifier l'interaction sin la fonction coupling par une somme de fonction périodiques, de façcon à avoir une fonction de couplage de la forme
Fixez $a = 0.5$ et $b = 1.0$, et simulez. Qu'observez-vous?
Prenez N = 81; et modifiez les conditions initiales des oscillateurs pour que les phases forment un X sur la grille de simulation.
% conditions initiales en X bruité
phi0 = [1 1 -1 -1 -1 -1 -1 1 1;
-1 1 1 -1 -1 -1 1 1 -1;
-1 -1 1 1 -1 1 1 -1 -1;
-1 -1 -1 1 1 1 -1 -1 -1;
-1 -1 -1 -1 1 -1 -1 -1 -1;
-1 -1 -1 1 1 1 -1 -1 -1;
-1 -1 1 1 -1 1 1 -1 -1;
-1 1 1 -1 -1 -1 1 1 -1;
1 1 -1 -1 -1 -1 -1 1 1];
phi0 = reshape(phi0,N,1) + 2*pi/20*randn(N,1);
Modifiez la fonction couplingStrength pour avoir une force de couplage constante.
function y = couplingStrength(t)
% COUPLINGSTRENGTH force de couplage variable
y = 0.3;
end
Simulez. Avec un peu de chance vous reconnaîtrez le X, qui devrait rester visible. Si ça ne marche pas, quels paramètres pouvez vous changer pour améliorer la stabilité?