The Method of Multiple Scales

This page details the dynamical systems on which the Method of Multiple Scales (MMS) is applied, how it is applied, how steady state responses can be derived, and how their stability can be evaluated.

Detailed descriptions of the MMS and its variants can be found in [1, 2, 3, 4]. Only the version used in the oscilate package is described below.

Dynamical systems of interest

Coupled oscillators system

Systems considered are typically composed of \(N\) coupled nonlinear equations of the form

\[\begin{split}\begin{cases} \ddot{x}_0 + \omega_0^2 x_0 & = f_0(\boldsymbol{x}, \dot{\boldsymbol{x}}, \ddot{\boldsymbol{x}}, t), \\ & \vdots \\ \ddot{x}_{N-1} + \omega_{N-1}^2 x_{N-1} & = f_{N-1}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \ddot{\boldsymbol{x}}, t). \end{cases}\end{split}\]

The \(x_i(t)\) (\(i=0,...,N-1\)) are the oscillators’ coordinates, and represent the degrees of freedom (dof) of the system.

\[\boldsymbol{x}(t)^\intercal = [x_0(t), x_1(t), \cdots, x_{N-1}(t)]\]

is the vector containing all the oscillators’ coordinates (\(^\intercal\) denotes the transpose), \(\omega_i\) are their natural frequencies, \(t\) is the time, \(\dot{(\bullet)} = \textrm{d}(\bullet)/\textrm{d}t\) denotes a time-derivative. \(f_i\) are functions which can contain:

Weak linear terms in \(x_i\), \(\dot{x}_i\), or \(\ddot{x}_i\).
Weak linear coupling terms involving \(x_j\), \(\dot{x}_j\), or \(\ddot{x}_j\) with \(j \neq i\).
Weak nonlinear terms. Taylor expansions are performed to approximate nonlinear terms as polynomial nonlinearities.
Forcing terms, which can be:
- Hard (appearing at leading order) or weak (small).
- Primarily harmonic, e.g., \(F \cos(\omega t)\), where \(F\) and \(\omega\) are the forcing amplitude and frequency, respectively.
- Modulated by any function (constant, linear, or nonlinear) to model parametric forcing (e.g., \(x_i(t) F \cos(\omega t)\)).

Internal resonance relations among oscillators can be specified in a second step by expressing the \(\omega_i\) as a function of a reference frequency. See Dynamical_system. Detuning can also be introduced during this step.

Nonlinear system. — Illustration of a system of coupled nonlinear oscillators that can be studied using `oscilate`. \(f_{\text{nl},i}(x_i, \dot{x}_i, \ddot{x}_i)\), \(g_{\text{nl},i}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \ddot{\boldsymbol{x}})\) and \(f_{\text{f},_i}(x_i, \dot{x}_i, \ddot{x}_i)\) represent the intrinsic nonlinear of oscillator \(i\), the nonlinear couplings with other oscillators, and the parametric forcing coeffiicent. A unitary mass is considered for each oscillator, and \(\omega_i^2\) and \(c_i\) are the stiffness (natural frequencies) and the (weak) viscous damping coefficients. Note that this illustration does not represent weak linear terms nor complex coupled parametric forcing, which can be accounted for in the `oscilate` package.

Complex form

It will later be shown that it may be convenient to rewrite the coupled \(2^{\text{nd}}\) order equations as coupled \(1^{\text{st}}\) order ones through the introduction of complex variables \(z_i(t)\) (see section Complex form of the MMS). These are defined through the transformation

\[\begin{split}\begin{cases} x_i & = z_i + \bar{z}_i, \\ \dot{x}_i & = \textrm{j} \omega_{i} (z_i - \bar{z}_i), \end{cases}\end{split}\]

where \(\bar{(\bullet)}\) denotes the complex conjugate, or equivalently,

\[\begin{split}\begin{cases} z_{i} & = \frac{1}{2} (x_{i} - \frac{\textrm{j}}{\omega_i} \dot{x}_{i}), \\ \bar{z}_{i} & = \frac{1}{2} (x_{i} + \frac{\textrm{j}}{\omega_i} \dot{x}_{i}). \end{cases}\end{split}\]

The \(z_i(t)\), \(i=0,...,N-1\), complex variables therefore encapsulate the state space of the system. Now, noting that

\[\ddot{x}_{i} = 2 \textrm{j} \omega_{i} \dot{z}_i + \omega_{i}^2 (z_{i} - \bar{z}_{i}),\]

one obtains the coupled, complex first order equations on the complex variables,

\[\begin{split}\begin{cases} \dot{z}_0 & = \textrm{j} \omega_0 z_0 + g_0(\boldsymbol{z}, \bar{\boldsymbol{z}}, \dot{\boldsymbol{z}}, t), \\ & \vdots \\ \dot{z}_{N-1} & = \textrm{j} \omega_{N-1} z_{N-1} + g_{N-1}(\boldsymbol{z}, \bar{\boldsymbol{z}}, \dot{\boldsymbol{z}}, t), \end{cases}\end{split}\]

where

\[\boldsymbol{z}(t)^\intercal = [z_0(t), z_1(t), \cdots, z_{N-1}(t)]\]

is the vector of complex coordinates, and functions \(g_i\) are related to \(f_i\) such that

\[g_i(\boldsymbol{z}, \bar{\boldsymbol{z}}, \dot{\boldsymbol{z}}, t) = - \frac{\textrm{j}}{2\omega_i} f_i(\boldsymbol{x}, \dot{\boldsymbol{x}}, \ddot{\boldsymbol{x}}, t).\]

The dynamical system is transformed from oscillator to complex form using complex_form().

Application of the MMS

Asymptotic series and time scales

The starting point is to introduce asymptotic series and multiple time scales in the initial dynamical system. The solution for oscillator \(i\) is sought as a series expansion up to order \(N_e\) (for a leading order term \(\epsilon^0 = 1\)). This expansion takes the form

\[x_i(t) = x_{i,0}(t) + \epsilon x_{i,1}(t) + \epsilon^2 x_{i,2}(t) + \cdots + \epsilon^{N_e} x_{i,N_e}(t) + \mathcal{O}(\epsilon^{N_e+1}).\]

Time scales are introduced as follows:

\[t_0 = t, \; t_1 = \epsilon t, \; t_2 = \epsilon^2 t, \cdots, t_{N_e} = \epsilon^{N_e} t,\]

where \(t_0\) is the fast time, i.e. the time used to describe the oscillations, while \(t_1, \; t_2,\; \cdots,\; t_{N_e}\) are slow times, associated to amplitude and phase variations of the solutions in time. In addition, the chain rule gives

\[\begin{split}\begin{aligned} \dfrac{\textrm{d}(\bullet)}{\textrm{d}t} & = \sum_{i=0}^{N_e} \epsilon^{i} \dfrac{\partial(\bullet)}{\partial t_i} + \mathcal{O}(\epsilon^{N_e+1}), \\ \dfrac{\textrm{d}^2(\bullet)}{\textrm{d}t^2} & = \sum_{j=0}^{N_e}\sum_{i=0}^{N_e} \epsilon^{i+j} \dfrac{\partial}{\partial t_j}\dfrac{\partial(\bullet)}{\partial t_i} + \mathcal{O}(\epsilon^{N_e+1}). \end{aligned}\end{split}\]

The introduction of asymptotic series and time scales are performed using asymptotic_series() and time_scales().

Scaling

The construction of the MMS system requires a scaling of the parameters. Most scalings are already passed to the MMS through the sub_scaling parameter. However, the natural frequencies also need to be scaled as they can contain both a leading order term and a detuning term. Natural frequencies \(\omega_i\) are defined as a function of the reference frequency \(\omega_{\textrm{ref}}\) through the ratio_omega_osc optional parameter, which is then used in oscillators_frequencies(). This allows to define internal resonance relations among the oscillators. If these internal resonances are not perfect, detunings can be introduced through the detunings optional parameter, which needs to be scaled and part of the sub_scaling parameter.

