# Continuation of quasi-periodic solutions with two-frequency Harmonic Balance Method

Louis Guillot, Pierre Vigué, Christophe Vergez, Bruno Cochelin

► To cite this version:

Louis Guillot, Pierre Vigué, Christophe Vergez, Bruno Cochelin. Continuation of quasi-periodic solutions with two-frequency Harmonic Balance Method. Journal of Sound and Vibration, 2017, 394, pp.434-450. hal-01413398

HAL Id: hal-01413398

https://hal.science/hal-01413398

Submitted on 12 Dec 2016

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 - ShareAlike 4.0 International License
# Continuation of quasi-periodic solutions with two-frequency Harmonic Balance Method

Louis Guillot<sup>a</sup>, Pierre Viguë<sup>b</sup>, Christophe Vergez<sup>b</sup>, Bruno Cochelin<sup>b</sup>

a : École Normale Supérieure de Cachan, France

b : Aix Marseille Univ, CNRS, Centrale Marseille, LMA, Marseille, France

louis.guillot@ens-cachan.fr, {vigue, vergez} @lma.cnrs-mrs.fr, bruno.cochelin@centrale-marseille.fr

## Abstract

The continuation of quasi-periodic solutions for autonomous or forced nonlinear systems is presented in this paper. The association of the Asymptotic Numerical Method, a robust continuation method, and a two-frequency Harmonic Balance Method, is performed thanks to a quadratic formalism. There is no need for a *priori* knowledge of the solution: the two pulsations can be unknown and can vary along the solution branch, and the double Fourier series are computed without needing a harmonic selection. A norm criterion on Fourier coefficients can confirm *a posteriori* the accuracy of the solution branch. On a forced system, frequency-locking regions are approximated, without blocking the continuation process. The continuation of these periodic solutions can be done independently. On an autonomous system an example of solution is shown where the number of Fourier coefficients is increased to improve the accuracy of the solution.

## 1. Introduction

Time-periodic solutions are commonly investigated in dynamical systems. They can be a source of nuisance, causing noise or potentially destructive motions (for example in rotating machinery) or they can be sought after (for example in many musical instruments relying on auto-oscillations: winds, bowed strings, etc). The well-known continuation of periodic solutions aims at producing bifurcation diagrams representing the existence, stability and other characteristics of such periodic solutions (amplitude, pulsation, etc) with respect to some parameters of the nonlinear system.

However, other solutions can arise, among which, quasi-periodic solutions. Our aim in this paper is the continuation of two-pulsation, quasi-periodic solutions. The direct computation of quasi-periodic solutions through numerical integration can be difficult. Because of dependence on initial conditions, some solutions may be overlooked (as in the periodic case). Moreover, performing integration on long intervals to get rid of transient solutions - a brute force, expensive tactic when studying periodic solutions - needs a new kind of stopping criterion since one cannot easily determine if the steady-state solution is reached when non periodic solutions are considered. These drawbacks have lead to different approaches to the study of quasi-periodicity.

Chua and Ushida [1] adapted the Harmonic Balance Method (HBM) to compute the steady-state response to a quasi-periodic forcing. Another early work based on a frequency-domain approach was the incremental HBM devised by Lau *et al.* [2]. Kaas-Petersen [3, 4] reformulated the problem in the time domain in terms of Poincaré maps. An interesting advantage to this approach is that it provides a stability analysis [5, 6].

The Alternating Frequency-Time method proposed by Cameron [7] was an important breakthrough, allowing the consideration of strong nonlinearities by computing them in the time domain. Several works used it to find quasi-periodic solutions [8] and coupled it with arc-length continuation [9, 10].

Finally, instead of computing a trajectory on an invariant tori, the continuation of these tori was performed independently by Rasmussen [11] and Schilder *et al.* [12, 13].

The method of obtaining and continuing quasi-periodic solutions presented in this paper is different and original. The aim is to create a resilient method that is capable of coping with either forced or autonomous systems, and is fast enough to tackle algebraic systems with a few thousand unknowns. The association of the Asymptotic Numerical Method (ANM), a continuation technique based on Taylor series, with the Harmonic Balance Method (HBM), proved its effectiveness to continue periodic solutions of forced or autonomous
systems [14]. This framework requires a quadratic formalism, which can be obtained for most problems thanks to auxiliary variables. The present paper extends the principle of coupling ANM and HBM to the quasi-periodic case. The ANM is briefly recalled in section 2. It is a very robust continuation method and its choice, differing from the widespread arc-length continuation, aims to deal appropriately with difficult situations, for instance frequency-locking.

Several aspects of the HBM are presented in section 3. After a general sketch of its association with the ANM, an original implementation of the periodic HBM is given in section 3.2. It is based on the complex representation of the Fourier basis. Then, its adaptation for quasi-periodic solutions is presented in section 3.3. The method is capable of dealing with periodically forced systems as well as autonomous systems, illustrated respectively in sections 4 and 5. In this last case, the two pulsations are both unknown. The method presented here does not need to select specific spectral contributions, and to the authors' knowledge it is able to compute solutions with more harmonics than existing methods.

## 2. Asymptotic Numerical Method (ANM)

Let us consider a vector-valued, smooth function:

$$ \mathbf{R}: \mathbb{R}^{N+1} \to \mathbb{R}^N, \quad (\mathbf{X}, \lambda) \mapsto \mathbf{R}(\mathbf{X}, \lambda) \quad (1) $$

Assuming that a first regular solution $(\mathbf{X}_0, \lambda_0)$ is known, the aim of the Asymptotic Numerical Method (ANM) is to follow the solution branch of the equation:

$$ \mathbf{R}(\mathbf{X}, \lambda) = 0 \quad (2) $$

with respect to $\lambda$. Without a great loss of generality $\mathbf{R}$ can be considered quadratic, meaning that there are $\mathbf{C}_0, \mathbf{C}_1 \in \mathbb{R}^N, \mathbf{L}_0, \mathbf{L}_1$ real $N \times N$ matrices, $\mathbf{Q}_0$ a bilinear application of $\mathbb{R}^N \times \mathbb{R}^N \to \mathbb{R}^N$, so that:

$$ \mathbf{R}(\mathbf{X}, \lambda) = \mathbf{C}_0 + \lambda \mathbf{C}_1 + \mathbf{L}_0 \mathbf{X} + \lambda \mathbf{L}_1 \mathbf{X} + \mathbf{Q}_0(\mathbf{X}, \mathbf{X}) \quad (3) $$

As noted when the method was designed, the transformation from a general function $\mathbf{R}$ to the quadratic formalism required in Eq. (3) is possible for usual smooth functions, adding auxiliary variables [15, 16, 17]. The continuation parameter is now treated as an unknown; thus, a vector of all unknowns is defined as $\mathbf{V} := (\mathbf{X}, \lambda)$, and the idea of the ANM is to develop the solution branch as a (truncated) power series of a path parameter $a$, in the vicinity of a known solution $\mathbf{V}_0$. Let $\mathbf{V}_1$ be a tangent vector at $\mathbf{V}_0$ and let $a$ be defined as:

$$ a := (\mathbf{V} - \mathbf{V}_0)^t \mathbf{V}_1 \quad (4) $$

Then $\mathbf{V}$ is developed to the order $n$ ($n = 15$ or 20 in applications):

$$ \mathbf{V}(a) = \mathbf{V}_0 + a\mathbf{V}_1 + a^2\mathbf{V}_2 + \dots + a^n\mathbf{V}_n \quad (5) $$

Operators are combined as follows:

$$ \mathbf{C} = \mathbf{C}_0, \quad \mathbf{L}(\mathbf{V}) = \lambda \mathbf{C}_1 + \mathbf{L}_0 \mathbf{X}, \quad \mathbf{Q}(\mathbf{V}, \mathbf{V}) = \lambda \mathbf{L}_1 \mathbf{X} + \mathbf{Q}_0(\mathbf{X}, \mathbf{X}) \quad (6) $$

so Eq. (3) now reads as:

$$ \mathbf{R}(\mathbf{V}) = \mathbf{Q}(\mathbf{V}, \mathbf{V}) + \mathbf{L}(\mathbf{V}) + \mathbf{C} \quad (7) $$

The series of Eq. (5) is substituted in Eq. (2), and powers of $a$ are collected to obtain a family of linear systems on all $\mathbf{V}_i$:

* Order 0: $\mathbf{Q}(\mathbf{V}_0, \mathbf{V}_0) + \mathbf{L}(\mathbf{V}_0) + \mathbf{C} = 0$ which is true since $\mathbf{V}_0$ is a solution of Eq. (2).

