3.1. Equilibria and bifurcations in A

We solve for equilibria $(\chi _{\text {eq}}, v_{\text {eq}})$ by setting the derivatives in the system (3.1) to zero. This immediately gives $v_{\text {eq}} = 0$ as expected on physical grounds, but the equation for $\chi _{\text {eq}}$ must be solved numerically. For illustrative purposes, in the following analysis, we show some results for $H = 10$ , as this is a reasonable value for a physical oscillator one could construct from permanent magnetsFootnote ¹. (The effect of varying both H and $A_0$ is discussed later in Section 3.2.) At $A_0 = 0$ , we obtain an equilibrium value $\chi _{\text {eq}} \approx 0.0248$ using Newton’s method.

We then construct a bifurcation diagram that tracks $\chi _{\text {eq}}$ as $A_0$ is increased from zero. The bifurcation diagram is generated exactly by employing a parametrization approach, which we briefly outline here. For a fixed value of H, we solve

(3.2) $$ \begin{align} \chi_{\text{eq}} = \frac{A_0 - H \chi_{\text{eq}}} {[ 1 + ( A_0 - H \chi_{\text{eq}} )^2 ]^{5/2}}. \end{align} $$

We seek $\chi _{\text {eq}}$ and $A_0$ as functions of a parameter $\theta $ . Thus, the change of variable ${H \chi _{\text {eq}} = A_0 - \tan \theta _{\text {eq}}}$ is made in (3.2). This gives rise to the algorithm:

• • give values of the parameter $\theta _{\text {eq}}$ : $0 < \theta _{\text {eq}} < \pi /2$ ;

• • calculate $\chi _{\text {eq}} = \sin \theta _{\text {eq}} \cos ^4 \theta _{\text {eq}}$ from (3.2);

• • create $A_0 =\tan \theta _{\text {eq}} + H \chi _{\text {eq}}$ from the change of variables.

The bifurcation curve is now able to be drawn in the parametric form $( A_0 ( \theta _{\text {eq}} ) , \chi _{\text {eq}} ( \theta _{\text {eq}} ) )$ . Note that the bounds imposed on $\theta _{\text {eq}}$ restrict $\chi _{\text {eq}}$ (and subsequently $A_0$ ) to physically relevant positive values.

Figure 2 displays a bifurcation diagram for the equilibrium location $\chi _{\text {eq}}$ of magnet 1, against the fixed location $A_0$ of magnet 2, for the illustrative case $H = 10$ , and similar curves are obtained with many other values of the force H. The two small squares on the plot indicate the locations of saddle-node bifurcations, at which the curve folds over and the equilibrium location $\chi _{\text {eq}}$ for magnet 1 becomes multi-valued. For $H = 10$ , these events occur at approximately $A_0 = 2.267$ and $A_0 = 3.392$ .

Figure 2 Bifurcation diagram for parameter $A_0$ , with $H = 10$ . The two small red squares on the curve denote the locations of saddle-node bifurcations (SN) at $A_0 \approx 2.267$ and $A_0 \approx 3.392$ . The dashed section indicates an unstable equilibrium.

The dashed portion of this curve signifies that the equilibrium is unstable. The stability of the equilibria in Figure 2 is obtained in the usual way, by linearization of the governing equations (3.1) near $\chi _{\text {eq}}$ . Further details are not provided here (although Section 3.3 discusses the particular case of pure centres). This then shows that, in the interval $0 \leq A_0 < 2.267$ , Figure 2 indicates there is a single equilibrium and that it is a (neutrally stable) centre, characterized by pure oscillations centred about the equilibrium. At $A_0 \approx 2.267$ , a second centre and a saddle are born from a saddle-node bifurcation. Then, in the interval $2.267 < A_0 < 3.392$ , there are two (neutrally stable) centres separated in the phase plane by a (unstable) saddle. We shall from hereon refer to this as the bistable region in parameter space. At $A_0 \approx 3.392$ , the saddle collides with the right-hand centre and both equilibria disappear in another saddle-node bifurcation, leaving the system with a single centre for $A_0> 3.392$ .

The bifurcation analysis presented above only describes the linear behaviour near the equilibria. To understand any effects of the nonlinearity in the system, phase-plane portraits are plotted for representative values of $A_0$ in Figures 3(a), 3(b) and 3(c) for the three values $A_0 = 1$ , $A_0 = 2.5$ and $A_0 = 4$ , respectively. To do this, the system (3.1) was integrated numerically from a set of initial conditions evenly spaced along the $\chi $ -axis, using a fourth-order Runge–Kutta method. These plots confirm that the linearized analysis does indeed give a reliable qualitative description of the true nonlinear dynamics near the equilibria. The system has two homoclinic orbits in the bi-stable region ( $2.267 < A_0 < 3.392$ ) and a single homoclinic orbit at each saddle-node bifurcation.

Figure 3 Phase portraits for $H = 10$ for the three stationary magnet 2 locations (a) $A_0 = 1$ , (b) $A_0 = 2.5$ and (c) $A_0 = 4$ , showing centre and saddle equilibria. Each separate curve represents a solution trajectory (with different colours).

