Courgeon Matthieu
email: matthieu.courgeon@univ-ubs.fr

Duhaut Dominique
email: dominique.duhaut@univ-ubs.fr

Artificial Companions as Personal Coach for Children: The Interactive Drums Teacher

Keywords: Artificial companions, human robot interaction, virtual humans, e-learning. ACM Classification H.1.2 User/Machine Systems, H.5.2 User-centered design, Prototyping, H.5.1 Animations, Artificial, augmented, and virtual realities, H.5.3 Asynchronous interaction

The MOCA Project that aims at designing and studying the interaction and relationship between artificial companions and children in everyday life at home activities. Artificial companions are digital embodied entities that can be either robotic or virtual. In this paper, we focus on a single activity, subpart of the whole project: a coaching application that uses two artificial companions to teach the basics of drums to children. One device is a Nao robot, the other is a virtual agent. This application offers a wide variety of interaction situation that will enable us to conduct several experiments. Future directions are discussed, in terms of user-computer interaction and in terms of coaching and e-learning applications.

INTRODUCTION

Over the last twenty years, an increasing number of digital devices were included as part of our daily life; market studies show that this dynamic will continue to accelerate. Yet most if not allof these devices do not have the capacity to establish a relationship with the user. In affective computing, this ability is considered as a key factor to ensure significant added value in everyday life. Where, in history, the tool is an extension of the hand of man, the collective of modern digital devices must form a true extension of himself in the world.

To do so, one must use this constantly evolving collective as artificial entities incarnated in each of the different devices. From the relationship created with these entities, emerge the concept of "artificial companions". Artificial companions have to be reliable, fair, unselfish, safe, fast, credible and empathetic. The MOCA project aims to develop and study this relationship through virtual agents, personal robots and communicating devices. In this article, we present a sample application showing the complementarity of Robotics and Virtual Agent: A virtual coach for children dedicated to interactive learning of drums.

RELATED WORK

Artificial companions are embodied entities that may use a robotic body or a virtual body. Embodied Conversational Agents (ECA) are commonly used in affective computing [START_REF] Picard | Affective computing[END_REF], as they display a human like form and expressive capacities (see [START_REF] Rist | A Review of the Development of Embodied Presentation Agents and Their Application[END_REF] for a review). However, in most of this system, the agent is only a service provider, like a virtual butler.

Using ECA as coaches have already been explored in several specific domains. For example, they have been used in assisting users to quit smoking [START_REF] Grolleman | Break the habit! designing an e-therapy intervention using a virtual coach in aid of smoking cessation[END_REF], training them for job interviews [START_REF] Hoque | Mach: My automated conversation coach[END_REF], and so on. Yet, few of these applications are building a personal and social link with the user. [START_REF] Fasola | A socially assistive robot exercise coach for the elderly[END_REF] created and studied a robot designed to engage elderly users in physical exercise. While it is not designed to create any social relationship with the user, the robot was effective. They also compared the real robot to its virtual counterpart displayed on a TV. Their results show 'a clear preference by older adults for the physically embodied robot coach over the virtual coach in terms of enjoyableness, helpfulness and social attraction'.

Fasola and Matarić

Lisetti et al. [START_REF] Lisetti | I can help you change! an empathic virtual agent delivers behavior change health interventions[END_REF] created a coach to help users to quit drinking that uses a dynamic empathy expressive system to establish a personal relationship with the user. They show that social behaviors improve users' motivation and engagement, and increase system's acceptability.

CALONIS [START_REF] Wilks | Calonis: An artificial companion for the care of cognitively impaired patients[END_REF] is one of the few attempt to model the social relationship with the user. CALONIS aims at providing cognitive support for cognitively disabled army veterans. Unlike previous works in this area, its design considers the relationship that emerges from interactions, mostly to ensure this relation will not interfere with the real career and his patient. The robot is being involved in a three-party relation: patient-career-robot.

Pesty and Duhaut [START_REF] Pesty | Artificial Companion: building an impacting relation[END_REF] defend the importance of building an "impacting relation" with artificial companions. They argue that this relation 'will enable believable interactive ECAs or robots to become believable impacting companions able to build a kind of lifelong relationship with the human'.

MOTIVATIONS

Artificial companions are expected to be more than simple butlers. Their interaction domain is wide and does not limit to providing service to the user. They are not specific to a task and can fulfill different type of roles. In the application that is presented in the remaining of this paper, our companions are coaches and teachers.

While real teachers cannotand must notbe replaced by virtual entities, we believe that artificial companions are appropriate as complementary learning coaches. For example, when a child is learning a musical instrument, he must practice alone at home between the lessons. In this context, having virtual companions to motivate him and give insightful feedback would be valuable.

Inducing motivation to practice any learning activity however require to create some form of relationship with the user. Unlike most previous researches using a single companion, we believe that a unique multipurpose companion is not adapted. For example, the companion dedicated to games should establish a friendship relation, but the companion dedicated to homework should not, if it is to maintain a minimum level of authority over the child. Furthermore, using several companions as a collective better allows for exploiting their complementarity.

