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]}];