# A harmonic-based method for computing the stability of periodic solutions of dynamical systems

Arnaud Lazarus, Olivier Thomas

► To cite this version:

Arnaud Lazarus, Olivier Thomas. A harmonic-based method for computing the stability of periodic solutions of dynamical systems. Comptes Rendus Mécanique, 2010, 338 (9), pp.510-517. 10.1016/j.crme.2010.07.020. hal-01452004

HAL Id: hal-01452004

https://hal.science/hal-01452004

Submitted on 25 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A harmonic-based method for computing the stability of periodic solutions of dynamical systems

## Une méthode fréquentielle pour le calcul de stabilité des solutions périodiques des systèmes dynamiques

Arnaud Lazarus *, Olivier Thomas

Structural Mechanics and Coupled Systems Laboratory, Cnam, 2, rue Conté, 75003 Paris, France

### ABSTRACT

In this Note, we present a harmonic-based numerical method to determine the local stability of periodic solutions of dynamical systems. Based on the Floquet theory and the Fourier series expansion (Hill method), we propose a simple strategy to sort the relevant physical eigenvalues among the expanded numerical spectrum of the linear periodic system governing the perturbed solution. By mixing the harmonic-balance method and asymptotic numerical method continuation technique with the developed Hill method, we obtain a purely-frequency based continuation tool able to compute the stability of the continued periodic solutions in a reduced computation time. To validate the general methodology, we investigate the dynamical behavior of the forced Duffing oscillator with the developed continuation technique.

**Keywords:**
Dynamical systems
Stability
Hill's method
Continuation procedure
Harmonic-balance method

### RÉSUMÉ

Dans cette Note, nous présentons une méthode numérique fréquentielle pour déterminer la stabilité des solutions périodiques d'un système dynamique. La méthode, basée sur la théorie de Floquet et le développement en série de Fourier (méthode de Hill), consiste à extraire les valeurs propres physiques de l'ensemble des valeurs propres numériques du système perturbé étendu dans le domaine fréquentiel. En combinant alors la méthode de l'équilibrage harmonique et la méthode asymptotique numérique avec la précédente méthode de Hill, on obtient un outils de continuation purement fréquentiel où le calcul de la stabilité des solutions suivies est quasiment immédiat. Afin de valider la méthode, nous appliquons la méthode de continuation à un oscillateur de Duffing forcé.

**Mots-clés:**
Systèmes dynamiques
Stabilité
Méthode de Hill
Méthode de continuation
Méthode de l'équilibrage harmonique

### 1. Introduction

Determining the local stability of a dynamical system periodic solution is of primary interest in an engineering context since only stable solutions are experimentally encountered. Moreover, a change in the stability can lead to significant qualitative, and possibly dramatic, changes in the system response. A commonly-used method consists in computing the

<footnote>* Corresponding author.
E-mail addresses: arnaud.lazarus@polytechnique.edu (A. Lazarus), olivier.thomas@cnam.fr (O. Thomas).</footnote>
monodromy matrix of the linear periodic system governing the evolution of the perturbed solution and compare its eigenvalues (Floquet multipliers) moduli to one. Two well-known numerical strategies are classically used. The first technique consists in a time-integration of the Jacobian matrix of the initial dynamical system over one period [1]. In the second technique, the monodromy matrix is simply a by-product of a shooting continuation method [2].

The efficiency of those methods has been proved for years, but, inherently, the eigenvalues accuracy depends on the chosen time-steps size, leading to possible long time computations. Furthermore, those time-domain approaches are not necessarily best suited for nonlinear periodic solution computed, for instance, with a harmonic-balance method. In the latter, eigenvalues of the linear operator can naturally be computed with Hill's method also based on the Fourier series expansion. This harmonic-based method deals with linear periodic systems and has been introduced by Hill one century ago for the determination of the lunar perigee [3]. Since then, it has been applied to a wide range of physical problems like the computation of energy eigenvalues of the Schrödinger equation [4,5], the stability of limit cycles in electrical circuits [6] or the vibratory behavior of flexible rotating machines [7,8]. However, although Hill's method often provides satisfactory results, the meaning of the computed eigenvalues is generally misunderstood and may lead to wrong results. Actually, this method is not of a trivial use since it requires to approximate the spectra of infinite-dimensional operators, i.e. to sort the most converged eigenvalues among all the numerical ones. In this Note, we propose a simple numerical strategy to sort these eigenvalues and therefore properly determine the stability of periodic solutions of dynamical systems.

