Analytic Prediction of Isolated Forced Response Curves from Spectral Submanifolds

TL;DR

This paper introduces a spectral submanifold-based analytical method to predict isolated forced responses in complex mechanical systems, surpassing numerical techniques especially in high-dimensional cases.

Contribution

It develops an explicit analytic criterion for isolated forced responses using spectral submanifold theory, applicable to high-dimensional nonlinear vibration problems.

Findings

Accurate predictions in high-dimensional systems

Explicit criterion for isolated responses

Refinable predictions without numerical solutions

Abstract

We show how spectral submanifold theory can be used to provide analytic predictions for the response of periodically forced multi-degree-of-freedom mechanical systems. These predictions include an explicit criterion for the existence of isolated forced responses that will generally be missed by numerical continuation techniques. Our analytic predictions can be refined to arbitrary precision via an algorithm that does not require the numerical solutions of the mechanical system. We illustrate all these results on low- and high-dimensional nonlinear vibration problems. We find that our SSM-based forced-response predictions remain accurate in high-dimensional systems, in which numerical continuation of the periodic response is no longer feasible.

Figures

Tables

Table 1. Table 1 : Parameter values for Example 1 .

Symbol	Value
$m$	$1\mathrm{kg}$
$c_1$	$0.03N\cdot s/m$
$c_2$	$\sqrt{3}\cdot 0.03N\cdot s/m$
$k$	$3N/m$
$\kappa$	$0.4N/m^3$
$\alpha$	$-0.6N\cdot {\left(s/m\right)}^3$
$P$	$3N$

Table 2. Table 2 : Notation used in the discretized beam example.

Symbol	Meaning (unit)
$L$	Length of beam $\left(\mathrm{mm}\right)$
$h$	Height of beam $\left(\mathrm{mm}\right)$
$b$	Width of beam $\left(\mathrm{mm}\right)$
$\rho$	Density $\left(\mathrm{kg}/{\mathrm{mm}}^3\right)$
$E$	Young’s Modulus $\left(\mathrm{kPa}\right)$
$I$	Area moment of inertia $\left({\mathrm{mm}}^4\right)$
$\kappa$	Coefficient cubic spring $\left(N/{\mathrm{mm}}^3\right)$
$\gamma$	Coefficient cubic damper $\left(N\cdot s/{\mathrm{mm}}^3\right)$
$A$	Cross-section of beam $\left({\mathrm{mm}}^2\right)$
$P$	External forcing amplitude $\left(N\right)$

Table 3. Table 3 : Geometric and material parameters for the Bernoulli beam.

Symbol	Value
$L$	$2700\mathrm{mm}$
$h$	$10\mathrm{mm}$
$b$	$10\mathrm{mm}$
$\rho$	$1780\cdot {10}^{-9}\mathrm{kg}/{\mathrm{mm}}^3$
$E$	$45\cdot {10}^6\mathrm{kPa}$
$\kappa$	$6N/{\mathrm{mm}}^3$
$\gamma$	$-0.02N\cdot s/{\mathrm{mm}}^3$
$\alpha$	$1.25\cdot {10}^{-4}$
$\beta$	$2.5\cdot {10}^{-4}$
$P$	$0.1N$

Table 4. Table 4 : Resonant terms in the non-autonomous part of the SSM for small damping and forcing frequency Ω ≈ Im ​ λ 1 Ω Im subscript 𝜆 1 \Omega\approx\text{Im}\lambda_{1} .

	$i=1$	$i=2$
$a_{i,𝐤}$	$k_1=k_2$	$k_1=k_2-2$
$b_{i,𝐤}$	$k_1=k_2+2$	$k_1=k_2$

Equations

$$\mathbf{M}\ddot{\mathbf{y}}+\mathbf{C}\dot{\mathbf{y}}+\mathbf{Ky}+\mathbf{g}(\mathbf{y},\dot{\mathbf{y}})=\epsilon \mathbf{f}(\Omega t),0\le \epsilon \ll 1,$$

$$\mathbf{g}(\mathbf{y},\dot{\mathbf{y}})=\mathcal{O}\left({\left|\mathbf{y}\right|}^2,\left|\mathbf{y}\right|\left|\dot{\mathbf{y}}\right|,{\left|\dot{\mathbf{y}}\right|}^2\right),$$

$$\dot{\mathbf{x}}$$

$$=\mathbf{Ax}+{\mathbf{G}}_{\text{p}}(\mathbf{x})+\epsilon {\mathbf{F}}_{\text{p}}(\Omega t).$$

$$\dot{\vec{x}}=\vec{A}\vec{x},$$

$$\text{Re}({\lambda }_{2n})\le \text{Re}({\lambda }_{2n-1})\le \dots \le \text{Re}({\lambda }_1)<0,$$

$$\dot{\mathbf{q}}=\underset{\vec{\Lambda }}{\underbrace{\text{diag}({\lambda }_1,{\lambda }_2\dots ,{\lambda }_{2n})}}\mathbf{q}+{\mathbf{G}}_{\text{m}}(\mathbf{q})+\epsilon {\mathbf{F}}_{\text{m}}(\mathbf{\varphi }).$$

$${\mathbb{C}}^{2n}=\mathcal{E}\oplus \mathcal{C}.$$

$$\vec{\Lambda}=\left[\begin{array}[]{cc}\vec{\Lambda}_{\mathcal{E}}&0\\ 0&\vec{\Lambda}_{\mathcal{C}}\end{array}\right],\quad\text{Spect}\left(\vec{\Lambda}_{\mathcal{E}}\right)=\left\{\lambda_{1},\lambda_{2}\right\},\quad\text{Spect}\left(\vec{\Lambda}_{\mathcal{C}}\right)=\left\{\lambda_{3},\ldots,\lambda_{2n}\right\},$$

$$\Sigma (\mathcal{E})=\text{Int}\left[\frac{{\min }_{\lambda \in \text{Spect}(\vec{\Lambda })}\text{Re}\lambda }{{\max }_{\lambda \in \text{Spect}({\vec{\Lambda }}_{\mathcal{E}})}\text{Re}\lambda }\right]\in {\mathbb{N}}^{+}.$$

$$a\text{Re}{\lambda }_1+b\text{Re}{\lambda }_2\ne \text{Re}{\lambda }_l,\forall {\lambda }_l\in \text{Spect}({\vec{\Lambda }}_{\mathcal{C}}),2\le a+b\le \Sigma (\mathcal{E}),a,b\in {\mathbb{N}}_0.$$

$$\vec{W}(\mathbf{s},\varphi ):\mathcal{U}\subset {\mathbb{C}}^2\times S^1\to {\mathbb{C}}^{2n}.$$