While small personal robots, such as Nao, are designed to be attractive to children, they are rather limited in its motion dynamics. Nao for example, walks slowly and it is not capable of holding objects, except for very small and light objects. Because of these limitations, we cannot use the Nao robot to demonstrate drums gestures.

On the other end, virtual agent might be less attractive than robots [START_REF] Lisetti | I can help you change! an empathic virtual agent delivers behavior change health interventions[END_REF], as they are only present in virtual worlds. However, they can move without limitation and even faster and wider than real humans. It's not difficult to animate an agent that will play drums with enough precision. We will not discuss the realism issue of gesture dynamics in this paper, but it is clearly a challenge to make a believable virtual drummer.

The combination of both technologiesrobots and agentsis in our view an efficient approach for a drums coaching system. That is why we designed the system that we describe in the following section.

SYSTEM

The drum coaching system is composed of a game drum kit (Xbox drumband), a Nao robot, and a virtual agent managed by the MARC toolkit [START_REF] Courgeon | MARC: a framework that features emotion models for facial animation during human-computer interaction[END_REF]. The child is sitting in front of the Xbox drum kit and interact with Nao and the virtual agent on a 48inch TV. (See fig 1)

Description of User-Side Drumming System

As the MOCA project aims at designing childrencomputer interaction, we try to restrain interaction devices we use to general public game controllers. In this case, we used an Xbox MadCatz "Drums Band" controller. It is composed of 4 drums and 3 cymbals, plus a central kick pedal. Drums and cymbals are pressure sensitive. Their precision is limited but they can differentiate a strong stroke from a light one. The kick pedal however is an on-off switch.

Sound is generated within the computer. The drums controller is a xInput device that only transmit events, but does not generate sound effect. To ensure the minimal latency, we use a sampler based engine using ASIO audio. ASIO is a very low latency (max 10ms) audio system that enable a very high responsivity. Using high level libraries, such as OpenAL for example, introduces a delay of about 200ms to 300ms. The sound has thus a temporal offset that is very disturbing when playing the drums. Sounds sampler are 44khz wav files of real drums elements. Their volume is modulated by the stroke intensity.

Description of Virtual-Side Drumming System

The real drum kit is facing a TV screen on which a virtual character is sitting in front of a virtual drum kit. The virtual drums kit is designed to react on stroke event by generating sound (also using ASIO) and by animating visual element. For example, virtual cymbals oscillate softly when they are stroked. It increases the audio-visual consistency.

To facilitate the animation of the virtual character, all gesture are dynamically computed using a rhythm description sheet. These descriptions are timed sequence of "body part / drums element" pair (i.e. "Right hand / Crash Cymbal" "Left hand / Snare drum", or "Right foot / Kick Pedal")

The character animation is performed by the MARC Toolkit [START_REF] Courgeon | MARC: a framework that features emotion models for facial animation during human-computer interaction[END_REF]. The animation of MARC is keyframe based. For each possible stroke combinations, we manually created three key poses: preparation, stroke, release. The animation is then created dynamically by combining the keyframes, body parts by body parts. To increased motion credibility, the interpolation between keyframes is not linear, but a squared factor for prestrike interpolation and square root factor for post stroke interpolation. It create an acceleration effect with a speed apex on the stroke (fig. 2). On the stroke apex, we added a 16ms delay. Our software is running at 60hz. The 16ms delay is meant to ensure that the stroke impact is always visible (at least one rendering frame will render this position)

Figure 2. Stroke dynamics based on 3 keyframes.

Nao Animation

Nao animation is also keyframe based. We have predefined a set of gestures that can be used during the interaction. We used a custom animation editor and Nao is remotely controlled using the Proxy included in Nao by Aldebaran.

Do we Need a Robot?

In the application presented here, Nao does not really take part in the tutoring. However, its physical presence is used as a motivation factor. By dancing, and commenting the interaction, it supposedly creates a more complex and more fun interaction system. However, a user study must be led to validate this hypothesis.

Interaction and Scenario

For now, interaction is based on a scripted scenario. The scenario consists of an introduction, in which the two companions (virtual and robotic) introduce the activity and joke about one another. They create a convivial situation to increase children engagement. During this phase, the user is mostly an observer and the two companions are looking alternately at the child and at each other.

Once the introduction is finished, the child is engaged in a series of drums exercise that start simple and increased in time. These exercises are each described by a rhythm description sheet. They can easily be modified or new exercise can be added.

Each exercise is done as follows: the virtual drummer shows the exercise once and then starts to play is on a loop softly. During this phase, the virtual point of view switch to a "from the top" view of the virtual drum kit, making it easier to see what the virtual drummer does. Then Nao give a few specific advices (included in the exercise description) and encourage the child to try to play. During the exercise, Nao is motivating the child, detecting mistakes and providing advices accordingly. For example, detecting rhythm errors, Nao will say "Your rhythm is not good, you can try to focus on what you hear more than on what you see." (Note:

focusing on reproducing the gesture rather than the music is a common mistake that we observed).

More adaptive advices will be added as we proceed to more interaction session with children.

One of the key challenges is to preserve child's engagement in the activity, even when he fails an exercise multiple time.