To write the MMS system it is convenient to introduce the leading order natural frequencies

\[\omega_{i,0} = r_i \omega_{\textrm{ref}},\]

where \(r_i\) stands for ratio_omega_osc[i].

The multiple scales system

Introducing the asymptotic series, the time scales and the scaled parameters in the initial dynamical system (see Dynamical_system) results in \(N_e+1\) dynamical systems, each one appearing at different orders of \(\epsilon\). Denoting time scales derivatives as

\[\textrm{D}_i(\bullet) = \partial (\bullet) / \partial t_i,\]

introducing the vector of time scales

\[\boldsymbol{t}^\intercal = [t_0, t_1, \cdots, t_{N_e}],\]

where \(^\intercal\) denotes the transpose, and the vectors of asymptotic coordinates

\[\boldsymbol{x}_i(\boldsymbol{t})^\intercal = [x_{0,i}(\boldsymbol{t}), x_{1,i}(\boldsymbol{t}), \cdots, x_{N-1, i}(\boldsymbol{t})],\]

which contains all the asymptotic terms of order \(i\), the MMS equations can be written as

\[\begin{split}\begin{aligned} & \epsilon^0 \rightarrow \; \begin{cases} \textrm{D}_0^2 x_{0,0} + \omega_{0,0}^2 x_{0,0} & = f_{0}^{(0)}(t_0, t_1), \\ & \vdots \\ \textrm{D}_0^2 x_{N-1,0} + \omega_{N-1,0}^2 x_{N-1,0} & = f_{N-1}^{(0)}(t_0, t_1), \end{cases} \\[15pt] & \epsilon^1 \rightarrow \; \begin{cases} \textrm{D}_0^2 x_{0,1} + \omega_{0,0}^2 x_{0,1} & = f_{0}^{(1)} (\boldsymbol{x}_0, \textrm{D}_0 \boldsymbol{x}_0, \textrm{D}_0^2 \boldsymbol{x}_0, \textrm{D}_1 \boldsymbol{x}_0, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_0, t_0, t_1), \\ & \vdots \\ \textrm{D}_0^2 x_{N-1,1} + \omega_{N-1,0}^2 x_{N-1,1} & = f_{N-1}^{(1)}(\boldsymbol{x}_0, \textrm{D}_0 \boldsymbol{x}_0, \textrm{D}_0^2 \boldsymbol{x}_0, \textrm{D}_1 \boldsymbol{x}_0, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_0, t_0, t_1), \end{cases} \\[15pt] & \epsilon^2 \rightarrow \; \begin{cases} \textrm{D}_0^2 x_{0,2} + \omega_{0,0}^2 x_{0,2} & = f_{0}^{(2)} (\boldsymbol{x}_0, \cdots, \textrm{D}_0 \textrm{D}_2 \boldsymbol{x}_0, \boldsymbol{x}_1, \cdots, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_1, t_0, t_1), \\ & \vdots \\ \textrm{D}_0^2 x_{N-1,2} + \omega_{N-1,0}^2 x_{N-1,2} & = f_{N-1}^{(2)}(\boldsymbol{x}_0, \cdots, \textrm{D}_0 \textrm{D}_2 \boldsymbol{x}_0, \boldsymbol{x}_1, \cdots, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_1, t_0, t_1), \end{cases} \\[10pt] & \hspace{3cm} \vdots \\[10pt] & \epsilon^{N_e} \rightarrow \; \begin{cases} \textrm{D}_0^2 x_{0,N_e} + \omega_{0,0}^2 x_{0,N_e} & = f_{0}^{(N_e)} (\boldsymbol{x}_0, \cdots, \textrm{D}_0 \textrm{D}_{N_e} \boldsymbol{x}_0, \cdots, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_{N_e-1}, t_0, t_1), \\ & \vdots \\ \textrm{D}_0^2 x_{N-1,N_e} + \omega_{N-1,0}^2 x_{N-1,N_e} & = f_{N-1}^{(N_e)}(\boldsymbol{x}_0, \cdots, \textrm{D}_0 \textrm{D}_{N_e} \boldsymbol{x}_0, \cdots, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_{N_e-1}, t_0, t_1). \end{cases} \end{aligned}\end{split}\]

Consider oscillator \(i\) at order \(j\):

The left-hand side term represents a harmonic oscillator of frequency \(\omega_{i,0}\) oscillating with respect to the fast time \(t_0\).
The right-hand side term \(f_{i}^{(j)}\) is analogous to a forcing generated by all combinations of terms that appear on oscillator \(i\)’s equation at order \(\epsilon^j\). This can involve lower order terms \(x_{i,\ell}, \; \ell \leq j\), coupling terms \(x_{k, \ell}, \; k \neq j,\; \ell \leq j\), their derivatives and cross-derivatives with respect to the time scales, and physical forcing terms. The later is responsible for the dependency on \(t_0,\; t_1\). The reason why slower time scales are not involved will be explained in the following.

Function \(f_{i}^{(j)}\) tends to get increasingly complex as the order increases because the initial equations generate more high order terms than low order ones.

This operation is performed using compute_EqO().

Note that internal resonance relations can be given through the ratio_omega_osc optional parameter, which is then used in oscillators_frequencies().

Frequency of interest

The response of \(x_i\) will be analysed at a frequency \(\omega\), defined as

\[\omega = \omega_{\textrm{MMS}} + \epsilon \sigma,\]

where \(\omega_{\textrm{MMS}}\) is the central MMS frequency, controlled through the ratio_omegaMMS optional parameter and expressed in terms of omega_ref, and \(\sigma\) is a detuning about that frequency. In case the forced response is studied, \(\omega\) corresponds to the forcing frequency. In case the free response is studied, \(\omega\) corresponds to the frequency of free oscillations, which generates the backbone curve of the forced response. Note that \(\omega t = \omega_{\textrm{MMS}} t_0 + \sigma t_1\). This is the reason why the forcing only involves these two time scales in the right-hand side functions of the MMS system.

Iteratively solving the MMS system

The multiple scales system can be solved iteratively by solving successively the systems of equations at each order.

Leading order solution

The leading order solution for oscillator \(i\) must satisfy

\[\textrm{D}_0^2 x_{i,0} + \omega_{i,0}^2 x_{i,0} = f_{i}^{(0)}(t_0, t_1).\]

It is sought as

\[x_{i,0}(\boldsymbol{t}) = x_{i,0}^\textrm{h}(\boldsymbol{t}) + x_{i,0}^\textrm{p}(t_0, t_1),\]

where \(x_{i,0}^\textrm{h}(\boldsymbol{t})\) and \(x_{i,0}^\textrm{p}(t_0, t_1)\) are the leading order homogeneous and particular sollutions, respectively.

It is now conveninent to introduce the slow times vector

\[\boldsymbol{t}_s^\intercal = [t_1, \cdots, t_{N_e}].\]

This way, one can express the leading order solutions as

\[\begin{split}\begin{cases} x_{i,0}^\textrm{h}(\boldsymbol{t}) & = A_i(\boldsymbol{t}_s) e^{\textrm{j} \omega_{i,0} t_0} + cc = |A_i(\boldsymbol{t}_s)| \cos(\omega_{i,0} t_0 + \arg{A_i(\boldsymbol{t}_s)}), \\ x_{i,0}^\textrm{p}(t_0, t_1) & = B_i e^{\textrm{j} \omega t} + cc = B_i e^{\textrm{j} (\omega_{\textrm{MMS}} t_0 + \sigma t_1)} + cc = |B_i| \cos(\omega_{\textrm{MMS}} t_0 + \sigma t_1 + \arg{B_i}), \end{cases}\end{split}\]