$$\mathbf{\Lambda W}(\mathbf{s},\varphi )+{\mathbf{G}}_{\text{m}}(\mathbf{W}(\mathbf{s},\varphi ))+\epsilon {\mathbf{F}}_{\text{m}}(\varphi )=D_{\mathbf{s}}\mathbf{W}(\mathbf{s},\varphi )\mathbf{R}(\mathbf{s},\varphi )+D_{\varphi }\mathbf{W}(\mathbf{s},\varphi )\Omega ,$$

$$\dot{\mathbf{s}}=\mathbf{R}(\mathbf{s},\varphi ).$$

$$\dot{\rho }$$

$$\dot{\psi }$$

$$a(\rho )$$

$$b(\rho )$$

$$f_1(\rho ,\Omega )$$

$$f_2(\rho ,\Omega )$$

$$g_1(\rho ,\Omega )$$

$$g_2(\rho ,\Omega )$$

$$\mathbf{F}(\mathbf{u})=\left[\begin{array}{c}{F}_1(\mathbf{u})\\ F_2(\mathbf{u})\end{array}\right]=\left[\begin{array}{c}a(\rho )+\epsilon \left(f_1(\rho ,\Omega )\cos (\psi )+f_2(\rho ,\Omega )\sin (\psi )\right)\\ (b(\rho )-\Omega )\rho +\epsilon \left(g_1(\rho ,\Omega )\cos (\psi )-g_2(\rho ,\Omega )\sin (\psi )\right)\end{array}\right]=0,$$

$$\mathbf{F}(\mathbf{u}):{\mathbb{R}}^3\to {\mathbb{R}}^2,\mathbf{u}=\left[\begin{array}{c}\rho \\ \Omega \\ \psi \end{array}\right].$$

$$G(\rho ;\Omega )=(b(\rho )-\Omega )\rho +\epsilon \left(g_1(\rho ,\Omega )\frac{1-K^{\pm }(\rho ;\Omega {)}^2}{1+K^{\pm }(\rho ;\Omega {)}^2}-g_2(\rho ,\Omega )\frac{2K^{\pm }(\rho ;\Omega )}{1+K^{\pm }(\rho ;\Omega {)}^2}\right)=0,$$

$$K^{\pm }(\rho ;\Omega )=\frac{-\epsilon f_2(\rho ,\Omega )\pm \sqrt{{\epsilon }^2\left(f_1(\rho ,\Omega {)}^2+f_2(\rho ,\Omega {)}^2\right)-a(\rho {)}^2}}{a(\rho )-\epsilon f_1(\rho ,\Omega )}.$$

$$a({\rho }_0)=0,{\partial }_{\rho }a({\rho }_0)\ne 0.$$

$$\tilde{\mathbf{F}}(\mathbf{u})=\left[\begin{array}{c}\underset{a(\rho )}{\underbrace{\text{Re}({\lambda }_1)\rho +\text{Re}({\gamma }_1){\rho }^3}}+\epsilon \left(\text{Re}(c_{1,0})\cos (\psi )+\text{Im}(c_{1,0})\sin (\psi )\right)\\ (\underset{b(\rho )}{\underbrace{\text{Im}({\lambda }_1)+\text{Im}({\gamma }_1){\rho }^2}}-\Omega )\rho +\epsilon \left(\text{Im}(c_{1,0})\cos (\psi )-\text{Re}(c_{1,0})\sin (\psi )\right)\end{array}\right]=0.$$

$$\epsilon <\frac{1}{\Vert c_{1,0}\Vert }\sqrt{\frac{4|\text{Re}({\lambda }_1){|}^3}{27\text{Re}({\gamma }_1)}}.$$

$${\epsilon }_{\text{m}}=\frac{1}{\Vert c_{1,0}\Vert }\sqrt{\frac{4|\text{Re}({\lambda }_1){|}^3}{27\text{Re}({\gamma }_1)}}.$$

$$\dot{\mathbf{x}}=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cccc}[c]0& 0& 1& 0\text{vspace}1mm\\ 0& 0& 0& 1\text{vspace}1mm\\ -\frac{2k}{m}& \frac{k}{m}& -\frac{c_1+c_2}{m}& \frac{c_2}{m}\text{vspace}1mm\\ \frac{k}{m}& -\frac{2k}{m}& \frac{c_2}{m}& -\frac{c_1+c_2}{m}\text{vspace}1mm\end{array}\right]}}\mathbf{x}+\underset{{\mathbf{G}}_{\text{p}}(\mathbf{x})}{\underbrace{\left[\begin{array}{c}[c]0\text{vspace}1mm\\ 0\text{vspace}1mm\\ -\frac{\kappa }{m}{x}_1^3-\frac{\alpha }{m}{x}_3^3\text{vspace}1mm\\ 0\text{vspace}1mm\end{array}\right]}}+\epsilon \underset{{\mathbf{F}}_{\text{p}}(\Omega t)}{\underbrace{\left[\begin{array}{c}0\text{vspace}1mm\\ 0\text{vspace}1mm\\ \frac{P}{m}\cos (\Omega t)\text{vspace}1mm\\ 0\text{vspace}1mm\end{array}\right]}},$$

$${\lambda }_{1,2}=\left(-{\zeta }_1\pm \mathrm{i}\sqrt{1-{\zeta }_1^2}\right){\omega }_1,{\zeta }_1=\frac{c_1}{2m{\omega }_1},{\omega }_1=\sqrt{\frac{k}{m}},$$

$${\lambda }_{3,4}=\left(-{\zeta }_2\pm \mathrm{i}\sqrt{1-{\zeta }_2^2}\right){\omega }_2,{\zeta }_2=\frac{c_1+2c_2}{2m{\omega }_2},{\omega }_2=\sqrt{\frac{3k}{m}},$$

$$\mathbf{T}=\left[\begin{array}{cccc}[r]1& 1& 1& 1\\ 1& 1& -1& -1\\ {\lambda }_1& {\bar{\lambda }}_1& {\lambda }_3& {\bar{\lambda }}_3\\ {\lambda }_1& {\bar{\lambda }}_1& -{\lambda }_3& -{\bar{\lambda }}_3\end{array}\right],$$