* Order 1:

$$ \mathbf{Q}(\mathbf{V}_0, \mathbf{V}_1) + \mathbf{Q}(\mathbf{V}_1, \mathbf{V}_0) + \mathbf{L}(\mathbf{V}_1) = 0 \quad (8a) $$

$$ \text{path parameter definition (4) leads to } \mathbf{V}_1^t \mathbf{V}_1 = 1 \quad (8b) $$
If $\mathbf{J}_{\mathbf{V}_0}$ denotes the Jacobian matrix of $\mathbf{R}$ evaluated at $\mathbf{V}_0$, i.e.

$$ \mathbf{J}_{\mathbf{V}_0} := \mathbf{Q}(\mathbf{V}_0, \cdot) + \mathbf{Q}(\cdot, \mathbf{V}_0) + \mathbf{L}(\cdot) \quad (9) $$

then Eq. (8a) can be written $\mathbf{J}_{\mathbf{V}_0} \mathbf{V}_1 = 0$, and forms with Eq. (8b) an invertible system.

• Order $2 \le p \le n$:

$$ \mathbf{J}_{\mathbf{V}_0} \mathbf{V}_p + \sum_{i=1}^{p-1} \mathbf{Q}(\mathbf{V}_i, \mathbf{V}_{p-i}) = 0 \quad (10a) $$

$$ \mathbf{V}_p^t \mathbf{V}_1 = 0 \quad (10b) $$

The $n$ linear systems (Eq. (8) and Eq. (10)) can all be solved once the Jacobian matrix $\mathbf{J}_{\mathbf{V}_0}$ is computed and inverted.

The series is computed from $\mathbf{V}_0$ to $\mathbf{V}(a_{max})$, from which another continuation step can start (after a correction process if the residual norm is too large) to continue iteratively the solution branch. The choice of $a_{max}$ has been discussed in other papers, for example [18]; given a tolerance criterion $\varepsilon_R$ on the residual norm, a simple estimate is

$$ a_{max} = \left( \frac{\varepsilon_R}{\left\| \sum_{i=1}^{n} \mathbf{Q}(\mathbf{V}_i, \mathbf{V}_{n+1-i}) \right\|} \right)^{\frac{1}{n+1}} \quad (11) $$

### 3. Harmonic Balance Method (HBM) coupled with the ANM

A sketch of the Harmonic Balance Method (HBM) coupled with the ANM is first given in section 3.1, in order to understand the approach and its different steps. Then an original implementation of the method is detailed in section 3.2 in the periodic case. Finally the extension that allows to deal with quasi-periodic solutions is presented in section 3.3.

#### 3.1. Coupling HBM and ANM: a sketch of the method

Let us consider an autonomous system with a real parameter $\lambda$ (continuation parameter):

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

where $\mathbf{x}: t \in \mathbb{R} \mapsto \mathbf{x}(t) \in \mathbb{R}^{N_d}$ and $f: \mathbb{R}^{N_d} \times \mathbb{R} \to \mathbb{R}^{N_d}$. The principle of the Harmonic Balance Method is to replace $\mathbf{x}$ by its truncated Fourier series in Eq. (12), in order to solve a new algebraic problem on the Fourier coefficients and the pulsation.

First, it is assumed that thanks to auxiliary variables, Eq. (12) is recast quadratically as:

$$ \dot{\mathbf{m}}\dot{\mathbf{u}} = \mathbf{c}_0 + \lambda\mathbf{c}_1 + \mathbf{l}_0\mathbf{u} + \lambda\mathbf{l}_1\mathbf{u} + \mathbf{q}(\mathbf{u}, \mathbf{u}) \quad (13) $$

where $\mathbf{u}$ is the vector of all variables (including auxiliary variables if required), of length $N_{eq}$; $\mathbf{c}_0, \mathbf{c}_1 \in \mathbb{R}^{N_{eq}}$, $\mathbf{m}, \mathbf{l}_0$ and $\mathbf{l}_1$ are $N_{eq} \times N_{eq}$ matrices, and $\mathbf{q}$ is a bilinear operator, from $\mathbb{R}^{N_{eq}} \times \mathbb{R}^{N_{eq}}$ to $\mathbb{R}^{N_{eq}}$.

Then, larger operators on the complex Fourier coefficients of $\mathbf{u}$ and the pulsation are built. Finally, a system on real and imaginary parts of the Fourier coefficients plus the pulsation is formed. When read as $\mathbf{R}(\mathbf{X}, \lambda)$, it verifies the assumptions of the ANM, namely:

• there are $N+1$ unknowns and $N$ equations;

• thanks to the underlying quadratic formalism, the Jacobian matrix of $\mathbf{R}$ evaluated at $(\mathbf{X}, \lambda)$, is computed exactly.

The next section describes the periodic HBM. Though the method is already known, and its efficient coupling with the ANM (thanks to a quadratic formulation) was already described in [14], this section describes an original implementation and serves as a model for the quasi-periodic HBM, for which little modifications will be made.
## 3.2. Description of the periodic HBM coupled with the ANM

Assuming the quadratic recast of Eq. (13), the truncated Fourier series of $\mathbf{u}$ is denoted as:

$$ \mathbf{u}(t) = \sum_{k=-H}^{H} \mathbf{U}_k e^{ik\omega t} \qquad (14) $$

The problem is now to find the complex vectors $\mathbf{U}_k$ and the pulsation $\omega$. Since $\mathbf{u}(t)$ is real,

$$ \mathbf{U}_0 \in \mathbb{R}, \quad \mathbf{U}_{-k} = \overline{\mathbf{U}_k} \qquad (15) $$

Thus only $2H + 1$ real vectors need to be stored: $\mathbf{U}_0$, and $\mathrm{Re}(\mathbf{U}_k)$, $\mathrm{Im}(\mathbf{U}_k)$ (real and imaginary parts) for $k \in [1, H]$. This represents $(2H + 1)N_{eq}$ real coefficients. The pulsation $\omega$ is also an unknown, so the total number of real unknowns $N_P$ is:

$$ N_P = (2H + 1)N_{eq} + 1 \qquad (16) $$

With $\mathbf{U}$ defined as $\mathbf{U} = [\mathbf{U}_0 \mathbf{U}_1 \dots \mathbf{U}_H]^t$, operators acting on $\mathbf{U}$, $\omega$, $\lambda$ are now defined to transform Eq. (13) into a quadratic system on these unknowns. First of all, since:

$$ \dot{\mathbf{u}}(t) = \sum_{k=-H}^{H} ik\omega \mathbf{U}_k e^{ik\omega t} = i\omega \sum_{k=-H}^{H} k \mathbf{U}_k e^{ik\omega t} \qquad (17) $$

the left-hand side of Eq. (13) is transformed into $i\omega\mathbf{M}\mathbf{U}$, where $\mathbf{M}$ is defined as a diagonal block matrix $\mathrm{Diag}(km, 0 \le k \le H)$. Secondly, operators $\mathbf{c}_0$ and $\mathbf{c}_1$ do not act on Fourier series, so $\mathbf{C}_0$ (respectively $\mathbf{C}_1$) is defined by concatenation of $\mathbf{c}_0$ (resp. $\mathbf{c}_1$) and $HN_{eq}$ zeros. Thirdly, linear operators $\mathbf{l}_0$ and $\mathbf{l}_1$ act linearly on Fourier series, so $\mathbf{L}_0$ (resp. $\mathbf{L}_1$) is defined as a diagonal block matrix $\mathrm{Diag}(l_0, 0 \le k \le H)$ (resp. $\mathbf{l}_1$), i.e. $H + 1$ times the same block. Finally, the Fourier series of a product of two variables can be obtained by convolution. This approach differs from [14], where the authors give the explicit expression for the Fourier coefficients of variables products. The advantage of the approach below is that it will require little adaptation for the quasi-periodic case. A function (denoted as **Conv**) is dedicated to computing the Fourier series of a product, and reads as follows:

### **Conv**

1. Entries: $\mathbf{U}^i, \mathbf{U}^j$ Fourier coefficients (from $k = 0$ to $H$) of the $i$-th and $j$-th variables of $\mathbf{u}$

2. Recreate the respective coefficients vectors $\mathbf{U}^i, \mathbf{U}^j$ by conjugation and symmetry (Eq. (15)).

3. Compute the central part of the convolution of $\mathbf{U}^i$ and $\mathbf{U}^j$: for two vectors $\mathbf{U}^i$ and $\mathbf{U}^j$ of length $2H + 1$, the convolution has length $4H + 1$. Only its central part of length $2H + 1$ is kept.