Based on the Floquet theory and the Fourier series expansion of both the Jacobian matrix and perturbed solution, we express the infinite-dimensional harmonic state operator of the linear periodic system governing the disturbance of the periodic solutions. We propose a simple criterion based on the energy distribution of the computed eigenvectors to extract the most converged eigenvalues (Floquet exponents) of the truncated problem and therefore determine the local stability of the solution. By mixing the standard harmonic-balance method with the so-called asymptotic numerical method continuation technique based on the quadratic recast of the dynamical system [9], it becomes a straightforward and efficient procedure to continue the periodic solutions as well as their associated Floquet exponents. To validate the general methodology, the developed harmonic-based continuation tool is applied to the forced Duffing oscillator. We confirm not only the convergence of the proposed Hill method but also its considerable running time benefit as compared to the time-domain monodromy matrix computation.

## 2. Stability of periodic solutions

### 2.1. General formulation

We consider continuous-time dynamical systems governed by the general equation

$$ \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), t, \lambda) \quad (1) $$

where $\mathbf{x}$ is an $N$-dimensional state vector and $\mathbf{f}$ is a nonlinear $N$-dimensional vector field that depends on a control parameter $\lambda$. In the following, $\mathbf{f}$ may explicitly depend on $t$ (the system is non-autonomous) or not (the system is autonomous).

We consider a periodic solution $\mathbf{x}_0(t)$ of (1), with minimal period $T$, at the particular control parameter value $\lambda = \lambda_0$. The stability of this periodic solution is studied by superimposing a small disturbance $\mathbf{y}(t)$:

$$ \mathbf{x}(t) = \mathbf{x}_0(t) + \mathbf{y}(t) \quad (2) $$

Substituting (2) into (1), assuming that $\mathbf{f}$ is at least twice continuously differentiable, expanding the result in a Taylor series about $\mathbf{x}_0$, and retaining only linear terms in the disturbance, we obtain:

$$ \dot{\mathbf{y}}(t) = \mathbf{J}(t)\mathbf{y}(t) \quad (3) $$

with

$$ \mathbf{J}(t) = \mathbf{J}(\mathbf{x}_0(t), t, \lambda_0) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x}_0(t), t, \lambda_0) \quad (4) $$

$\mathbf{J}(t)$ is the $N \times N$ Jacobian matrix of $\mathbf{f}$, at $\mathbf{x} = \mathbf{x}_0$ and $\lambda = \lambda_0$. The stability study of the periodic solution $\mathbf{x}_0$ consists in finding if the disturbance $\mathbf{y}(t)$, solution of (3), fades away or is amplified as $t$ is increased. Since $\mathbf{x}_0(t)$ is $T$-periodic in time, $\mathbf{J}(t)$ is also $T$-periodic in time considering its definition (4). Consequently, the system (3) is a linear system with periodic coefficients. The Floquet theory [10,11], used in the following, specifically deals with this kind of dynamical systems.

### 2.2. Floquet theory

The following developments are classical and can be found, e.g., in the textbooks [1,12]. The $N$-dimensional linear system (3) has $N$ linearly independent solutions $\mathbf{y}_n(t)$, so that any solution $\mathbf{y}(t)$ of (3) can be written:

$$ \mathbf{y}(t) = \sum_{n=1}^{N} c_n \mathbf{y}_n(t) \quad (5) $$
where $c_n$ are N constants that depend upon the initial conditions. For all $n = 1, \dots, N$, every fundamental solution $\mathbf{y}_n(t)$ can be expressed in the Floquet form

$$ \mathbf{y}_n(t) = \mathbf{p}_n(t) e^{\alpha_n t} \qquad (6) $$

where $\mathbf{p}_n(t + T) = \mathbf{p}_n(t)$ is a T-periodic N-dimensional complex vector and $\alpha_n$ is a complex number. Then, it follows from (6) and from the T-periodicity of $\mathbf{p}_n$ that

$$ \mathbf{y}_n(t + T) = \mathbf{p}_n(t + T) e^{\alpha_n(t+T)} = \mathbf{y}_n(t) e^{\alpha_n T} \qquad (7) $$

The $\alpha_n$ represent the so-called Floquet exponents. Replacing $\alpha_n$ by $\alpha_n + 2ik\pi/T$ in (7) shows that the $\alpha_n$ are unique to within an additive integer multiple of $2i\pi/T$ [11].

Considering (6) or (7), the values of the Floquet exponents $\alpha_n$ can be used to determine the stability of periodic solutions.

