Stabilization of direct numerical simulation for finite truncations of circular cumulant expansions

TL;DR

This paper investigates numerical instabilities in direct simulations of truncated circular cumulant equations for coupled oscillators and proposes stabilization methods effective near the Ott-Antonsen manifold.

Contribution

It introduces two stabilization approaches for direct numerical simulations of truncated circular cumulant equations and evaluates their effectiveness.

Findings

Stabilization techniques work well near the Ott-Antonsen manifold.

Effectiveness decreases as deviation from the Ott-Antonsen manifold increases.

Methods enable stable simulations with up to 20 cumulants.

Abstract

We study a numerical instability of direct simulations with truncated equation chains for the "circular cumulant" representation and two approaches to its suppression. The approaches are tested for a chimera-bearing hierarchical population of coupled oscillators. The stabilization techniques can be efficiently applied without significant effect on the natural system dynamics within a finite vicinity of the Ott-Antonsen manifold for direct numerical simulations with up to 20 cumulants; with increasing deviation from the Ott-Antonsen manifold the stabilization becomes more problematic.

Equations

$${\dot{\phi }}_k=\Omega (t)+\mathrm{Im}(2h(t)e^{-i{\phi }_k}),k=1,...,N,$$

$${\dot{\phi }}_k={\Omega }_k+\mathrm{Im}(2h(t)e^{-i{\phi }_k})+\sigma {\xi }_k(t),k=1,...,N$$

$${\dot{\varkappa }}_j=j(i{\Omega }_0-\gamma ){\varkappa }_j+h{\delta }_{1j}-h^{\ast }(j^2{\varkappa }_{j+1}+j\sum _{m=1}^j{\varkappa }_{j-m+1}{\varkappa }_m)-{\sigma }^2(j^2{\varkappa }_j+j\sum _{m=1}^{j-1}{\varkappa }_{j-m}{\varkappa }_m),$$

$${\dot{\phi }}_k=\Omega +\frac{1+A}{2N}\sum _{j=1}^N\sin ({\phi }_j-{\phi }_k-\alpha )+\frac{1-A}{2N}\sum _{j=1}^N\sin ({\psi }_j-{\phi }_k-\alpha )+\sigma {\xi }_k(t),$$

$${\dot{\psi }}_k=\Omega +\frac{1+A}{2N}\sum _{j=1}^N\sin ({\psi }_j-{\psi }_k-\alpha )+\frac{1-A}{2N}\sum _{j=1}^N\sin ({\phi }_j-{\psi }_k-\alpha )+\sigma {\eta }_k(t),$$

$${\dot{\varkappa }}_j=j(i{\Omega }_0-\gamma ){\varkappa }_j+h_{\phi }{\delta }_{1j}-h_{\phi }^{\ast }(j^2{\varkappa }_{j+1}+j\sum _{m=1}^j{\varkappa }_{j-m+1}{\varkappa }_m)-{\sigma }^2(j^2{\varkappa }_j+j\sum _{m=1}^{j-1}{\varkappa }_{j-m}{\varkappa }_m),$$

$${\dot{\kappa }}_j=j(i{\Omega }_0-\gamma ){\kappa }_j+h_{\psi }{\delta }_{1j}-h_{\psi }^{\ast }(j^2{\kappa }_{j+1}+j\sum _{m=1}^j{\kappa }_{j-m+1}{\kappa }_m)-{\sigma }^2(j^2{\kappa }_j+j\sum _{m=1}^{j-1}{\kappa }_{j-m}{\kappa }_m),$$

$${\dot{\varkappa }}_j=j(i{\Omega }_0-\gamma ){\varkappa }_j+h{\delta }_{j1}-h^{\ast }\left(j^2{\varkappa }_{j+1}+j{\sum }_{m=1}^j{\varkappa }_{j-m+1}{\varkappa }_m\right)$$