4. Return the vector starting from the central one (constant coefficient)

The third step of the function **Conv** may use existing convolution procedures efficiently. For instance, in Matlab, it is performed through the instruction `conv(U^i, U^j, 'same')`. Then an operator $\mathbf{Q}_0$ is defined, giving the Fourier series of each variables product specified in $\mathbf{q}(\mathbf{u}, \mathbf{u})$ thanks to **Conv**.

The system can be autonomous or forced. If it is autonomous, when $\mathbf{u}$ is solution then $\mathbf{u}(\cdot + \tau)$ is also solution. Another equation is then added to determine a unique solution, called phase condition. Examples of such conditions include: a variable is null at $t = 0$; the first sine coefficient of a variable is null. They often read as a simple linear combination of complex Fourier coefficients, for instance $\mathrm{Im}(\mathbf{U}_1^1) = 0$. Given a solution $\mathbf{u}$, since:

$$ \mathbf{u}(t + \varphi) = \sum_{k=-H}^{H} \mathbf{U}_k e^{ik\omega(t+\varphi)} = \sum_{k=-H}^{H} e^{ik\omega\varphi} \mathbf{U}_k e^{ik\omega t} \qquad (18) $$
a phase condition of the form Im($\mathbf{U}_1^1$) = 0 is satisfied exactly with the right choice of $\varphi$. If a solution has been obtained, its Fourier series can be time-shifted exactly in order to satisfy the required phase condition. This will not be the case with quasi-periodic solutions (section 3.3.4).

If the system is forced, there are two ways to deal with it : including the forcing term in the constant operator $c_0$ since it does not depend on the variables in $\mathbf{u}$, or adding a nonlinear differential system that admits the forcing term as a unique solution. More details are given in section 4.

Either for autonomous or forced systems, the complex operators are combined:

$$ \mathbf{Z} := (\mathbf{U}, \omega, \lambda) \qquad (19a) $$

$$ \mathbf{C}_{\mathcal{C}} := \mathbf{C}_0 \qquad (19b) $$

$$ \mathbf{L}_{\mathcal{C}}(\mathbf{Z}) := \lambda \mathbf{C}_1 + \mathbf{L}_0 \mathbf{U} \qquad (19c) $$

$$ \mathbf{Q}_{\mathcal{C}}(\mathbf{Z}, \mathbf{Z}) := \mathbf{Q}_0(\mathbf{U}, \mathbf{U}) - i\omega\mathbf{M}\mathbf{U} \qquad (19d) $$

$$ \mathbf{R}_{\mathcal{C}}(\mathbf{Z}) := \mathbf{C}_{\mathcal{C}} + \mathbf{L}_{\mathcal{C}}(\mathbf{Z}) + \mathbf{Q}_{\mathcal{C}}(\mathbf{Z}, \mathbf{Z}) \qquad (19e) $$

and so Eq. (13) now reads for periodic solutions $\mathbf{R}_{\mathcal{C}}(\mathbf{Z}) = 0$. Let $\mathbf{X}$ the vector of real unknowns be defined as:

$$ \mathbf{X} = [\mathbf{U}_0 \operatorname{Re}(\mathbf{U}_1) \dots \operatorname{Re}(\mathbf{U}_H) \operatorname{Im}(\mathbf{U}_1) \dots \operatorname{Im}(\mathbf{U}_H) \omega]^t \in \mathbb{R}^{N_P}, \qquad (20) $$

the real quadratic system $\mathbf{R}(\mathbf{X}, \lambda) = 0$ for the ANM is defined as the collection of equations on constant, cosine and sine Fourier coefficients plus the phase condition:

$$ \mathbf{R}(\mathbf{X}, \lambda) = \begin{pmatrix} \text{lines 1 to } N_{eq} \text{ of Re } (\mathbf{R}_{\mathcal{C}}(\mathbf{Z})) & \text{Size :} \\ \text{lines } N_{eq} + 1 \text{ to } (H+1)N_{eq} \text{ of Re } (\mathbf{R}_{\mathcal{C}}(\mathbf{Z})) & N_{eq} \\ \text{lines } N_{eq} + 1 \text{ to } (H+1)N_{eq} \text{ of Im } (\mathbf{R}_{\mathcal{C}}(\mathbf{Z})) & HN_{eq} \\ \text{phase condition} & HN_{eq} \\ \multicolumn{2}{c}{1} \end{pmatrix} \qquad (21) \\ \text{total size : } (2H+1)N_{eq} + 1 = N_P $$

Similarly, the Jacobian matrix $\mathbf{J}_V$ of $\mathbf{R}$ evaluated at $\mathbf{V} := (\mathbf{X}, \lambda)$ is computed in two steps, first in its complex form, then sorting real and imaginary parts.

### 3.3. Quasi-periodic HBM

Similarly to the periodic HBM, the differential system is assumed to fit in the formalism of Eq. (13), recalled here :

$$ \mathbf{m}\dot{\mathbf{u}} = \mathbf{c}_0 + \lambda\mathbf{c}_1 + \mathbf{l}_0\mathbf{u} + \lambda\mathbf{l}_1\mathbf{u} + \mathbf{q}(\mathbf{u}, \mathbf{u}) \qquad (13) $$

and the quasi-periodic HBM looks for the solution as a double Fourier sum:

$$ \mathbf{u}(t) = \sum_{k_1=-H}^{H} \sum_{k_2=-H}^{H} \mathbf{U}_{k_1, k_2} e^{i(k_1\omega_1 + k_2\omega_2)t} \qquad (22) $$

In this development, since $\mathbf{u}(t)$ is real,

$$ \mathbf{U}_{0,0} \in \mathbb{R}, \quad \mathbf{U}_{-k_1, -k_2} = \overline{\mathbf{U}_{k_1, k_2}} \qquad (23) $$

Several steps of implementation are now presented as they differ from the periodic case : the enumeration of Fourier coefficients, the derivation operator, the computation of variables products, and the phase conditions.

#### 3.3.1. Enumeration of coefficients

Thanks to Eq. (23), only half of the coefficients plus the central one need to be stored. More precisely, in the indexes table $k_1, k_2 \in [-H, H]$, half of the table can be found applying complex conjugation and central symmetry (with respect to the coefficient $0,0$) to the other half. Then, let the letter $C$ represent coefficients that are left out ; the right half of the table is kept, starting from the central coefficient, enumerating column after column:
<table><tr><td>C</td><td>...</td><td>C</td><td>1, -H</td><td>H, -H</td></tr><tr><td>⋮</td><td></td><td>⋮</td><td>↓</td><td>↓</td></tr><tr><td></td><td></td><td>C</td><td></td><td></td></tr><tr><td>C</td><td></td><td><b>0,0</b></td><td>⋮</td><td>⋮</td></tr><tr><td>⋮</td><td></td><td>0, 1</td><td></td><td></td></tr><tr><td></td><td></td><td>↓</td><td>↓</td><td>↓</td></tr><tr><td>C</td><td></td><td>0, H</td><td>1, H</td><td>H, H</td></tr></table>

Then $\mathbf{U}$ is defined as

$$ \mathbf{U} = [ \underbrace{\mathbf{U}_{0,0}}_{\text{one block of size } N_{eq}} \underbrace{\mathbf{U}_{0,1} \dots \mathbf{U}_{0,H}}_{H \text{ blocks of size } N_{eq}} \underbrace{(\mathbf{U}_{k_1,k_2})_{1 \le k_1 \le H \atop -H \le k_2 \le H}}_{H(2H+1) \text{ blocks of size } N_{eq}} ]^t \quad (24) $$

The vector of all real unknowns, denoted $\mathbf{X}$, is

$$ \mathbf{X} = [\mathbf{U}_{0,0} \ \mathrm{Re}(\mathbf{U}_{0,k_2})_{1 \le k_2 \le H} \ \mathrm{Re}(\mathbf{U}_{k_1,k_2})_{1 \le k_1 \le H \atop -H \le k_2 \le H} \ \mathrm{Im}(\mathbf{U}_{0,k_2}) \ \mathrm{Im}(\mathbf{U}_{k_1,k_2}) \ \omega_1 \ \omega_2]^t \in \mathbb{R}^N, \quad (25) $$

$$ \text{where } N = \underbrace{N_{eq}}_{0,0} + 2 \left( \underbrace{H N_{eq} + \underbrace{H(2H+1) N_{eq}}_{k_1, k_2}}_{\text{pulsations}} + \underbrace{2}_{\text{pulsations}} \right) = (4(H+1)H+1)N_{eq} + 2 \sim 4N_{eq}H^2 \quad (26) $$