* If $\Re(\alpha_n) < 0$ for all $n$, all fundamental solutions $\mathbf{y}_n$ converge toward zero as $t$ is increased, so does disturbance $\mathbf{y}(t)$. The periodic solution is said to be asymptotically stable.

* If a subscript $n$ exists such that $\Re(\alpha_n) > 0$, the corresponding fundamental solutions increase exponentially, so does disturbance $\mathbf{y}(t)$. The periodic solution is unstable in this case.

The above statement stands for non-autonomous initial systems (1). If the initial system (1) is autonomous, one of the Floquet exponent is always $\alpha_n = 2ik\pi/T$. In this circumstances, only the set of remaining eigenvalues can provide information on the periodic solution stability.

## 2.3. Hill's method

Hill's method is a harmonic-based numerical approach used to determine the solutions of linear periodic systems like (3), which is used here to calculate the Floquet exponents $\alpha_n$. This method has been widely used together with perturbation methods to analytically obtain the stability of periodic solutions, for instance in the case of the Mathieu equation (see e.g. [13]).

The unknown periodic functions $\mathbf{p}_n(t)$ introduced in (6) are expressed by the general Fourier series

$$ \mathbf{p}_n(t) = \sum_{k=-\infty}^{+\infty} \mathbf{p}_n^k e^{ik\omega t} \qquad (8) $$

where the fundamental frequency $\omega$ reads $\omega = 2\pi/T$ and $\mathbf{p}_n^k$ are N-dimensional complex vectors. By replacing (8) in the Floquet form (6), any fundamental solution $\mathbf{y}_n(t)$ is an infinite sum of harmonic contributions

$$ \mathbf{y}_n(t) = \sum_{k=-\infty}^{+\infty} \mathbf{p}_n^k e^{(\alpha_n + ik\omega)t} \qquad (9) $$

Being T-periodic as well, the Jacobian matrix $\mathbf{J}(t)$ is also expanded in the infinite Fourier series

$$ \mathbf{J}(t) = \sum_{h=-\infty}^{+\infty} \mathbf{J}^h e^{ih\omega t} \qquad (10) $$

where the $\mathbf{J}^h$ are $N \times N$ matrices. By replacing the solution $\mathbf{y}_n(t)$ and the Jacobian $\mathbf{J}(t)$ by their Fourier series in the linear periodic system (3), one obtains the following vector equation, for all $n = 1, \dots, N$:

$$ \sum_{k=-\infty}^{+\infty} (\alpha_n + ik\omega) \mathbf{p}_n^k e^{(\alpha_n + ik\omega)t} = \sum_{k=-\infty}^{+\infty} \sum_{h=-\infty}^{+\infty} \mathbf{J}^h \mathbf{p}_n^k e^{[\alpha_n + i(k+h)\omega]t} \qquad (11) $$

Since the sums have an infinite number of terms, replacing superscript $k$ by $k-h$ in the left-hand side term of the above equation does not changes it. It leads to, for all $n = 1, \dots, N$:

$$ \sum_{k=-\infty}^{+\infty} \left[ \sum_{h=-\infty}^{+\infty} \mathbf{J}^h \mathbf{p}_n^{k-h} - (\alpha_n + ik\omega) \mathbf{p}_n^k \right] e^{(\alpha_n + ik\omega)t} = \mathbf{0} \qquad (12) $$

By applying the harmonic-balance method to the above equation, i.e. by separately equating to zero each harmonic of (12) (for each value of $k$), the above equation can be rewritten in the following infinite-dimensional eigenproblem

$$ (\mathbf{H} - s\mathbf{I})\mathbf{q} = \mathbf{0} \qquad (13) $$
where **H** is the infinite-dimensional Hill matrix

$$ \mathbf{H} = \begin{bmatrix} \mathbf{J}^0 + i\omega\mathbf{I} & \mathbf{J}^{-1} & \mathbf{J}^{-2} & \cdots \\ \mathbf{J}^1 & \mathbf{J}^0 & \mathbf{J}^{-1} & \cdots \\ \mathbf{J}^2 & \mathbf{J}^1 & \mathbf{J}^0 - i\omega\mathbf{I} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \tag{14} $$

$s$ is a complex number, $\mathbf{q}$ is an infinite-dimensional vector and $\mathbf{I}$ is the identity matrix of appropriate size.

By comparing (13) and (12), the eigenvalues and eigenvectors $s$ and $\mathbf{q}$ are related to $\alpha_n$ and $\mathbf{p}_n$ by the following relations:

$$ s_n^l = \alpha_n + il\omega \quad \text{and} \quad \mathbf{q}_n^l = [\dots \mathbf{p}_n^{-1-l} \mathbf{p}_n^{0-l} \mathbf{p}_n^{1-l} \dots] \tag{15} $$

valid for all $l = 0, \pm 1, \dots, \pm \infty$ and $n = 1, \dots, N$. In the above equation, $\bullet^T$ means the vector transpose of $\bullet$. The above relations show that each Floquet exponent $\alpha_n$, for a given $n$, is associated to an infinite set of eigenvalues $s_n^l$, for all $l = 0, \pm 1, \dots, \pm \infty$. Consequently, knowing the infinite set of $s_n^l$ would enable to calculate the $\alpha_n$ and then to assess the periodic solution stability.

However, for evident practical purposes and numerical computations, the infinite-dimensional problem (12)–(14) is truncated to a finite dimension. Namely, matrix $\mathbf{H}$ is truncated to an $N(2H+1) \times N(2H+1)$ dimension, so that $l = 0, \pm 1, \dots, \pm H$ in (15). Let $\hat{s}_n^l$ and $\hat{\mathbf{q}}_n^l$ denote the $N(2H+1)$ approximate eigensolutions of Eq. (13). If the $\alpha_n$ are to be computed, two questions need to be solved: (1) do the $(\hat{s}_n^l, \hat{\mathbf{q}}_n^l)$ tend to $(s_n^l, \mathbf{q}_n^l)$ as $H$ tends to infinity and (2) for a given $n$, among all the computed $\hat{s}_n^l$, which ones, as $l$ is varied, can be used to compute $\alpha_n$.

Firstly investigated by Poincaré [14] more than one century ago, those questions were still argued twenty years ago, in particular for the computation of the Schrödinger equation [4]. Recently, papers based on the Fredholm theory provide rigorous proof of the convergence [15,16].

However, a delicate issue remains the choice of the eigensolutions of the truncated problem that can be used in a numerical procedure to evaluate the $\alpha_n$, $n = 1, \dots, N$. A simple idea consists in analysing the truncated problem eigenvectors distribution $\hat{\mathbf{q}}_n^l = [\mathbf{p}_n^{-H-l} \dots \mathbf{p}_n^{0-l} \dots \mathbf{p}_n^{H-l}]^T$ and considering the one associated to $l = 0$:

$$ \hat{s}_n^0 = \alpha_n \quad \text{and} \quad \hat{\mathbf{q}}_n^0 = [\mathbf{p}_n^{-H} \dots \mathbf{p}_n^{-1} \mathbf{p}_n^0 \mathbf{p}_n^1 \dots \mathbf{p}_n^H] \tag{16} $$

Since all the above developments are based on the assumption that the Fourier series (8) and (10) are convergent ($\mathbf{p}_n^k$ tends to zero as $k$ tends to infinity), among all the computed eigenvectors, the components of the $N$ ones associated with $l = 0$ (the $\hat{\mathbf{q}}_n^0$) have almost symmetric shapes since they involve only the $\mathbf{p}_n^k$ with $|k| < H$. For this reason, we assume that the $\hat{\mathbf{q}}_n^0$ converge faster to $\mathbf{q}_n^0$ than the others, for $l \neq 0$.

With the above assumption, the numerical procedure to compute the $N$ Floquet exponents $\alpha_n$ is the following. The Hill matrix $\mathbf{H}$, truncated to order $H$ is assembled and its $N(2H+1)$ eigensolutions $(\hat{s}_n^l, \hat{\mathbf{q}}_n^l)$ computed. Then, the $N$ ones associated to $l = 0$ are numerically selected by considering that the $N$ eigenvectors $\hat{\mathbf{q}}_n^0$ are the $N$ ones with the most symmetric shapes. The $N$ associated eigenvalues are then the $N$ Floquet exponents: $\alpha_n = \hat{s}_n^0, n = 1, \dots, N$. Their real parts are then compared to 0 to assess the stability of the periodic solution.

### 3. The HBM-ANM-Hill method

One of the crucial points in the numerical method mentioned above is the expansion of the Jacobian $J(t)$ in the Fourier series (Eq. (10)). Indeed, in the most general situation, the nonlinear evolution equation (1) shows nonlinearities of any kinds and the same is true for the Jacobian dependence on the periodic solution. Therefore, computation of its Fourier coefficients $J^h$ can be very cumbersome.