$${\mathbf{T}}^{-1}=\left[\begin{array}{cccc}[r]-\frac{{\bar{\lambda }}_1}{2({\lambda }_1-{\bar{\lambda }}_1)}& -\frac{{\bar{\lambda }}_1}{2({\lambda }_1-{\bar{\lambda }}_1)}& \frac{1}{2({\lambda }_1-{\bar{\lambda }}_1)}& \frac{1}{2({\lambda }_1-{\bar{\lambda }}_1)}\text{vspace}1mm\\ \frac{{\lambda }_1}{2({\lambda }_1-{\bar{\lambda }}_1)}& \frac{{\lambda }_1}{2({\lambda }_1-{\bar{\lambda }}_1)}& -\frac{1}{2({\lambda }_1-{\bar{\lambda }}_1)}& -\frac{1}{2({\lambda }_1-{\bar{\lambda }}_1)}\text{vspace}1mm\\ -\frac{{\bar{\lambda }}_3}{2({\lambda }_3-{\bar{\lambda }}_3)}& \frac{{\bar{\lambda }}_3}{2({\lambda }_3-{\bar{\lambda }}_3)}& \frac{1}{2({\lambda }_3-{\bar{\lambda }}_3)}& -\frac{1}{2({\lambda }_3-{\bar{\lambda }}_3)}\text{vspace}1mm\\ \frac{{\lambda }_3}{2({\lambda }_3-{\bar{\lambda }}_3)}& -\frac{{\lambda }_3}{2({\lambda }_3-{\bar{\lambda }}_3)}& -\frac{1}{2({\lambda }_3-{\bar{\lambda }}_3)}& \frac{1}{2({\lambda }_3-{\bar{\lambda }}_3)}\end{array}\right].$$

Full text

1 Introduction

For an $n$-degree-of-freedom, periodically forced, nonlinear mechanical system, the forced response curve (FRC) gives the amplitude of the periodic response of the system as a function of the frequency of the periodic forcing. The FRC may contain isolated branches of periodic solutions, also known as isolas, that are detached from the main FRC. A small change in the forcing amplitude might result in the merger of the isola with the main branch of the FRC (cf. Detroux et al. Detroux2018 and Noël et al. noel2015isolated ), which can lead to an unexpected and significant increase in the response amplitude.

The existence of isolated branches of periodic solutions in the frequency response of nonlinear oscillatory systems has been known since the 1950s abramson1955response . For an extensive review of the subject, we refer the reader to Habib et al. Habib2017 . It is broadly agreed that the identification of isolas is difficult, because numerical continuation techniques are generally initiated on a non-isolated solution branch and will therefore miss any isolated branch. Similarly, a frequency sweep of the full system will generally not capture an isolated response unless the sweep is initialized on one.

The detection of isolas and the prediction of their behavior under changing system parameters can be critical in practice because, their merger with the main FRC may lead to a dramatic shift in the resonance frequency and response amplitude. Habib et al. Habib2017 use singularity theory in combination with averaging for the prediction and identification of isolas in a specific single-degree-of-freedom mechanical system with nonlinear damping. The averaging method they use (cf. Sanders et al. sanders2007averaging ), however, requires both the forcing amplitude and the nonlinear damping coefficients to be small. Hill et al. hill2016analytical use a second-order normal form technique to obtain analytical expressions for the autonomous conservative backbone curves (i.e. amplitude-frequency plots of nonlinear periodic orbits) of a specific two-degree-of-freedom mechanical system. They give leading-order criteria for the intersection of this backbone curve with the forced response curve and postulate this location to be a potential starting point for an isola, which is to be constructed numerically in a separate effort. This procedure also relies on the smallness of the nonlinear and damping coefficients, as well as on the absence of quadratic nonlinearities.

In summary, while the significance of isolas is broadly recognized, their existence has only been studied in specific, low-dimensional examples under restrictions on the nonlinearities. A conclusive analytical criterion for predicting isolas in multi-degree of freedom systems without costly numerical simulations, therefore, has been unavailable.

In this work, we seek to fill this gap by developing a generally applicable methodology for the prediction of isolas in multi-degree-of-freedom, forced mechanical systems. Our approach is based on the mathematically rigorous theory of spectral submanifolds (SSMs) that are the unique, smoothest, nonlinear continuations of spectral subspaces of the linearized, unforced limit of a mechanical system (cf. Haller and Ponsioen Haller2016 ). The reduced dynamics on a two-dimensional SSM serves as an exact, single-degree-of-freedom reduced-order model that can be constructed for each vibration mode of the full nonlinear system (cf. ponsioen2018automated ; Jain2018 ; Breunung2017 ; Szalai2017 ; kogelbauer2018rigorous ).

By construction, these rigorously, simplified two-dimensional reduced models will capture all isolas that are remnants of periodic orbit families of the conservative limit of the system. As we show for a cubic-order approximation, the reduced SSM dynamics gives a closed form first-order prediction for isolas that can even be calculated by hand in simple examples. Higher-order refinements to this analytic formula can be recursively constructed and have been implemented in the publicly available MATLAB script SSMtool111SSMtool is available at: www.georgehaller.com. We show the use of the analytic formula as well as its numerical refinements on simple and more complicated examples.

2 System set-up

We consider $n$-degree-of-freedom, periodically forced mechanical systems of the form

[TABLE]

where $\mathbf{y}\in\mathbb{R}^{n}$ is the generalized position vector; $\mathbf{M}=\mathbf{M}^{T}\in\mathbb{R}^{n\times n}$ is the positive definite mass matrix; $\mathbf{C}=\mathbf{C}^{T}\in\mathbb{R}^{n\times n}$ is the damping matrix; $\mathbf{K}=\mathbf{K}^{T}\in\mathbb{R}^{n\times n}$ is the stiffness matrix and $\mathbf{g}(\mathbf{y},\dot{\mathbf{y}})$ denotes all the nonlinear terms in the system. These nonlinearities are assumed to be analytic for simplicity. The external forcing $\varepsilon\mathbf{f}(\Omega t)$ does not depend on the positions and velocities.

System (1) can be transformed into a set of $2n$ first-order ordinary differential equations by introducing the change of variables $\mathbf{x}_{1}=\mathbf{y}$, $\mathbf{x}_{2}=\dot{\mathbf{y}}$, with $\mathbf{x}=(\mathbf{x}_{1},\mathbf{x}_{2})\in\mathbb{R}^{2n}$, which gives

[TABLE]

The transformed first-order system (9) has a fixed point at $\mathbf{x}=\mathbf{0}$ when the system is unforced ($\varepsilon=0$); $\mathbf{A}\in\mathbb{R}^{2n\times 2n}$ is a constant matrix and $\mathbf{G}_{\text{p}}(\mathbf{x})$ is an analytic function containing all the nonlinearities. Note that $\mathbf{M}^{-1}$ is well-defined because $\mathbf{M}$ is assumed positive definite.

The linearized part of system (9) is

[TABLE]

where the matrix $\vec{A}$ has $2n$ eigenvalues $\lambda_{k}\in\mathbb{C}$ for $k=1,\ldots,2n$. Counting multiplicities, we sort these eigenvalues based on their real parts in the decreasing order

[TABLE]