where \(A_i\) is a slow time-dependent complex amplitude to be determined while \(B_i\) is a time-independent function of the forcing parameters. \(cc\) denotes the complex conjugate. Note that in most situations, forcing does not appear at leading order (i.e. forcing is weak), so \(B_i=0\).

In the following it will be convenient to use the notations

\[\begin{split}\begin{split} \boldsymbol{A}(\boldsymbol{t}_s)^\intercal & = [A_0(\boldsymbol{t}_s), A_1(\boldsymbol{t}_s), \cdots, A_{N-1}(\boldsymbol{t}_s)], \\ \boldsymbol{B}^\intercal & = [B_0, B_1, \cdots, B_{N-1}]. \end{split}\end{split}\]

The leading order solutions are defined in sol_order_0().

Higher order solutions

Once the leading order solutions are computed, they can be injected in the \(1^\textrm{st}\) higher order equations, where they (and their derivatives) appear as forcing terms, potentially together with physical forcing. The \(1^\textrm{st}\) higher order equation for oscillator \(i\) is

\[\textrm{D}_0 x_{i,1} + \omega_{i,0}^2 x_{i,1} = f_{i}^{(1)}(\boldsymbol{x}_0, \textrm{D}_0 \boldsymbol{x}_0, \textrm{D}_0^2 \boldsymbol{x}_0, \textrm{D}_1 \boldsymbol{x}_0, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_0, t_0, t_1).\]

The forcing terms that involve oscillations at \(\omega_{i,0}\) would force the oscillator on its natural frequency. Moreover, damping is always weak in the MMS, so damping terms of the form \(c \textrm{D}_0 x_{i,1}\) do not appear at this order. The aforementioned forcing terms would thus lead to unbounded solutions, which is unphysical. These forcing terms, called secular terms, must therefore be eliminated. For instance, the \(1^\textrm{st}\) higher order equation for oscillator \(i\) with the secular terms cancelled is

\[\textrm{D}_0 x_{i,1} + \omega_{i,0}^2 x_{i,1} = \breve{f}_{i}^{(1)}(\boldsymbol{x}_0, \textrm{D}_0 \boldsymbol{x}_0, \textrm{D}_0^2 \boldsymbol{x}_0, \textrm{D}_1 \boldsymbol{x}_0, \textrm{D}_0\textrm{D}_1 \boldsymbol{x}_0, t_0, t_1),\]

where \(\breve{f}_{i}^{(1)}\) is \(f_{i}^{(1)}\) with the secular terms cancelled, i.e. without terms oscillating as \(\omega_{i,0}\). After cancelation of the secular terms, each oscillator’s equation can be solved as a forced harmonic oscillator with the independent variable \(t_0\).

Note that only the particular solutions are considered when solving higher order terms, i.e.

\[x_{i,1}(\boldsymbol{t}) = \underbrace{x_{i,1}^\textrm{h}(\boldsymbol{t})}_{=0} + x_{i,1}^\textrm{p}(t_0, t_1).\]

This choice can be justified if one assumes that initial conditions are of leading order. Though this is questionable, it is assumed here. A more advanced treatment of the initial conditions of higher order solutions is detailed in [5] and section Complex form of the MMS.

The higher order solutions \(x_{i,1}(\boldsymbol{t})\) are expressed as a function of the leading order unknown amplitudes \(\boldsymbol{A}(\boldsymbol{t}_s)\), their slow time derivatives \(\textrm{D}_i\boldsymbol{A}(\boldsymbol{t}_s), \; i=1, ..., N_e\), and forcing terms if any (including the hard forcing amplitudes \(\boldsymbol{B}\)).

This process is repeated successively at each order, i.e. the computed solutions are introduced in the next higher order system of equations, secular terms are cancelled and the next higher order solutions are computed.

The secular terms are identified in secular_analysis() and the leading order solutions are computed in sol_higher_order(). Note that sol_higher_order() is applied on equations with only \(t_0\) as the independent variable so as to allow the use of dsolve(). This is enforced using system_t0(), which temporarily ignores the dependency of \(\boldsymbol{A}(\boldsymbol{t}_s)\) on the slow time scales.

Secular analysis

At this stage, the solutions are all expressed in terms of the unknown amplitudes \(\boldsymbol{A}(\boldsymbol{t}_s)\) and their slow time derivatives \(\textrm{D}_1 A_i(\boldsymbol{t}_s),\; \cdots,\; \textrm{D}_{N_e} A_i(\boldsymbol{t}_s)\). These can be obtained from the elimination of the secular terms (called secular conditions), as described below.

The \(i^\textrm{th}\) MMS equation at \(1^\textrm{st}\) higher order involves the slow time derivative \(\textrm{D}_1 A_i(\boldsymbol{t}_s)\), which appears in the secular term. It is coming from the chain rule

\[\begin{split}\begin{cases} \dfrac{\textrm{d} x_i(t)}{\textrm{d}t} & = \dfrac{\partial x_{i,0}(\boldsymbol{t})}{\partial t_0} + \epsilon \dfrac{\partial x_{i,1}(\boldsymbol{t})}{\partial t_0} + \epsilon \dfrac{\partial x_{i,0}(\boldsymbol{t})}{\partial t_1} + \mathcal{O}(\epsilon^2), \\ \dfrac{\textrm{d}^2 x_i(t)}{\textrm{d}t^2} & = \dfrac{\partial^2 x_{i,0}(\boldsymbol{t})}{\partial t_0^2} + \epsilon \dfrac{\partial^2 x_{i,1}(\boldsymbol{t})}{\partial t_0^2} + 2 \epsilon \dfrac{\partial^2 x_{i,0}(\boldsymbol{t})}{\partial t_0 \partial t_1} + \mathcal{O}(\epsilon^2). \end{cases}\end{split}\]

In addition, the \(\textrm{D}_1 A_j(\boldsymbol{t}_s),\; j\neq i\) do not appear in the \(i^\textrm{th}\) MMS equation as couplings among oscillators are weak. It is thus possible to use the secular conditions in order to express the \(\textrm{D}_1 A_i(\boldsymbol{t}_s)\) as a function of \(\boldsymbol{A}(\boldsymbol{t}_s)\).

This process can be done successively at each order to obtain the system of complex modulation (or slow time) equations

\[\begin{split}\begin{cases} \textrm{D}_1 A_i(\boldsymbol{t}_s) & = f_{A_i}^{(1)}(\boldsymbol{A}, t_1), \\ & \vdots \\ \textrm{D}_{N_e} A_i(\boldsymbol{t}_s) & = f_{A_i}^{(N_e)}(\boldsymbol{A}, t_1). \end{cases}\end{split}\]

\(f_{A_i}^{(j)}(\boldsymbol{A}, t_1)\) are functions governing the modulation of \(A_i\) with respect to the slow time \(t_j\). Note the dependency of \(f_{A_i}^{(j)}\) on \(t_1\) due to the possible presence of forcing. The complex modulation equations are derived in secular_analysis().

The above system of \(1^\textrm{st}\) order PDE can theoretically be solved to obtain the complex amplitudes \(\boldsymbol{A}\). However, this approach is not the prefered one as

It is more convenient to deal with real variables than complex ones to get a physical meaning from the analysis,
It is more convenient to deal with autonomous systems (without the explicit \(t_1\)-dependency) than nonautonomous ones,
The PDEs are complex.

The first two points can be achieved introducing new coordinates, as described thereafter.

Modulation equations

Polar coordinates

As discussed previously, it is more convenient to deal with real variables than complex ones. This can be done introducing the polar coordinates \(a_i\) and \(\phi_i\) for oscillator \(i\) such that

\[A_i(\boldsymbol{t}_s) = \dfrac{1}{2} a_i(\boldsymbol{t}_s) e^{\textrm{j} \phi_i(\boldsymbol{t}_s)}.\]