To this end, we included several mechanisms. Frist, Nao is always positive. It congratulates the child and give positive comments. Then, each time an exercise is failed the expectations are lowered. The coach is thus more and more tolerant to mistakes. Most children we eventually get the exercise right and move on to the next one. If the child keeps on failing, Nao eventually stops it and the interaction moves on to the next level. On one of the exercises, Nao also start singing and dancing.

Empirical Tests with Children

We conducted a set of first empirical tests by bringing the system to classes of children at different ages. The first test was conducted on 8 and 9 y-o children. While they were really engaged in the task at first and quite fond of the Nao robot, the drumming experience was somehow difficult for them. Except for the few ones already learning music as a hobby, they could not manage to play the exercises right, which led to a quick disengagement from the task.

Our second informal test was conducted on older children, 11 and 12 y-o (fig. 3). This time, most of the child were strongly engaged in the task for the entire duration and all exercises were successfully finished. Surprisingly, even when they had already seen the scripted interaction between the two companions several times, they seemed to rather enjoy it. We informally asked a few questions. Every single one of them were interested in having such companions at home. This is, in our view, a strong informal result that validates our approach. Further experimentation must be led to fully validate our system. We had a lot of questions about what other tasks the companions could do and interact with, showing a strong interest of children for the technology.

Future Improvements

One of the main improvement that must be made is the ability to speak with the companions in natural language. The current system features a set of vocal commands, for example to restart or skip an exercise, but the vocal interaction is limited. Natural language processing is a significant challenge but could help to establish a relationship. Our experiments will thus include a Wizard of Oz approach to enable such interactions.

Children have different abilities regarding music and rhythm. We are planning on adding an auto adjustment system to the sequence of exercise, skipping exercises when the child is outperforming, and adding sub-steps when the child is underperforming.

One other approach would be to include our coaching system as part of real drumming learning. Ideally, the children would have real, usual, lessons with a teacher, and use our coaching system at home, when they practice between lessons. In order to create a meaning full system, the teacher should be able to configure the exercise remotely and get feedback on the child activity between lessons.

Designing such a system would allow us to test the real impact of our system on children learning. We could compare progress two groups of children, with and without the virtual coaches, and over a long period of time, at least several weeks.

FUTURE EXPERIMENTATION

Preliminary studies with children have shown us that they are attracted to our virtual drum teacher. However, we need to validate how this system could be beneficial on a long term basis.

We intend to conduct an experiment with 3 conditions. 1) Using only the robot, the agent being replaced by charts indicating where to play, in which order. 2) Using only the virtual drummer, showing the exercises, with minimal explanations and no motivational components, and finally 3) the whole setup. Our main hypothesis is that the combination of both is more engaging and more motivating. Consequently, it should lead to more progress in learning the drums.

We could also test the physical presence effect. By using a virtual Nao robot versus a real robot, we could compare the engagement and motivation level.

Finally, this drums coaching system will be part of a bigger MOCA experiment involving several kinds of activity (from games to homeworks) and including more companion. This experimentation will aim at evaluating the global acceptance of artificial companions by children left alone for some time in a living lab. Companions will be endowed with personality, style, and will be designed to build a relationship with the user.

Figure 1 .

 1 Figure 1. The interactive system: a Xbox Drums Band controller, a Nao robot and one of the MARC toolkit virtual agent on a virtual drum.

Figure 3 -

 3 Figure 3-Photo of a class of 12 y-o children discovering the virtual drum teacher (faces blurred for anonymity). In this class, we used a projector instead of a TV.

ACKNOWLEDGMENTS

The work described in this paper was funded by the French National Research Agency (ANR): project MOCA (ANR-12-CORD-019) and by the Cap Digital Cluster.


Louis Guillot
email: louis.guillot@ens-cachan.fr

Pierre Vigué
email: vigue@lma.cnrs-mrs.fr

Christophe Vergez
email: vergez@lma.cnrs-mrs.fr

Bruno Cochelin
email: bruno.cochelin@centrale-marseille.fr

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

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.

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 [START_REF] Chua | Algorithms for computing almost periodic steady-state response of nonlinear systems to multiple input frequencies[END_REF] 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. [START_REF] Lau | Incremental Harmonic Balance Method with multiple time scales for aperiodic vibration of nonlinear systems[END_REF]. Kaas-Petersen [START_REF] Kaas-Petersen | Computation of quasi-periodic solutions of forced dissipative systems[END_REF][START_REF] Kaas-Petersen | Computation, continuation, and bifurcation of torus solutions for dissipative maps and ordinary differential equations[END_REF] 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 [START_REF] Choi | Response and stability analysis of piecewise-linear oscillators under multi-forcing frequencies[END_REF][START_REF] Kim | Quasi-periodic response and stability analysis for non-linear systems: a general approach[END_REF].