Since the number of unknowns grows as $H^2$, increasing the number of harmonics (to improve the accuracy of the solution, for example) is much more expensive in computation time than in the periodic case (where $N_P \sim 2N_{eq}H$).

### 3.3.2. Derivation operator

To define a derivation operator $\mathbf{M}$ so that $\mathbf{M}\mathbf{U}$ is the Fourier series of $\mathbf{m}\dot{\mathbf{u}}$, one recalls that:

$$ \dot{\mathbf{u}}(t) = \sum_{k_1=-H}^{H} \sum_{k_2=-H}^{H} i(k_1\omega_1 + k_2\omega_2) \mathbf{U}_{k_1,k_2} e^{i(k_1\omega_1 + k_2\omega_2)t} \quad (27) $$

thus the following diagonal block matrix $\mathbf{M}$ is defined as follows:

$$ \mathbf{M} = i \mathrm{Diag} (0\mathbf{m} \ \omega_2\mathbf{m} \ \dots \ H\omega_2\mathbf{m} \ (\omega_1 - H\omega_2)\mathbf{m} \ \dots \ (\omega_1 + H\omega_2)\mathbf{m} \ \dots \ (H\omega_1 + H\omega_2)\mathbf{m}) \quad (28) $$

One notices that $\mathbf{M}$ depends on $\omega_1$ and $\omega_2$. Like in the periodic case, where $\mathbf{m}\dot{\mathbf{u}}$ is turned into $i\omega\mathbf{M}\mathbf{U}$, it is more efficient to leave the unknowns $\omega_1$ and $\omega_2$ out of constant matrices that can thus be computed once and combined when needed. That is why two constant, diagonal block matrices $\mathbf{M}_1, \mathbf{M}_2$ are defined as follows:

$$ \mathbf{M}_1 = \mathrm{Diag}(\underbrace{0\mathbf{m} \dots 0\mathbf{m}}_{H+1 \text{ blocks of size } N_{eq}} \underbrace{1\mathbf{m} \dots 1\mathbf{m}}_{2H+1 \text{ blocks}} \dots \underbrace{H\mathbf{m} \dots H\mathbf{m}}_{2H+1 \text{ blocks}}) \quad (29a) $$

$$ \mathbf{M}_2 = \mathrm{Diag}(\underbrace{0\mathbf{m} \dots H\mathbf{m}}_{H+1 \text{ blocks}} \underbrace{-H\mathbf{m} \dots H\mathbf{m}}_{2H+1 \text{ blocks}} \dots \underbrace{-H\mathbf{m} \dots H\mathbf{m}}_{2H+1 \text{ blocks}}) \quad (29b) $$

where the sequence of $2H+1$ blocks, $-H\mathbf{m} \dots H\mathbf{m}$, is repeated $H$ times

so that the Fourier series of $\mathbf{m}\dot{\mathbf{u}}$ is $(i\omega_1\mathbf{M}_1 + i\omega_2\mathbf{M}_2)\mathbf{U}$.
### 3.3.3. Computation of variables products

Similarly to the periodic case, the double Fourier series of a product of two variables can be computed through a two-dimensional convolution of the matrices of their Fourier coefficients. The function **Conv2** dedicated to computing the Fourier series of a product reads almost the same as in the periodic case:

**Conv2**

1. Entries: $U^i, U^j$ vector of Fourier coefficients of the i-th and j-th variables of $\mathbf{u}$

2. Recreate the respective coefficients matrices $U^i, U^j$ by conjugation and symmetry (Eq. (23))

3. Compute the central part of the convolution of $U^i$ and $U^j$

4. Return the vector version of this list of coefficients, starting from the central one (Eq. (24))

Once again, the third step may use existing convolution procedures efficiently. For instance, the 2D convolution of matrices is performed in Matlab through the instruction `conv2(U^i, U^j, 'same')`, and it is fast (compared, for instance, to the linear algebra computations involved for continuation, on $N \times N$ systems).

### 3.3.4. Phase conditions

In the periodic case, for an autonomous system, the pulsation of the solution is unknown, and one equation, called phase condition, is added to the system. As explained above, a phase shift of a solution gives a solution verifying exactly the phase condition. For a quasi-periodic solution, the two incommensurable pulsations $\omega_1$ and $\omega_2$ are unknown, so two equations must be added. These are supposed to affect the j-th and k-th variables and to be for instance of the form: $\mathrm{Im}(\mathbf{U}_{1,0}^j) = 0, \mathrm{Im}(\mathbf{U}_{0,1}^k) = 0$. One could consider the real parts instead, or any linear combination of real and imaginary parts. Let us assume a solution $\mathbf{u}$ is known. A time shift now reads:

$$ \mathbf{u}(t + \varphi) = \sum_{k_1=-H}^{H} \sum_{k_2=-H}^{H} \mathbf{U}_{k_1,k_2} e^{i(k_1\omega_1 + k_2\omega_2)(t+\varphi)} = \sum_{k_1=-H}^{H} \sum_{k_2=-H}^{H} e^{i(k_1\omega_1 + k_2\omega_2)\varphi} \mathbf{U}_{k_1,k_2} e^{i(k_1\omega_1 + k_2\omega_2)t} \quad (30) $$

The coefficient $\mathbf{U}_{1,0}$ is multiplied by $e^{i\omega_1\varphi}$, and $\mathbf{U}_{0,1}$ is multiplied by $e^{i\omega_2\varphi}$. When the two pulsations are assumed incommensurable two remarks can be proved:

1. for all $\varepsilon > 0$, there exists a phase shift $\varphi$ so that

$$ |\mathrm{Im}(e^{i\omega_1\varphi}\mathbf{U}_{1,0}^j)| < \varepsilon, \quad |\mathrm{Im}(e^{i\omega_2\varphi}\mathbf{U}_{0,1}^k)| < \varepsilon \quad (31) $$

2. almost surely, there is no phase shift that implies the exact phase conditions.

*Demonstration.* Let $\theta_{1,0}$ (resp. $\theta_{0,1}$) denote the argument of $\mathbf{U}_{1,0}^j$ (resp. $\mathbf{U}_{0,1}^k$).

$$ \exists \varphi \in \mathbb{R}, \quad \begin{cases} \mathrm{Im}(\mathbf{U}_{1,0}^j e^{i\omega_1\varphi}) = 0 \\ \mathrm{Im}(\mathbf{U}_{0,1}^k e^{i\omega_2\varphi}) = 0 \end{cases} \Leftrightarrow \exists \varphi \in \mathbb{R}, \quad \begin{cases} \theta_{1,0} + \omega_1\varphi \equiv 0[\pi] \\ \theta_{0,1} + \omega_2\varphi \equiv 0[\pi] \end{cases} \quad (32) $$

$$ \Leftrightarrow \exists n_{1,0}, n_{0,1} \in \mathbb{Z}, \quad \varphi = \frac{n_{1,0}\pi - \theta_{1,0}}{\omega_1} = \frac{n_{0,1}\pi - \theta_{0,1}}{\omega_2} \quad (33) $$

$$ \Leftrightarrow \exists n_{1,0}, n_{0,1} \in \mathbb{Z}, \quad n_{1,0}\pi\omega_2 - \theta_{1,0}\omega_2 = n_{0,1}\pi\omega_1 - \theta_{0,1}\omega_1 \quad (34) $$

$$ \Leftrightarrow (\theta_{1,0}; \theta_{0,1}) \in \mathcal{A}, \text{ where} \quad (35) $$

$$ \mathcal{A} = \left\{ \left( x ; \frac{\omega_2}{\omega_1}x + \frac{n_{0,1}\pi\omega_1 - n_{1,0}\pi\omega_2}{\omega_1} \right), x \in \mathbb{R}, n_{1,0}, n_{0,1} \in \mathbb{Z} \right\} \subset \mathbb{R}^2 $$

The set $\mathcal{A}$ is a dense, null subset of $\mathbb{R}^2$. This can be seen more clearly using the linear map $(x, y) \mapsto (x, y - \frac{\omega_2}{\omega_1}x)$. This map transforms $\mathcal{A}$ into the product of $\mathbb{R}$ and a dense, countable subset of $\mathbb{R}$. Therefore,
Eq. (32) has no solution, almost surely. However, for any $(\theta_{1,0}; \theta_{0,1}) \in \mathbb{R}^2$, there are elements in $\mathcal{A}$ arbitrarily close to them.