\(a_i\) and \(\phi_i\) correspond to the amplitude and phase of the leading solution for oscillator \(i\), respectively. With these new coordinates and introducing

\[\begin{split}\begin{aligned} \boldsymbol{a}(\boldsymbol{t}_s)^\intercal & = [a_0(\boldsymbol{t}_s), a_1(\boldsymbol{t}_s), \cdots, a_{N-1}(\boldsymbol{t}_s)], \\ \boldsymbol{\phi}(\boldsymbol{t}_s)^\intercal & = [\phi_0(\boldsymbol{t}_s), \phi_1(\boldsymbol{t}_s), \cdots, \phi_{N-1}(\boldsymbol{t}_s)], \end{aligned}\end{split}\]

the modulation equations on the \(A_i(\boldsymbol{t}_s)\) can be split into real and imaginary terms, leading to

\[\begin{split}\begin{aligned} & \epsilon^1 \rightarrow \; \begin{cases} \textrm{D}_1 a_i & = \hat{f}_{a_i}^{(1)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1), \\ a_i \textrm{D}_1 \phi_i & = \hat{f}_{\phi_i}^{(1)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1), \end{cases} \\[5pt] & \hspace{3cm} \vdots \\[5pt] & \epsilon^{N_e} \rightarrow \; \begin{cases} \textrm{D}_{N_e} a_i & = \hat{f}_{a_i}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1), \\ a_i \textrm{D}_{N_e} \phi_i & = \hat{f}_{\phi_i}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1). \end{cases} \end{aligned}\end{split}\]

The \(\epsilon^j \rightarrow\) indicate a system of 2 equations originating from the secular analysis at order \(j\). As one has

\[\textrm{D}_j A_i = \textrm{D}_j \dfrac{1}{2} a_i e^{\textrm{j} \phi_i} = \dfrac{1}{2} \textrm{D}_j \left( a_i \right) e^{\textrm{j} \phi_i} + \textrm{j} \dfrac{1}{2} a_i e^{\textrm{j} \phi_i} \textrm{D}_j \left( \phi_i \right),\]

it is convenient to pre-multiply the modulation equations on the \(A_i(\boldsymbol{t}_s)\) by \(e^{-\textrm{j} \phi_i(\boldsymbol{t}_s)}\) or even \(\gamma e^{-\textrm{j} \phi_i(\boldsymbol{t}_s)}\) with, for instance, \(\gamma = 2\). This avoids the presence of \(\cos(\phi_i)\), \(\sin(\phi_i)\) and many \(1/2\) terms in the modulation equations. The modulation functions on polar coordinates are therefore defined as

\[\begin{split}\begin{cases} \hat{f}_{a_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1) & = \Re\left[ 2 e^{-\textrm{j} \phi_i(\boldsymbol{t}_s)} f_{A_i}^{(j)}(\boldsymbol{A}, t_1) \right], \\ \hat{f}_{\phi_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\phi}, t_1) & = \Im\left[ 2 e^{-\textrm{j} \phi_i(\boldsymbol{t}_s)} f_{A_i}^{(j)}(\boldsymbol{A}, t_1) \right]. \end{cases}\end{split}\]

The modulation equations system on the polar coordinates involves only real variables, but functions \(\hat{f}_{a_i}^{(j)}\) and \(\hat{f}_{\phi_i}^{(j)}\) are still nonautonomous due to the explicit dependency on \(t_1\).

The polar coordinates are introduced in polar_coordinates(). The real modulation equations are only computed for the autonomous system, as described below.

Autonomous phase coordinates

The presence of nonautonomous terms stems from forcing terms, which involve \(\cos(\sigma t_1 - \phi_i), \; \sin(\sigma t_1 - \phi_i)\) in the polar modulation functions of oscillator \(i\). A change of phase coordinates is required to make this autonomous. Moreover, the change of phase coordinate is necessary even in the absence of forcing for a convenient representation of the leading order solution. Indeed, the solution for \(x^{\textrm{h}}_{i,0}(\boldsymbol{t})\) written in terms of the current polar coordinates is

\[x^{\textrm{h}}_{i,0}(\boldsymbol{t}) = a_i(\boldsymbol{t}_s) \cos(\omega_{i,0} t_0 + \phi_i(\boldsymbol{t}_s)).\]

However, one would eventually like to express the oscillations of oscillator \(i\) in terms of the frequency \(\omega\). To force its appearance, we recall that

\[\omega_{i,0} = r_i \omega_{\textrm{ref}}, \quad \omega_{\textrm{ref}} = \frac{1}{r_{\textrm{MMS}}} \omega_{\textrm{MMS}}, \quad \textrm{and} \quad \omega_{\textrm{MMS}} = \omega - \epsilon \sigma.\]

Introducing this in the leading order solution leads to

\[x^{\textrm{h}}_{i,0}(\boldsymbol{t}) = a_i(\boldsymbol{t}_s) \cos\left( \frac{r_i}{r_{\textrm{MMS}}} \omega t_0 - \frac{r_i}{r_{\textrm{MMS}}} \sigma t_1 + \phi_i(\boldsymbol{t}_s)\right).\]

It therefore appears convenient to introduce the new phase coordinate \(\beta_i(\boldsymbol{t}_s)\) as

\[\beta_i(\boldsymbol{t}_s) = \frac{r_i}{r_{\textrm{MMS}}} \sigma t_1 - \phi_i(\boldsymbol{t}_s),\]

which allows to write the leading order solution as

\[x^{\textrm{h}}_{i,0}(\boldsymbol{t}) = a_i(\boldsymbol{t}_s) \cos\left( \frac{r_i}{r_{\textrm{MMS}}} \omega t_0 - \beta_i(\boldsymbol{t}_s)\right).\]

In addition, and as discussed previously, the introduction of these new phase coordinates removes the explicit dependency of the modulation functions on \(t_1\), which was due to terms \(\cos(\sigma t_1 - \phi_i), \; \sin(\sigma t_1 - \phi_i)\). The forcing phase being zero (i.e. reference phase), the \(\beta_i(\boldsymbol{t}_s)\) can be seen as the phases relative to the forcing.

Introducing the notation

\[\boldsymbol{\beta}(\boldsymbol{t}_s)^\intercal = [\beta_1(\boldsymbol{t}_s), \beta_2(\boldsymbol{t}_s), \dots, \beta_{N_e}(\boldsymbol{t}_s)],\]

the modulation equations can be rewritten as