The system (3.1) is an undamped ( $D=0$ ), autonomous, Hamiltonian system. Therefore, its solution trajectories are closed curves (here, either homoclinic or periodic) given by the level sets of the Hamiltonian (that is, by contours of constant energy). As a result, the system has a first integral, so that an algebraic expression in terms of $\chi $ and v may be found for every solution trajectory, although this is not developed further here. Further, we see from the phase portrait of the bi-stable case in Figure 3(b) that the system resembles an undamped particle in a double-well potential with two local energy minima, whereas the phase portraits shown in Figures 3(a) and 3(c) only have one local minimum at the centre.

3.2. Codimension-1 and 2 bifurcations in A and H

The bifurcation diagram in Figure 2 for the single parameter $H = 10$ showed two saddle-node bifurcations, at which the curve folded to produce multiple equilibria. The parametrization technique used to produce this graph was outlined in Section 3.1 and we generalize it to track the locations of the saddle-node bifurcations for all values of H. Since these saddle-node bifurcation points occur at equilibria, the algebraic conditions in Section 3.1 still hold and a further condition must be added to them, to seek those points at which the parametrized bifurcation curve $( A_0 ( \theta _{\text {eq}} ) , \chi _{\text {eq}} ( \theta _{\text {eq}} ) )$ becomes vertical in Figure 2. This occurs when

$$ \begin{align*} \frac{\mathit{d} A_0}{\mathit{d} \theta_{\text{eq}}} = \frac{\mathit{d}}{\mathit{d} \theta_{\text{eq}}}(\tan\theta_{\text{eq}} + H \chi_{\text{eq}}) = 0. \end{align*} $$

Differentiation of the parametric equations in the algorithm given in Section 3.1 yields

$$ \begin{align*} \frac{\mathit{d} A_0}{\mathit{d} \theta_{\text{eq}}} = \sec^2 \theta_{\text{eq}} + H \cos^3 \theta_{\text{eq}} ( 5\cos^2 \theta_{\text{eq}} - 4 ) = 0. \end{align*} $$

This equation is rearranged and solved for H in terms of the parameter $\theta _{\text {eq}}$ . The location of the saddle-node bifurcation, for primary bifurcation constant $A_0$ and unfolding constant H, is calculated from the pair of equations

(3.3) $$ \begin{align} \begin{aligned} &A_0 = \tan\theta_{\text{eq}} \bigg[ 1 + \frac{1}{( 4 - 5\cos^2 \theta_{\text{eq}} )} \bigg],\\ &H = \frac{1}{\cos^5 \theta_{\text{eq}} ( 4 - 5\cos^2 \theta_{\text{eq}} )}, \end{aligned} \end{align} $$

for $0 < \theta _{\text {eq}} < \pi / 2$ .

The saddle-node bifurcation curve is plotted in Figure 4 and was calculated in parametric form from (3.3). This diagram shows the system behaviour for all positions of magnet 2 ( $A_0$ ) and force ratios (H). The bistable region (for which we have two centres and a saddle equilibrium) lies between the two saddle-node bifurcation curves, labelled SN in the diagram. A cusp forms in the saddle-node curve at approximately the point $(A_0 ,H) = (1.6238,3.5449)$ , and this is indicated on the figure. If $A_0 < 1.6238$ , then the system is stable (with a single centre equilibrium) for all $H> 0$ and the same is true when $H < 3.5449$ as we vary $A_0$ .

Figure 4 Stability diagram showing the saddle-node bifurcation curve (labelled SN) in the ( $A_0$ , H) parameter plane. At the indicated cusp point, there is a codimension-2 cusp bifurcation.

The location of the cusp point is revealed from the parametric description of the saddle-node curve, along which the condition $\mathit {d} A_0 / \mathit {d} \theta _{\text {eq}} = 0$ already applies for this codimension-1 bifurcation. To obtain the location of the cusp, which is a codimension-2 condition, we append to this the additional requirement $\mathit {d} H / \mathit {d} \theta _{\text {eq}} = 0$ . This is obtained by differentiating the expression for H in the SN condition (3.3), which yields

$$ \begin{align*} \frac{\mathit{d} H}{\mathit{d} \theta_{\text{eq}}} = \frac{5 \sin\theta_{\text{eq}} ( 4 - 7\cos^2 \theta_{\text{eq}} )} {\cos^6 \theta_{\text{eq}} ( 4 - 5\cos^2 \theta_{\text{eq}} )^2} = 0. \end{align*} $$

The root of interest to this equation is

$$ \begin{align*} \cos^2 \theta_{\text{eq}} = 4/7, \end{align*} $$

and this is substituted into the parametric conditions (3.3) for the saddle-node curve. This shows that the cusp occurs at the point

$$ \begin{align*} &A_0 = 15 \sqrt{3} / 16 \approx 1.6238,\\ &H = 7^3 \sqrt{7} / 4^4 \approx 3.5449 , \end{align*} $$

and this point is indicated on the diagram in Figure 4.

3.3. Natural frequencies