These two remarks highlight a theoretical difference between the periodic and the quasi-periodic cases. However, it is not entirely relevant to numerical applications, where the ratio of the two pulsations, $\frac{\omega_1}{\omega_2}$, is bound to be rational (since the pulsations are rational too). More details will be given when necessary (see for instance section 4.2). In examples studied here, to find an initial solution verifying phase conditions, the authors found that a solution can be shifted to verify them approximately, in the sense of Eq. (31), provided that $\varepsilon$ is not chosen too small. Corrections (e.g. Newton-Raphson corrections) can then be applied to provide a solution, at any given tolerance on the residual norm.

# 4. Forced Van der Pol oscillator

## 4.1. Quadratic formulations, continuation of periodic solution

An example of a forced system that exhibits periodic and quasi-periodic solutions is a forced Van der Pol oscillator:

$$ \ddot{x} - \mu_1 \dot{x} + \mu_2 x \dot{x} + \mu_3 x^2 \dot{x} + a_1 x = \cos(\lambda t) \quad (36) $$

Values of dimensionless parameters are chosen as: $\mu_1 = \mu_2 = 0.1$, $\mu_3 = a_1 = 1$. A first quadratic formulation, emphasizing the operators of Eq. (13), is:

$$ \begin{array}{l@{\qquad}l@{\qquad}l} \dot{x} = 0 & + y & + 0 & (37a) \\ \dot{y} = \cos(\lambda t) & + \mu_1 y - a_1 x & - \mu_2 xy - \mu_3 yz & (37b) \\ \underbrace{0}_{\mathbf{m}\mathbf{u}} = \underbrace{0}_{\mathbf{c}_0} & \underbrace{+ z}_{\mathbf{l}_0 \mathbf{u}} & \underbrace{- x^2}_{\mathbf{q}(\mathbf{u},\mathbf{u})} & (37c) \end{array} $$

Operators $c_1$ and $l_1$ are null. Similarly to example 4 in [14], the forcing term $\cos(\lambda t)$ is placed in the second line of the constant operator. It means that its constant amplitude (1) is added to the second equation on $\cos(\omega t)$ (line $N_{eq} + 2$ in the system of Eq. (21)). The additional equation is then simply: $\omega = \lambda$. The drawback is that if the system, forced at a pulsation $\lambda$, is able to adopt a pulsation $\omega \neq \lambda$, the additional condition must be changed manually to reflect that: an *a priori* knowledge of the solution becomes necessary.

Another way (cited in [19]) of including a forcing term is to consider the following nonlinear differential system:

$$ \dot{c} = c + \lambda s - c(c^2 + s^2) \quad (38a) $$

$$ \dot{s} = s - \lambda c - s(c^2 + s^2) \quad (38b) $$

with the phase equation $s(t = 0) = 0$, admitting $c = \cos(\lambda t)$, $s = -\sin(\lambda t)$ as a unique stable solution in the steady state regime. Hence, a second possible quadratic formulation is:

$$ \begin{array}{l@{\qquad}l@{\qquad}l} \dot{x} = y & + 0 & + 0 & (39a) \\ \dot{y} = c + \mu_1 y - a_1 x & + 0 & - \mu_2 xy - \mu_3 yz & (39b) \\ 0 = z & + 0 & - x^2 & (39c) \\ \dot{c} = c & + \lambda s & - cr & (39d) \\ \dot{s} = s & - \lambda c & - sr & (39e) \\ \underbrace{0}_{\mathbf{m}\mathbf{u}} = \underbrace{-r}_{\mathbf{l}_0 \mathbf{u}} & \underbrace{+0}_{\lambda \mathbf{l}_1 \mathbf{u}} & \underbrace{+c^2 + s^2}_{\mathbf{q}(\mathbf{u},\mathbf{u})} & (39f) \end{array} $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Forced Van der Pol system (Eq. (36)), continuation (with $H = 25$) of the periodic solution that has the same pulsation as the forcing pulsation, peak-to-peak amplitude of $x$ versus the forcing pulsation. Solid line indicates a stable solution, dashed line, unstable. The dot with letters NS indicates a Neimark-Sacker bifurcation.</figcaption></figure>

In this second formulation, constant operators $c_0$ and $c_1$ are null. Both formulations (Eq. (37), (39)) can be used for the continuation of periodic or quasi-periodic solutions. However, in the second one, variables $c$, $s$ and $r$ must be developed in Fourier series although they are expected to have only one non-zero coefficient.

The periodic solution at pulsation $\omega = \lambda$ is continued, and the bifurcation diagram is presented in fig. 1. The stability is computed thanks to Hill's method [20]. The solution is stable for low forcing pulsations, then a Neimark-Sacker bifurcation occurs at $\lambda \approx 1.798$, the solution becomes unstable and there is an emerging quasi-periodic solution branch that will be continued.

The quadratic formulation of Eq. (37) is used to continue the quasi-periodic solution, now with the following additional equations:

$$ \lambda = \omega_1, \quad \mathrm{Im}(x_{0,1}) = 0 \qquad (40) $$

The first condition is the same forcing condition as before, while the second condition means: the coefficient of $\sin((0\omega_1 + 1\omega_2)t)$ for the variable $x$ is null. In fig. 2, the coefficient of $\cos(0\omega_1 + 1\omega_2)t$ of the variable $x$ (i.e. $\mathrm{Re}(x_{0,1})$) is plotted versus the forcing pulsation $\omega_1$. Following the concept of a secondary Hopf bifurcation, in the vicinity of the Neimark-Sacker bifurcation, the new (quasi-periodic) branch diverges from the previous periodic branch with a first Fourier coefficient growing as a square root.

A criterion proves useful to verify the accuracy of the solution branch avoiding costly comparisons to time integration. One can use an energy criterion: how many terms in the Fourier series must be retained to find an arbitrary percentage, e.g. 99.5 %, of the 2-norm of a reference solution? For a variable, e.g. $x$, let us define the 2-norms of the partial sums:

$$ \forall k \in [0, H[, \quad \|x\|_k^2 = \sum_{k_1=-k}^{k} \sum_{k_2=-k}^{k} |x_{k_1,k_2}|^2, \quad \|x\| = \sum_{k_1=-H}^{H} \sum_{k_2=-H}^{H} |x_{k_1,k_2}|^2 \qquad (41) $$

and the threshold $H_x \in [1, H]$ so that:

$$ \|x\|_{H_x}^2 \ge \frac{99.5}{100} \|x\|^2, \quad \|x\|_{H_x-1}^2 < \frac{99.5}{100} \|x\|^2 \qquad (42) $$

Then, a handful of behaviors are observed. If the number $H$ is large enough, the Fourier decomposition has numerically converged, and $H_x$, small before $H$, represents the number of harmonics needed to produce a (likely) relevant solution. For the forced Van der Pol, this criterion is applied to the solution branch computed
<figure><img src="image_3.png" /><figcaption>Figure 2: Forced Van der Pol (Eq. (36)), continuation of the quasi-periodic solution with $H = 6$, coefficient of $\cos(\omega_2 t)$ for the variable $x$ versus the continuation parameter $\lambda$, that is, the forcing pulsation $\omega_1$. Solid blue line: solution branch; blue dots: ends of continuation steps. Letters NS indicate a Neimark-Sacker bifurcation.</figcaption></figure>

with $H = 6$. It gives $H_x = 2$ for all three variables $(x, y, z)$, suggesting that an approximation of the quasi-periodic solution can be computed with $H = 2$. As an illustration, fig. 3 presents solutions computed with $H = 2$ and compares them with time integration (using Matlab ODE suite [21]) on a long time interval. By long, the authors mean with respect to the meaningful time scales here, the two periods $\frac{2\pi}{\omega_1}$ and $\frac{2\pi}{\omega_2}$. It is observed that the solution with $H = 2$ may be slightly incorrect pointwise, which means that each curve (HBM and time integration) can be distinguished in the phase diagram: they leave slightly different areas empty; they have different self-intersections and so, different moiré patterns. But the approximate shape of the solution is obtained (enveloppe in the phase diagram, time between two extrema).

If $H_x$ is close to $H$ it is observed that:

* either $H$ should be increased to compute accurately the solution;

* or the branch is going through a strong frequency-locking region, which will be studied in section 4.2.

## 4.2. Frequency-locking regions

As the branch of quasi-periodic solution is continued, the pulsation $\omega_2$ varies slightly but stays close to the natural pulsation $a_1 = 1$. The pulsation ratio $\frac{\omega_1}{\omega_2}$ is thus roughly equal to $\omega_1 = \lambda$, and so a difficulty arises: what happens when the pulsation ratio is rational? It is a major assumption of the method that $\omega_1$ and $\omega_2$ are incommensurable, so one could expect the method to fail. Indeed, if there are $p, q \in \mathbb{N}^*$ such that $\frac{\omega_1}{\omega_2} = \frac{p}{q}$, it can introduce redundancy in the Fourier decomposition, meaning:

$$ \exists k_1, k_2, h_1, h_2 \in [-H, H], \quad (k_1, k_2) \neq (h_1, h_2), \quad k_1\omega_1 + k_2\omega_2 = h_1\omega_1 + h_2\omega_2 \quad (43) $$

and then the problem is ill-posed: at least two coefficients in the Fourier series are assigned to the same frequency. For a given ratio $\frac{\omega_1}{\omega_2} = \frac{p}{q}$, Eq. (43) admitting solutions is equivalent to:
<figure><img src="image_4.png" /><figcaption>Figure 3: Forced Van der Pol, two examples of solutions from the quasi-periodic solution branch computed with $H = 2$. From a) to d) : for $\lambda \simeq 4.29$, quasi-periodic HBM ($H = 2$) in blue, time integration (ODE suite from Matlab) in red. a, b) Phase diagram ($x, \dot{x} = y$). c, d) Time comparison of the quasi-periodic HBM ($H = 2$) and time integration. From e to g) : similarly for $\lambda \simeq 2.75$. e, f) Phase diagram. g, h) Time comparison.</figcaption></figure>