\[\begin{split}\begin{aligned} & \epsilon^1 \rightarrow \; \begin{cases} \textrm{D}_1 a_0(\boldsymbol{t}_s) & = f_{a_0}^{(1)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_0 \textrm{D}_1 \beta_0(\boldsymbol{t}_s) & = f_{\beta_0}^{(1)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ & \vdots \\ \textrm{D}_1 a_{N-1}(\boldsymbol{t}_s) & = f_{a_{N-1}}^{(1)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_{N-1} \textrm{D}_1 \beta_{N-1}(\boldsymbol{t}_s) & = f_{\beta_{N-1}}^{(1)}(\boldsymbol{a}, \boldsymbol{\beta}), \end{cases} \\[5pt] & \hspace{3cm} \vdots \\[5pt] & \epsilon^{N_e} \rightarrow \; \begin{cases} \textrm{D}_{N_e} a_0(\boldsymbol{t}_s) & = f_{a_0}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_0 \textrm{D}_{N_e} \beta_0(\boldsymbol{t}_s) & = f_{\beta_0}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ & \vdots \\ \textrm{D}_{N_e} a_{N-1}(\boldsymbol{t}_s) & = f_{a_{N-1}}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_{N-1} \textrm{D}_{N_e} \beta_{N-1}(\boldsymbol{t}_s) & = f_{\beta_{N-1}}^{(N_e)}(\boldsymbol{a}, \boldsymbol{\beta}). \end{cases} \end{aligned}\end{split}\]

The above system is the key result of the application of the MMS as it governs the modulation of leading order amplitudes and phases, on which all higher order solutions depend. It can be solved numerically or analytically, if analytical solutions exist, though it is generally not the case. An analytical resolution is shown in Van der Pol oscillator. The modulation system can also be rewritten in a more compact form as discussed in the following.

The autonomous phase coordinates are introduced in autonomous_phases() and the modulation equations are computed in modulation_equations().

All solutions previously computed using the complex amplitudes \(\boldsymbol{A}(\boldsymbol{t}_s)\) can be rewritten in terms of the polar coordinates \(\boldsymbol{a}(\boldsymbol{t}_s),\; \boldsymbol{\beta}(\boldsymbol{t}_s)\) using sol_x_polar().

Reconstituted equations

The reconstitution method consists in combining the modulation equations at each slow time scales in order to obtain modulation equations in the physical time. The approach used in the oscilate package corresponds to the Complete Inconsistent Method discussed in [3], in opposition to consistent approaches which treat each order individually.

Recalling the chain rule

\[\dfrac{\textrm{d}(\bullet)}{\textrm{d}t} = \sum_{i=0}^{N_e} \epsilon^{i} \textrm{D}_i(\bullet),\]

and reintroducing the physical time, the systems at each order can be summed up to write

\[\begin{split}\begin{cases} \dfrac{\textrm{d}}{dt} a_0(t) & = f_{a_0}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_0 \dfrac{\textrm{d}}{dt} \beta_0(t) & = f_{\beta_0}(\boldsymbol{a}, \boldsymbol{\beta}), \\ & \vdots \\ \dfrac{\textrm{d}}{dt} a_{N-1}(t) & = f_{a_{N-1}}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_{N-1} \dfrac{\textrm{d}}{dt} \beta_{N-1}(t) & = f_{\beta_{N-1}}(\boldsymbol{a}, \boldsymbol{\beta}), \end{cases}\end{split}\]

where

\[\begin{split}\begin{cases} f_{a_i}(\boldsymbol{a}, \boldsymbol{\beta}) & = \sum_{j=1}^{N_e} \epsilon^{j} f_{a_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta}) & = \sum_{j=1}^{N_e} \epsilon^{j} f_{\beta_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\beta}). \end{cases}\end{split}\]

The MMS system then obtained represents a nonlinear autonomous system of \(2N\) coupled \(1^{\textrm{st}}\) order PDEs, with the physical time \(t\) as the independent variable. Like the time scales-dependent modulation equations, they can be solved numerically or analytically, if analytical solutions exist.

These physical time-dependent modulation equations are also computed in modulation_equations().

Note that after reintroduction of the physical time, the leading order, homogeneous solution for oscillator \(i\) takes the form

\[x^{\textrm{h}}_{i,0}(t) = a_i(t) \cos\left( \frac{r_i}{r_{\textrm{MMS}}} \omega t - \beta_i(t)\right).\]

Complex form of the MMS

The application of the MMS to the system’s equations in oscillator form was described above. Alternatively, the MMS can be applied to the system’s equations in complex form (see section Complex form). The procedure is very similar, but the fact that the complex variables \(\boldsymbol{z}\) embed the state space \((\boldsymbol{x}, \dot{\boldsymbol{x}})\) allows to extract more information from the system’s dynamics, resulting in more accurate solutions and modulation equations. The differences between the oscillator and complex forms of the equations are stressed thereafter.

The multiple scales system

Series expansions of the \(z_i(t)\) variables are performed such that

\[z_i(t) = z_{i,0}(t) + \epsilon z_{i,1}(t) + \epsilon^2 z_{i,2}(t) + \cdots + \epsilon^{N_e} z_{i,N_e}(t) + \mathcal{O}(\epsilon^{N_e+1}),\]

time scales are introduced in the same way as before, the chain rule is used (only the \(1^{\text{st}}\) order time derivatives are needed), parameters are scaled, oscillators’ natural frequencies and the frequency of interest are expressed in terms of a reference frequency, possibly with detunings. This leads to the complex MMS system

\[\begin{split}\begin{aligned} & \epsilon^0 \rightarrow \; \begin{cases} \textrm{D}_0 z_{0,0} & = \textrm{j} \omega_{0,0} z_{0,0} + g_{0}^{(0)}(t_0, t_1), \\ & \vdots \\ \textrm{D}_0 z_{N-1,0} & = \textrm{j} \omega_{N-1,0} z_{N-1,0} + g_{N-1}^{(0)}(t_0, t_1), \end{cases} \\[15pt] & \epsilon^1 \rightarrow \; \begin{cases} \textrm{D}_0 z_{0,1} & = \textrm{j} \omega_{0,0} z_{0,1} + g_{0}^{(1)} (\boldsymbol{z}_0, \bar{\boldsymbol{z}}_0, \textrm{D}_0 \boldsymbol{z}_0, \textrm{D}_1 \boldsymbol{z}_0, t_0, t_1), \\ & \vdots \\ \textrm{D}_0 z_{N-1,1} & = \textrm{j} \omega_{0,0} z_{N-1,1} + g_{N-1}^{(1)}(\boldsymbol{z}_0, \bar{\boldsymbol{z}}_0, \textrm{D}_0 \boldsymbol{z}_0, \textrm{D}_1 \boldsymbol{z}_0, t_0, t_1), \end{cases} \\[10pt] & \hspace{3cm} \vdots \\[10pt] & \epsilon^{N_e} \rightarrow \; \begin{cases} \textrm{D}_0 z_{0,N_e} & = \textrm{j} \omega_{0,0} z_{0,N_e} + g_{0}^{(N_e)} (\boldsymbol{z}_0, \cdots, \textrm{D}_{N_e} \boldsymbol{z}_0, \cdots, \textrm{D}_1 \boldsymbol{z}_{N_e-1}, t_0, t_1), \\ & \vdots \\ \textrm{D}_0 z_{N-1,N_e} & = \textrm{j} \omega_{N-1,0} z_{N-1,N_e} + g_{N-1}^{(N_e)}(\boldsymbol{z}_0, \cdots, \textrm{D}_{N_e} \boldsymbol{z}_0, \cdots, \textrm{D}_1 \boldsymbol{z}_{N_e-1}, t_0, t_1). \end{cases} \end{aligned}\end{split}\]

Like with the oscillator form, these equations represent harmonic oscillators of natural frequencies \(+ \omega_{i,0}\) in the fast time, and subjected to forcing-like terms.

Leading order solution

In the same way as with the oscillator form, the leading order solution for oscillator \(i\) satisfies

\[\textrm{D}_0 z_{i,0} = \textrm{j} \omega_{i,0} z_{i,0} + g_{i}^{(0)}(t_0, t_1),\]

such that it takes the form

\[z_{i,0}(\boldsymbol{t}) = z_{i,0}^\textrm{h}(\boldsymbol{t}) + z_{i,0}^\textrm{p}(t_0, t_1).\]

These homogeneous and particular solutions can be expressed

\[\begin{split}\begin{cases} z_{i,0}^\textrm{h}(\boldsymbol{t}) & = A_i(\boldsymbol{t}_s) e^{\textrm{j} \omega_{i,0} t_0}, \\ z_{i,0}^\textrm{p}(t_0, t_1) & = B_{i,+} e^{\textrm{j} \omega t} + B_{i,-} e^{-\textrm{j} \omega t} = B_{i,+} e^{\textrm{j} (\omega_{\textrm{MMS}} t_0 + \sigma t_1)} + B_{i,-} e^{-\textrm{j} (\omega_{\textrm{MMS}} t_0 + \sigma t_1)}, \end{cases}\end{split}\]

where \(A_i\) is a slow time-dependent complex amplitude to be determined while \(B_{i,\pm}\) are time-independent functions of the forcing parameters. Note that in most situations, forcing does not appear at leading order (i.e. forcing is weak), so \(B_{i,\pm}=0\).

In the following it will be convenient to use the notations

\[\begin{split}\begin{split} \boldsymbol{A}(\boldsymbol{t}_s)^\intercal & = [A_0(\boldsymbol{t}_s), A_1(\boldsymbol{t}_s), \cdots, A_{N-1}(\boldsymbol{t}_s)], \\ \boldsymbol{B}_\pm^\intercal & = [B_{0,\pm}, B_{1,\pm}, \cdots, B_{N-1,\pm}], \\ \boldsymbol{B}^\intercal & = [\boldsymbol{B}_+^\intercal, \boldsymbol{B}_-^\intercal]. \end{split}\end{split}\]

The leading order solutions are defined in sol_order_0().

Higher order solutions

The higher order solutions can be derived in the same way as with the oscillator form, noticing that the terms in \(e^{+\textrm{j}\omega_{i,0} t_0}\) in the equation of oscillator \(i\) are secular and must therefore be eliminated. The \(1^{\textrm{st}}\) higher order solution of oscillator \(i\) is hence governed by the equation

\[\textrm{D}_0 z_{i,1} = \textrm{j} \omega_{i,0} z_{i,1} + \breve{g}_{i}^{(1)}(\boldsymbol{z}_0, \bar{\boldsymbol{z}}_0, \textrm{D}_0 \boldsymbol{z}_0, \textrm{D}_1 \boldsymbol{z}_0, t_0, t_1),\]

where \(\breve{g}_{i}^{(1)}\) is \(g_{i}^{(1)}\) with the secular terms cancelled.

It is interesting to stress that this cancellation of secular terms is equivalent to fulfilling the condition that the particular solution of a linear problem must be orthogonal to the kernel of the adjoint problem. To apply this, it is convenient to introduce the differential operator \(\mathscr{L}_i\) such that

\[\mathscr{L}_i z_{i,1} = (\textrm{D}_0 - \textrm{j} \omega_{i,0}) z_{i,1}.\]

This way, the initial problem writes

\[\mathscr{L}_i z_{i,1} = \breve{g}_{i}^{(1)},\]

and the associated adjoint problem is

\[\mathscr{L}^{\dagger}_i \breve{g}_{i}^{(1)} = z_{i,1},\]

where the adjoint operator \(\mathscr{L}^{\dagger}_i\) is defined such that

\[\langle \breve{g}_{i}^{(1)}, \mathscr{L}_i z_{i,1} \rangle = \langle \mathscr{L}^{\dagger}_i \breve{g}_{i}^{(1)}, z_{i,1} \rangle.\]

Given the \(\omega_{i,0}\)-periodicity of the homogeneous solution, it is relevant to use the metric associated to the \(L^2\) inner product

\[\langle f, g \rangle = \frac{\omega_{i,0}}{2\pi} \int_0^{2 \pi / \omega_{i,0}} \bar{f} g \textrm{d} t_0.\]

From there, one can show that the operator is anti self-adjoint, i.e. \(\mathscr{L}^{\dagger}_i = - \mathscr{L}_i\), and the adjoint kernel is

\[\textrm{ker}(\mathscr{L}^{\dagger}_i) = \{ e^{\textrm{j} \omega_{i,0} t_0} \}.\]

The condition to fulfil is therefore

\[\langle e^{\textrm{j} \omega_{i,0} t_0}, {g}_{i}^{(1)} \rangle = \frac{\omega_{i,0}}{2\pi} \int_0^{2 \pi / \omega_{i,0}} e^{-\textrm{j} \omega_{i,0} t_0} {g}_{i}^{(1)} \textrm{d} t_0 = 0,\]

such that the adjoint kernel orthogonality condition is equivalent to imposing a zero \(e^{\textrm{j} \omega_{i,0} t_0}\) component on \({g}_{i}^{(1)}\), i.e. cancelling the secular terms.

After the cancellation of secular terms, the higher order solutions \(z_{i,1}(\boldsymbol{t})\) can be expressed as a function of the leading order unknown amplitudes \(\boldsymbol{A}(\boldsymbol{t}_s)\), their slow time derivatives \(\textrm{D}_i\boldsymbol{A}(\boldsymbol{t}_s), \; i=1, ..., N_e\), and forcing terms if any (including the hard forcing amplitudes \(\boldsymbol{B}\)).

At this stage, there is a key difference to highlight between the present approach and the application of the MMS to the system in oscillator form. In the complex form, the secular terms cancelled are those oscillating at \(e^{+\textrm{j} \omega_{i,0} t_0}\), but not \(e^{-\textrm{j} \omega_{i,0} t_0}\). Therefore, when reconstructing the oscillator’s motion \(x_{i} = z_{i} + \bar{z}_{i}\), there can be higher order terms contributing to the term in \(\cos(\omega_{i,0} t_0), \; \sin(\omega_{i,0} t_0)\). In other words, the higher order particular solutions enrich the main oscillations term. On the contrary, with the oscillator form there are no contributions for higher order particular solutions. Indeed, the cancellation of secular terms is more stringent, imposing no terms in \(\cos(\omega_{i,0} t_0)\) nor \(\sin(\omega_{i,0} t_0)\), therefore setting to zero both \(e^{+\textrm{j} \omega_{i,0} t_0}\) and \(e^{-\textrm{j} \omega_{i,0} t_0}\) terms. Contributions to the main oscillations term can therefore only arise from the homogeneous higher order solutions, but there are no conditions imposing a proper choice of higher order homogeneous amplitudes. These conditions can only be derived from a formulation of the system’s dynamics that embeds its state space, as shown in [5].

There is a choice left for the higher order homogeneous solutions. With the present complex approach, it is relevant to set these to zero so as to fully capture oscillations at \(e^{+\textrm{j} \omega_{i,0} t_0}\) through the leading order term. Their amplitude in the full solutions \(z_i\) are thus simply the complex amplitudes \(A_i\).

As with the oscillator form, the process described above is repeated successively at each order, i.e. the computed solutions are introduced in the next higher order system of equations, secular terms are cancelled and the next higher order solutions are computed.

The secular terms are identified in secular_analysis() and the leading order solutions are computed in sol_higher_order().

Secular analysis and modulation equations

The secular analysis is carried out in the same way as with the oscillator form. A small, practical difference is that the secular terms are multiplied by \(\textrm{j}\) to retrieve the system of complex modulation equations

\[\begin{split}\begin{cases} \textrm{D}_1 A_i(\boldsymbol{t}_s) & = f_{A_i}^{(1)}(\boldsymbol{A}, t_1), \\ & \vdots \\ \textrm{D}_{N_e} A_i(\boldsymbol{t}_s) & = f_{A_i}^{(N_e)}(\boldsymbol{A}, t_1), \end{cases}\end{split}\]

with a structure similar to that of the oscillator form. From there on, the derivation of modulation equations is unchanged. Still, one must keep in mind that these modulation equations will govern the amplitude and phase evolution of the complex coordinates, only indirectly that of the oscillators’ motion.

The complex modulation equations are derived in secular_analysis().

Reconstruction of the oscillators’ motion

The oscillators’ motion is simply reconstructed at each order from the definition of the complex coordinates. For oscillator \(i\) and at order \(j\), this results in

\[x_{i,j}(\boldsymbol{t}) = z_{i,j}(\boldsymbol{t}) + \bar{z}_{i,j}(\boldsymbol{t}).\]

At this stage, the complex form approach branches back with the oscillator form one.

The oscillators’ motion is retrieved using sol_xO().

Evaluation at steady state

Steady state solutions

Steady state conditions

At steady state, the solutions’ amplitudes and phases are time-independent. One therefore has, for each oscillator \(i=1,...,N\),

\[\begin{split}\begin{cases} \dfrac{\textrm{d}}{dt} a_i & = 0, \\ \dfrac{\textrm{d}}{dt} \beta_i & = 0, \\ \end{cases}\end{split}\]

and the homogeneous steady state solutions take the form

\[x^{\textrm{h}}_{i,0}(t) = a_i \cos\left( \frac{r_i}{r_{\textrm{MMS}}} \omega t - \beta_i \right).\]

MMS modulation equations at steady state

The steady state amplitudes \(\boldsymbol{a}\) and phases \(\boldsymbol{\beta}\) are governed by the modulation equations evaluated at steady state, which take the form

\[\begin{split}\begin{cases} f_{a_0}(\boldsymbol{a}, \boldsymbol{\beta}) & = 0, \\ f_{\beta_0}(\boldsymbol{a}, \boldsymbol{\beta}) & = 0, \\ & \vdots \\ f_{a_{N-1}}(\boldsymbol{a}, \boldsymbol{\beta}) & = 0, \\ f_{\beta_{N-1}}(\boldsymbol{a}, \boldsymbol{\beta}) & = 0. \end{cases}\end{split}\]

This is now an algebraic system of equations.

MMS solutions at steady state

Solving the modulation equations at steady state yields the steady state solutions \(\boldsymbol{a}, \boldsymbol{\beta}\). The system can be solved directly for \(\boldsymbol{a}\) and \(\boldsymbol{\beta}\), yielding explicit analytical solutions, but this is often complex as the system is nonlinear.

A possibility that is sometimes available to obtain analytical solutions is to rearrange the equations to isolate the phase terms \(\cos(f(\boldsymbol{\beta})),\; \sin(f(\boldsymbol{\beta}))\) where \(f(\boldsymbol{\beta})\) is a linear function of the \(\beta_i\). Then the equations can be squared and summed up to obtain an equation on \(a_i\) only and/or \(a_j\) as a function of \(a_i\). The \(a_j\) can be expressed as a function of \(a_i\), leading to a polynomial equation on \(a_i\). The resulting polynomial equation can rarely be solved directly as the polynomial involved are often of high order. However, the polynomial is often quadratic in the detuning \(\sigma\) and forcing amplitude \(F\). It can therefore be solved for \(\sigma\) and \(F\) with \(a_i\) seen as a parameter. This yields an implicit solution for \(a_i(\sigma) \Rightarrow a_i(\omega)\) and \(a_i(F)\), from which one can deduce the other amplitudes \(a_j\) and phases \(\boldsymbol{\beta}\), thus reconstructing the oscillators’ solutions.

The processus described above is not always feasible. The blocking points are typically to

get rid of phase terms \(\cos(f(\boldsymbol{\beta})),\; \sin(f(\boldsymbol{\beta}))\) in the equations,
express every amplitude \(a_j\) in terms of a single one, \(a_i\),
end up with a polynomial of order 2 in \(\sigma\) and \(F\)

These difficulties become more pronounced when the system involves several oscillators, in which case the amplitudes and phases may only be computed numerically.

To facilitate the derivation of an analytical solution, it is possible to consider the backbone curve (bbc) of the forced solution rather than the forced solution itself. This bbc is computed in the absence of damping and forcing, therefore simplifying the system. Typically, this reduces the number of phase terms appearing. The solving procedure is then the same as that described previously.

Stability analysis

The stability analysis is described in details in Stability analysis.

Stability analysis

Stability information

Consider a steady state solution \((\hat{\boldsymbol{a}} , \hat{\boldsymbol{\beta}})\) such that, for \(i=1,...,N\),

\[\begin{split}\begin{cases} f_{a_i}(\hat{\boldsymbol{a}}, \hat{\boldsymbol{\beta}}) & = 0, \\ f_{\beta_i}(\hat{\boldsymbol{a}}, \hat{\boldsymbol{\beta}}) & = 0. \end{cases}\end{split}\]

The aim is to determine the stability state of that steady solution, which corresponds to a fixed point in the phase space.

Jacobian matrix

Let’s first introduce the vector of polar coordinates and polar modulation functions

\[\begin{split}\boldsymbol{x}^{(\textrm{p})\intercal} & = [a_0, \beta_0, \cdots, a_{N-1}, \beta_{N-1}], \\ \boldsymbol{f}^{(\textrm{p})\intercal} & = [f_{a_0}, f_{\beta_0}^*, \cdots, f_{a_{N-1}}, f_{\beta_{N-1}}^*].\end{split}\]

Note that the appearance of \(f_{\beta_i}^*\) in \(\boldsymbol{f}^{(\textrm{p})}\) requires \(a_i \neq 0\), which strongly constraints the type of steady state solutions that can be considered in the approach described below. To relax this constraint, one can use a change of coordinates from polar to cartesian ones. This will be discussed in following sections, after the description of this polar approach.

Using the vectors of polar coordinates and modulation functions, one can write the modulation equations system as

\[\dfrac{\textrm{d} \boldsymbol{x}^{(\textrm{p})}}{\textrm{d}t} = \textrm{J}^{(\textrm{p})} \boldsymbol{x}^{(\textrm{p})},\]

where we introduced the Jacobian matrix

\[\begin{split}\textrm{J}^{(\textrm{p})} = \dfrac{\partial \boldsymbol{f}^{(\textrm{p})} }{ \partial \boldsymbol{x}^{(\textrm{p})} } = \begin{bmatrix} \frac{\partial f_{a_0}}{\partial a_0} & \frac{\partial f_{a_0}}{\partial \beta_0} & \cdots & \frac{\partial f_{a_0}}{\partial \beta_{N-1}} \\ \frac{\partial f_{\beta_0}^*}{\partial a_0} & \frac{\partial f_{\beta_0}^*}{\partial \beta_0} & \cdots & \frac{\partial f_{\beta_0}^*}{\partial \beta_{N-1}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_{\beta_{N-1}}^*}{\partial a_0} & \frac{\partial f_{\beta_{N-1}}^*}{\partial \beta_0} & \cdots & \frac{\partial f_{\beta_{N-1}}^*}{\partial \beta_{N-1}} \end{bmatrix}.\end{split}\]

Perturbation of the steady state solution

Let us now consider a small perturbation \(\tilde{\boldsymbol{x}}^{(\textrm{p})}\) of the steady state solution such that

\[\boldsymbol{x}^{(\textrm{p})} = \hat{\boldsymbol{x}}^{(\textrm{p})} + \tilde{\boldsymbol{x}}^{(\textrm{p})} \quad \Leftrightarrow \quad \tilde{\boldsymbol{x}}^{(\textrm{p})} = \boldsymbol{x}^{(\textrm{p})} - \hat{\boldsymbol{x}}^{(\textrm{p})}.\]

Using a first order Taylor expansion for \(\textrm{d} \tilde{\boldsymbol{x}}^{(\textrm{p})} / \textrm{d} t\), one can write

\[\dfrac{\textrm{d} \tilde{\boldsymbol{x}}^{(\textrm{p})}}{\textrm{d}t} = \left.\textrm{J}^{(\textrm{p})}\right|_{\hat{\boldsymbol{x}}^{(\textrm{p})}} \tilde{\boldsymbol{x}}^{(\textrm{p})} + \mathcal{O}(||\tilde{\boldsymbol{x}}^{(\textrm{p})}||^2),\]

where \(\left.\textrm{J}^{(\textrm{p})}\right|_{\hat{\boldsymbol{x}}^{(\textrm{p})}}\) denotes the Jacobian matrix evaluated on the steady state solution. The perturbation solution takes the form

\[\tilde{\boldsymbol{x}}^{(\textrm{p})} = \sum_{i=1}^{2N} C_i \boldsymbol{\psi}_i e^{\lambda_i t},\]

where \((\lambda_i, \boldsymbol{\psi}_i),\; i=1, ..., 2N\) are the eigensolutions of the Jacobian (evaluated on \(\hat{\boldsymbol{x}}^{(\textrm{p})}\)).

Stability condition

The steady state solution \(\hat{\boldsymbol{x}}^{(\textrm{p})}\) is considered stable if a small perturbation \(\tilde{\boldsymbol{x}}^{(\textrm{p})}\) vanishes in time, such that solutions close from \(\hat{\boldsymbol{x}}^{(\textrm{p})}\) are converging towards it. This condition is fulfilled if

\[\Re[\lambda_i] < 0, \quad \forall i,\]

meaning that all eigenvalues of the Jacobian evaluated on the steady state solution must have negative real parts. If this condition is not met, the system is either quasi stable or unstable.

Bifurcations

A bifurcation occurs when the stability state of the system changes.

Simple bifurcations occur when at least one eigenvalue crosses the imaginary axis through 0.
Such bifurcations include saddle node and pitchfork bifurcations, which cause jumps of the response and the appearance of lower symmetry solutions, respectively.
Neimark-Sacker bifurcations occur when a pair of complex conjugate eigenvalues with nonzero real parts cross the imaginary axis.
These bifurcations lead to non periodic solutions.

Simple bifurcations can be detected by evaluating the sign of the Jacobian, as

\[\det \left.\textrm{J}^{(\textrm{p})}\right|_{\hat{\boldsymbol{x}}^{(\textrm{p})}} = \prod_{i=1}^{2N} \lambda_i,\]

thereby making \(\det \left.\textrm{J}^{(\textrm{p})}\right|_{\hat{\boldsymbol{x}}^{(\textrm{p})}}\) an important stability indicator. Neimark-Sacker bifurcations are more difficult to detect. Information from the trace of \(\left.\textrm{J}^{(\textrm{p})}\right|_{\hat{\boldsymbol{x}}^{(\textrm{p})}}\) can be considered, or the Routh-Hurwitz criterion can be used. This is not detailed here.

Bifurcation curves

Bifurcation curves are curves constructed by evaluating the coordinates of bifurcation points when varying one or more parameters. The stability state of a solution changes when the response curve crosses a bifurcation curve.

Stability analysis in cartesian coordinates

Limitations of polar coordinates

As mentioned previously, the approach described above fails when one of the oscillator’s leading order amplitude is 0. Indeed, the Jacobian is constructed using the modulation functions

\[f_{\beta_i}^*(\boldsymbol{a}, \boldsymbol{\beta}) = \frac{f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta})}{a_i},\]