Here, we compute the natural frequency of each periodic orbit in the unforced and undamped system. This was done to derive a resonance condition for the forced system considered later. This will enable the forcing frequency to be related to the natural frequency of the corresponding unforced orbit.

The eigenvalues of the Jacobian matrix at a centre equilibrium give an analytic expression for the natural frequency of nearby orbits. The Jacobian matrix for our general system takes the form

$$ \begin{align*} \begin{bmatrix} 0 & 1 \\ - \mu_N^2 & -D \end{bmatrix} \quad\textrm{with } \mu_N^2 = 1 + \frac{H [ 1 - 4( A - H \chi)^2 ]}{[ 1 + (A - H \chi)^2 ]^{7/2}}. \end{align*} $$

In the absence of damping ( $D=0$ ), this matrix has purely imaginary eigenvalues $\pm \mu _N i$ for the centres, in which case, $\mu _N^2> 0$ . Close to a centre in the $( \chi , v )$ phase plane, the natural frequency associated with that centre at $\chi _{\text {eq}}$ is $\mu _N$ with corresponding period $2\pi / \mu _N$ .

To obtain natural frequencies of orbits far from a centre in the phase plane, we return to numerical methods and take the dominant frequency from each orbit’s frequency spectrum. Natural frequencies for $H=10$ and $A=2.4$ are illustrated in the phase plane in Figure 5(a) and along the $\chi $ -axis in Figure 5(b).

Figure 5 (a) Phase plane coloured by natural frequency of orbits for $H=10$ and $A=2.4$ . (b) Natural frequencies along the $\chi $ -axis.

The natural frequency curve has local maxima at the centre equilibria and, as expected, approaches zero when crossing homoclinic orbits. For large displacements $\chi $ , the magnetic force becomes negligible compared with the spring force, so our natural frequency approaches that of the spring, which is $1$ in our dimensionless variables.

In the forced case, to be considered later in Section 5, we will vary A between $2.4$ and $2.8$ for the case $H=10$ . For reference, we therefore show in Figure 6 the change in the natural frequency curve as A is varied in increments of $0.05$ . The movement of the sharp cusps, indicative of the homoclinic orbits, will affect the qualitative behaviour of the forced system. We use this to guide the choice of initial conditions in the subsequent section.

Figure 6 Natural frequency curves for varying A at $H=10$ .

4.1. Linearized forced solution

The location $A(\tau )$ of the lower magnet consists of an average position $A_0$ plus a periodically forced component of amplitude $\epsilon $ , as described in (2.4). If $\epsilon $ is assumed to be a small parameter, then a linearized solution for the location of upper magnet 1 is obtained in the form

(4.1) $$ \begin{align} \chi(\tau) = \chi_{\text{eq}} + \epsilon X_1(\tau) + \mathcal{O} ( \epsilon^2 ) , \end{align} $$

in which $\chi _{\text {eq}}$ is an equilibrium position for the magnet, calculated from (3.2). This approximate form (4.1) is substituted into the nonlinear system (2.6) and only terms of first order in $\epsilon $ are retained. This eventually leads to the linearized differential equation

(4.2) $$ \begin{align} \frac{\mathit{d}^2 X_1}{\mathit{d} \tau^2} + D \frac{\mathit{d} X_1}{\mathit{d} \tau} + ( 1 + H \Gamma_0 ) X_1 = \Gamma_0 \cos ( \Omega \tau ) \end{align} $$

for the linearized perturbation function $X_1 (\tau )$ in (4.1), where we define

(4.3) $$ \begin{align} \Gamma_0 = \frac{1 - 4 ( A_0 - H \chi_{\text{eq}} )^2} {[ 1 + ( A_0 - H \chi_{\text{eq}} )^2 ]^{7/2}}. \end{align} $$

The solution to the linearized equation (4.2) consists of transients, plus a periodic forced response which we write as

(4.4) $$ \begin{align} X_1 (\tau) = \frac{\Gamma_0 [ \mathcal{T}_1 \cos (\Omega \tau) + \mathcal{T}_2 \sin (\Omega \tau) ]} {( \mathcal{T}_1 )^2 + ( \mathcal{T}_2 )^2}. \end{align} $$

This expression involves the two constants

$$ \begin{align*} \mathcal{T}_1 = 1 + H \Gamma_0 - \Omega^2 \quad \textrm{and} \quad \mathcal{T}_2 = D \Omega. \end{align*} $$

For later use, we note that (4.4) describes a sinusoidal oscillation with amplitude

(4.5) $$ \begin{align} \mathcal{A}_1 = \frac{\Gamma_0}{\sqrt{(\mathcal{T}_1 )^2 + ( \mathcal{T}_2 )^2}}. \end{align} $$

In the absence of damping, $D=0$ , a (linearized) primary resonance occurs when the forcing frequency has the value

(4.6) $$ \begin{align} \Omega_{\text{res}} = \sqrt{1 + H \Gamma_0}. \end{align} $$

4.2. Weakly nonlinear solution near resonance

Further understanding of the nonlinear structure of the primary resonance in (4.4) is gained through an asymptotic analysis of the nonlinear equation (2.6) to third order, near the resonant frequency (4.6). To this end, we write