$$ \exists k_1, k_2, h_1, h_2 \in [-H, H], \quad (k_1, k_2) \neq (h_1, h_2), \quad (k_1 - h_1)\omega_1 + (k_2 - h_2)\omega_2 = 0 \quad (44) $$

$$ \Leftrightarrow \exists k_1, k_2, h_1, h_2 \in [-H, H], \quad (k_1, k_2) \neq (h_1, h_2), \quad (k_1 - h_1)p\omega_1 + (k_2 - h_2)q\omega_2 = 0 \quad (45) $$

$$ \Leftrightarrow \exists k_1 \neq h_1, k_2 \neq h_2 \in [-H, H], \quad \frac{h_2 - k_2}{k_1 - h_1} = \frac{p}{q} \quad (46) $$

$$ \Leftrightarrow \exists j_1, j_2 \in [-2H, 2H] \setminus \{0\}, \quad \frac{j_1}{j_2} = \frac{p}{q} \quad (47) $$

Thus, for a given $H$, there is a finite number of potentially problematic values $\frac{j_1}{j_2}, j_1, j_2 \in [-2H, 2H]$. Displayed on the same graph as the pulsation ratio (fig. 4 with $H = 6$), they enlighten the continuation process : in their vicinity, the ANM reduces the step size to deal appropriately with ill-conditioned systems. We observe in this example that the major part of the issue is caused by the rationals $\frac{j_1}{j_2}$ with small integers $j_1$ and $j_2$, namely, in $[-H, H]$ instead of $[-2H, 2H]$. Only two other values, $\frac{7}{2}$ and $\frac{7}{3}$, called here additional values, engender a shortening of the continuation steps. However, it is less important, and solutions in these regions agree with time integration (not shown here). The choice of $H$ is essential: for a given ratio $\frac{p}{q}$, choosing a smaller number $H$ can remove solutions in Eq. (43). The energy criterion (Eq. (42)) can then be used to confirm the relevance of this choice. For instance, it was shown in a previous section that instead of $H = 6$, the solution branch can be computed with $H = 2$. Then, as shown in fig. 5, continuation steps are significantly shortened only around $\lambda = 2$ and 3.

The most important cases of frequency-locking observed on this system are subharmonic resonances, i.e.

$$ \exists q \in [1, 2H], \quad \frac{\omega_1}{\omega_2} = q \quad \text{i.e.} \quad \omega_2 = \frac{\omega_1}{q} \quad (48) $$
a)

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

b)

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

c)

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

Figure 4: Forced Van der Pol, quasi-periodic solution branch computed with $H = 6$. a) Pulsation ratio $\omega_1/\omega_2$ versus the forcing pulsation $\omega_1$. Solid blue line: solution branch; red dots: ends of continuation steps; dashed black lines: rational values $\frac{j_1}{j_2}, j_1$, $j_2 \in [-H, H]$; mixed black lines: additional values, $\frac{7}{2}$ and $\frac{7}{3}$. b) Zoom window: around the subharmonic resonances $\omega_2 = \frac{\omega_1}{2}$, pulsation ratio computed by the quasi-periodic HBM (blue, solid) and by the periodic HBM (red, dashed). c) Zoom window: similar to b), around $\omega_2 = \frac{\omega_1}{3}$.
<figure><img src="image_8.png" /><figcaption>Figure 5: Forced Van der Pol, quasi-periodic solution branch computed with $H = 2$. Pulsation ratio $\omega_1/\omega_2$ versus the forcing pulsation $\omega_1$. Solid blue line: solution branch; red dots: ends of continuation steps; dashed black lines emphasize the two rational values (integers 2 and 3) where a reduction in continuation step size is observed.</figcaption></figure>

Then any combination $k_1\omega_1 + k_2\omega_2$ can be read as $(qk_1 + k_2)\omega_2$, and so the development in double Fourier series is likely to be highly redundant. Thanks to the robustness of the ANM, the continuation of the branch through these values ($q = 2, 3, 4$ for $H = 6$) is possible. For a given integer $q$, the size of the continuation steps, though small, stops decreasing and there is a continuation step crossing the boundary, the pulsation ratio being either superior or inferior to the rational number at the endpoints of the step. New periodic solutions that were not displayed on the bifurcation diagram of the periodic solution $\omega = \lambda$ (fig. 1) can now be found easily and investigated.

Figure 6 a, shows that in the vicinity of the subharmonic resonance $\lambda = 2$, the quasi-periodic HBM gives qualitatively good results, with an approximate shape of the solution in the phase diagram $(x, \dot{x})$. However, time integration starting from this solution converges to a slightly different periodic solution. A comparison is made between three versions of the solution: as obtained with the quasi-periodic method; with the periodic method; steady-state of time integration. The two last ones are superposed on the figure.

In case the first quadratic formulation was chosen (Eq. (37)) the additional condition has to be modified into $\omega = \frac{\lambda}{2}$. No change is needed if the second quadratic formulation (Eq. (39)) is used. A short isolated solution branch is continued. Presented in fig. 6 the heart-shaped bifurcation diagram exhibits two solutions, one stable, one unstable, for $\lambda$ comprised between 1.982 and 2.016. Therefore, an ideal continuation of the quasi-periodic branch should exhibit a plateau $\frac{\omega_1}{\omega_2} = 2$, for this interval of $\lambda$. Although this cannot be achieved by definition (redundancy in Fourier decomposition), it is not an issue. Once an “approximate plateau” has been observed it is easy to go back to a continuation algorithm dedicated to periodic solutions in order to compute it accurately. Fig. 4 displays a (smoothed) example of a pattern known as devil’s staircase<fnref n="1" />: a monotonous function exhibiting steps of finite length for all rational values.

In a similar manner, isolated branches of periodic solutions are found when the quasi-periodic solution branch passes through $\lambda = 3$ and 4. The bifurcation diagrams of these periodic solutions are presented in figure 7. Once more, the stable periodic solution has a narrow range of existence. The case $\omega_2 = \frac{\omega_1}{4}$ is different from the previous two cases. The branch spreads over a shorter interval (a few $10^{-4}$ instead of a few $10^{-2}$), which is consistent with the theory<fnref n="2" />: larger (coprime) integers $p, q$ in a rational pulsation ratio are associated with smaller frequency-locking regions. One can notice that the existence interval does not include

<footnote n="1">¹See for instance [22], section 4.2.</footnote>

<footnote n="2">²See for instance [22], section 4.2.</footnote>
a)

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

b)

<figure><img src="image_10.png" /><figcaption>Figure 6: Forced Van der Pol. a) Phase diagram ($x, \dot{x} = y$) in the vicinity of the subharmonic resonance $\omega_2 = \frac{\omega_1}{2}$. Comparison of the quasi-periodic HBM (blue), periodic HBM (green) and time integration (red). b) Bifurcation diagram of the periodic solution $2\omega = \lambda$ (solid line: stable solution, dashed: unstable, B: bifurcation).</figcaption></figure>

the integer value $\lambda = 4$.

a)

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

b)