which are defined only if \(a_i=0\). This prevents evaluating the stability of

Trivial solutions, for which no oscillator responds,
1 mode solutions, whose stability can be affected under perturbation from another mode.

These limitations can be overcome using a change of coordinates.

Cartesian coordinates

The leading order homogeneous solution for oscillator \(i\) can be written as

\[\begin{split}\begin{split} x_{i,0}^{\textrm{h}}(t) & = a_i \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t - \beta_i \right), \\ & = a_i \cos(\beta_i) \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right) + a_i \sin(\beta_i) \sin\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right). \end{split}\end{split}\]

It then appears natural to introduce the cartesian coordinates

\[\begin{split}\begin{cases} p_i = a_i \cos(\beta_i), \\ q_i = a_i \sin(\beta_i), \end{cases}\end{split}\]

in order to rewrite the solution as

\[x_{i,0}^{\textrm{h}}(t) = p_i \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right) + q_i \sin\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right).\]

In the following it will be convenient to use the cartesian coordinates vectors

\[\begin{split}\begin{aligned} \boldsymbol{p}(t)^\intercal & = [p_0(t), p_1(t), \cdots, p_{N-1}(t)], \\ \boldsymbol{q}(t)^\intercal & = [q_0(t), q_1(t), \cdots, q_{N-1}(t)]. \end{aligned}\end{split}\]