$$ \begin{align*} \begin{aligned} \chi(\tau) &= \chi_{\text{eq}} + X(\tau) , \\ A(\tau) &= A_0 + A_F(\tau) \quad \textrm{with } A_F(\tau) = \epsilon \cos (\Omega\tau). \end{aligned} \end{align*} $$

To simplify the analysis, we only consider the undamped case $D=0$ here. We further define the variable

$$ \begin{align*} Y(\tau) = A_F (\tau) - H X(\tau), \end{align*} $$

and carry out a Taylor-series expansion of the nonlinear equation (2.6) to third order in this variable.

After very considerable algebra, we obtain the weakly nonlinear approximate differential equation

$$ \begin{align*} \frac{\mathit{d}^2 Y}{\mathit{d} \tau^2} + ( 1 + H \Gamma_0 ) Y - H \Delta_0 Y^2 - H \Lambda_0 Y^3 = \epsilon ( 1 - \Omega^2 ) \cos (\Omega\tau ). \end{align*} $$

In this expression, the constant $\Gamma _0$ is as defined in (4.3) and there are two additional constants

$$ \begin{align*} \begin{aligned} \Delta_0 & = \frac{5 Y_{\text{eq}} [ 3 - 4 ( Y_{\text{eq}} )^2 ]} {2 [ 1 + ( Y_{\text{eq}} )^2 ]^{9/2}}, \nonumber \\ \Lambda_0 & = \frac{5 [ 1 - 12 ( Y_{\text{eq}} )^2 + 8 ( Y_{\text{eq}} )^4 ]} {2 [ 1 + ( Y_{\text{eq}} )^2 ]^{11/2}} , \end{aligned} \end{align*} $$

in which the constant $Y_{\text {eq}} = A_0 - H \chi _{\text {eq}} $ is also defined.

Near the primary resonance (4.6), we seek an approximate solution in the form

(4.7) $$ \begin{align} Y(\tau) = \mathcal{A}_W \cos( \Omega\tau) + \mathcal{B}_W \sin( \Omega\tau) ,\end{align} $$

and retain only terms involving the forcing frequency $\Omega $ . A substantial amount of algebra gives rise to the two algebraic equations

(4.8) $$ \begin{align} \begin{aligned} \mathcal{A}_W [ 1 + H \Gamma_0 - \Omega^2 - \tfrac{3}{4} H \Lambda_0 ( \mathcal{A}_W^2 + \mathcal{B}_W^2 ) ] & = \epsilon ( 1 - \Omega^2 ), \\ \mathcal{B}_W [ 1 + H \Gamma_0 - \Omega^2 - \tfrac{3}{4} H \Lambda_0 ( \mathcal{A}_W^2 + \mathcal{B}_W^2 ) ] & = 0 \end{aligned} \end{align} $$

for the two amplitude constants $\mathcal {A}_W$ and $\mathcal {B}_W$ in the assumed form (4.7) of the weakly nonlinear solution.

In the absence of damping, $D=0$ , the appropriate solution in the system (4.8) is simply $\mathcal {B}_W = 0$ , from which it follows that $\mathcal {A}_W$ satisfies the cubic equation

(4.9) $$ \begin{align} \mathcal{A}_W^3 + \frac{4[ \Omega^2 - 1 - H \Gamma_0 ]}{3 H \Lambda_0} \mathcal{A}_W + \frac{4\epsilon ( 1 - \Omega^2 )}{3 H \Lambda_0} = 0. \end{align} $$

Since this cubic has either one or three real roots for $\mathcal {A}_W$ , it follows from (4.9) that, near the primary resonance, there is either one or else three forced orbits, and this is discussed more fully in Section 5.6. To solve this equation numerically for a given $\Omega $ , we cast (4.9) in terms of a companion matrix and then used MATLAB’s eig routine to calculate the roots $\mathcal {A}_W$ , keeping only those with zero imaginary part.

4.3. Computing forced periodic orbits

As a check on the weakly nonlinear results of Section 4.2, we also compute periodic forced oscillations of the fully nonlinear model (2.4). We express the displacement of the upper mass in the form

(4.10) $$ \begin{align} \chi(\tau) = C_0 + \sum_{n=1}^N [ C_n \cos(n\Omega\tau) + S_n \sin(n\Omega\tau) ], \end{align} $$

and use Newton’s method to calculate the vector of coefficients

$$ \begin{align*} \mathbf{u} = [ C_0 , C_1 , \ldots , C_N; S_1 , \ldots , S_N ]. \end{align*} $$

The number N of Fourier modes in (4.10) is taken to be as large as practicable and we use $N=151$ . The Newton–Galerkin method we employ involves multiplying the nonlinear ODE (2.6) successively by the basis functions $\cos (\ell \Omega \tau )$ , $\ell =0,1,2, \ldots ,N$ and $\sin (\ell \Omega \tau )$ , $\ell = 1,2, \ldots ,N$ , and integrating over one period $0 < \tau < 2\pi / \Omega $ . This creates the error vector