$$-{\sigma }^2\left(j^2{\varkappa }_j+j{\sum }_{m=1}^{j-1}{\varkappa }_{j-m}{\varkappa }_m\right)-G(\gamma +{\sigma }^2)(j-1)(j-2)2^j{\varkappa }_j,$$

$$a_j=q_1(b_1{)}^j+q_2(b_2{)}^j,$$

$$b_{1,2}=(R_0\pm \delta R)e^{i({\Phi }_0\pm \delta \varphi )}$$

$$\displaystyle\textstyle\dot{\varkappa}_{j}=j(i\Omega_{0}-\gamma)\varkappa_{j}+h\delta_{j1}-h^{\ast}\Big{(}j^{2}\varkappa_{j+1}e^{-\lambda(j-1)}\qquad\qquad\qquad$$

$$\displaystyle\textstyle{}+j\sum_{m=1}^{j}\varkappa_{j-m+1}\varkappa_{m}\Big{)}-{\sigma^{2}}\Big{(}j^{2}{\varkappa}_{j}+j\sum_{m=1}^{j-1}\varkappa_{j-m}\varkappa_{m}\Big{)}.$$

Full text

1 Introduction

In studies of collective phenomena, especially interesting are the cases, where a network of elements exhibits nontrivial dynamics, while the inherent dynamics of elements is simple and coupling between elements or/and an external driving of them are small. This is the case of populations of self-sustained periodic oscillators with a weak mutual coupling or/and subject to a weak common forcing. For description of such system a phase reduction proved [1] to be applicable and the paradigmatic Kuramoto model [2] was found to be relevant for many physical, chemical and biological systems [3, 4]. The examples include both continuous distributed parameter systems and networks of discrete elements.

Later on, Kuramoto system as well as a broad class of paradigmatic phase systems of form

[TABLE]

where $\Omega(t)$ and $h(t)$ are arbitrary real and complex-valued functions of time, were found to be partially integrable. Watanabe–Strogatz [5, 6, 7, 8] and Ott–Antonsen theories [12] were developed for such systems. While WS theory allowed studying the dynamics of finite populations [9, 10, 11] and made a solid foundation for OA theory, the latter one, constructed for the thermodynamic limit of large populations, and especially its extension to the case of nonidentical elements [13, 14, 15], found broad applications and helped to elucidate many sophisticated collective phenomena [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27].

Simultaneously, the peculiarity of the systems of type (1) necessitated construction of a perturbation theory for the OA approach, since real system can be close to (1) but never perfectly possess such a form. In [28], a formalism of “circular cumulants” formally corresponding to the Kuramoto–Daido order parameters $a_{n}=\langle(e^{i\varphi})^{n}\rangle=N^{-1}\sum_{k=1}^{N}e^{in\varphi_{k}}$ was suggested. Specifically, with a moment-generating function $F(\zeta)=\sum_{j=0}^{\infty}a_{j}\zeta^{j}/j!$, one can introduce circular cumulants $\varkappa_{j}$ via power series of the generating function $\zeta\frac{\partial}{\partial\zeta}\ln{F(\zeta)}=\sum_{j=1}^{\infty}\varkappa_{j}\zeta^{j}$. For instance, the first three cumulants are: $\varkappa_{1}=a_{1}$, $\varkappa_{2}=a_{2}-a_{1}^{2}$, $\varkappa_{3}=(a_{3}-3a_{2}a_{1}+2a_{1}^{3})/2$ . This formalism turned out to be fruitful for generalization of the Ott–Antonsen theory [28, 29, 30].

In particular, an important case of violation of the applicability conditions of OA theory is the one of population on phase elements

[TABLE]

with natural frequencies $\Omega_{k}$, which are distributed around the mean value $\Omega_{0}$ according to the Lorentzian distribution of half-width $\gamma$, and individual white gaussian noises of strength $\sigma$: $\langle\xi_{k}(t)\rangle=0$, $\langle\xi_{k}(t)\xi_{m}(t^{\prime})\rangle=2\delta_{km}\delta(t-t^{\prime})$. The collective dynamics of population (2) for large $N$ is governed by the following system of equations for circular cumulants [28]:

[TABLE]

where $\delta_{1j}=1$ for $j=1$ and $\delta_{1j}=0$ for $j\neq 1$. For $\sigma=0$, equation system (3) admits a particular solution $\varkappa_{1}=a_{1}$, $\varkappa_{j\geq 2}=0$; this case reproduces the OA theory.

While equations (3) are useful for analytical studies and for construction of perturbation theories [28, 29, 30], direct numerical simulations with them require truncations of the infinite system. Truncations with more than the first two cumulants kept were revealed to be able to cause a numerical instability. This issue turned out to be especially important for chimera-bearing systems [7, 16, 19], where neutrally stable or weakly stable degrees of freedom appear to be rather snesitive to the numeric instabilities. On the other hand, chimera regimes are in focus of many research works (e.g., [31, 32]) also due to similarities with turbulence intermittency in fluid dynamics [33, 34].

In this paper we suggest two approaches to suppression of the numerical instability of direct simulations with truncated chains (3) and test them for a chimera-bering system [7, 16] described in section 2. In section 3, the stabilization by means of a dissipation term for the high-order cumulants is studied. In section 4, we consider the approach based on the suppression of terms causing instability. In section 5, the findings are summarized.

2 Hierarchical population of coupled oscillators

Abrams et al. [16] studied the collective dynamics of an ensemble of hierarchically coupled oscillators, composed by two subpopulations of equal size $N$; the strength of a Kuramoto–Sakaguchi coupling between oscillators [2, 35, 36] within each subpopulation, $(1+A)/2$, is larger than the coupling strength $(1-A)/2$ between subpopulations:

[TABLE]

where $\alpha$ is the coupling phase shift, $\sigma$ is the intrinsic noise strength, $\xi_{k}$ and $\eta_{k}$ are independent white Gaussian noises. For vanishing noise, the system satisfies the conditions of Watanabe–Strogatz and Ott–Antonsen theories [5, 6, 12], which provided a substantial progress in understanding its dynamics [7, 16]. The chimera states, where one of subpopulations is fully synchronized $\psi_{1}=\psi_{2}=...=\Psi$, while the other one is in a partial synchrony state, were attracting the main interest of researchers [16].

However, the intrinsic noise violates these properties, and the circular cumulant approach proved to be useful for the description of the noisy system [28, 29]. In the thermodynamic limit of large $N$, one can introduce the circular cumulant description, with cumulants $\varkappa_{j}$ for subpopulation $\{\varphi_{k}\}$ and $\kappa_{j}$ for $\{\psi_{k}\}$ (cf. equation (3)):

[TABLE]

where $h_{\varphi}=0.25[(1+A)\varkappa_{1}+(1-A)\kappa_{1}]e^{-i\alpha}$ and $h_{\psi}=0.25[(1+A)\kappa_{1}+(1-A)\varkappa_{1}]e^{-i\alpha}$. In the classical “Abrams’ system” [16] the oscillator natural frequencies are identical: $\gamma=0$. In this paper we use the Abrams’ system (4) as a model for testing the approaches to stabilizing numerical simulation with truncated equation chains (5), (6).

3 Introducing auxiliary dissipation

From numerical viewpoint, eqs. (3) is a system of ODEs, which, if properly truncated, can be solved with the Runge-Kutta method. However, for some systems, truncation of the system can lead to a numerical instability. In this section we introduce additional dissipation-like terms into the original equations,

[TABLE]

and study their effect on the numerical stability of calculations. The choice of such a form of the stabilizing term is guided by the following reasons:

$\bullet$ The $G$-term introduces additional dissipation into the system dynamics. Similar to natural dissipation terms, it depends on parameters $\gamma$ and $\sigma$, but is rather strong for large $j$.