<figure><img src="image_12.png" /><figcaption>Figure 7: Forced Van der Pol, bifurcation diagrams for two periodic subharmonic resonances. Solid line: stable solution, dashed: unstable, B: bifurcation. a) Subharmonic resonance $3\omega = \lambda$. b) $4\omega = \lambda$.</figcaption></figure>

On this example, when the number $H$ is large enough to engender issues in the vicinity of subharmonic resonances (for instance, $H \ge 2$ around the subharmonic resonance $\omega_2 = \frac{\omega_1}{3}$), the norm of the Fourier coefficients exhibits important variation, with hills or valleys around the resonances. Figure 8 presents the 2-norm and 1-norm of $x$, around the subharmonic resonance 3, for $H = 2$ and $H = 6$. It is observed that the hills or valleys tend to be greater for $H = 6$ than for $H = 2$.

The existence of redundancy in the Fourier decomposition can also be observed directly on the matrix of coefficients. Figure 9 presents the modulus of Fourier coefficients for the variable $x$, around the subharmonic resonances 2 and 3, for $H = 6$. A block surrounding the central coefficient ($k_1 = k_2 = 0$) is almost replicated, following a simple pattern since $\omega_2 = \frac{\omega_1}{q} \Leftrightarrow \forall k_1, k_2, (k_1 + 1)\omega_1 + (k_2 - q)\omega_2 = k_1\omega_1 + k_2\omega_2$. For instance, in the vicinity of the resonance 2, the pulsations 0, $\omega_1 - 2\omega_2$, $2\omega_1 - 4\omega_2$ and $3\omega_1 - 6\omega_2$ are almost equal. The central coefficient ($k_1, k_2$) = (0, 0) and its replications $(1, -2), (2, -4)$ and $(3, -6)$ are marked with black squares. Similarly for other coefficients, we observe that a right-shift of 2 columns and up-shift of 1 row links similarly close coefficients.
<figure><img src="image_13.png" /><figcaption>Figure 8: Forced Van der Pol, norm of x versus the forcing pulsation for the quasi-periodic HBM, around the subharmonic resonance $\omega_2 = \frac{\omega_1}{3}$, for $H = 2$ (solid blue line) and $H = 6$ (dashed red line). a) 2-norm. b) 1-norm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9: Forced Van der Pol, modulus of the coefficients for the variable x, for $H = 6$, log color scale. a) Around the resonance $\omega_2 = \frac{\omega_1}{2}$. Black squares highlight the central coefficient ($k_1, k_2$) = (0,0) and its replications (1, -2), (2, -4) and (3, -6). b) Around the resonance $\omega_2 = \frac{\omega_1}{3}$. This time, black squares highlight the central coefficient (0,0) and its replications (1, -3) and (2, -6). A translucent mask is added on the left halves to remind that they are only the complex conjugate of the right halves, according to section 3.3.1. Only the right half is computed. The left half illustrates that in the double Fourier series, indices $k_1$ and $k_2$ span $[-H, H]$.</figcaption></figure>

## 5. Coupled Van der Pol oscillators

Two coupled Van der Pol oscillators constitute an autonomous system exhibiting quasi-periodicity. The equations are:

$$ \ddot{x}_1 + a_1 \dot{x}_1 + \Omega_1^2 x_1 = a_2 \lambda (\dot{x}_1 + \dot{x}_2) - a_3 \lambda (\dot{x}_1 + \dot{x}_2)(x_1 + x_2) - a_4 \lambda (\dot{x}_1 + \dot{x}_2)(x_1 + x_2)^2 \quad (49a) $$

$$ \ddot{x}_2 + b_1 \dot{x}_2 + \Omega_2^2 x_2 = b_2 \lambda (\dot{x}_1 + \dot{x}_2) - b_3 \lambda (\dot{x}_1 + \dot{x}_2)(x_1 + x_2) - b_4 \lambda (\dot{x}_1 + \dot{x}_2)(x_1 + x_2)^2 \quad (49b) $$

Values of dimensionless parameters are chosen as: $\Omega_1 = 1$, $a_1 = 0.01$, $a_2 = 0.5$, $a_3 = a_4 = 2$; $\Omega_2 = 2.5$, $b_1 = 0.025$, $b_2 = 1$, $b_3 = b_4 = 4$. A quadratic formulation, emphasizing the operators of Eq. (13), is:
a)

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

b)

<figure><img src="image_16.png" /><figcaption>Figure 10: Coupled Van der Pol oscillators, periodic bifurcation diagram. Solid line: stable solution; dashed: unstable; NS: Neimark-Sacker bifurcation. The first periodic branch, with a pulsation close to $\Omega_1$ and a great amplitude for the first oscillator $x_1$, is plotted in red. The second periodic branch, with a pulsation close to $\Omega_2$ and a great amplitude for the second oscillator $x_2$, is plotted in blue. a) Peak-to-peak amplitude of $x_1$ versus the continuation parameter $\lambda$. b) Amplitude of $x_2$ versus $\lambda$.</figcaption></figure>

$$ \dot{x}_1 = y_1 \qquad +0 \qquad +0 \qquad (50a) $$

$$ \dot{y}_1 = -a_1 y_1 - \Omega_1^2 x_1 \qquad +a_2 \lambda (y_1 + y_2) - a_3 \lambda w - a_4 \lambda v \qquad +0 \qquad (50b) $$

$$ \dot{x}_2 = y_2 \qquad +0 \qquad +0 \qquad (50c) $$

$$ \dot{y}_2 = -b_1 y_2 - \Omega_2^2 x_2 \qquad +b_2 \lambda (y_1 + y_2) - b_3 \lambda w - b_4 \lambda v \qquad +0 \qquad (50d) $$

$$ 0 = r \qquad +0 \qquad -(x_1 + x_2)^2 \qquad (50e) $$

$$ 0 = w \qquad +0 \qquad -(x_1 + x_2)(y_1 + y_2) \qquad (50f) $$

$$ \underbrace{0}_{\text{m}\bar{\text{u}}} = \underbrace{v}_{\text{l}_0 \text{u}} \qquad +0 \underbrace{\lambda l_1 \text{u}}_{\text{q(u,u)}} \qquad -(r(y_1 + y_2)) \qquad (50g) $$

In this formulation, constant operators $c_0$ and $c_1$ are null. The bifurcation diagram of the periodic solution is shown in fig. 10. The first periodic branch, with a pulsation close to $\Omega_1$ and a great amplitude for the first oscillator $x_1$, is stable for small values of $\lambda$. It has a NS bifurcation (called NS1) at $\lambda \simeq 0.35$ where it becomes unstable. The second periodic branch, with a pulsation close to $\Omega_2$ and a great amplitude for the second oscillator $x_2$, is unstable at first, between a Hopf bifurcation (at $\lambda = 0.025$) and a NS bifurcation (at $\lambda \simeq 0.03$, called NS2), where it becomes stable.

The quasi-periodic solution branch is continued using the quadratic formulation (Eq. (50)), now with the following phase conditions:

$$ \operatorname{Im}(x_{1,(0,1)}) = 0, \quad \operatorname{Im}(x_{2,(1,0)}) = 0 \qquad (51) $$

These conditions mean that the coefficient of $\sin(\omega_2 t)$ for $x_1$ and the coefficient of $\sin(\omega_1 t)$ for $x_2$ are null. Therefore, in the vicinity of the NS1 (respectively NS2) bifurcation, the rise of the quasi-periodic branch, starting from the periodic branch, is lead by a $\cos(\omega_2 t)$ for $x_1$, and by $\cos(\omega_1 t)$ for $x_2$. As shown in fig. 11, similarly to fig. 2, one notices square-root growths of these coefficients around Neimark-Sacker bifurcations. Note that around NS1 (resp. NS2) the coefficient of $\cos(\omega_1 t)$ for the variable $x_2$ (resp. the coefficient of $\cos(\omega_2 t)$ for the variable $x_1$) is not null: the periodic branch of $x_2$ (resp. $x_1$) affected by the NS bifurcation is of pulsation close to $\Omega_1$ (resp. $\Omega_2$).

Near NS2, the pulsation ratio $\frac{\omega_1}{\omega_2}$ (fig. 12) is slightly inferior to the ratio $\frac{\Omega_1}{\Omega_2} = \frac{1}{2.5}$ and starts increasing when $\lambda$ increases. The 2-norm criterion (Eq. (42)), tested with $H = 7$, suggests that $H = 3$ is sufficient to continue the branch at least from NS2 to $\lambda = 0.2$, where the number $H$ can be increased. Thus, there is no problematic frequency-locking phenomenon when the pulsation ratio $\frac{\omega_1}{\omega_2}$ goes through the rational value $\frac{2}{5}$.
a)

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