$$ \begin{align*} \mathbf{E} = [ E_0 , E_1 , \ldots , E_N; E_{N+1} , \ldots , E_{2N+1} ]. \end{align*} $$

This has components

$$ \begin{align*} \begin{aligned} E_0 & = C_0 - \frac{\Omega}{2\pi} \mathcal{R}^{(C)}_0, \\ E_{\ell} & = C_{\ell} [ 1 - \ell^2 \Omega^2 ] + S_{\ell} ( D\ell\Omega ) - \frac{\Omega}{\pi} \mathcal{R}^{(C)}_{\ell} , \quad \ell = 1 , \ldots , N, \\ E_{N + \ell} & = - C_{\ell} ( D\ell\Omega ) + S_{\ell} [ 1 - \ell^2 \Omega^2 ] - \frac{\Omega}{\pi} \mathcal{R}^{(S)}_{\ell} , \quad \ell = 1 , \ldots , N , \end{aligned} \end{align*} $$

in which we define

$$ \begin{align*} \begin{aligned} \mathcal{R}^{(C)}_{\ell} & = \int_0^{2\pi / \Omega} \frac{ [ A(\tau) - H \chi ] \cos(\ell\Omega\tau)} {( 1 + [ A(\tau) - H \chi ]^2 )^{5/2}} \, \mathit{d} \tau , \quad \ell = 0, 1 , 2 , \ldots , N, \\ \mathcal{R}^{(S)}_{\ell} & = \int_0^{2\pi / \Omega} \frac{ [ A(\tau) - H \chi ] \sin(\ell\Omega\tau)} {( 1 + [ A(\tau) - H \chi ]^2 )^{5/2}} \, \mathit{d} \tau , \quad \ell = 1 , 2 , \ldots , N. \end{aligned} \end{align*} $$

Newton’s method is used to solve the system of algebraic equations $\mathbf {E} ( \mathbf {u} ) = \mathbf {0}$ . This creates the coefficients needed in the Fourier representation (4.10) of the displacement $\chi (\tau )$ . The amplitude of these nonlinear periodic solutions calculated from the numerical solution is then

(4.11) $$ \begin{align} \mathcal{A}_N = \tfrac{1}{2} ( \max \{ \chi(\tau) \} - \min \{ \chi(\tau) \} ). \end{align} $$

5.1. Solution behaviour: quasiperiodic orbits

We first explore quasiperiodic orbits, which appear in areas of phase space sufficiently far from the homoclinic connections of the unforced system. Example trajectories of this type in phase space for $\Omega = 0.1$ are shown in Figure 7. Because this is a forced system, solution trajectories are now objects in a higher-dimensional phase space and so cannot be characterized purely using the $(\chi , v)$ phase plane as for the previous Section 3. Nevertheless, the quasiperiodic orbits in Figure 7 are presented here as projections onto the $\chi , v$ axes, which contain some memory of the unforced system. Thus, there is a moderate-amplitude quasiperiodic orbit centred on the left-hand unforced equilibrium (drawn in red) and a second one centred on the right-hand equilibrium (drawn in yellow). We denote these orbits as $\Gamma _r$ and $\Gamma _y$ , respectively. A quasiperiodic orbit of much larger amplitude (in blue) surrounds all three equilibria. This is denoted as $\Gamma _b$ . These three solution trajectories were obtained simply by varying the initial value of $\chi $ , so that $\chi _0 = 0.06$ , $0.2$ and $-0.15$ . The initial condition of v was set as $v_0=0$ for all three.

Figure 7 Examples of quasiperiodic orbits for $\Omega = 0.1$ . The initial conditions used were $v_0 = 0$ and $\chi _0 = -0.15$ for the large orbit ( $\Gamma _b$ in blue), $\chi _0 = 0.06$ for the left-hand smaller orbit ( $\Gamma _r$ in red) and $\chi _0 = 0.2$ for the smaller right-hand orbit ( $\Gamma _y$ in yellow).

To explore further the indicators of quasiperiodicity, we focus on $\Gamma _r$ in Figure 7 (obtained with $\chi _0 = 0.06$ ). Its time series $\chi (\tau )$ is drawn in Figure 8(a) and shows the almost-periodic nature of the orbit. We also plot the frequency spectrum of its time series in Figure 8(b). This was generated via a discrete Fourier transform of the time-series data. For a quasiperiodic orbit, we expect the solution to consist of a finite number of distinct frequencies which, importantly, are not all rational multiples of each other. This is not immediately discernible from the frequency spectrum so we turn to additional indicators of quasiperiodicity.

Figure 8 (a) Example time series $\chi (\tau )$ for $\Gamma _r$ in Figure 7 and (b) the corresponding frequency spectrum showing a discrete set of frequencies.

We construct a stroboscopic map for the example orbit by sampling our trajectory (starting at time $\tau = 0$ ) with a “strobing period” equal to our parametric forcing period; that is, we let $T_{\text {strobe} }= 2 \pi / \Omega $ . For quasiperiodic orbits, the stroboscopic map shows a set of densely spaced points on a closed loop. This is illustrated in Figure 9.