However, in the special case where the harmonic-balance method (HBM) combined with the asymptotic numerical method (ANM) is used to compute the Fourier coefficients of the periodic solution $\mathbf{x}_0$ (as explained in [9]), the computation of the $J^h$ is rendered very efficiently. The leading idea, also used for the ANM, is to systematically recast the dynamical system (1) in a quadratic polynomial form. From (1), the new quadratic system can be written as follows

$$ m(\dot{\mathbf{z}}) = c(\lambda, t) + l(\mathbf{z}, t, \lambda) + q(\mathbf{z}, \mathbf{z}, t, \lambda) \tag{17} $$

where $c, l(\cdot), q(\cdot, \cdot)$ are constant, linear and quadratic operators in $\mathbf{z}$, respectively, and $m(\cdot)$ is a linear operator. The unknown vector $\mathbf{z}$ contains the original components of the state vector $\mathbf{x}$ and some new variables added to get the quadratic form. This procedure is fully explained in [9], which also shows that this quadratic recast can be applied to a large class of smooth systems with a few algebraic manipulations and a few additions of auxiliary variables.

A periodic solution of (17), of minimal period $T = 2\pi/\omega$, is expanded into the Fourier series
$$ z_0(t) = \sum_{p=-H}^{H} z_0^p e^{ip\omega t} \qquad (18) $$

The HBM combined to the ANM continuation method enables to compute the $z_0^p$, $p = -H, \dots, H$, and $\omega$ for various contiguous values of $\lambda$, to detect possible bifurcations and compute different branches of solutions. For externally driven systems, the control parameter is often the forcing frequency $\lambda = \omega$. Assuming the quadratic recast (17), for any $\lambda = \lambda_0$, the Jacobian $J(t)$ defined by (4) may be expressed in the form

$$ J(t) = J_C(\lambda_0, t) + J_L(z_0(t), t, \lambda_0) + J_Q(z_0(t), z_0(t), t, \lambda_0) \qquad (19) $$

where $J_C, J_L(\cdot, \cdot)$ and $J_Q(\cdot, \cdot)$ are constant, linear and quadratic $N \times N$ matrices with respect to $z_0$. From then, it becomes straightforward to find the Fourier coefficients $J^h$ of $J(t)$ according to the ones of $z_0(t)$:

$$ J^0 = J_C(\lambda_0, t) + J_L(z_0^p, t, \lambda_0) + \sum_{p=-H}^{H} J_Q(z_0^p, z_0^{-p}, t, \lambda_0) \qquad (20) $$

$$ J^h = J_L(z_0^h, t, \lambda_0) + \sum_{p=-H}^{H} J_Q(z_0^p, z_0^{-p}, t, \lambda_0) \qquad (21) $$

Then, for each parameter $\lambda_0$, the Hill matrix $H$ is expressed in terms of the already computed $z_0^p$ and $\omega$ and the stability of the periodic solution $x_0(t)$ is simply obtained following the previous Hill's method.

The selection of the $N$ eigenvectors $q_n^0$ among the $N(2H+1)$ eigenvectors of (13) is done by computing for each $q_n^l$, $l = -H, \dots, H$, $n = 1, \dots, N$, the median value of its components moduli. If $q_i$ is the $i$-th component of vector $q$ of dimension $Q = N(2H+1)$, the median is $\sum_i |q_i| / \sum_i |q_i|$. The selected $q_n^0$ are the $N$ ones with the median value closest to $Q/2$, which are supposed to have the most symmetrical shape. Knowing those eigenvectors, the $N$ corresponding eigenvalues are the $\alpha_n$, used for the stability study.

## 4. Validation example

To give more insight in the stability study technique explained above as well as to give a validation on a well-known example, the stability of the periodic solutions of a forced hardening Duffing oscillator is presented here in detail. The studied Duffing equation is written:

$$ \ddot{u} + \mu \dot{u} + \omega_0^2 u + \Gamma u^3 = F \cos(\omega t) \qquad (22) $$

The particular values of parameters used here are: $\omega_0 = 1$, $\Gamma = 1$, $F = 1$, $\mu = 0.05$. The periodic solution $u(t)$ of (22), for various values of $\omega$ around the primary resonance, are computed with the HBM combined with the ANM used to continue the solution branches, as explained in [9]. The first step is to write (22) in the form of a first-order dynamical system of dimension $N=2$:

$$ \begin{cases} \dot{u} = v \\ \dot{v} = F \cos(\omega t) - \mu v - \omega_0^2 u - \Gamma u^3 \end{cases} \qquad (23) $$