assuming that the real part of each eigenvalue is less than zero and hence the fixed point of Eq. (10) is asymptotically stable. We further assume that the constant matrix $\vec{A}$ is semisimple, i.e., the algebraic multiplicity of each $\lambda_{k}$ is equal to its geometric multiplicity: $\text{alg}(\lambda_{k})=\text{geo}(\lambda_{k})$. We can, therefore, identify $2n$ linearly independent eigenvectors $\vec{v}_{k}\in\mathbb{C}^{2n}$, with $k=1,\ldots,2n$, each spanning a real eigenspace $E_{k}\subset\mathbb{R}^{2n}$ with $\text{dim}(E_{k})=2\times\text{alg}(\lambda_{k})$ in case $\text{Im}(\lambda_{k})\neq 0$, or $\text{dim}(E_{k})=\text{alg}(\lambda_{k})$ in case $\text{Im}(\lambda_{k})=0$.

3 Non-autonomous spectral submanifolds and their reduced dynamics

Because $\mathbf{A}$ is semisimple, the linear part of system (9) is diagonalized by a linear change of coordinates $\mathbf{x}=\mathbf{T}\mathbf{q}$, with $\mathbf{T}=\left[\mathbf{v}_{1},\mathbf{v}_{2},\ldots,\mathbf{v}_{2n}\right]\in\mathbb{C}^{2n\times 2n}$ and $\mathbf{q}\in\mathbb{C}^{2n}$, yielding

[TABLE]

We now consider the two-dimensional modal subspace $\mathcal{E}=\text{span}\left\{\vec{v}_{1},\vec{v}_{2}\right\}\subset\mathbb{C}^{2n}$ with $\vec{v}_{2}=\bar{\vec{v}}_{1}.$ The remaining linearly independent eigenvectors $\vec{v}_{3},\ldots,\vec{v}_{2n}$ span a complex subspace $\mathcal{C}\subset\mathbb{C}^{2n}$ such that the full phase space of (12) can be expressed as the direct sum

[TABLE]

We write the diagonal matrix $\vec{\Lambda}$ as

[TABLE]

with $\vec{\Lambda}_{\mathcal{E}}=\text{diag}(\lambda_{1},\lambda_{2})$ and $\vec{\Lambda}_{\mathcal{C}}=\text{diag}(\lambda_{3},\ldots,\lambda_{2n})$.

Following Haller and Ponsioen Haller2016 , we now define a non-autonomous spectral submanifold (SSM), $\mathcal{W}(\mathcal{E})$, corresponding to the spectral subspace $\mathcal{E}$ of $\vec{\Lambda}$ as a two-dimensional invariant manifold of the dynamical system (12) that:

• (i)

Perturbs smoothly from $\mathcal{E}$ at the trivial fixed point $\mathbf{q}=\mathbf{0}$ under the addition of the $\mathcal{O}(\varepsilon)$ terms in Eq. (12).

• (ii)

Is strictly smoother than any other invariant manifold with the same properties.

We also define the absolute spectral quotient $\Sigma(\mathcal{E})$ of $\mathcal{E}$ as the positive integer

[TABLE]

Additionally, we introduce the non-resonance conditions

[TABLE]

We now restate the following result from Haller and Ponsioen Haller2016 , for the existence of an SMM in system (12).

Theorem 3.1

Under the non-resonance conditions (16), the following hold for system (12):

• (i)

There exists a unique two-dimensional, time-periodic, analytic SSM, $\mathcal{W}(\mathcal{E})$.

• (ii)

$\mathcal{W}(\mathcal{E})$* can be viewed as an embedding of an open set $\mathcal{U}$ into the phase space of system (12) via the map*

[TABLE]

We can approximate $\mathbf{W}(\mathbf{s},\phi)$ in a neighborhood of the origin using a Taylor expansion in the parameterization coordinates $\mathbf{s}$, having coefficients that depend periodically on the phase variable $\phi$.

• (iii)

There exists a polynomial function $\vec{R}(\mathbf{s},\phi):\mathcal{U}\rightarrow\mathcal{U}$ satisfying the invariance relationship

[TABLE]

such that the reduced dynamics on the SSM can be expressed as

[TABLE]

Proof

: We have simply restated the main theorem by Haller and Ponsioen Haller2016 , which is based on the more abstract results of Cabré et al. Cabre2003 ; Cabre2003b ; Cabre2005 for mappings on Banach spaces. ∎

We give an illustration of a time-periodic SSM in Fig. 1. We have assumed a case in which the SSM has three limit cycles for a given forcing frequency, with two of these limit cycles contained in an isola. The SSM approach can be viewed as a refinement and extension of the seminal work of Shaw and Pierre Shaw1993 , who envision nonlinear normal modes as invariant manifolds that are locally graphs over two-dimensional modal subspaces of the linearized system. We explain in detail how to construct time-periodic SSMs in Appendix A.

Our next result concerns the dynamics on the SSM described in Theorem 3.1

Theorem 3.2

The dynamics on the two-dimensional SSM given in Theorem 3.1 can approximately be written in polar coordinates $(\rho,\psi)$ as

[TABLE]

where

[TABLE]

with $2M+1$ denoting the order of the expansion.

Proof

: We derive this result in Appendix B. ∎

In the unforced limit ($\varepsilon=0$), the reduced system (20)-(21) can have fixed points but no nontrivial periodic orbits. This is because (20) decouples from (21), representing a one-dimensional ordinary differential equation that cannot have non-constant periodic solutions. By construction, the trivial fixed point of (20)-(21) is asymptotically stable and will persist for $\varepsilon>0$. These persisting fixed points satisfy the system of equations

[TABLE]

where

[TABLE]

As we show in Appendix C, under appropriate non-degeneracy conditions, the zeros of (28) form a one-dimensional manifold, which, after a projection onto the amplitude-frequency space, will represent the FRC. The stability of these fixed points (which correspond to periodic solutions of the full mechanical system) is determined by the real parts of the eigenvalues of the Jacobian of $\mathbf{F}(\mathbf{u})$.

Theorem 3.3

The amplitude $\rho$ of the $T$-periodic orbits of the reduced dynamics (20)-(21) are given by the zeros of the function

[TABLE]

where

[TABLE]

Proof

We derive this result in Appendix D.∎

The zero-level set of Eq. (29) yields the forced-response curve in the $(\Omega,\rho)$-space. This curve will consist of two segments obtained from $K^{+}$ and $K^{-}$ in Eq. (30). The two segments meet exactly at the point where the square root term in the definition of $K^{\pm}(\rho;\Omega)$ is equal to zero. We sketch this for a damped, nonlinear, periodically forced mechanical system in Fig. 2.

Because the isolated branches of periodic solutions are also a part of the FRC, the zero-level set of $G(\rho;\Omega)$ can predict isolas as well. In contrast, detecting isolas by numerical continuation requires one to start on the isola and hence assumes a priori knowledge of an isolated branch of periodic solutions. We will discuss this in section 4.