Figure 9 Stroboscope map for $\Gamma _r$ in Figure 7, showing densely spaced points on a closed loop.

If we display the position of magnet 2 $A(\tau )$ as a third dimension of our system, then quasiperiodic trajectories $(\chi (\tau ), v(\tau ), A(\tau ))$ lie on the surface of an “invariant” torus for all time $\tau $ Footnote ¹. Figure 10 illustrates this torus for $\Gamma _r$ from Figure 7, which was obtained with initial condition $\chi _0 = 0.06$ .

Figure 10 Invariant torus for $\Gamma _r$ ( $\chi _0 = 0.06$ ) shown in Figure 7.

5.2. Solution behaviour: chaotic orbits

A second type of behaviour occurs in these forced orbits when a trajectory jumps between quasiperiodic orbits that enclose one or both centres. An illustration of such behaviour for this system is given in Figure 11. This diagram is again produced with forcing frequency $\Omega = 0.1$ . The solution was started from rest ( $v_0 = 0$ ), with initial displacement $\chi _0 = -0.08$ for magnet 1. In terms of the earlier analogy in which the unforced system was likened to a particle in a double well, the forced particle here can be thought of as transitioning between small oscillations in either well and large oscillations covering both wells when it gains enough energy. Such jumps occur intermittently – without a fixed period – so the solution is genuinely chaotic. These intermittent jumps may be seen more clearly from the time series for the trajectory, which is presented in Figure 12(a). The frequency spectrum for this solution shows a continuous distribution of frequencies, as is expected for chaos. This is illustrated in Figure 12(b).

Figure 11 Example of chaotic orbit for $\Omega = 0.1, \chi _0 = -0.08, v_0 = 0$ .

Figure 12 (a) Example time series $\chi (\tau )$ of chaotic orbit with $\Omega = 0.1$ , $\chi _0 = - 0.08$ , $v_0 = 0$ showing intermittent jumps. (b) The corresponding frequency spectrum showing a continuous set of frequencies.

We also computed Lyapunov exponents for this orbit, using the QR decomposition method [Reference Dieci, Russell and Van Vleck3] with 250 000 samples from $\tau =0$ to $\tau =5000$ . This resulted in $\lambda _{1} = 0.0239$ and $\lambda _2 = -0.0239$ . As $\lambda _1$ is a nonnegligible positive value, this implies that nearby trajectories diverge exponentially, thus providing further evidence that the solution presented in Figure 11 is indeed chaotic.

With the forcing function A represented as an extra dimension, the chaotic trajectory $(\chi (\tau ), v(\tau ), A(\tau ))$ illustrated above now lies within the volume of a strange attractor in the phase space. Iterates of the corresponding stroboscopic map now fill out a 2D region in the $(\chi , v)$ -plane, displayed in Figure 13 Footnote ¹.

Figure 13 Stroboscope map for the chaotic orbit filling out a 2D region of the $(\chi , v)$ plane at $A=2.8$ .

5.3. Solution behaviour: periodic orbits

A third type of behaviour in this forced system consists of periodic orbits, which enclose either a single centre or jump between them in a fixed pattern. Three such periodic orbits are shown in Figure 14. There are two small orbits centred on the left and right equilibria of the unforced system, and a third large orbit that encloses all three of the equilibria.

Figure 14 Three periodic orbits for the forced system with $\Omega = 0.1$ . The initial conditions used were $v_0 = 0$ and $\chi _0 = -0.0603, 0.0144, 0.2486$ for the large orbit (in blue), smaller left orbit (in red) and the smaller right orbit (in yellow), respectively.

It is instructive to consider the large-amplitude forced periodic orbit (in blue) in Figure 14 in greater detail. This solution was generated from initial condition ${\chi _0 = -0.0603}$ . The time series $\chi (\tau )$ for this solution is shown in Figure 15(a), and clearly illustrates the periodic jumps between orbits centred around the left and right centre equilibria of the unforced system. The periodic nature of this solution is underscored by its frequency spectrum presented in Figure 15(b), since this consists of spikes only at integer multiples of the driving frequency (which is also the frequency of the orbit). This is further supported by the stroboscopic map, which contains only a single point and is not shown here in the interest of space.

Figure 15 (a) Example time series $\chi (\tau )$ of periodic orbit with $\Omega = 0.1, \chi _0 = -0.0603, v_0 = 0$ showing regular jumps between centres. (b) The corresponding frequency spectrum contains only integer multiples of $\Omega = 0.1$ .

5.4. Solution behaviour: resonant orbits