$\bullet$ Since the stabilizing term must make a stronger impact on the higher-order cumulants, i.e. grow with $j$ much faster than the $\sigma^{2}j^{2}\varkappa_{j}$-term, we adopt the $j$-dependence of the form $\propto j^{2}s^{j}\varkappa_{j}$. As simulations with only the first or two first cumulants are always stable, it is reasonable to amend the $j^{2}$-multiplier by changing it to $(j-1)(j-2)$. The numerical simulations revealed that the choice of $s=2$ is optimal.

We perform numerical simulations with the following initial conditions: the second subpopulation is synchronized $\kappa_{j}=\delta_{1j}$ and the elements of the first subpopulation form a two-bunch distribution, which is a superposition of two wrapped Cauchy distributions with complex order parameters $b_{1}$ and $b_{2}$. In terms of moments $a_{j}=\langle{e^{ij\varphi}}\rangle$, the latter means

[TABLE]

where $q_{1}$ and $q_{2}=1-q_{1}$ are the fractions of elements in the first and second bunch, respectively. Terms $(b_{1})^{j}$ and $(b_{2})^{j}$ represent the probability density distributions of $\varphi$ corresponding to OA solutions with complex order parameters $b_{1}$ and $b_{2}$, respectively. The distance between $b_{1}$ and $b_{2}$ is a convenient parameter for controlling the proximity to the OA manifold; the case $b_{1}=b_{2}$ yields an OA solution.

Let us first consider the case of equipartition between two bunches: $q_{1}=q_{2}=0.5$. In this case, simulations can be performed accurately, due to the fact that all odd cumulants, except the first one, turn to zero [30]. Indeed, one can perform simulations for $\varkappa_{1}$ and $\varkappa_{2n}$ with terms $\varkappa_{2n+1}=0$ wiped-out from equations for $n\geq 1$; the numerical calculations are stable for any order of truncation in this case. The trajectory of the system is shown in figure 2.

Figure 2 provides the results which can guide the choice of the values of coefficient $G$ for numerical simulations. In the graph, one can see the dependence of the critical value $G$, above which the numerical calculations are stable, on the number of cumulants kept in a truncated chain, $N_{\mathrm{trunc}}$. The choice of value of $G$, moderately exceeding the critical value, allows achieving compromise: the simulations are free of numerical instability, while the system trajectory is practically unaffected by the stabilizing term.

In figure 3(a) we demonstrate the explosion of numerical simulation due to numerical instability below the stabilization threshold; one can see that the blow-up is driven by the highest-order cumulant, the equation for evolution of which is affected by truncation of the cumulant chain. Above the threshold, in figure 3(b), the evolution is regular, fast fluctuations of cumulants decay with time, and the system forms a hierarchy of smallness for cumulants as predicted theoretically [28, 29].

Stabilization can be achieved with a minor distortion of the original system dynamics for not more than 20 cumulants and a non-large distance between the two bunches in the unsynchronized subpopulation. When using more than 20 cumulants, the value of $G$ required to stabilize the numerical calculations grows very quickly with the number of cumulants: the stabilizing term becomes significantly distorting the ‘true’ dynamics of the original system.

Let us consider the case of noise strength $\sigma=0.01$. In the unsynchronized subpopulation, the oscillators are equally distributed between two bunches, i.e. $q_{1}=q_{2}=0.5$. In figure 2 one can see that the critical value of $G$ for the noisy case (red circles) is lower than for the noise-free case (blue triangles). In figure 4(a) the example of system trajectory is plotted. One can notice that the trajectory for a marginal value of $G$, although being significantly disturbed by fluctuations, does not diverge with the trajectory calculated for a larger value ​​of $G$, for which the computational fluctuations are suppressed.

Testing with the equipartition case $q_{1}=q_{2}$ at $\sigma=0$ is important, since in this case the integration of the modified equation system can be accurately performed without a numerical instability; one can check the results of numerical integration with stabilization terms against the numerical solutions for the unaltered system. However, the specificity of the equipartition solutions can also effect their numerical stability properties. Hence, an additional treatment for the case of unequal partition between two bunches is required. The numerical results for $\sigma=0$, $q_{1}=0.4$ and $q_{2}=0.6$ are presented in figures 2 (green squares) and 4(b). In figure 2 one can see that the characteristic critical values ​​of $G$ for equal and unequal partitions between two bunches are close. For a noisy system, smaller values ​​of the coefficient $G$ are required, since noise makes a stabilizing effect.