4 Analytic criterion for isolas

We will now give an analytic criterion for the emergence of isolas in terms of the function $a(\rho)$ defined in equation (22). An essential question is how approximate zeros obtained from finite-order Taylor series expansions persist as $M\rightarrow\infty$. Jentzsch jentzsch1916untersuchungen proved that in the limit of the order of the Taylor series expansion going to infinity, the non-persistent spurious zeros come arbitrarily close to the boundary of the domain of convergence. Hurwitz hurwitz1888ueber showed that in the same limit, the uniform convergence of the Taylor series polynomial leads to a good approximation of the genuine zeros inside the circle of convergence. Christiansen and MadsenChristiansen2006 numerically verified this behavior on several examples. We will call such a zero $\rho_{0}$ of $a(\rho)$, non-spurious if it converges to a genuine zero of $a(\rho)$ in the limit of $M\rightarrow\infty$.

Theorem 4.1

Assume that $\rho_{0}\neq 0$ is a non-spurious transverse zero of $a(\rho)$, i.e.,

[TABLE]

Then, for $\varepsilon>0$ small enough, system (28) has an isola that perturbs from the unforced damped backbone curve $\Omega=b(\rho)$ near the amplitude value $\rho_{0}$.

Proof

We derive this result in Appendix E. ∎

In order to verify if a non-trivial zero of the Taylor series expansion of $a(\rho)$ is also non-spurious, we compute the (generally complex) zeros of the function $a(\rho)$ for increasing order of approximation $M$. As we discussed before Theorem 4.1, spurious zeros will converge to the circle defining the radius of convergence of $a(\rho)$, whereas non-spurious zeros stay bounded away from that circle and hence converge to the genuine zeros of $a(\rho)$.

In Fig. 3, we sketch qualitatively the statement of Theorem 4.1: a non-spurious, transverse zero $\rho_{0}$ of $a(\rho)$ indicates a nearby isola.

4.1 Leading-order analytic formula for isolas

For higher-order approximations of $a(\rho)$, we can determine the roots of $a(\rho)$ numerically. Restricting ourselves to a $3^{\text{rd}}$-order approximation of the SSM, we can, however, extract even an analytic criterion for the existence of an isola. Following the work of Breunung and Haller Breunung2017 , we truncate the parameterization $\mathbf{W}(\mathbf{s},\phi)$ and the reduced dynamics $\mathbf{R}(\mathbf{s},\phi)$ at $\mathcal{O}(\varepsilon|\mathbf{s}|,\varepsilon^{2})$, which they justify by introducing the scaling $\mathbf{s}\rightarrow\varepsilon^{\frac{1}{4}}\mathbf{s}$, such that the zero problem of the reduced system can be written as

[TABLE]

Here we have $f_{1}=g_{2}=\text{Re}(c_{1,\mathbf{0}})$ and $f_{2}=g_{1}=\text{Im}(c_{1,\mathbf{0}})$. We now show in Theorem 4.2 that this approximation gives an analytically computable condition for the existence of an isola.

Theorem 4.2

Assume that for system (32), $\textup{Re}(\gamma_{1})>0$ is satisfied and the cubic-order zero $\rho_{1}=\sqrt{|\textup{Re}(\lambda_{1})|/\textup{Re}(\gamma_{1})}$ of $a(\rho)$ is non-spurious. Then the following holds:

• (i)

For $\varepsilon>0$ small enough, an isola of the type described in Theorem 4.1 exists near the point $(\Omega,\rho)=(b(\rho_{0}),\rho_{0})$ of the damped backbone curve.

• (ii)

The isola will be disconnected from the main FRC for $\varepsilon>0$ values satisfying

[TABLE]

• (iii)

The isola will merge with the main FRC approximately at the $\varepsilon$ value

[TABLE]

Proof

We derive this result in Appendix F.∎

5 Numerical Examples

5.1 The modified Shaw–Pierre example

As a typical benchmark, we first consider a modified version of the example of Shaw and Pierre Shaw1994 , in which an additional cubic nonlinear damper is added, as in the single degree-of-freedom example of Habib et al. Habib2017 . The equations of motion of our two-degree-of-freedom system in first-order form are given by

[TABLE]

where $\mathbf{x}=[x_{1},x_{2},x_{3},x_{4}]^{\top}=[y_{1},y_{2},\dot{y}_{1},\dot{y}_{2}]^{\top}$. The matrix $\mathbf{A}$ has the eigenvalue pairs

[TABLE]

assume that both modes are underdamped, i.e., $0<\zeta_{1}<1$ and $0<\zeta_{2}<1$. The matrix $\mathbf{T}$ that transforms our system to complex modal coordinates is composed of the eigenvectors of our system, i.e.,

[TABLE]

with the inverse of $\mathbf{T}$ given by

[TABLE]

In this example, we can compute the cubic coefficient of $a(\rho)$ explicitly. Specifically, we have

[TABLE]

Therefore, for $\alpha<0$, the reduced dynamics on the third-order autonomous SSM will have a nontrivial zero. If, additionally, this zero is non-spurious, then Theorem 4.2 guarantees the existence of an isola. Using Eq. (34), the isola will merge with the main FRC for

[TABLE]

We verify this analytic prediction numerically in Example 1 below.

Example 1

We choose the parameter values listed in Table 1 and compute the forced response curve for system (35). Note that for this choice of damping parameters, the non-resonance conditions (16) are satisfied.

Plugging in the parameter values of Table 1 into Eq. (42), we observe that the third-order coefficient of the autonomous part of the SSM is

[TABLE]

We now numerically verify that this transverse zero is non-spurious by computing the (complex) roots of $a(\rho)$ for an increasing order $M$ of expansion in formula (22) using SSMtool. In Fig. 4, we show these roots, up to $50^{\text{th}}$ order, with lighter colors indicating higher orders of approximation. Eq. (44) and Fig. 4 allow us to conclude from statement (i) of Theorem 4.2 the existence of an isola near the amplitude value $|\rho_{1}^{\pm}|$, where $\rho_{1}^{\pm}$ are the two nontrivial, non-spurious zeros of $a(\rho)$ seen in Fig. 4.

By Eq. (43), the isola will merge with the main branch of the forced response curve approximately at the parameter value

[TABLE]

We now verify this analytic prediction for the isola merger numerically. In Fig. 5, we show in red the leading-order forced response curves for $\varepsilon=0.0027$ and $\varepsilon=0.0029$. Also shown in black are the forced response curves of the full system obtained via the periodic-orbit toolbox of coco Dankowicz2013 . We conclude that the FRC obtained from our two-dimensional, SSM-reduced system perfectly predicts the behavior of the full system.

Example 2