A fourth behaviour type is a resonant quasiperiodic or periodic orbit. The quasiperiodic case presents itself on the stroboscope map as a number of disjoint closed loops and an example of these is presented in Figure 16, obtained for frequency $\Omega = 0.9$ and with initial conditions $\chi _0 = -0.01$ , $v_0 = 0$ . These structures are a consequence of the Kolmogorov–Arnol’d–Moser theory in resonant Hamiltonian systems (see [Reference Guckenheimer and Holmes5, pages 219–223]) and are referred to as “KAM island chains” in some modern literature such as that by Kroetz et al. [Reference Kroetz, Portela and Viana9] and Nieto et al. [Reference Nieto, Seoane, Barrio and Sanjuán13]. A chain of s islands implies the existence of a resonant periodic orbit with frequency $\Omega / s$ that traverses the island centres (centres of closed loops in Figure 16) in a given order. The corresponding trajectory for this $s = 7$ quasiperiodic resonant orbit is illustrated in Figure 17. We find that stroboscope maps, such as in Figure 16, are most convenient for identifying resonances, because the phase plane trajectory (Figure 17), the frequency spectrum and the Lyapunov exponents (here, $\lambda _{\chi } = 0.00093, \lambda _v = -0.00093$ ); otherwise, all appear similar to those of a nonresonant orbit.

Figure 16 Stroboscope map for an $s=7$ resonant orbit obtained with $\Omega = 0.9, \chi _0 = -0.01, v_0 = 0$ , showing a chain of 7 KAM islands.

Figure 17 Example of the trajectory for the $s=7$ resonant orbit obtained with $\Omega = 0.9, \chi _0 = -0.01, v_0 = 0$ .

We again extend the phase space by including as an additional axis the forcing function $A(\tau )$ . This creates an invariant torus for the $s=7$ resonant solution and the resulting structure is plotted in Figure 18. This gives a long and twisted object in the phase space, compared with a nonresonant quasiperiodic orbit, and the structure in Figure 18 intersects the $A=2.8$ plane $s=7$ times.

Figure 18 Invariant torus for the $s=7$ resonant orbit.

5.5. Damped solution behaviour

When a sufficiently small amount of damping is introduced ( $D<10^{-3}$ ), periodic and resonant periodic forced orbits of the undamped system are preserved. Damping causes nearby trajectories to be drawn towards them. To illustrate this, we return to the $s=7$ resonance from Section 5.4. We sample initial conditions within the leftmost KAM island in Figure 16 (for which $-0.01 < \chi < -0.007$ ) and apply damping $D=0.0005$ . The resulting layered stroboscope map is shown in Figure 19. Trajectories that start from initial conditions that are sufficiently close to the island centre (periodic resonance) are drawn towards the resonance as a result of the damping, while others move towards a different nonresonant periodic orbit. As damping is increased to $D=0.001$ , the $s=7$ resonance is destroyed and initial conditions in Figure 19 lead to trajectories that are all damped towards the nonresonant periodic orbit.

Figure 19 Layered stroboscope maps of initial conditions close to the $s=7$ resonance with damping. Initial conditions sufficiently close to the island centre are damped towards the resonance, but others are damped towards some nonresonant periodic orbit. Colour indicates the time along each trajectory.

Figure 20 The response amplitudes for three different approximations near the primary resonance of the lowest equilibrium point, with $H=10$ , $A_0=2.6$ and forcing amplitude $\epsilon =0.2$ . The solid blue line represents the linearized forced solution amplitude $\epsilon \mathcal {A}_1$ , the dashed black lines are the predictions of the weakly nonlinear amplitude $\mathcal {A}_W$ and the red circles are fully nonlinear amplitudes $\mathcal {A}_N$ .

5.6. Observed behaviours for varying driving frequency $\Omega $

We now explore the possible system behaviours for a range of parametric driving frequencies $\Omega $ . We return briefly to a discussion of the solution near primary resonance, examined in Section 4. For definiteness, we again consider separation distance $H=10$ and average location $A_0=2.6$ of the lower forcing magnet, illustrated in Figures 2–6. Figure 4 shows that the unforced system (3.1) possesses three different equilibria $\chi _{\text {eq}}$ at these parameter values; when subject to forcing, each of these equilibria develops its own complicated resonance structure, resulting in solution behaviour of extraordinary complexity at certain forcing frequencies $\Omega $ .

For brevity, we consider only the forced periodic behaviour of the lowest of these three equilibria at $H=10$ and $A_0=2.6$ , for which $\chi _{\text {eq}} = 0.020275$ (see Figure 2). This is subject to forcing with amplitude $\epsilon =0.2$ and over a variety of driving frequencies $\Omega $ , and the results are illustrated in Figure 20. The solid blue lines correspond to values of the linearized forcing amplitude $\epsilon \mathcal {A}_1$ calculated from (4.5) in Section 4.1. There is a primary (linearized) resonance at the frequency $\Omega _{\text {res}} \approx 0.8511$ given in (4.6) at which the linearized amplitude $\epsilon \mathcal {A}_1$ is predicted to increase without bound. The predictions of the weakly nonlinear solution of Section 4.2 are also shown in Figure 20 and are drawn using dashed black lines. These are the weakly nonlinear amplitudes $\mathcal {A}_W$ obtained from the cubic equation (4.9), and they show that nonlinear effects are responsible for bending the resonance peak to the left of the diagram, to produce three simultaneous periodic forced solutions for frequencies $\Omega $ below that of the primary resonance. This prediction is confirmed, at least qualitatively, by the numerical periodic solution of Section 4.3, for which the nonlinear solution amplitudes $\mathcal {A}_N$ at different frequencies $\Omega $ were computed using (4.11). These are shown by small red circles. Two distinct forced nonlinear solutions are computed numerically, and while a third is suspected, it was not able to be obtained. In addition, the nonlinear results also show the emergence of the first superharmonic resonance at approximately $\Omega =0.4$ and this is evident in Figure 20.