Cartesian modulation equations

The cartesian modulation equations can be obtained from the polar ones. To do so, one can write

\[\begin{split}\begin{cases} \dfrac{\textrm{d} p_i}{\textrm{d}t} & = \dfrac{\textrm{d} a_i}{\textrm{d}t} \cos(\beta_i) - a_i \sin(\beta_i) \dfrac{\textrm{d} \beta_i}{\textrm{d}t}, \\ \dfrac{\textrm{d} q_i}{\textrm{d}t} & = \dfrac{\textrm{d} a_i}{\textrm{d}t} \sin(\beta_i) + a_i \cos(\beta_i) \dfrac{\textrm{d} \beta_i}{\textrm{d}t}. \end{cases}\end{split}\]

Then, by identification, one necessarily has

\[\begin{split}\begin{cases} f_{p_i}(\boldsymbol{p}, \boldsymbol{q}) & = f_{a_i}(\boldsymbol{a}, \boldsymbol{\beta}) \cos(\beta_i) - f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta}) \sin(\beta_i), \\ f_{q_i}(\boldsymbol{p}, \boldsymbol{q}) & = f_{a_i}(\boldsymbol{a}, \boldsymbol{\beta}) \sin(\beta_i) + f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta}) \cos(\beta_i), \end{cases}\end{split}\]