In this example, we add a quintic nonlinear damper to system (35), which yields the modified equations of motion in first-order form

[TABLE]

We again use the parameter values in Table 1 and additionally select the quintic damping coefficient $\beta=1.2\,\mathrm{N\cdot(s/m)^{5}}$. We use SSMtool to calculate the functions included in the reduced dynamics (20)-(21) up to $5^{\text{th}}$ order in $\rho$. The function $a(\rho)$ now has two positive, non-spurious zeros located at $\rho^{+}_{1}=0.13$ and $\rho^{+}_{2}=0.17$. Therefore, Theorem 4.1 implies the existence of two separate isolas bifurcating from the damped backbone curve under periodic forcing.

We show the extracted forced response curves for three different values of $\varepsilon$ in Fig. 7. As we have predicted above, two isolas are born out of the non-trivial transverse zeros of $a(\rho)$ along the autonomous backbone curve. The isola with lower amplitudes is unstable, whereas the isola with higher amplitudes is partially stable. If we increase the forcing amplitude $\varepsilon$, the two isolas merge. Increasing $\varepsilon$ further will make the merged isolas merge with the lower FRC branch. The branches of the forced response curve, extracted using SSM theory, are again verified using the periodic-orbit toolbox of coco Dankowicz2013 . In order to initialize the continuation algorithm, we integrate the full system to provide an initial solution guess that is used to start the continuation process. For higher amplitude values, our $5^{\text{th}}$-order approximation slightly deviates from the numerical continuation results, as expected. For lower amplitudes, however, our SSM-based prediction perfectly matches the numerical result.

5.2 A discretized, forced Bernoulli beam with a cubic spring

Inspired by a similar example of Peeters et al. Peeters2009 , we now construct a reduced-order model for a discretized, cantilevered Bernoulli beam with a cubic spring and damper attached to the free end of the beam. We obtain the reduced model by computing the dynamics on the slowest, two-dimensional, time-periodic SSM of the system. This example is complex enough so that numerical continuation for obtaining the forced response is no longer feasible numerically, whereas our SSM-based predictions can still be carried out and will be seen to be accurate.

We consider a square beam of length $L$, with cross-section $A$, situated in a Cartesian coordinate system of $(x,y,z)$ and basis $(\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})$. We list the relevant beam parameters in Table 2.

The line of points coinciding with the $x$-axis is called the beam’s neutral axis. The Bernoulli hypothesis states that initially straight material lines, normal to the neutral axis, remain (a) straight and (b) inextensible, and (c) rotate as rigid lines to remain perpendicular to the beam’s neutral axis after deformation. These kinematic assumptions are satisfied the displacement field,

[TABLE]

where $(u_{x},u_{y},u_{z})$ are the components of the displacement field $\mathbf{u}(x,y,z,t)$ of a material point located at $(x,y,z)$. The transverse displacement of a material point with initial coordinates on the beam’s neutral axis at $z=0$ is denoted by $w(x)$. The rotation angle of a transverse normal line about the $y$-axis is given by $-\partial_{x}w(x)$. Using the Green-Lagrange strain tensor, we can express the relevant strains as

[TABLE]

We assume an isotropic, linearly elastic constitutive relation between the stresses and strains, i.e.

[TABLE]

which finally leads to the equation of motion of the beam

[TABLE]

We assume that the thickness of the beam is small compared to its length, i.e., $h\ll L$, and hence we can neglect the mixed partial derivative term in Eq. (52) (cf. Reddy and Mahaffey Reddy2013 ).

After discretization of (52), we obtain the set of ordinary differential equations

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{2m}$, and $m$ is the number of elements used in the discretization. Each node of the beam has two coordinates related to the transverse displacement $w(x)$ and the rotation angle $-\partial_{x}w(x)$ of the cross section. We assume structural damping by considering the damping matrix

[TABLE]

with parameters $\alpha$ and $\beta$. We apply cosinusoidal external forcing on the transverse displacement coordinate at the free end of the beam with forcing frequency $\Omega$ and forcing amplitude $\varepsilon P$. Additionally, we add a cubic spring and damper along this coordinate, with coefficients $\kappa$ and $\gamma$, respectively. As a result, the equations of motion of the beam can be written as

[TABLE]

We illustrate the kinematics, the forcing and the cubic spring and damper in Fig. 8.

We select $m=25$ for the number of discretized elements, which gives $\mathbf{x}\in\mathbb{R}^{50}$, resulting in a 100-dimensional phase space. We list the geometric and material parameter values in Table 3.

For these parameter values, the eigenvalues corresponding to the slowest eigenspace are

[TABLE]

As earlier, introducing the scaling $\mathbf{s}\rightarrow\varepsilon^{\frac{1}{4}}\mathbf{s}$, we obtain the approximations

[TABLE]

The function $a(\rho)$ in Eq. (57) has a non-trivial, transverse, positive zero at $\rho_{1}^{+}=0.413$. Fig. 9 shows this zero to be non-spurious. Therefore, by Theorem 4.2, an isola will perturb from the point $(\Omega=b(\rho_{1}^{+}),\rho_{1}^{+})$ of the autonomous backbone curve. Also by Theorem 4.2, the isola will merge with the main branch of the FRC approximately for

[TABLE]

To verify our predicted merger amplitude in (60), we perform a discrete numerical sweep of the full system, in which we force the system at different forcing frequencies and plot the resulting maximum absolute value of the transverse displacement of the tip of the beam, while keeping the forcing amplitude fixed (see Fig. 10).

While numerical continuation remains a powerful tool in verifying our analytic predictions in lower dimensions, its use becomes unfeasible in higher dimensions. For this reason, Fig. 10 only shows a discrete set of periodic responses computed from a point-wise, long-term numerical integration leading to a steady state, as opposed to a continuous FRC obtained from numerical continuation.

6 Conclusion

We have used the exact reduced dynamics on two-dimensional time-periodic spectral submanifolds (SSMs) to extract forced-response curves (FRCs) and predict isolas in arbitrary multi-degree-of-freedom mechanical systems without performing costly numerical simulations. We showed that for a cubic-order approximation, the reduced dynamics on the SSM gives an analytic prediction for the isolas, valid for any mode of a multi-degree-of-freedom oscillatory system. For simple examples, these predictions can explicitly be expressed as functions of the system parameters. Our lower-order predictions can be refined to higher-orders using the publicly available MATLAB script SSMtool222SSMtool is available at: www.georgehaller.com.