b)

<figure><img src="image_18.png" /><figcaption>Figure 11: Coupled Van der Pol oscillators, quasi-periodic solution branch. Some significant Fourier coefficients exhibit a square-root growth starting from Neimark-Sacker bifurcations. a) Coefficient of $cos(\omega_2 t)$ for $x_1$. b) Coefficient of $cos(\omega_1 t)$ for $x_2$.</figcaption></figure>

a)

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

b)

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

c)

<figure><img src="image_21.png" /><figcaption>Figure 12: Coupled Van der Pol oscillators, quasi-periodic solution branch computed with $H = 4$. a) Pulsation ratio $\omega_1/\omega_2$ versus the continuation parameter $\lambda$. Solid blue line: solution branch; red dots: ends of continuation steps. b) Modulus of the coefficients for the variable $x_1$, log color scale (see also fig. 9), at the Neimark-Sacker bifurcation NS1: close to the periodic branch, the quasi-periodic solution has few coefficients with significant amplitude, and almost all have indices of the form $(k_1, 0)$ near NS1. c) Similar to b), modulus of coefficients for $x_1$ near NS2. Significant coefficients have indices of the form $(0, k_2)$.</figcaption></figure>

Comparisons with time integration, between NS2 and $\lambda \simeq 0.07$, confirm that this portion of the branch is computed accurately with $H = 3$. For $\lambda$ between 0.07 and 0.2, comparison with time integration is not meaningful on long intervals, as the quasi-periodic solution seems to be unstable<fnref n="3" />: starting from the quasi-periodic solution, time integration eventually converges to one or the other periodic solutions found above,

<footnote n="3">³The authors prefer not to assert that as a fact, since no stability analysis has been performed.</footnote>
a)

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

b)

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

c)

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

d)

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

e)

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

f)

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

g)

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

h)

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

i)

<figure><img src="image_30.png" /><figcaption>Figure 13: Coupled Van der Pol oscillators, $\lambda = 0.36$, an example of quasi-periodic solution, comparison of different computations. a) Phase diagram ($x_1, \dot{x}_1 = y_1$), for the quasi-periodic HBM with $H = 4$. b) Diagram ($x_1, \dot{x}_1 = y_1$), $H = 10$. c, d) Comparison between the quasi-periodic HBM ($H=4$, blue solid line) and time integration (red dashed line). e, f) Comparison between $H = 10$ and time integration. g to i) Modulus of the coefficients for the variable $x_1$, log color scale. g) $H = 4$. h) $H = 10$. i) $H = 20$.</figcaption></figure>

and this is still true with a more accurate quasi-periodic solution (i.e. with $H$ greater than 3).

As the continuation process keeps moving towards increasing values of $\lambda$ (and near NS1), the 2-norm criterion and comparisons with time integration both suggest an increase of the number $H$. A low number provides only rough approximations of the solutions : in the phase diagram ($x_1, \dot{x}_1 = y_1$) (fig. 13, a) the envelope matches the one from time integration but there is an important pointwise divergence as time increases (fig. 13, c, d). Adding terms in the Fourier series quickly increases the size of the systems ($N \sim 4N_{eq}H^2$) and therefore it slows computations down. However, it is possible for this system to continue the branch with $H = 10$ on a standard laptop. Agreement between the quasi-periodic HBM and time integration is then excellent (fig. 13, e, f).
# 6. Conclusion

An original implementation of the HBM, coupled with the ANM continuation, is presented and extended to a two-frequency HBM. The continuation of branches of quasi-periodic solutions is thus made possible. As in the periodic case, the quasi-periodic HBM coupled with the ANM is able to compute and continue solutions without need for prior knowledge, i.e. the two pulsations can be unknown and can vary along the solution branch, and there is no optimized harmonic selection. The method is applied successfully on a forced system and on an autonomous system, and the results are confirmed by time integration. The accuracy of the solution can be estimated a posteriori by computing 2-norms of Fourier coefficients.

The ANM is very robust and allows the branch to be computed even when the pulsation ratio encounters simple rational values. If necessary, the periodic solution that it represents can be continued on its own. Around these points the number of terms in Fourier series can be reduced in order to improve the behavior of the algorithm.

# Acknowledgements

This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the *Investissements d'Avenir* French Government program managed by the French National Research Agency (ANR).

# References

[1] L. O. Chua, A. Ushida. Algorithms for computing almost periodic steady-state response of nonlinear systems to multiple input frequencies, IEEE Transactions on Circuits and Systems, vol. CAS-28, **10**, p. 953-971, 1981.

[2] S. L. Lau, Y. K. Cheung, S. Y. Wu. Incremental Harmonic Balance Method with multiple time scales for aperiodic vibration of nonlinear systems, Journal of Applied Mechanics, **50**, p. 871-876, 1983.

[3] C. Kaas-Petersen. Computation of quasi-periodic solutions of forced dissipative systems, Journal of Computational Physics, **58**, p. 395-408, 1985.

[4] C. Kaas-Petersen. Computation, continuation, and bifurcation of torus solutions for dissipative maps and ordinary differential equations, Physica **25D**, p. 288-306, 1987.

[5] S.-K. Choi, S. T. Noah. Response and stability analysis of piecewise-linear oscillators under multi-forcing frequencies, Nonlinear Dynamics, **3**, p. 105-121, 1992.

[6] Y.-B. Kim. Quasi-periodic response and stability analysis for non-linear systems: a general approach, Journal of Sound and Vibration, **192** (4), p. 821-833, 1996.

[7] T. M. Cameron, J. H. Griffin. An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems, Journal of Applied Mechanics, **56**, p. 149-154, 1989.

[8] N. Coudeyras, S. Nacivet, J.-J. Sinou. Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal, Journal of Sound and Vibration **328**, p. 520-540, 2009.

[9] M. Guskov, J.-J. Sinou, F. Thouverez. Multi-dimensional harmonic balance applied to rotor dynamics, Mechanics Research Communications, **35**, p. 537-545, 2008.

[10] L. Peletan, S. Baguet, M. Torkhani, G. Jacquet-Richardet. Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor-stator dynamics. Nonlinear Dynamics, **78** (**4**), p. 2501-2515, 2014.

[11] B. Rasmussen. Numerical Methods for the Continuation of Invariant Tori, PhD thesis, Georgia Institute of Technology, 2003.

[12] F. Schilder, H. M. Osinga, W. Vogt. Continuation of Quasi-periodic Invariant Tori, SIAM Journal of Applied Dynamical Systems, **4** (3), p. 459-488, 2005.
[13] F. Schilder, W. Vogt, S. Schreiber, H. M. Osinga. Fourier methods for quasi-periodic oscillations, International Journal for Numerical Methods in Engineering, **67**, p. 629-671, 2006.

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

[15] L. Azrar, B. Cochelin, N. Damil, M. Potier-Ferry. An asymptotic numerical method to compute the post-buckling behavior of elastic plates and shells, International Journal for Numerical Methods in Engineering, **36**, p. 1251-1277, 1993.

[16] B. Cochelin, N. Damil, M. Potier-Ferry. The Asymptotic Numerical Method : an efficient perturbation technique for nonlinear structural mechanics, Revue Européenne des Éléments Finis, **3**, p. 281-297, 1994.

[17] N. Damil, M. Potier-Ferry. A new method to compute perturbed bifurcations : application to the buckling of imperfect elastic structures, International Journal of Engineering Sciences N9, **28**, p. 943-957, 1990.

[18] B. Cochelin, F. Pérignon. Résolution de petits systèmes algébriques par la MAN sous Matlab. Revue Européenne des Éléments, **13** (**1-2**), p. 79-96, 2004.

[19] B. Krauskopf, H. M. Osinga, J. Galan-Vioque. Numerical Continuation Methods for Dynamical Systems: Path following and boundary value problems, Springer Netherlands jointly published with Canopus Publishing Limited, Bristol (UK), 2007.

[20] A. Lazarus, O. Thomas, A harmonic-based method for computing the stability of periodic solutions of dynamical systems. Comptes Rendus Mecanique, **338**, 510-517, 2010.

[21] L.F. Shampine, M.W. Reichelt. The Matlab ODE suite, SIAM Journal on Scientific Computing, **18**, p. 1-22, 1997.

[22] A. H. Nayfeh, B. Balachandran. Applied nonlinear dynamics : analytical, computational and experimental methods. John Wiley & Sons, Weinheim (Germany), 2008.