and then to transform it into a new one with quadratic nonlinear terms only:

$$ \begin{pmatrix} \dot{u} \\ \dot{v} \\ 0 \end{pmatrix} = \begin{pmatrix} 0 \\ F \cos(\omega t) \\ 0 \end{pmatrix} + \begin{pmatrix} v \\ -\mu v - \omega_0^2 u \\ w \end{pmatrix} + \begin{pmatrix} 0 \\ -\Gamma u w \\ -u^2 \end{pmatrix} \qquad (24) $$

The new unknown variables $v = \dot{u}$ and $w = u^2$ have been added. If $z$ is defined as $z(t) = [u(t) \ v(t) \ w(t)]^T$ and $\omega = \lambda$ is the control parameter, the above dynamical system has the form of (17). The Jacobian associated to the system (23) reads:

$$ J(t) = \begin{bmatrix} 0 & 1 \\ -\omega_0^2 & -\mu \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ -3\Gamma u^2(t) & 0 \end{bmatrix} \qquad (25) $$

It is naturally in the form of Eq. (19). The HBM-ANM method enables to compute all Fourier coefficients of $u, v$ and $w$ up to harmonics of order $H$. Those Fourier coefficients are then used to compute the ones of $J(t)$ using Eqs. (20) and (25), the Hill matrix (14) is assembled, the eigenvalue problem is solved and the stability of periodic solution $u(t)$ is determined by comparing $\Re(\alpha_n)$ to zero.

Fig. 1 shows the results of the HBM-ANM-Hill method applied to the Duffing oscillator. The first three harmonics of $u(t)$ are shown as a function of the excitation frequency, around the primary resonance. Details around the superharmonics resonances of order 2 ($\omega \simeq \omega_1/2$) and 3 ($\omega \simeq \omega_1/3$) are also shown. The shape of the branches qualitatively conforms to the one obtained by perturbation methods (see e.g. [17] for the special shape of the third harmonics, Fig. 1 (right)). The
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Resonance curves of a forced Duffing oscillator with $\omega_0 = 1$, $\Gamma = 1$, $F = 1$, $\mu = 0.05$: amplitude of the first three harmonics as a function of forcing frequency. Simulations with HBM-ANM, stability with the Hill method, with $H = 12$ harmonics. '-': stable branches; '...': unstable branches. 'SN': saddle node bifurcations; 'SB': symmetry breaking bifurcations. 'H<sub>1</sub>b<sub>j</sub>': i-th harmonic, j-th branch. (left): general view; (middle), (right): details of superharmonic resonance of order 2 and 3.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The $N = 2$ Floquet multipliers moduli $|\rho_n|$ as functions of the excitation frequency $\omega$, corresponding to the simulation of Fig. 1. 'o': $\rho_n$ computed with the Hill's method; '•': $\rho_n$ computed with the monodromy matrix, for various numbers $N_s$ of time samples per period $T$. (left): general view; (middle), (right): details of superharmonic resonance of order 2 and 3. Gray region corresponds to unstable periodic solutions.</figcaption></figure>

stability has been determined by the Hill method. The saddle node (SN) bifurcations associated to the primary resonance as well as the secondary resonance of order 3 are precisely located at the turning points (the points with a vertical tangent), thus validating our method [17]. The particular symmetry breaking bifurcation (SB) around the secondary resonance of order 2 is also well recovered [17]. The primary branch (denoted as $b_1$ in Fig. 1), for which $u(t)$ is symmetric since all even harmonics are zero, is found unstable for $\omega \in [0.695 \ 0.735]$, giving birth to two stable branches (branches $b_2$ and $b_3$) with non-zero even harmonics. Those two branches have opposite constant component ($H_0b_2$, $H_0b_3$), identical to the first harmonics ($H_1b_2$, $H_1b_3$) and opposite to the second harmonics (with identical amplitude and phase difference of $\pi$, $H_2b_2$, $H_2b_3$). All simulations have been obtained by truncating the Fourier series to $H = 12$ harmonics.

Fig. 2 enables to compare our stability results to those obtained with a classical time-domain method. The latter consists in integrating in time over one period the linear transformation given in Eq. (7). By gathering the fundamental solutions
<caption>**Table 1**</caption>

Relative and absolute time to compute the diagrams of Fig. 2 with the HBM-ANM continuation method on a standard desk PC (with an Intel core 2 quad CPU, 2.83 GHz) with matlab R2007-b, including the stability computation with the Hill method and with the time-integration method to compute the monodromy matrix with $N_s$ sample points over one period.