For mechanical systems of high degrees of freedom, numerical continuation techniques for forced response curves become computationally expensive. Instead, using the non-autonomous SSM and the corresponding reduced dynamics on the SSM, we are able to approximate all possible FRCs for different forcing amplitudes, as our expressions depend symbolically on the forcing amplitude $\varepsilon$. An additional advantage of the results derived here is that the isolas are uncovered by the transverse intersection of the zero-level sets of our two reduced equations. The isolas will generally be missed by numerical continuation techniques, which require starting on an isolated solution branch. As we have shown, our predictions for the main FRC branches, as well as for isolas, remain valid and computable in high dimensional problems in which numerical continuation is no longer a viable alternative for constructing these curves. Using the general results of Haller and Ponsioen Haller2016 , one can extend the periodic approach to detect quasi-periodic responses and isolas under quasi-periodic forcing.

Appendix A Coefficient equations for the the non-autonomous SSMs

As the SSM perturbs smoothly with respect to the parameter $\varepsilon$, we can expand $\mathbf{W}(\mathbf{s},\phi)$ and $\mathbf{R}(\mathbf{s},\phi)$ in $\varepsilon$, obtaining the expressions

[TABLE]

The forcing terms in system (12) are of order $\varepsilon$, therefore, the leading-order terms in Eqs. (61)-(62) are not functions of $\phi$. We now substitute Eqs. (61)-(62) into the invariance Eq. (18) and collect terms of equal order in $\varepsilon$. Given that $\mathbf{G}_{\text{m}}(\mathbf{q})=\mathcal{O}(\left|\mathbf{q}\right|^{2})$, we can Taylor expand $\mathbf{G}_{\text{m}}(\mathbf{W}(\mathbf{s},\phi))$ around $\varepsilon=0$, to obtain

[TABLE]

Collecting terms of $\mathcal{O}(1)$, we obtain

[TABLE]

and, subsequently, collecting terms of $\mathcal{O}(\varepsilon)$ we find that

[TABLE]

Eq. (64) is derived and solved in Ponsioen et al. ponsioen2018automated to compute autonomous SSMs using the Kronecker product. The construction of the autonomous SSMs has been automated and cast into the MATLAB-based computational tool SSMtool. In the following sections, the $\mathcal{O}(\varepsilon)$ terms of the SSM will simply be referred to as the non-autonomous part of the SSM. Note that in order to solve Eq. (65), Eq. (64) needs to be solved first for the autonomous SSM and for the reduced dynamics, since the respective terms obtained from Eq. (64) are needed in Eq. (65).

The autonomous and non-autonomous parts of the SSM and the reduced dynamics, which have previously been derived from an expansion in $\varepsilon$, are in turn Taylor expanded in the parameterization coordinates $\mathbf{s}$, which we explicitly express as

[TABLE]

Here we have made us of the multi-index notation $\mathbf{m}\in\mathbb{N}_{0}^{2}$. We now assume that Eq. (64) has already been solved for $\mathbf{W}_{0}(\mathbf{s})$ and $\mathbf{R}_{0}(\mathbf{s})$ and therefore their coefficients $W_{i,\mathbf{m}}^{0}$ and $R_{i,\mathbf{m}}^{0}$ are known.

Theorem A.1

For $\phi\in S^{1}$, the coefficient equation related to the $\textbf{k}^{\text{th}}$-power term of the $i^{\text{th}}$ row of the non-autonomous invariance Eq. (65), is equal to

[TABLE]

where $P_{i,\mathbf{k}}(\phi)$ can be written as

[TABLE]

Proof

Assuming that $\phi\in S^{1}$, we obtain that for the $i^{\text{th}}$ row, the $\mathbf{k}^{\text{th}}$-power terms on the right-hand side of Eq. (65) can be expressed as

[TABLE]

where we have made use of the multi-index notation

[TABLE]

with $\mathbf{e}_{j}$ denoting a unit vector.

The $\mathbf{k}^{\text{th}}$-power terms on the left-hand side of the $i^{\text{th}}$ row of Eq. (65) can be written as

[TABLE]

Therefore, the coefficient equation related to the $\textbf{k}^{\text{th}}$-power term of the $i^{\text{th}}$ row of the non-autonomous invariance Eq. (65) is

[TABLE]

where

[TABLE]

which concludes the proof of Theorem A.1.∎

A.1 Solving the non-autonomous invariance equation for $|\mathbf{k}|=0$

For $|\mathbf{k}|=\mathbf{0}$, Eq. (70) reduces to

[TABLE]

which is equivalent to the $0^{\text{th}}$-order expansion of Breunung and Haller Breunung2017 . Assuming that the forcing term $F_{i,\mathbf{0}}(\phi)$ can be written as

[TABLE]

we express $W_{i,\mathbf{0}}^{1}(\phi)$ and $R_{i,\mathbf{0}}^{1}(\phi)$ in the following form

[TABLE]

As the SSM is constructed over a two-dimensional spectral subpace $\mathcal{E}$, corresponding to the eigenvalues $\lambda_{1}=\text{Re}\lambda_{1}+\mathrm{i}\text{Im}\lambda_{1}$ and $\lambda_{2}=\bar{\lambda}_{1}=\text{Re}\lambda_{1}-\mathrm{i}\text{Im}\lambda_{1}$, we can write the solution of Eq. (81) as

[TABLE]

Our main goal is to obtain the forced response curve around the frequency $\text{Im}\lambda_{1}$, which corresponds to the spectral subspace $\mathcal{E}$. For lightly damped systems where $\text{Re}\lambda_{1}$ is small, we obtain small denominators in Eq. (84) if the forcing frequency $\Omega$ is approximately equal to $\text{Im}\lambda_{1}$. We, therefore, intend to remove this near-resonance by setting

[TABLE]

A.2 Solving the non-autonomous invariance equation for $|\mathbf{k}|>0$

For $|\mathbf{k}|>0$, the solution to the non-autonomous invariance Eq. (70) is given by

[TABLE]

where we have let $P_{i,\mathbf{k}}=\alpha_{i,\mathbf{k}}\mathrm{e}^{\mathrm{i}\phi}+\beta_{i,\mathbf{k}}\mathrm{e}^{-\mathrm{i}\phi}$. Using the same reasoning as in section A.1, we want to choose $c_{i,\mathbf{k}}$ and $d_{i,\mathbf{k}}$ in Eq. (86) in such a way that the coefficients $a_{i,\mathbf{k}}$ and $b_{i,\mathbf{k}}$, which become large when the damping is small and the forcing frequency $\Omega$ is close to $\text{Im}\lambda_{1}$, are completely removed. These terms are the so-called resonant terms. In Table 4, we list the resonant terms $a_{i,\mathbf{k}}$ and $b_{i,\mathbf{k}}$ for small damping and for $\Omega\approx\text{Im}\lambda_{1}$.

The terms listed in Table 4 will be removed from the expressions of $\mathbf{W}_{1}(\mathbf{s},\phi)$ and included into the non-autonomous part of the reduced dynamics $\mathbf{R}_{1}(\mathbf{s},\phi)$, in order to avoid small denominators in $\mathbf{W}_{1}(\mathbf{s},\phi)$.