The resonance structure shown in Figure 20 pertains only to the lowest of three equilibria possible at the parameter values $H=10$ , $A_0=2.6$ . Additional resonances are possible from the other two equilibria, so that many more periodic forced orbits than shown in Figure 20 are expected. In addition, quasi-periodic and chaotic solutions are also possible, so that the nonlinear solutions in Figure 20 are likely to be embedded in an extremely complex solution environment. Their stability is therefore of interest.

We first illustrate two different nonlinear periodic solutions at the forcing frequency $\Omega =0.1$ in Figure 21(a). These two solutions correspond to the two red circles at the left-most value of $\Omega $ in Figure 20 and were obtained using the algorithm described in Section 4.3. These show the displacement $\chi $ as a function of time $\tau $ over a single period $0 < \tau < 2\pi / \Omega $ . The smaller-amplitude solution sketched with red dashed lines is the lower of the two points in Figure 20 and is approximated very closely by the linearized solution (4.4). The larger solution in Figure 21(a), sketched using a solid blue line, comes from the upper branch of the primary resonance in Figure 20 and its behaviour with time is clearly highly nonsinusoidal.

Figure 21 (a) Two solutions computed by integrating forward in time from the two periodic solutions at $\Omega =0.1$ taken from Figure 20 ( $0 < \tau < 2000$ ). (b) The smaller-amplitude solution (red line) remains periodic, but the large-amplitude solution (solid blue lines) is drawn into a chaotic strange attractor.

The stability of the two solutions shown in Figure 21(a) has been examined here simply by integrating the nonlinear equations (2.6) forward in time for $0 < \tau < 2000$ , starting with initial values taken from these two periodic forced solutions at ${\Omega =0.1}$ . The results are shown as trajectories in the $(\chi ,v)$ -plane in Figure 21(b). The smaller-amplitude solution remains periodic in time and so it appears simply as a small-amplitude ellipse in Figure 21(b) (in red). This small-amplitude periodic solution is therefore stable. By contrast, the large-amplitude solution is rapidly drawn away from the periodic configuration in Figure 21(a) and onto the strange attractor of a genuinely chaotic orbit. This is also shown in Figure 21(b) (in blue), where the initial point $\chi (0)$ is sketched with an asterisk, after which the trajectory produces the chaotic orbit indicated.

The forced periodic orbits in Figures 20–21(b) exist within a highly complex solution environment and so we conclude this presentation of results by examining how this changes with increasing forcing frequency $\Omega $ . This is done by layering stroboscope maps computed for a large sample of initial conditions (taken along the $\chi $ -axis and with $A(0) = 2.8$ ). The resulting maps give an overview of which initial conditions lead to which behaviour type. Areas densely filled with points are regions of chaos, corresponding to a slice through the strange attractor in the $(\chi ,v,A)$ phase space. Single closed loops correspond to quasiperiodic orbits (QPOs), while multiple disjoint closed loops correspond to resonant orbits. Notice that, since QPOs and resonances only fill out curves and not areas in the $(\chi , v)$ -plane, we will often fail to observe such orbits with our sample of initial conditions. Hence, it is likely that any empty “holes” in the chaotic region may contain such orbits. We also note that some loops corresponding to a QPO may not appear closed (that is, they may have a “dashed” appearance), but this is simply an artefact of insufficient integration time. The evolution of the strange attractor slice as $\Omega $ increases is explored in Figures 22–25.

Figure 22 Layered stroboscope maps for $\Omega = 0.1$ .

The layered stroboscope map for $\Omega = 0.1$ , shown in Figure 22, indicates a chaotic region containing two large holes. These are regions corresponding to QPOs and resonant orbits. The outer region surrounding the attractor also contains QPOs and resonant orbits, but these are not easily seen in Figure 22 due to the small sample of initial conditions. In Figure 23, the forcing frequency is increased to the value $\Omega = 0.9$ and, as a result, the chaotic attractor region changes shape. In particular, the holes in its structure are reduced in size.

Figure 23 Layered stroboscope maps for $\Omega = 0.9$ .

Figure 24 Layered stroboscope maps for $\Omega = 1$ , showing only quasiperiodic orbits.

Figure 25 Layered stroboscope maps for $\Omega = 20$ .

In Figure 24, the forcing frequency is increased to the value $\Omega =1$ . Surprisingly, for this case, chaotic trajectories are not observed in the small region of the $(\chi , v)$ solution space illustrated (although they do occur at much larger amplitudes). Consequently, the trajectories shown consist of QPOs only and there even do not appear to be any resonant orbits.

For $\Omega = 20$ (the highest frequency tested), the attractor appears to break apart into multiple “bands” of chaos, separated by regions of QPOs and resonant orbits. This results in the beautiful and elaborate layered stroboscopic map shown in Figure 25.