From the provided results on the stability boundary of the numerical simulations, one can conclude that, for the considered initial conditions, the stabilizing term is not required if less than 8 cumulants are used. For a simulation with larger number of cumulants, a stabilizing term is required.

The efficiency of stabilization procedure depends on initial conditions. For initial conditions, the distance between the complex order parameters

[TABLE]

in amplitude and phase ($\delta R$, $\delta\phi$) of two bunches is required to be non-large, otherwise the suppression of the numerical instability becomes impossible (simulation with two first cumulants is always stable). For these situations, the required value of the coefficient $G$ becomes too large and affects the natural dynamics of the system. In figure 5, the dependence of the critical value of $G$ on the phase difference $\delta\phi$ is plotted for a fixed distance in amplitude between two bunches.

4 Suppression of the terms responsible for instability

An alternative approach to stabilization is the suppression of the terms which are primarily responsible for numerical instability in the original system. We modify one of the terms of the original equation by an additional exponential factor $e^{-\lambda(j-1)}$:

[TABLE]

This modification is guided by the following reason:

$\bullet$  The numerical instability is related to the $\varkappa_{j+1}$-term in $\dot{\varkappa}_{j}$. We modify it in such a way that for large $j$ it vanishes, while for small $j$ it is nearly unchanged, that is the modification does not affect the accuracy of calculations significantly.

We apply this stabilization procedure to the same examples as in the previous section. The second subpopulation is nearly perfectly synchronized, in the first subpopulation oscillators are partitioned between two bunches. The cases to be considered:

1) Oscillators are equally distributed between two bunches, i.e. $q_{1}=q_{2}=0.5$; no noise, $\sigma=0$.

2. Oscillators are equally distributed between two bunches, i.e. $q_{1}=q_{2}=0.5$, but noise is added to the system $\sigma=0.01$.

3. Oscillators are unequally distributed between two bunches, $q_{1}=0.4$ and $q_{2}=0.6$; $\sigma=0$.

The minimal exponential coefficient $\lambda$ required for stabilization as a function of the number of cumulants $N_{\mathrm{trunc}}$ is presented in figure 7. Initial conditions are the same as in the previous section (see the caption for figure 2).

We also calculate the dependence of the coefficient $\lambda$ on the phase difference $\delta\phi$ between two bunches. In figure 7, one can see such a dependence for the case of unequal partition of oscillators between two bunches, $q_{1}=0.4$ and $q_{2}=0.6$.

Similarly to the previous stabilization method, this method can be implemented without affecting the natural system dynamics for direct numerical simulation with no more than $20$ cumulants. The strength of the required stabilizing correction depends on the proximity to the Ott–Antonsen manifold; too far away from the OA manifold (for a large distance between two bunches), numerical simulations with large number of cumulants become impossible.

5 Conclusion

We have considered two approaches to suppress numerical instability of direct simulations of truncated systems of circular cumulant equations (3). With the first approach, an additional dissipation term for high-order cumulants is introduced (see equation 7). With the second approach, the the coupling term term $\propto\varkappa_{j+1}$, responsible for instability of truncated chains, is exponentially suppressed for high-order cumulants (see equation 8). The stability domains have been constructed for these stabilization methods; the dependence of stability on the initial conditions has been studied as well.

One can conclude, that the stabilization can be practically achieved without significant effect on the natural system dynamics for no more than 20 cumulants. Stabilization efficiently works in a small but finite vicinity of the Ott–Antonsen manifold. For a larger deviation from the OA manifold, direct numerical simulations with fewer number of cumulants can be stabilized.

\ack

The work has been financially supported by the Russian Science Foundation (Grant no. 19-42-04120).