Appendix B Proof of Theorem 3.2

The $\mathcal{O}(\varepsilon)$ approximation of the reduced dynamics for $\mathbf{s}$ can be written as

[TABLE]

where the first row of Eq. (87) takes the form

[TABLE]

and the second row of Eq. (87) is, by construction, the complex conjugate of the first row. Introducing a change to polar coordinates, $s_{1}=\rho\mathrm{e}^{\mathrm{i}\theta}$, $s_{2}=\bar{s}_{1}=\rho\mathrm{e}^{-\mathrm{i}\theta}$, dividing by $\mathrm{e}^{\mathrm{i}\theta}$ and introducing the new phase coordinate $\psi=\theta-\phi$, we obtain

[TABLE]

We obtain the result listed in Theorem 3.2 by splitting Eq. (89) into its real and imaginary part. ∎

Appendix C Extracting the forced response curve

For convenience, we restate our zero problem (28),

[TABLE]

where

[TABLE]

If there exists a regular point $\mathbf{p}=(\rho,\Omega,\psi)$, such that $\mathbf{F}(\mathbf{p})=\mathbf{0}$ in (90) and the Jacobian of $\mathbf{F}$ evaluated at $\mathbf{p}$ is surjective, then by the implicit function theorem, locally there exists a one-dimensional submanifold of $\mathbb{R}^{3}$ which will represent the forced response curve when projected onto the $(\Omega,\rho)$-space. Equivalently, the zero-level sets of $F_{1}(\mathbf{u})$ and $F_{2}(\mathbf{u})$ in (90), which we will denote by $\mathcal{M}_{1}^{\mathbf{p}}$ and $\mathcal{M}_{2}^{\mathbf{p}}$, will be two two-dimensional submanifolds in the $(\rho,\psi,\Omega)$-space that, locally around $\mathbf{p}$, intersect each other transversely, yielding the forced response curve. We illustrate this concept in Fig. 11, which is a typical picture for a damped non-linear periodically forced mechanical system with a hardening nonlinearity.

Appendix D Proof of Theorem 3.3

Let $\mathbf{u}_{0}$ be a regular point of the map $\mathbf{F}(\mathbf{u})$ (28) such that $\mathbf{F}(\mathbf{u}_{0})=\mathbf{0}$, and $D_{\mathbf{u}}\mathbf{F}(\mathbf{u}_{0})$ is surjective. Then, by the implicit function theorem, locally there exists a one-dimensional submanifold of $\mathbb{R}^{3}$ around $\mathbf{u}_{0}$. We express $\psi$ as a function of $\rho$ and $\Omega$, using the tangent half-angle substitution and the trigonometric identities,

[TABLE]

Setting Eq. (20) equal to zero and substituting the identities into Eq. (20), we obtain a quadratic equation in $K$,

[TABLE]

which has the solution,

[TABLE]

Substituting Eq. (93), together with the trigonometric identities in (91), into Eq. (21), we obtain the result stated in Theorem 3.3. ∎

Appendix E Proof of Theorem 4.1

We now consider $\varepsilon$ to be a variable in our zero problem (28), i.e.

[TABLE]

If there exists a non-spurious non-trivial transverse zero $\rho_{0}:a(\rho_{0})=0$ and $\partial_{\rho}a(\rho_{0})\neq 0$, then by restricting ourselves to the autonomous backbone curve (see Ponsioen et al. ponsioen2018automated ), i.e.

[TABLE]

we have found a solution

[TABLE]

The Jacobian of $\mathbf{F}$ with respect to $\rho$ and $\Omega$, evaluated at the solution $(\mathbf{u}_{0},0)$, is given by the square matrix

[TABLE]

which is invertible. Therefore, by the implicit function theorem, we can continue our solution as a two-dimensional submanifold of $\mathbb{R}^{4}$. Locally, we can express $\rho$ and $\Omega$ as a function of $\psi$ and $\varepsilon$. For $\varepsilon>0$, an isola is born out of the non-trivial transverse zero on the autonomous backbone curve located at $(\Omega,\rho)=(b(\rho_{0}),\rho_{0})$. For a fixed forcing ampltide $\varepsilon$, the isola is parameterized by $\psi$ (as illustrated in Fig. 3). ∎

Appendix F Proof of Theorem 4.2

In the setting of (32), our implicit function (29) will reduce to

[TABLE]

Any zero $\rho_{0}$ that makes the square root term in Eq. (96) vanish, will also be a zero of (96) itself by setting $\Omega=b(\rho_{0})$, and therefore will be on the forced response curve and on the autonomous backbone curve. Additionally, at this point, two segments of the FRC will meet and create a fold over the $\Omega$ direction. We set the argument inside the square root in Eq. (96) equal to zero and rewrite it as

[TABLE]

Restricting $\rho\in\mathbb{R}^{+}_{0}$, then for $\text{Re}(\gamma_{1})>0$, the third-order autonomous function $a(\rho)$ will have a trivial transverse zero and a non-trivial transverse zero located at

[TABLE]

such that $a(\rho_{0})=0$, $a(\rho_{1})=0$, $\partial_{\rho}a(\rho_{0})\neq 0$ and $\partial_{\rho}a(\rho_{1})\neq 0$. Under the assumption that the cubic order zero $\rho_{1}$ is a non-spurious zero for the function $a(\rho)$, then, using the same type of argument as in the proof of Theorem 4.1, an isola will be born out of this non-trivial transverse zero for system (32).

We note that between $\rho_{0}$ and $\rho_{1}$, $a(\rho)$ will have a local minimum at,

[TABLE]

Therefore, for $\varepsilon>0$ small enough, the function $a(\rho)$ will have three intersections with the constant curves $\pm\varepsilon\left\lVert c_{1,\mathbf{0}}\right\rVert$, meaning that we have found three zeros, $0<\rho_{0}^{\text{a}}<\rho_{1}^{\text{a}}<\rho_{1}^{\text{b}}$ of Eq. (97), that correspond to three folding points of the FRC over the $\Omega$ direction. In this setting, $\rho_{0}^{\text{a}}$ corresponds to the maximum amplitude of the main FRC branch, $\rho_{1}^{\text{a}}$ corresponds to the minimum amplitude of the isola, whereas $\rho_{1}^{\text{b}}$ will be the maximum amplitude of the isola.

We can increase $\varepsilon$ such that $\rho_{0}^{\text{a}}=\rho_{1}^{\text{a}}$, which merges the maximum amplitude of the main FRC branch with the minimum amplitude of the isola, which is exactly at

[TABLE]

whereas for

[TABLE]

the isola will be disconnected from the main FRC, proving Theorem 4.2. ∎