Solving the mean field Kuramoto model with NDSolve

Question

I would like to numerically solve a system of differential equations that describes the dynamics of $N$ coupled oscillatory units (Kuramoto model) via their phase variables $\phi_i$:

$\frac{\partial\phi_i}{\partial t} = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin(\phi_j-\phi_i)$

A simple implementation with NDSolve works nicely:

tfinal = 20;
stdϕ = 2.5;
n = 100;
k = 1.2;

SeedRandom[0];
ϕ0s = Mod[RandomVariate[NormalDistribution[0, stdϕ], n], 2 π];
ω = RandomVariate[NormalDistribution[1, 0.1], n];

trajs = NDSolve[{
  Table[{
    ϕ[i]'[t] == ω[[i]] + k/n Sum[Sin[ϕ[j][t] - ϕ[i][t]], {j, n}],
    ϕ[i][0] == ϕ0s[[i]]
  }, {i, n}]
  },
  Table[ϕ[i], {i, n}], {t, 0, tfinal}][[1]];

For suitable coupling strengths, the oscillators synchronize as shown in the following animation:

It is possible to simplify the differential equations via a mean field interaction.

Mean field:

$R e^{\text{i}\psi} = \frac{1}{N} \sum_{j=1}^N e^{\text{i}\phi_j}$

Mean field coupling:

$\frac{\partial\phi_i}{\partial t} = \omega_i + R K \sin(\psi-\phi_i)$

The corresponding mathematica code for this differential-algebraic equation is:

trajs = NDSolve[{
  Table[{
    ϕ[i]'[t] == ω[[i]] +  Abs[r[t]] k Sin[Arg[r[t]] - ϕ[i][t]],
    ϕ[i][0] == ϕ0s[[i]]
   }, {i, n}],
   r[t] == Sum[Exp[I ϕ[i][t]], {i, n}]/n
   },
   Join[Table[ϕ[i], {i, n}], {r}], {t, 0, tfinal}][[1]];

Problem: Instead of a solution that is computed in less time, NDSolve throws an error:

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

I also tried adding initial conditions for both function values and derivatives, but that did not help. Different "EquationSimplification" options (as listed here) also did not help.

How can I solve this?

---

Edit

Here is the code for the animation:

animationFrames = Table[
   Show[
    ParametricPlot[
     Evaluate@ReIm[Sum[Exp[I ϕ[i][t] /. trajs], {i, n}]/n], {t, 
      0, tt}, PlotStyle -> Lighter@ColorData[97, 1], 
     PlotRange -> 1.1 {{-1, 1}, {-1, 1}}, Axes -> None, 
     ImageSize -> {Automatic, 300}],
    Graphics[{
      Gray, AbsoluteThickness[2], Circle[],
      Orange, 
      Disk[#, 0.03] & /@ 
       Table[{Cos[#], Sin[#]} &@ϕ[i][tt] /. trajs, {i, n}],
      ColorData[97, 1], 
      Disk[ReIm[
        Sum[Exp[I Evaluate[ϕ[i][tt] /. trajs]], {i, n}]/n], 0.1]
      }]
    ]
   , {tt, 0.001, tfinal, 0.01 tfinal}];


Export["kuramoto_sync.gif", animationFrames, 
 AnimationRepetitions -> Infinity]

Bonus: Is there a way to "precompute" the trajectory of the centroid (blue bobble) and store it in an interpolating function, so it's not required to calculate ReIm[Sum[Exp[I ϕ[i][t] /. trajs], {i, n}]/n] in each frame? Of course, one option is to sample the centroid trajectory in advance like:

{rx, ry} = 
  Table[ListInterpolation[ri, {0, tfinal}], {ri, 
    Table[ReIm[
       N@Sum[Exp[I Evaluate[ϕ[i][t] /. trajs]], {i, n}]/n], {t, 
       0, tfinal, 0.001 tfinal}]\[Transpose]}];

Answer 1

The DAE solver is generally weaker than the ODE solver of NDSolve, at least now, so the simplest way to resolve your problem is avoiding DAE system. This can be achieved by substituting definition of r[t] into the rest part:

Table[{ϕ[i]'[t] == ω[[i]] + 
     k Abs[r[t]] Sin[Arg[r[t]] - ϕ[i][t]], ϕ[i][0] == ϕ0s[[i]]}, {i, 
   n}] /. r[t] -> Sum[Exp[I ϕ[i][t]], {i, n}]/n

This isn't the end, because Mathematica will quit in v9.0.1, or spit out ndinid warning in v11.2 and v11.3 and fail when one tries to solve the new system with NDSolve.

So there probably exists a bug in NDSolve, but in which part? Can we circumvent it? Luckily a hint is given in the warning message:

The Context of this variable is discretevariableNDSolve`, so the bug may be something related to DiscontinuityProcessing in NDSolve. (I associate the variable name with this option because we already found several bugs related to it in this site. ) Let's try turning it off:

tfinal = 20;
stdϕ = 2.5;
n = 10;
k = 1.2;

SeedRandom[0];
ϕ0s = Mod[RandomVariate[NormalDistribution[0, stdϕ], n], 2 π];
ω = RandomVariate[NormalDistribution[1, 0.1], n];

trajs2 = NDSolveValue[
   Table[{ϕ[i]'[t] == ω[[i]] + 
        k Abs[r[t]] Sin[Arg[r[t]] - ϕ[i][t]], ϕ[i][0] == ϕ0s[[i]]}, {i, 
      n}] /. r[t] -> Sum[Exp[I ϕ[i][t]], {i, n}]/n, ϕ /@ Range@n, {t, 0, 
    tfinal}, Method -> {DiscontinuityProcessing -> False}];

trajs2 // ListLinePlot

OK, now we confirm another bug of DiscontinuityProcessing .

---

Update

We can make NDSolve calculate r[t] for us by differentiating definition of r[t] to transform it into an ODE. Notice we no longer need /. r[t] -> … in this approach:

sum = Sum[Exp[I ϕ[i][t]], {i, n}]/n;

trajs3 = NDSolveValue[
 {Table[{ϕ[i]'[t] == ω[[i]] + k Abs[r[t]] Sin[Arg[r[t]] - ϕ[i][t]], 
         ϕ[i][0] == ϕ0s[[i]]}, {i, n}], 

   r'[t] == D[sum, t], 
   r[0] == sum /. ϕ[i_][t] :> ϕ0s[[i]]}, 

 {ϕ /@ Range@n, r}, {t, 0, tfinal}, Method -> DiscontinuityProcessing -> False]

If[$VersionNumber < 10, ReIm = {Re@#, Im@#} &];
ParametricPlot[ReIm@trajs3[[2]][t], {t, 0, tfinal}]