The Alternating Frequency-Time method proposed by Cameron [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] 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 [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] and coupled it with arc-length continuation [START_REF] Guskov | Multi-dimensional harmonic balance applied to rotor dynamics[END_REF][START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor-stator dynamics[END_REF].

Finally, instead of computing a trajectory on an invariant tori, the continuation of these tori was performed independently by Rasmussen [START_REF] Rasmussen | Numerical Methods for the Continuation of Invariant Tori[END_REF] and Schilder et al. [START_REF] Schilder | Continuation of Quasi-periodic Invariant Tori[END_REF][START_REF] Schilder | Fourier methods for quasi-periodic oscillations[END_REF].

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 Asympotic 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 [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. 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.

Asymptotic Numerical Method (ANM)

Let us consider a vector-valued, smooth function :

R : R N +1 -→ R N , (X, λ) → R(X, λ) (1)
Assuming that a first regular solution (X 0 , λ 0 ) is known, the aim of the Asymptotic Numerical Method (ANM) is to follow the solution branch of the equation:

R(X, λ) = 0 (2)
with respect to λ. Without a great loss of generality R can be considered quadratic, meaning that there are

C 0 , C 1 ∈ R N , L 0 , L 1 real N × N matrices, Q 0 a bilinear application of R N × R N → R N , so that: R(X, λ) = C 0 + λC 1 + L 0 X + λL 1 X + Q 0 (X, X) (3)
As noted when the method was designed, the transformation from a general function R to the quadratic formalism required in Eq. ( 3) is possible for usual smooth functions, adding auxiliary variables [START_REF] Azrar | An asymptotic numerical method to compute the postbuckling behavior of elastic plates and shells[END_REF][START_REF] Cochelin | The Asymptotic Numerical Method : an efficient perturbation technique for nonlinear structural mechanics[END_REF][START_REF] Damil | A new method to compute perturbed bifurcations : application to the buckling of imperfect elastic structures[END_REF]. The continuation parameter is now treated as an unknown; thus, a vector of all unknowns is defined as V := (X, λ), 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 V 0 . Let V 1 be a tangent vector at V 0 and let a be defined as:

a := (V -V 0 ) t V 1 ( 4
)
Then V is developed to the order n (n = 15 or 20 in applications):

V(a) = V 0 + aV 1 + a 2 V 2 + ... + a n V n (5)
Operators are combined as follows:

C = C 0 , L(V) = λC 1 + L 0 X, Q(V, V) = λL 1 X + Q 0 (X, X) (6)
so Eq. ( 3) now reads as:

R(V) = Q(V, V) + L(V) + C (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 V i :

• Order 0 : Q(V 0 , V 0 ) + L(V 0 ) + C = 0 which is true since V 0 is a solution of Eq. (2).
• Order 1 :

Q(V 0 , V 1 ) + Q(V 1 , V 0 ) + L(V 1 ) = 0 (8a) path parameter definition (4) leads to V t 1 V 1 = 1 (8b)
If J V0 denotes the Jacobian matrix of R evaluated at V 0 , i.e.

J V0 := Q(V 0 , •) + Q(•, V 0 ) + L(•) (9)
then Eq. (8a) can be written J V0 V 1 = 0, and forms with Eq. (8b) an invertible system.

• Order 2 ≤ p ≤ n :

J V0 V p + p-1 i=1 Q(V i , V p-i ) = 0 (10a) V t p V 1 = 0 (10b)
The n linear systems (Eq. ( 8) and Eq. ( 10)) can all be solved once the Jacobian matrix J V0 is computed and inverted.

The series is computed from V 0 to 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 [START_REF] Cochelin | Résolution de petits systèmes algébriques par la MAN sous Matlab[END_REF]; given a tolerance criterion ε R on the residual norm, a simple estimate is

a max = ε R n i=1 Q(V i , V n+1-i ) 1 n+1 (11)

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.

Coupling HBM and ANM : a sketch of the method

Let us consider an autonomous system with a real parameter λ (continuation parameter):

ẋ = f (x, λ) (12)
where

x : t ∈ R → x(t) ∈ R N d and f : R N d × R → R N d .
The principle of the Harmonic Balance Method is to replace x by its truncated Fourier series in Eq. [START_REF] Schilder | Continuation of Quasi-periodic Invariant Tori[END_REF], 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:

m u = c 0 + λc 1 + l 0 u + λl 1 u + q(u, u) (13)
where u is the vector of all variables (including auxiliary variables if required), of length N eq ; c 0 , c 1 ∈ R Neq , m, l 0 and l 1 are N eq × N eq matrices, and q is a bilinear operator, from R Neq × R Neq to R Neq .

Then, larger operators on the complex Fourier coefficients of 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 R(X, λ), 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 R evaluated at (X, λ), 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 [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], this section describes an original implementation and serves as a model for the quasi-periodic HBM, for which little modifications will be made.

Description of the periodic HBM coupled with the ANM

Assuming the quadratic recast of Eq. ( 13), the truncated Fourier series of u is denoted as:

u(t) = H k=-H U k e ikωt ( 14
)
The problem is now to find the complex vectors U k and the pulsation ω. Since u(t) is real,

U 0 ∈ R, U -k = U k ( 15
)
Thus only 2H + 1 real vectors need to be stored : U 0 , and Re(U k ), Im(U k ) (real and imaginary parts) for k ∈ 1, H . This represents (2H + 1)N eq real coefficients. The pulsation ω is also an unknown, so the total number of real unknowns N P is:

N P = (2H + 1)N eq + 1 (16)
With U defined as U = [U 0 U 1 ...U H ] t , operators acting on U, ω, λ are now defined to transform Eq. ( 13) into a quadratic system on these unknowns. First of all, since:

u(t) = H k=-H ikωU k e ikωt = iω H k=-H kU k e ikωt ( 17
)
the left-hand side of Eq. ( 13) is transformed into iωMU, where M is defined as a diagonal block matrix Diag (km, 0 k H). Secondly, operators c 0 and c 1 do not act on Fourier series, so C 0 (respectively C 1 ) is defined by concatenation of c 0 (resp. c 1 ) and HN eq zeros. Thirdly, linear operators l 0 and l 1 act linearly on Fourier series, so L 0 (resp. L 1 ) is defined as a diagonal block matrix Diag (l 0 , 0 k H) (resp. 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 [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], 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: U i , U j Fourier coefficients (from k = 0 to H) of the i-th and j-th variables of u 2. Recreate the respective coefficients vectors U i , U j by conjugation and symmetry (Eq. ( 15)).

3. Compute the central part of the convolution of U i and U j : for two vectors U i and 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 Q 0 is defined, giving the Fourier series of each variables product specified in q(u, u) thanks to Conv.

The system can be autonomous or forced. If it is autonomous, when u is solution then u(• + τ ) 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 Im(U 1 1 ) = 0. Given a solution u, since:

u(t + ϕ) = H k=-H U k e ikω(t+ϕ) = H k=-H e ikωϕ U k e ikωt ( 18
)
a phase condition of the form Im(U 1 1 ) = 0 is satisfied exactly with the right choice of ϕ. 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 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:

Z := (U, ω, λ) (19a) C C := C 0 (19b) L C (Z) := λC 1 + L 0 U (19c) Q C (Z, Z) := Q 0 (U, U) -iωMU (19d) R C (Z) := C C + L C (Z) + Q C (Z, Z) (19e)
and so Eq. ( 13) now reads for periodic solutions R C (Z) = 0. Let X the vector of real unknowns be defined as:

X = [U 0 Re(U 1 ) . . . Re(U H ) Im(U 1 ) . . . Im(U H ) ω] t ∈ R N P , ( 20
)
the real quadratic system R(X, λ) = 0 for the ANM is defined as the collection of equations on constant, cosine and sine Fourier coefficients plus the phase condition:

Size :

R(X, λ) =     lines 1 to N eq of Re (R C (Z)) lines N eq + 1 to (H + 1)N eq of Re (R C (Z)) lines N eq + 1 to (H + 1)N eq of Im (R C (Z)) phase condition    
N eq HN eq HN eq 1 total size : (2H + 1)N eq + 1 = N P [START_REF] Shampine | The Matlab ODE suite[END_REF] Similarly, the Jacobian matrix J V of R evaluated at V := (X, λ) is computed in two steps, first in its complex form, then sorting real and imaginary parts.

Quasi-periodic HBM

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

m u = c 0 + λc 1 + l 0 u + λl 1 u + q(u, u) (13)
and the quasi-periodic HBM looks for the solution as a double Fourier sum:

u(t) = H k1=-H H k2=-H U k1,k2 e i(k1ω1+k2ω2)t (22)
In this development, since u(t) is real,

U 0,0 ∈ R, U -k1,-k2 = U k1,k2 (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.

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 ∈ -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:

C . . . C 1, -H H, -H . . . . . . ↓ ↓ C C 0,0 . . . . . . . . . 0, 1 ↓ ↓ ↓ C 0, H 1, H H, H
Then U is defined as

U = [ U 0,0 U 0,1 . . . U 0,H U 1,-H . . . U 1,H . . . U H,H ] t (24) = [ U 0,0 one block of size Neq (U 0,k2 ) 1 k2 H H blocks of size Neq (U k1,k2 ) 1 k1 H -H k2 H H(2H+1) blocks of size Neq ] t
The vector of all real unknowns, denoted X, is

X = [U 0,0 Re(U 0,k2 ) 1 k2 H Re(U k1,k2 ) 1 k1 H -H k2 H Im(U 0,k2 ) Im(U k1,k2 ) ω 1 ω 2 ] t ∈ R N , ( 25
)
where

N = N eq 0,0 +2 ( H 0,k2 N eq + H(2H + 1) k1,k2 N eq ) + 2 pulsations = (4(H + 1)H + 1)N eq + 2 ∼ 4N eq H 2 (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 ∼ 2N eq H).

Derivation operator

To define a derivation operator M so that MU is the Fourier series of m u, one recalls that:

u(t) = H k1=-H H k2=-H i(k 1 ω 1 + k 2 ω 2 )U k1,k2 e i(k1ω1+k2ω2)t (27)
thus the following diagonal block matrix M is defined as follows:

M = iDiag (0m ω 2 m . . . Hω 2 m (ω 1 -Hω 2 )m . . . (ω 1 + Hω 2 )m . . . (Hω 1 + Hω 2 )m) (28)
One notices that M depends on ω 1 and ω 2 . Like in the periodic case, where m u is turned into iωMU, it is more efficient to leave the unknowns ω 1 and ω 2 out of constant matrices that can thus be computed once and combined when needed. That is why two constant, diagonal block matrices M 1 , M 2 are defined as follows:

M 1 = Diag(0m . . . 0m
) (29b)
where the sequence of 2H + 1 blocks, -Hm . . . Hm, is repeated H times so that the Fourier series of m u is (iω

1 M 1 + iω 2 M 2 )U.

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 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 ×N systems).

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 ω 1 and ω 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 : Im(U j 1,0 ) = 0, Im(U k 0,1 ) = 0. One could consider the real parts instead, or any linear combination of real and imaginary parts. Let us assume a solution u is known. A time shift now reads:

u(t + ϕ) = H k1=-H H k2=-H U k1,k2 e i(k1ω1+k2ω2)(t+ϕ) = H k1=-H H k2=-H e i(k1ω1+k2ω2)ϕ U k1,k2 e i(k1ω1+k2ω2)t (30)
The coefficient U 1,0 is multiplied by e iω1ϕ , and U 0,1 is multiplied by e iω2ϕ . When the two pulsations are assumed incommensurable two remarks can be proved:

1. for all ε > 0, there exists a phase shift ϕ so that

|Im(e iω1ϕ U j 1,0 )| < ε, |Im(e iω2ϕ U k 0,1 )| < ε (31)
2. almost surely, there is no phase shift that implies the exact phase conditions.

Demonstration. Let θ 1,0 (resp. θ 0,1 ) denote the argument of U j 1,0 (resp. U k 0,1 ).

∃ ϕ ∈ R, Im(U j 1,0 e iω1ϕ ) = 0 Im(U k 0,1 e iω2ϕ ) = 0 ⇔ ∃ ϕ ∈ R, θ 1,0 + ω 1 ϕ ≡ 0[π] θ 0,1 + ω 2 ϕ ≡ 0[π] (32) ⇔ ∃ n 1,0 , n 0,1 ∈ Z, ϕ = n 1,0 π -θ 1,0 ω 1 = n 0,1 π -θ 0,1 ω 2 (33) ⇔ ∃ n 1,0 , n 0,1 ∈ Z, n 1,0 πω 2 -θ 1,0 ω 2 = n 0,1 πω 1 -θ 0,1 ω 1 (34) ⇔ (θ 1,0 ; θ 0,1 ) ∈ A, where (35)
A = x ; ω 2 ω 1 x + n 0,1 πω 1 -n 1,0 πω 2 ω 1 , x ∈ R, n 1,0 , n 0,1 ∈ Z ⊂ R 2
The set A is a dense, null subset of R 2 . This can be seen more clearly using the linear map (x, y) →

x, y -ω 2 ω 1

x . This map transforms A into the product of R and a dense, countable subset of R. Therefore, Eq. (32) has no solution, almost surely. However, for any (θ 1,0 ; θ 0,1 ) ∈ R 2 , there are elements in 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, ω1 ω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 ε 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.

Forced Van der Pol oscillator

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:

ẍ -µ 1 ẋ + µ 2 x ẋ + µ 3 x 2 ẋ + a 1 x = cos(λt) (36)
Values of dimensionless parameters are chosen as:

µ 1 = µ 2 = 0.1, µ 3 = a 1 = 1.
A first quadratic formulation, emphasizing the operators of Eq. ( 13), is:

ẋ = 0 + y + 0 (37a) ẏ = cos(λt) + µ 1 y -a 1 x -µ 2 xy -µ 3 yz (37b) 0 m u = c0 0 l0u + z q(u,u) -x 2 (37c)
Operators c 1 and l 1 are null. Similarly to example 4 in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], the forcing term cos(λ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(ωt) (line N eq + 2 in the system of Eq. ( 21)). The additional equation is then simply : ω = λ. The drawback is that if the system, forced at a pulsation λ, is able to adopt a pulsation ω = λ, the additional condition must be changed manually to reflect that: an a priori knowledge of the solution becomes necessary.

Another way (cited in [START_REF] Krauskopf | Numerical Continuation Methods for Dynamical Systems: Path following and boundary value problems[END_REF]) of including a forcing term is to consider the following nonlinear differential system:

ċ = c + λs -c(c 2 + s 2 ) (38a) ṡ = s -λc -s(c 2 + s 2 ) (38b)
with the phase equation s(t = 0) = 0, admitting c = cos(λt), s = -sin(λt) as a unique stable solution in the steady state regime. Hence, a second possible quadratic formulation is: 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.

ẋ = y + 0 + 0 (39a) ẏ = c + µ 1 y -a 1 x + 0 -µ 2 xy -µ 3 yz (39b) 0 = z + 0 -x 2 (39c) ċ = c + λs -cr (39d) ṡ = s -λc -sr (39e) 0 m u = l0u -r λl1u + 0 q(u,u) + c 2 + s 2 (39f)
The periodic solution at pulsation ω = λ is continued, and the bifurcation diagram is presented in fig. 1. The stability is computed thanks to Hill's method [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. The solution is stable for low forcing pulsations, then a Neimark-Sacker bifurcation occurs at λ 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:

λ = ω 1 , Im(x 0,1 ) = 0 (40)
The first condition is the same forcing condition as before, while the second condition means : the coefficient of sin((0ω 1 + 1ω 2 )t) for the variable x is null. In fig. 2, the coefficient of cos(0ω 1 + 1ω 2 )t of the variable x (i.e. Re(x 0,1 )) is plotted versus the forcing pulsation ω 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:

∀ k ∈ 0, H , x 2 k = k k1=-k k k2=-k |x k1,k2 | 2 , x 2 = H k1=-H H k2=-H |x k1,k2 | 2 (41)
and the threshold H x ∈ 1, H so that:

x 2 Hx 99.5 100 x 2 , x 2 Hx-1 < 99.5 100 x 2 (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 36)), continuation of the quasi-periodic solution with H = 6, coefficient of cos(ω 2 t) for the variable x versus the continuation parameter λ, that is, the forcing pulsation ω 1 . Solid blue line : solution branch; blue dots : ends of continuation steps. Letters NS indicate a Neimark-Sacker bifurcation.

with H = 6. It gives H x = 2 for all three variables (x, y, z), suggesting that an approximation of the quasiperiodic 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 [START_REF] Shampine | The Matlab ODE suite[END_REF]) on a long time interval. By long, the authors mean with respect to the meaningful time scales here, the two periods 2π ω 1 and 2π ω 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 distinguised 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.

Frequency-locking regions

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

∃ k 1 , k 2 , h 1 , h 2 ∈ -H, H , (k 1 , k 2 ) = (h 1 , h 2 ), k 1 ω 1 + k 2 ω 2 = h 1 ω 1 + h 2 ω 2 (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 ω 1 ω 2 = p q , Eq. ( 43) admitting solutions is equivalent to:

∃ k 1 , k 2 , h 1 , h 2 ∈ -H, H , (k 1 , k 2 ) = (h 1 , h 2 ), (k 1 -h 1 )ω 1 + (k 2 -h 2 )ω 2 = 0 (44) ⇔ ∃ k 1 , k 2 , h 1 , h 2 ∈ -H, H , (k 1 , k 2 ) = (h 1 , h 2 ), (k 1 -h 1 )pω 1 + (k 2 -h 2 )qω 1 = 0 (45) ⇔ ∃ k 1 = h 1 , k 2 = h 2 ∈ -H, H , h 2 -k 2 k 1 -h 1 = p q (46) ⇔ ∃ j 1 , j 2 ∈ -2H, 2H \ {0}, j 1 j 2 = p q (47)
Thus, for a given H, there is a finite number of potentially problematic values

j 1 j 2 , j 1 , j 2 ∈ -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 j 1 j 2 with small integers j 1 and j 2 , namely, in -H, H instead of -2H, 2H . Only two other values, 7 2 and 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 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 λ = 2 and 3.

The most important cases of frequency-locking observed on this system are subharmonic resonances, i.e. Then any combination k 1 ω 1 + k 2 ω 2 can be read as (qk 1 + k 2 )ω 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 ω = λ (fig. 1) can now be found easily and investigated.

∃ q ∈ 1, 2H , ω 1 ω 2 = q i.e. ω 2 = ω 1 q (48) a)
Figure 6 a, shows that in the vicinity of the subharmonic resonance λ = 2, the quasi-periodic HBM gives qualitatively good results, with an approximate shape of the solution in the phase diagram (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 ω = λ 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 λ comprised between 1.982 and 2.016. Therefore, an ideal continuation of the quasi-periodic branch should exhibit a plateau ω1 ω2 = 2, for this interval of λ. 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 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 λ = 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 ω 2 = ω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 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 1 See for instance [START_REF] Nayfeh | Applied nonlinear dynamics : analytical, computational and experimental methods[END_REF], section 4.2.

2 See for instance [START_REF] Nayfeh | Applied nonlinear dynamics : analytical, computational and experimental methods[END_REF], section 4. On this example, when the number H is large enough to engender issues in the vicinity of subharmonic resonances (for instance, H 2 around the subharmonic resonance ω 2 = ω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 ω 2 = ω1 q ⇔ ∀k 1 , k 2 , (k 1 + 1)ω 1 + (k 2 -q)ω 2 = k 1 ω 1 + k 2 ω 2 .
For instance, in the vicinity of the resonance 2, the pulsations 0, ω 1 -2ω 2 , 2ω 1 -4ω 2 and 3ω 1 -6ω 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. a) 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 .

Coupled Van der Pol oscillators

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

ẍ1 + a 1 ẋ1 + Ω 2 1 x 1 = a 2 λ( ẋ1 + ẋ2 ) -a 3 λ( ẋ1 + ẋ2 )(x 1 + x 2 ) -a 4 λ( ẋ1 + ẋ2 )(x 1 + x 2 ) 2 (49a) ẍ2 + b 1 ẋ2 + Ω 2 2 x 2 = b 2 λ( ẋ1 + ẋ2 ) -b 3 λ( ẋ1 + ẋ2 )(x 1 + x 2 ) -b 4 λ( ẋ1 + ẋ2 )(x 1 + x 2 ) 2 (49b)
Values of dimensionless parameters are chosen as: Ω 1 = 1, a 1 = 0.01, a 2 = 0.5, a 3 = a 4 = 2 ; Ω 2 = 2. Comparisons with time integration, between NS2 and λ 0.07, confirm that this portion of the branch is computed accurately with H = 3. For λ between 0.07 and 0.2, comparison with time integration is not meaningful on long intervals, as the quasi-periodic solution seems to be unstable3 : starting from the quasiperiodic solution, time integration eventually converges to one or the other periodic solutions found above, 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 λ (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 , ẋ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 ∼ 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).

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 sytem, 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.

M 2 =

 2 Diag(0m . . . Hm H+1 blocks -Hm . . . Hm 2H+1 blocks . . . -Hm . . . Hm 2H+1 blocks

Figure 1 :

 1 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.

Figure 3 :

 3 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 λ 4.29, quasi-periodic HBM (H = 2) in blue, time integration (ODE suite from Matlab) in red. a, b) Phase diagram (x, ẋ = y). c, d) Time comparison of the quasi-periodic HBM (H = 2) and time integration. From e to g) : similarly for λ 2.75. e, f) Phase diagram. g, h) Time comparison.