in order to write the modulation equations in cartesian coordinates

\[\begin{split}\begin{cases} \dfrac{\textrm{d}}{dt} p_0(t) & = f_{p_0}(\boldsymbol{p}, \boldsymbol{q}), \\ \dfrac{\textrm{d}}{dt} q_0(t) & = f_{q_0}(\boldsymbol{p}, \boldsymbol{q}), \\ & \vdots \\ \dfrac{\textrm{d}}{dt} p_{N-1}(t) & = f_{p_{N-1}}(\boldsymbol{p}, \boldsymbol{q}), \\ \dfrac{\textrm{d}}{dt} q_{N-1}(t) & = f_{q_{N-1}}(\boldsymbol{p}, \boldsymbol{q}). \end{cases}\end{split}\]

Jacobian matrix

As done previously for polar coordinates, let’s introduce the vectors of cartesian coordinates and modulation equations as

\[\begin{split}\boldsymbol{x}^{(\textrm{c})\intercal} & = [p_0, q_0, \cdots, p_{N-1}, q_{N-1}], \\ \boldsymbol{f}^{(\textrm{c})\intercal} & = [f_{p_0}, f_{q_0}, \cdots, f_{p_{N-1}}, f_{q_{N-1}}].\end{split}\]

Then one can write

\[\dfrac{\textrm{d} \boldsymbol{x}^{(\textrm{c})}}{\textrm{d}t} = \textrm{J}^{(\textrm{c})} \boldsymbol{x}^{(\textrm{c})},\]

where we introduced the Jacobian matrix

\[\begin{split}\textrm{J}^{(\textrm{c})} = \dfrac{\partial \boldsymbol{f}^{(\textrm{c})} }{ \partial \boldsymbol{x}^{(\textrm{c})} } = \begin{bmatrix} \frac{\partial f_{p_0}}{\partial p_0} & \frac{\partial f_{p_0}}{\partial q_0} & \cdots & \frac{\partial f_{p_0}}{\partial q_{N-1}} \\ \frac{\partial f_{q_0}}{\partial p_0} & \frac{\partial f_{q_0}}{\partial q_0} & \cdots & \frac{\partial f_{q_0}}{\partial q_{N-1}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_{q_{N-1}}}{\partial p_0} & \frac{\partial f_{q_{N-1}}}{\partial q_0} & \cdots & \frac{\partial f_{q_{N-1}}}{\partial q_{N-1}} \end{bmatrix}.\end{split}\]

Note that there are no constraints related to an oscillator’s amplitude being 0 here. This cartesian coordinates approach therefore allows to investigate how the stability of a steady state solution is affected by a perturbation from an oscillator who’s amplitude is 0 in that steady state solution.

The stability analysis with \(\textrm{J}^{(\textrm{c})}\) is carried out as described previously with \(\textrm{J}^{(\textrm{p})}\).