Finding a single equilibrium point of the Kuramoto model

Question

The Kuramoto model is used to describe a large variety of synchronization phenomena. It describes the time evolution of a group of n oscillators with $\omega_{i} \in\mathbb{R}, \frac{K}{n} > 0$: $$ \dot{\theta_{i}} = \omega_{i} -\frac{K}{n}\sum_{j=1}^{n}\sin (\theta_{i} - \theta_{j}), \,\,\,i \in \{ 1, \dots , n \}, $$ Is there a simple way to find a single equilibrium point of this system?

There are existing methods of solving for all equilibrium points using numerical algebraic geometry (see the linked question), but I really just need to find a particular solution. Of course, an analytical solution would be nice, but numerical methods are also fine as long as they can perform well when N is about 100.

For the numerical algebraic geometry method, see this linked question here. I am not familiar with numerical algebraic geometry, so I don't exactly know how to apply the method in the paper. Also, when $\omega_i$ are all the same, this paper provides an solution to find a particular solution via finding eigenvectors

Answer 1

I'll rewrite your system as $$ \sum_j \sin(\theta_i-\theta_j)=w_i $$ with $w_i=n\omega_i/K$. Note that $\sum_i w_i=0$ is a trivial necessary condition for solvability. If it fails, you have no chance to get any equilibria at all. Another trivial necessary condition is $|w_i|\le n$.

The following algorithm (dampened Newton with the initial guess of all $\theta_i=0$) is uncannily efficient (I cannot explain why, unfortunately) in the sense that I always got a solution satisfying the equations essentially with the machine precision when I knew that some solution existed, i.e., I set $w_i=\sum_j\sin(\Theta_i-\Theta_j)$ with some $\Theta_i\in[0,2\pi]$ (though, of course, it was not guaranteed to be the same I knew) and I got it more often than not when $w_i$ were chosen randomly but $|w_i|$ stayed reasonably below $n$ (<$0.7 n$ was fine but $0.8 n$ could give me trouble; on the other hand, then I just did not know whether it was the failure of the algorithm or just the absence of solutions). So, you are welcome to try it and see if it works for you. The language is Asymptote and the code should be more or less self-explanatory.

srand(seconds()); //seeds the randomness

int n=100; //the number of variables


real[] TH; //the suggested values
for(int k=0;k<n;++k) TH[k]=2*pi*unitrand(); 
//Try to change to anything you want.
//If you do pi/2*unitrand instead, it seems like you usually get th=TH (up to shift),
//otherwise it is another solution more often than not


real[] w;

//////////// Assigning the RHS to be random with 0 sum ///////////
real a=0.7*n;
for(int k=0;k<n;++k) w[k]=a*(2*unitrand()-1);
real s=sum(w)/n;
for(int k=0;k<n;++k) w[k]-=s;
//////////////////////////////////////////////////////////////////


/////////// Assigning the RHS to be the values at TH[i] ///////////
for(int k=0;k<n;++k)
{
w[k]=0;
for(int j=0;j<n;++j) w[k]+=sin(TH[k]-TH[j]);
}
////////////////////////////////////////////////////////////////////


/////////// Disable one of the above for testing //////////////////


real[] th; //the solution

for(int k=0;k<n;++k) th[k]=0; //initial guess of all 0


int M=100; //maximal number of iterations in dampened Newton

real[][] A; //The differential + dampening matrix 
real[][] B; //The discrepancy vector augmented with zeroes

for(int m=0;m<M;++m)
{
for(int i=0;i<n;++i)
{
B[i]=new real[]{-w[i]};
for(int j=0;j<n;++j) B[i][0]+=sin(th[i]-th[j]); //computing the discrepancy

A[i]=new real[]; 
for(int j=0;j<n;++j) {A[i][j]=-cos(th[i]-th[j]); }
A[i][i]=0;
A[i][i]=-sum(A[i]); 
}
//computing the differential part  

real s=0; 
for(int k=0;k<n;++k) s+=(B[k][0])^2/n;
s=sqrt(s);
//computing average discrepancy


real ss=(s/n)^(1-max(m-M/5,0)/M/0.8); 
for(int i=0; i<n;++i) 
{
A[n+i]=new real[];
for(int j=0;j<n;++j) A[n+i][j]=0;
A[i+n][i]=n*ss; 
B[n+i]=new real[]{0};} 
//computing the dampening part 
//this way of dampening is useful only when you finish near a positive local minimum
//of the sum of squares of the discrepancies (near-equilibrioum). If you are looking only
//for an exact equilibrium or nothing, you can discard it


real[][] I;
for(int i=0;i<n;++i)
{
I[i]=new real[];
for(int j=0; j<n;++j) I[i][j]=(i==j?n^2/10.0^14:0);
}
//small multiple of the identity matrix to avoid possible degeneracy

real[][] AA=inverse(I+transpose(A)*A)*transpose(A); 
real[][] V=AA*B;
//solving the (dampened) Newton equations by least squares

real sh=0;
for(int i=0;i<n;++i) {th[i]-=V[i][0]; sh+=V[i][0]^2/n;} //updating thetas
sh=sqrt(sh); //average size of the update shifts

//write(m,s,sh); //enable this line if you want to trace the iterations

if(sh<1/10.0^15) m=M; //stopping when the update shifts are below the machine precision
}

real S=0;
for(int i=0;i<n;++i) S+=B[i][0]^2/n;
S=sqrt(S);
write(S, sin(th[0]-th[1]), sin(TH[0]-TH[1])); pause(); 
//reporting the average discrepancy in the equations and 
//comparing th to TH (if TH initialization is used).
//You can change it to anything you want to know in the end.