2 Figure 4 :j 2 , j 1 , 2 Figure 5 :

 242125 Figure 4: Forced Van der Pol, quasi-periodic solution branch computed with H = 6. a) Pulsation ratio ω 1 /ω 2 versus the forcing pulsation ω 1 . Solid blue line : solution branch; red dots : ends of continuation steps; dashed black lines : rational values j 1 j 2 , j 1 , j 2 ∈ -H, H ; mixed black lines : additional values, 7 2 and 7 3 . b) Zoom window : around the subharmonic resonances ω 2 = ω 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 ω 2 = ω 1 3 .

Figure 6 :

 6 Figure 6: Forced Van der Pol. a) Phase diagram (x, ẋ = y) in the vicinity of the subharmonic resonance ω 2 = ω 1 2 . Comparison of the quasi-periodic HBM (blue), periodic HBM (green) and time integration (red). b) Bifurcation diagram of the periodic solution 2ω = λ (solid line : stable solution, dashed : unstable, B : bifurcation).

Figure 7 :

 7 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ω = λ. b) 4ω = λ.

Figure 8 :Figure 9 :

 89 Figure 8: Forced Van der Pol, norm of x versus the forcing pulsation for the quasi-periodic HBM, around the subharmonic resonance ω 2 = ω 1 3 , for H = 2 (solid blue line) and H = 6 (dashed red line). a) 2-norm. b) 1-norm. a) b)