<table><thead><tr><th rowspan="2">HBM-ANM procedure</th><th rowspan="2">Without stability</th><th rowspan="2">Hill</th><th colspan="3">Time-domain</th></tr><tr><th>N<sub>s</sub> = 100</th><th>N<sub>s</sub> = 500</th><th>N<sub>s</sub> = 1000</th></tr></thead><tbody><tr><td>Computation time</td><td>1 (158 s)</td><td>1.38 (218 s)</td><td>1.39 (219 s)</td><td>2.76 (436 s)</td><td>4.46 (705 s)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Eigenvalues of $\mathbf{H}$ as a function of the excitation frequency $\omega$, corresponding to the simulation of Fig. 1. '•': all $\hat{s}_n^k$, for $n = 1, 2$ and $k = 0, \pm 1, \dots, \pm H$; '$o$: $\alpha_n$, $n = 1, 2$ computed with the Hill's method. (left): real part; (right): imaginary part. Gray region corresponds to unstable periodic solutions.</figcaption></figure>

$\mathbf{y}_n(t)$ of Eq. (5) into the $N \times N$ matrix $\mathbf{Y}(t) = [\mathbf{y}_1(t) \ \mathbf{y}_2(t) \ \dots \ \mathbf{y}_N(t)]$, there exists an $N \times N$ constant matrix $\mathbf{\Phi}$, called the monodromy matrix [1], so that:

$$ \mathbf{Y}(t+T) = \mathbf{Y}(t)\mathbf{\Phi} \qquad (26) $$

Specifying the initial conditions $\mathbf{Y}(0) = \mathbf{I}$, where $\mathbf{I}$ is the $N \times N$ identity matrix, the monodromy matrix is simply $\mathbf{\Phi} = \mathbf{Y}(T)$ and the stability is obtained by comparing the Floquet multipliers $\rho_n$ moduli (its eigenvalues) to one. The time-integration is here performed by a fourth-order Runge-Kutta algorithm. The above equation, compared to Eq. (7), proves that the $\rho_n$ are linked to the Floquet exponents by, for all $n = 1, \dots, N$:

$$ \rho_n = e^{\alpha_n T} \Leftrightarrow \alpha_n = \frac{1}{T} \ln \rho_n + \frac{2\pi k}{T} i, \quad k = 0, \pm 1, \pm 2, \dots \qquad (27) $$

As compared to the Hill method, one has to set an additional parameter: the number $N_s$ of sample points over the computed period for the numerical integration. Fig. 2 shows $|\rho_n|$, $n = 1, 2$, as a function of $\omega$, computed with the Hill method through Eq. (27) as well as with the above described time-domain method, with various values of $N_s$. It shows that the Hill method appears as a limit to the time domain method as $N_s$ is increased. At least $N_s = 1000$ time integration points over the computed period are necessary to recover the right stability results and especially to place the saddle-node bifurcation points at the turning points of the solution branches. This large amount of integration points for the time-domain method is also related to an increased computational cost, as compared to the Hill method. Table 1 shows the needed running time to compute the continuation diagrams of Fig. 2. The frequency method is by far the more efficient approach to obtain the continued stability map of the dynamical system.

Fig. 3 highlights the relevance of our strategy to sort the most converged eigenvalues among all the computed $\hat{s}_n^l$. As expected, the computed imaginary parts $\Im(\hat{s}_n^l)$ are in qualitative agreement with the relation $\Im(\hat{s}_n^l) = \Im(\alpha_n) + il\omega$ (Eq. (15)) for the small values of $l$ visualized in Fig. 3(right). Due to the truncation error, this will not be the case when considering eigenvalues associated with high values of $|l|$. Since for a given $n$, the real parts of the $s_n^l$ are equal ($\Re(s_n^l) = \Re(\alpha_n)$ for all $l = 0, \pm 1, \dots, \pm \infty$), it would be theoretically possible to use any of them to determine the periodic solution stability. However, among the truncated spectrum $\hat{s}_n^l$, there exists a set of eigenvalues with a real part that leads to erroneous stability predictions, for values of $|l|$ close to $H$. Fig. 3(left) clearly shows the wrong instability branches computed for the forced Duffing oscillator with $H = 12$. Thus, a stability criterion based on the consideration of real parts of all eigenvalues would have led to erroneous results. Sorting the most symmetric eigenvectors $\hat{q}_n^l$ seems definitely a simple and correct method to follow the $N$ more convergent physical eigenvalues among all the numerical ones.