5 , b 1 Figure 11 :

 5111 Figure 11: Coupled Van der Pol oscillators, quasi-periodic solution branch. Some significant Fourier coefficients exhibit a squareroot growth starting from Neimark-Sacker bifurcations. a) Coefficient of cos(ω 2 t) for x 1 . b) Coefficient of cos(ω 1 t) for x 2 .a)

Figure 12 :

 12 Figure 12: Coupled Van der Pol oscillators, quasi-periodic solution branch computed with H = 4. a) Pulsation ratio ω 1 /ω 2 versus the continuation parameter λ. 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 ).

Figure 13 :

 13 Figure 13: Coupled Van der Pol oscillators, λ = 0.36, an example of quasi-periodic solution, comparison of different computations. a) Phase diagram (x 1 , ẋ1 = y 1 ), for the quasi-periodic HBM with H = 4. b) Diagram (x 1 , ẋ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.

Preprint submitted to ElsevierDecember 9, 2016

The authors prefer not to assert that as a fact, since no stability analysis has been performed.

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).

Neimark-Sacker bifurcation. The first periodic branch, with a pulsation close to Ω 1 and a great amplitude for the first oscillator x 1 , is plotted in red. The second periodic branch, with a pulsation close to Ω 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 λ. b) Amplitude of x 2 versus λ.

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 Ω 1 and a great amplitude for the first oscillator x 1 , is stable for small values of λ. It has a NS bifurcation (called NS1) at λ 0.35 where it becomes unstable. The second periodic branch, with a pulsation close to Ω 2 and a great amplitude for the second oscillator x 2 , is unstable at first, between a Hopf bifurcation (at λ = 0.025) and a NS bifurcation (at λ 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:

Im(x 1,(0,1) ) = 0, Im(x 2,(1,0) ) = 0 (51)

These conditions mean that the coefficient of sin(ω 2 t) for x 1 and the coefficient of sin(ω 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(ω 2 t) for x 1 , and by cos(ω 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(ω 1 t) for the variable x 2 (resp. the coefficient of cos(ω 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 Ω 1 (resp. Ω 2 ).

Near NS2, the pulsation ratio ω1 ω2 (fig. 12) is slightly inferior to the ratio Ω1 Ω2 = 1 2.5 and starts increasing when λ 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 λ = 0.2, where the number H can be increased. Thus, there is no problematic frequency-locking phenomenon when the pulsation ratio ω1 ω2 goes through the rational value 2 5 .