According to the expanded Floquet form (9), the convergents $\Im(\hat{s}_n^l)$ of Fig. 3(right) represents the location of the frequency spectrum of the disturbance $\mathbf{y}(t)$ as a function of the excitation frequency $\omega$. It is interesting to note that the
destabilization process is associated to a lock-in of the disturbance frequencies ($\tilde{\omega}(\alpha_1) = \tilde{\omega}(\alpha_2)$) on $\omega$-ranges embedding the unstable regions. This phenomenon is well known, especially when dealing with Flutter instabilities [18].

## 5. Conclusion

In this article, a purely frequency-based numerical method to determine the local stability of periodic solutions of dynamical systems, based on computing the eigenvalues of the problem Hill matrix, has been presented. Its efficiency, in term of accuracy and computation time, has been shown when associated to a harmonic balance based continuation method. While the monodromy matrix is simply a by-product of a shooting continuation method in the time-domain approach, here, the Hill matrix is simply a by-product of a harmonic-balance continuation method in the frequency domain. The main interest of this method, which explains its efficiency, is that there is no need to switch from one domain to another for computing the stability.

## Acknowledgements

The authors want to thank Bruno Cochelin, who is at the origin of the ANM continuation tool used in this article: our discussions and meetings, sometimes quick, have materialized! The second author (O.T.) also warmly thanks his 10 years companion and friend Cyril Touzé, for the uncountable and fruitful discussions on the nonlinear world and more specifically for the TGV monodromy matrix computation algorithm.

## References

[1] A.H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods, John Wiley & Sons, Inc., 1995.

[2] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, J.-C. Golinval, Nonlinear normal modes. Part II: Toward a practical computation using numerical continuation techniques, Mechanical Systems and Signal Processing 23 (2009) 195-216.

[3] G.W. Hill, On the part of the motion of the lunar perigee which is a function of the mean motions of the sun and moon, Acta Mathematica 8 (1) (1886) 1-36.

[4] J. Killingbeck, On the Hill determinant method, Journal of Physics A: Mathematical and General 19 (1986) 2903-2906.

[5] B. Deconinck, J.N. Kutz, Computing spectra of linear operators using the Floquet-Fourier-Hill method, Journal of Computational Physics 219 (2006) 296-321.

[6] F. Bonani, M. Gilli, Analysis of stability and bifurcations of limit cycles in Chua's circuit through the harmonic-balance approach, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 46 (1999) 881-890.

[7] C. Villa, J.-J. Sinou, F. Thouverez, Stability and vibration analysis of a complex flexible rotor bearing system, Communications in Nonlinear Science and Numerical Simulation 13 (2008) 804-821.

[8] A. Lazarus, D. Combescure, B. Prabel, A 3D finite element model for the vibration analysis of asymmetric rotating machines, Journal of Sound and Vibration, doi:10.1016/j.jsv.2010.03.029.

[9] B. Cochelin, C. Vergez, A high order purely-based harmonic balance formulation for continuation of periodic solutions, Journal of Sound and Vibration 324 (2009) 243-262.

[10] G. Floquet, Sur la théorie des équations différentielles, Annales Scientifiques de l'École Normale Supérieure 8 (1879) 3-132.

[11] P. Hartman, Ordinary Differential Equations, Society for Industrial and Applied Mathematics, 2002.

[12] J. Argyris, G. Faust, M. Haase, An Exploration of Chaos, Texts on Computational Mechanics, vol. VII, North-Holland, 1994.

[13] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, Wiley, 1979.

[14] H. Poincaré, Sur les déterminants d'ordre infini, Bulletin de la Société Mathématique de France 14 (1886) 77-90.

[15] J. Zhou, T. Hagiwara, M. Araki, Spectral analysis and eigenvalues computation of the harmonic state operators in continuous-time periodic systems, Systems & Control Letters 53 (2004) 141-155.

[16] C.W. Curtis, B. Deconinck, On the convergence of Hill's method, Mathematics of Computation 79 (2010) 169-187.

[17] W. Szemplinska-Stupnicka, The Behavior of Nonlinear Vibrating Systems, vol. 1: Fundamental Concepts and Methods: Applications to Single-Degree-of-Freedom Systems, Kluwer Academic Publishers, 1990.

[18] E. de Langre, Frequency lock-in is caused by coupled-mode flutter, Journal of Fluids and Structures 22 (2006) 783-791.