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

HAL
open science

# Comparison of two algorithms for Harmonic Balance and path continuation

Lukas Woiwode, Nidish Narayanaa Balaji, Jonas Kappauf, Fabia Tubita, Louis Guillot, Christophe Vergez, Bruno Cochelin, Aurélien Grolet, Malte Krack

► To cite this version:

Lukas Woiwode, Nidish Narayanaa Balaji, Jonas Kappauf, Fabia Tubita, Louis Guillot, et al. Comparison of two algorithms for Harmonic Balance and path continuation. Mechanical Systems and Signal Processing, 2020, 136, pp.106503. 10.1016/j.ymssp.2019.106503. hal-02424746

HAL Id: hal-02424746

https://hal.science/hal-02424746

Submitted on 28 Dec 2019

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.
Comparison of two algorithms for Harmonic Balance and path continuation

Lukas Woiwode<sup>d</sup>, Nidish Narayanaa Balaji<sup>e</sup>, Jonas Kappauf<sup>c</sup>, Fabia Tubita<sup>f</sup>, Louis Guillot<sup>b</sup>, Christophe

Vergez<sup>b</sup>, Bruno Cochelin<sup>b</sup>, Aurélien Grolet<sup>a</sup>, Malte Krack<sup>d</sup>

<sup>a</sup>Ecole Nationale Supérieure d'Arts et Métiers Lille - LISPEN, Lille, FRANCE

<sup>b</sup>Aix Marseille University, CNRS, Centrale Marseille, LMA, Marseille, FRANCE

<sup>c</sup>University of Kassel, Kassel, GERMANY

<sup>d</sup>University of Stuttgart, Stuttgart, GERMANY

<sup>e</sup>Rice University, Houston TX, USA

<sup>f</sup>Ecole Centrale de Lyon, Lyon, FRANCE

**Abstract**

In this work we apply and compare two algorithms for setting up Harmonic Balance equations and numerical continuation of the solution path for harmonically driven mechanical systems. The first algorithm relies on a predictor-corrector scheme and an Alternating Frequency-Time approach (AFT-PreCo). The second algorithm relies on a high-order Taylor series expansion of the solution path (asymptotic numerical method) and classical Harmonic Balance formulated entirely in the frequency domain (cHB-ANM). We conclude that the cHB-ANM is suited for smooth nonlinearities, for instance geometrically nonlinear finite element models. Here, cHB-ANM avoids aliasing errors and convinces with a numerically robust adjustment of the continuation step length and a continuous representation of the solution path. For non-smooth nonlinearities such as stick-slip friction or unilateral constraints, AFT-PreCo is better suited, and convinces with high numerical robustness and efficiency.

**Keywords:** Harmonic Balance Method, non-linear vibration, periodic solution, continuation method, friction damping, geometric nonlinearity
## Symbols and abbreviations

### Latin letters

$A_{\text{RMS}}$ root-mean-square value of the amplitude

$C$ vector of constant functions

$D$ coefficient matrix of velocity-proportional forces

$E_{HN}$ forward discrete Fourier transform matrix (N samples, H harmonics)

$E_{NH}^{-1}$ inverse discrete Fourier transform matrix (N samples, H harmonics)

$f$ generic function

$f_{\text{ex}}$ vector of excitation forces

$f_{\text{nl}}$ vector of nonlinear forces

$f_{\text{u}}$ unilateral spring force

$f_{\rho}$ friction force

$g$ coordinate of Coulomb slider

$H$ harmonic truncation order

$i$ imaginary unit $\sqrt{-1}$

$J_{X_0}$ Jacobian at expansion point $X_0$

$K$ coefficient matrix of coordinate-proportional forces

$L$ vector of linear functions

$M$ coefficient matrix of acceleration-proportional forces

$N$ number of time samples per period

$N_{\text{DOF}}$ number of degrees of freedom

$N_{\text{emp}}$ empirical rule for N

$N_{\text{eq}}$ number of equations

$N_{\text{fact}}$ number of Jacobian factorizations

$N_{\text{newt}}$ number of Newton iterations

$N_{\text{pt}}$ number of continuation steps

$P$ Taylor series order

$q$ vector of generalized coordinates

$Q$ vector of bilinear functions

$r$ residual of differential equation system

$R$ residual of algebraic equation system

$R_{\text{aux}}$ residual of auxiliary equations

$R_{\text{tot}}$ Concatenation of R and $R_{\text{aux}}$

$S$ dynamic stiffness matrix

$t_{\text{comp}}$ computation time

$T$ period of oscillation

$x$ vector of unknowns

$x_{\text{aux}}$ vector of auxiliary unknowns

$X$ $x$ extended by free parameter (excitation frequency)

$X_0$ Taylor series expansion point

$X_p$ Taylor series coefficients

### Greek letters

$\alpha$ step length

$\alpha_{\text{nom}}$ nominal step length (AFT-PreCo)

$\delta$ Dirac delta distribution

$\varepsilon_{\alpha}$ ANM threshold

$\varepsilon_{\omega}$ frequency error at resonance peak

$\varepsilon_{A}$ amplitude error at resonance peak

$\varepsilon_{\text{reg}}$ regularization parameter

$\varepsilon_{\text{tol}}$ error tolerance (algebraic equations)

$\eta$ normalized angular frequency

$\kappa$ linear stiffness of Coulomb slider

$\rho$ friction limit force

$\tau$ normalized time

$w_j$ weight functions

$\Omega$ fundamental angular osc. frequency

### Sub-, superscripts, operators

$|\square|$ magnitude

$\|\square\|$ norm (Euclidian if not specified)

$\square$ averaged value

$\tilde{\square}$ vector of time samples

$\dot{\square}$ derivative with respect to time

$\square'$ derivative with respect to time

$\hat{\square}_H$ normalized time

$\hat{\square}_H$ vector of Fourier coefficients (H harmonics)

$\hat{\square}(0)$ 0-th Fourier coefficient

$\hat{\square}_c(k)$ $k \ge 1$, k-th Fourier coefficient of cosine function

$\hat{\square}_s(k)$ $k \ge 1$, k-th Fourier coefficient of sine function

### Abbreviations

$\text{AFT}$ Alternating Frequency-Time (scheme)

$\text{ANM}$ asymptotic numerical method

$\text{cHB}$ classical Harmonic Balance (method)

$\text{DFT}$ discrete Fourier transform

$\text{HB}$ Harmonic Balance (method)

$\text{iFFT/FFT}$ inverse-/Fast Fourier Transform

$\text{PreCo}$ predictor-corrector (method)

$\text{SDOF}$ single degree of freedom (system)

$\text{MDOF}$ multiple degree of freedom (system)
# 1. Introduction

Harmonic Balance (HB) [1–3] is a method for approximating periodic solutions of nonlinear ordinary differential and differential-algebraic equations. HB can be interpreted as a Fourier-Galerkin method, where the (time-)dependent variables are represented by a truncated Fourier series, and the residual is required to be orthogonal to the Fourier basis functions [4]. This leads to an algebraic equation system for the sought Fourier coefficients of the approximation. HB avoids the costly computation of possibly long transients, and reasonable accuracy is often achieved already for a low truncation order. This permits to reduce the computational effort often by several orders of magnitude compared to numerical integration of the initial value problem [3]. Consequently, HB greatly gained importance in the field of (nonlinear) mechanical systems. For example, vibration problems of jointed structures with frictional contact interactions were thoroughly analyzed with HB, e. g. [5–9]. Besides, there are various review articles addressing the application of HB to a particular realm, e. g. on rotor dynamics [10] and brake squeal [11]. It is, however, not limited to vibration problems in mechanical systems, HB is also applied in other disciplines such as fluid dynamics [12] and electrical circuits [13].

A crucial task within HB is to determine the Fourier coefficients of the nonlinear terms. In some cases, in particular in the case of polynomial terms, an algebraic expansion is possible using trigonometric identities, or the continuous integral for the Fourier coefficients can be expressed in closed form. This is called *classical HB (cHB)*. A popular and more versatile alternative is to approximate the continuous integral using the discrete Fourier transform, which leads to the *Alternating Frequency-Time (AFT)* HB [3, 14]. An important parameter of this method is the number of time samples per oscillation period. In the special case of polynomial nonlinearities, the method is exact beyond a certain number of samples. Otherwise aliasing errors occur.

It is often of interest to determine how the solution evolves with a free parameter. For this task, numerical path continuation is commonly applied. The task of numerical path continuation methods is to recover the solution branch accurately and overcome potential turning points with respect to the free parameter. Particularly popular is the *predictor-corrector (PreCo)* method. From the current solution point, a prediction is made, for instance by going a certain step length in the direction of the tangent to the solution branch. A Newton-type method is then usually applied to iteratively correct the prediction and return to the solution branch. An interesting alternative is the *asymptotic numerical method (ANM)* [15–18]. Here, the solution path is expanded into a high-order Taylor series with respect to an arc length parameter. By remaining within the range of utility of the truncated power series, correction steps can be avoided. A challenge within the ANM is to efficiently calculate the high-order Taylor series coefficients. For this, the algebraic equation system is commonly recast to nonlinear terms of only quadratic order, by introducing auxiliary variables and additional equations.

The ANM can only be applied when the derivatives for the Taylor series expansion exist, i. e., it is only applicable to sufficiently smooth nonlinearities. However, many mechanical systems are modeled with non-smooth force-deformation relations, such as systems with stick-slip friction, unilateral constraints or materials undergoing phase transformations. A regularization is then usually required to approximate the non-smooth relations by smooth ones in order to apply the ANM.

In this work, we focus on two algorithms, AFT-PreCo and cHB-ANM. The goal is to identify and understand the individual strengths and weaknesses of these algorithms, both with regard to computational efficiency and robustness, as well as required level of experience to properly select critical method parameters. We analyzed various benchmark problems of harmonically driven mechanical systems, with either smooth or non-smooth nonlinearities, and different numbers of degrees of freedom. In Section 2, we recap the computational methods and the theoretical background required to understand their similarities and differences. In Section 3, we present selected results of representative benchmark problems. In Section 4, we develop guidelines towards the recommended application range of either method.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
## 2. Recap of Harmonic Balance and continuation methods

### 2.1. Harmonic Balance

In this work, we consider mechanical systems described by the equation of motion

$$r(t) = M\ddot{q}(t) + D\dot{q}(t) + Kq(t) + f_{nl}(q(t)) - f_{ex}(t) = 0, \quad (1)$$

where $q \in \mathbb{R}^{N_{DOF} \times 1}$ is the vector of generalized coordinates, overdots denote differentiation with respect to time $t$, $M, D, K \in \mathbb{R}^{N_{DOF} \times N_{DOF}}$ are coefficient matrices of acceleration-, velocity-, displacement-proportional forces, and $f_{nl}, f_{ex} \in \mathbb{R}^{N_{DOF} \times 1}$ are nonlinear and external forces, respectively. For ease of notation, we assume that $f_{nl}$ depends only on $q$; it is straight-forward to extend the following developments to the case where it depends also on $\dot{q}$ (and higher-order derivatives with respect to time [3]). We consider $T$-periodic external forces, $f_{ex}(t+T) = f_{ex}(t)$, and seek $T$-periodic responses $q(t)$. The requirement that $q(T) = q(0)$ where $q(t)$ has to satisfy Eq. (1) corresponds to a periodic boundary value problem on the domain $t \in [0, T]$.

#### 2.1.1. Method of mean weighted residuals, Fourier methods

The methods we apply in this work to solve the periodic boundary value problem can be interpreted as methods of mean weighted residuals, with a truncated Fourier series as ansatz,

$$q(t) = \hat{q}(0) + \sum_{k=1}^{H} \hat{q}_c(k) \cos(k\Omega t) + \hat{q}_s(k) \sin(k\Omega t), \quad (2)$$

with the angular frequency $\Omega = 2\pi/T$ and the harmonic truncation order $H$. Substitution of the ansatz (2) into Eq. (1) generally produces a non-zero residual $r(t) \neq 0$. The unknown Fourier coefficients are determined as the solution of the algebraic equation system

$$\frac{1}{T} \int_{0}^{T} r(t) w_{j}(t) dt = 0 \quad j = 1, \dots, 2H + 1, \quad (3)$$

which requires that the residual vanishes in a weighted average sense, with the weight functions $w_j(t)$; i.e., the residual is orthogonal to the truncated Fourier basis. When the $2H + 1$ Fourier basis functions $\{1, \cos(\Omega t), \dots, \sin(H\Omega t)\}$ are used as weights, Eq. (3) requires that the Fourier coefficients of the residual associated with the sine-/cosine vanish up to order $H$; i.e., $\hat{r}(0) = \hat{r}_c(1) = \dots = \hat{r}_s(H) = 0$. This algorithm was already introduced in 1965 by Urabe [4] and is now commonly known as Harmonic Balance (HB). HB can thus be interpreted as a Galerkin method where Fourier basis functions are used both as ansatz and weight functions. When Dirac delta distributions $\delta(t - t_j)$ with the shifts $t_j = jT/N$, $j = 0, \dots, N-1$ are used as weights, Eq. (3) requires that the residual vanishes at these collocation points; i.e., $r(t_0) = \dots = r(t_{N-1}) = 0$. Together with the Fourier ansatz in Eq. (2), this is known as Trigonometric Collocation. When we stack all Fourier coefficients of $q, r$ in vectors $\hat{q}_H, \hat{r}_H$, we can summarize the HB equations as

$$\hat{r}_H(\hat{q}_H) = S(\Omega)\hat{q}_H + \hat{f}_{nl,H} - \hat{f}_{ex,H} = 0, \quad (4)$$

where $S(\Omega)\hat{q}_H$ accounts for the linear terms in Eq. (1) with the dynamic stiffness matrix $S$, and $\hat{f}_{nl,H}$, $\hat{f}_{ex,H}$ are the vectors of Fourier coefficients of the nonlinear and external forces, respectively.

#### 2.1.2. Classical and Alternating-Frequency-Time Harmonic Balance

A challenge within HB is to determine the Fourier coefficients of the nonlinear terms. In some cases, for example when $f_{nl}$ is a polynomial in $q$ (and more generally also in $\dot{q}$), closed-form expressions can be established, either by algebraic expansion of $f_{nl}$ into a Fourier series (using trigonometric identities), or by analytical integration of the Fourier coefficients (e. g. using discrete convolutions). In general, the discrete
Fourier transform (DFT) can be useful to approximate the continuous integrals for the Fourier coefficients, e. g.,

$$ \hat{\mathbf{f}}_{\text{nl,c}}(k) = \frac{1}{T} \int_{0}^{T} \mathbf{f}_{\text{nl}}(\mathbf{q}(t)) \cos(k\Omega t) dt \approx \frac{1}{N} \sum_{j=0}^{N-1} \mathbf{f}_{\text{nl}}(\mathbf{q}(t_j)) \cos(k\Omega t_j), \quad (5) $$

where $\mathbf{q}(t_j)$ is evaluated in accordance with Eq. (2) (sampled) at the equidistant time instants $t_j = jT/N$. When we stack all samples of $\mathbf{q}(t_j)$ in a vector $\tilde{\mathbf{q}}_N$, we can summarize this as

$$ \hat{\mathbf{f}}_{\text{nl},H} \approx \mathbf{E}_{HN} \hat{\mathbf{f}}_{\text{nl},N} (\mathbf{E}_{HN}^{-1} \hat{\mathbf{q}}_H), \quad (6) $$

where $\mathbf{E}_{HN}$ is the DFT matrix (for approximation of the Fourier coefficients from the samples), and $\mathbf{E}_{HN}^{-1}$ is the inverse DFT matrix (for interpolation of the trigonometric polynomial). This procedure was also already proposed by Urabe [19] and is now commonly known as Alternating-Frequency-Time scheme. For a given frequency-domain representation of the approximation, $\hat{\mathbf{q}}_H$, the time samples $\mathbf{q}(t_j)$ are determined first. The nonlinear forces $\mathbf{f}_{nl}$ are then evaluated in the time domain. Finally, the Fourier coefficients $\hat{\mathbf{f}}_{nl,H}$ are approximated using the DFT.

In some cases, the sampling procedure in Eq. (6) (and Eq. (5)) is exact; i. e., the DFT yields the exact value of the continuous integral for the Fourier coefficients. This is the case if $\mathbf{f}_{nl}$, evaluated at the Fourier ansatz (2), is a Fourier series of finite truncation order $H_f$ and $N \ge 2H_f + 1$. This, in turn, happens if $\mathbf{f}_{nl}$ is a polynomial in $\mathbf{q}$ (and more generally also in $\hat{\mathbf{q}}$). This is shown in detail in Appendix A. The polynomial can be broken down to sums of products. A product of two Fourier series, truncated to orders $H_1$ and $H_2$, respectively, is again a Fourier series, with truncation order $H_1 + H_2$. Moreover, any linear operation and differentiation with respect to time yields a Fourier series with the same truncation order. The Fast Fourier Transform (FFT) is a computationally efficient DFT algorithm. Using the AFT scheme can in fact be the fastest and still exact opportunity to calculate $\hat{\mathbf{f}}_{nl,H}$ even in the case of polynomials [20]. If the number of samples, $N$, is not sufficiently high ($N < 2H_f + 1$), aliasing errors occur. Aliasing is inevitable if the considered periodic function has an infinite sequence of non-zero Fourier coefficients (no finite truncation order), see e. g. [3].

To summarize, the equation systems of classical HB (cHB) and Alternating-Frequency-Time (AFT) HB are:

$$ \text{cHB:} \qquad \hat{\mathbf{r}}_H(\hat{\mathbf{q}}_H) \qquad = \mathbf{0}, \quad (7) $$

$$ \text{AFT HB:} \qquad \mathbf{E}_{HN} \tilde{\mathbf{r}}_N (\mathbf{E}_{HN}^{-1} \hat{\mathbf{q}}_H) \qquad = \mathbf{0}. \quad (8) $$

We here utilized that the sampling procedure is exact for the linear terms in Eq. (1). For large $N$, the aliasing error typically becomes negligible compared to other errors (e. g. truncation error, tolerances of numerical root finding methods), and the solution of Eq. (8) should converge to the solution of Eq. (7).

### 2.1.3. Another interpretation of AFT HB

We can split Eq. (1) into the differential-algebraic equation system,

$$ \mathbf{0} = \mathbf{M}\ddot{\mathbf{q}} + \mathbf{D}\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} + \mathbf{v} - \mathbf{f}_{\text{ex}}, \quad (9) $$

$$ \mathbf{0} = \mathbf{v} - \mathbf{f}_{\text{nl}}(\mathbf{q}), \quad (10) $$

by introducing the auxiliary variables $\mathbf{v}$. Suppose that we apply cHB with truncation order $H_q$ only to the differential part (9), and Trigonometric Collocation with $N = 2H_v+1$ collocation points to the algebraic part (10). cHB applied to Eq. (9) simply leads to $\mathbf{S}(\Omega)\hat{\mathbf{q}}_{H_q} + \hat{\mathbf{v}}_{H_q} - \hat{\mathbf{f}}_{\text{ex},H_q}$. Trigonometric Collocation applied to Eq. (10) leads to $\mathbf{E}_{H_v N}^{-1} \hat{\mathbf{v}}_{H_v} = \tilde{\mathbf{f}}_{nl,N} (\mathbf{E}_{H_q N}^{-1} \hat{\mathbf{q}}_{H_q})$. From this, we can follow $\hat{\mathbf{v}}_{H_v} = \mathbf{E}_{H_v N} \tilde{\mathbf{f}}_{nl,N} (\mathbf{E}_{H_q N}^{-1} \hat{\mathbf{q}}_{H_q})$. For $H_v \ge H_q$ one can truncate this to $\hat{\mathbf{v}}_{H_q} = \mathbf{E}_{H_q N} \tilde{\mathbf{f}}_{nl,N} (\mathbf{E}_{H_q N}^{-1} \hat{\mathbf{q}}_{H_q})$, which is fully equivalent to Eq. (8). Thus, AFT HB can be interpreted as applying classical Harmonic Balance to the linear part of the equation of motion and Trigonometric Collocation to the nonlinear terms, but generally distinct truncation orders $H_q \ne H_v$ for the two sets of variables.
## 2.2. Numerical path continuation

It is often of interest to determine how the periodic solution evolves under variation of a free parameter. In this work, we focus on the excitation frequency $\Omega$ as free parameter, i. e., we analyze the frequency response. Thus, we are interested in the solutions $\mathbf{X} = [\mathbf{x}^T, \Omega]^T$ of the algebraic equation system

$$ \mathbf{R}(\mathbf{X}) = \mathbf{0}, \qquad (11) $$

where $\mathbf{R}$, $\mathbf{x} \in \mathbb{R}^{N_{eq} \times 1}$ are the residual vector and the vector of unknowns, $N_{eq}$ is the number of equations (and unknowns). The solutions form continuous branches under variation of the free parameter. We can directly consider the equation systems (7)-(8) and set $\mathbf{x} = \hat{\mathbf{q}}_H$ and e. g. $\mathbf{R} = \hat{\mathbf{r}}_H$ for cHB. The task of numerical path continuation methods is to recover the solution branch accurately and overcome potential turning points with respect to the free parameter (cf. Fig. 1a).

### 2.2.1. Predictor-Corrector (PreCo) method

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of continuation methods: (a) solution branch, (b) predictor-corrector (PreCo) method, (c) asymptotic numerical method (ANM)</figcaption></figure>

Suppose we already know a solution point $\mathbf{X}_0$. To obtain a next point on the solution branch, we can make a prediction,

$$ \mathbf{X}(\alpha) \approx \mathbf{X}_0 + \alpha \mathbf{X}_1, \qquad (12) $$

in the direction of the tangent to the solution branch, $\mathbf{X}_1$, with the step length $\alpha$. The tangent with unit length is defined by

$$ \underbrace{\frac{\partial \mathbf{R}}{\partial \mathbf{X}}}_{\mathbf{X}_0} \mid_{\mathbf{X}_1 = \mathbf{0}, \quad \mathbf{X}_1^T \mathbf{X}_1 = 1.} \qquad (13) $$

$\mathbf{X}_1$ is well-defined if the Jacobian $\mathbf{J}_{\mathbf{X}_0}$ has rank $N_{eq}$; then $\mathbf{X}_0$ is called a regular solution point. We usually have $\mathbf{R}(\mathbf{X}_0 + \alpha \mathbf{X}_1) \neq \mathbf{0}$, and a correction is needed to decrease the residual to a tolerable level and obtain an acceptable next solution point. Once the next solution point is found, the current continuation step is completed, and the process repeats from this point, until a termination criterion is met (e. g. the parameter
end value is reached). It is common to introduce an additional equation, which defines where along the solution branch the next point is located. A popular choice is [21–23]

$$ \mathbf{R}_{\text{tot}}(\mathbf{X}) = \left[ \frac{\mathbf{R}(\mathbf{X})}{(\mathbf{X} - \mathbf{X}_0)^{\text{T}}(\mathbf{X} - \mathbf{X}_0)} - \alpha^2 \right] = \mathbf{0}. \quad (14) $$

Here, the last equation requires that the next solution point lies on the hypersphere of radius $\alpha$ around the previous solution point $\mathbf{X}_0$ (cf. Fig. 1b). The number of equations is now equal to that of the unknowns, and a Newton-type procedure is commonly applied, which corresponds to the iterative solution of the equations

$$ \left. \frac{\partial \mathbf{R}_{\text{tot}}}{\partial \mathbf{X}} \right|_{\mathbf{X}^{(j)}} (\mathbf{X}^{(j+1)} - \mathbf{X}^{(j)}) = - \mathbf{R}_{\text{tot}}(\mathbf{X}^{(j)}), \quad (15) $$

starting from the initial guess $\mathbf{X}^{(0)} = \mathbf{X}_0 + \alpha \mathbf{X}_1$, until $\|\mathbf{R}_{\text{tot}}(\mathbf{X}^{(j)})\| \le \varepsilon_{\text{tol}}$. For well-behaved problems, the continuation step length $\alpha$ can be selected such that typically only 3 to 10 Newton iterations are needed. The required derivatives $\partial \mathbf{R}_{\text{tot}}/\partial \mathbf{X}$ are commonly calculated analytically. This is straight-forward with cHB where a closed-form expression is available for the residual as function of the unknowns. This is also easy with the AFT scheme thanks to the linearity of forward and inverse DFT [24]. Various other predictors and path parametrization strategies are possible, e. g. normal flow, the interested reader is referred to [3, 23], and also [25] where the performance of different parametrizations is compared.

### 2.2.2. Asymptotic numerical method (ANM)

The predictor in Eq. (12) of the PreCo method can be viewed as first-order part of a Taylor series expansion around $\mathbf{X}_0$. The ANM extends this to a higher-order predictor (order $P$),

$$ \mathbf{X}(\alpha) \approx \sum_{p=0}^{P} \alpha^p \mathbf{X}_p, \quad (16) $$

where $\mathbf{X}_p$ are the series coefficients. The domain of utility $[0, \alpha_{\text{max}}]$ of the expansion in Eq. (16) can be determined. To this end, the residual $\mathbf{R}$ is also expanded in a Taylor series. It is then assumed that the residual is dominated by the $(P+1)$-th order term, which yields

$$ \alpha_{\text{max}} = \left( \frac{\varepsilon_{\alpha}}{\| \mathbf{R}_{P+1} \|} \right)^{\frac{1}{P+1}}, \quad (17) $$

where $\varepsilon_{\alpha}$ is a user-defined tolerance (ANM threshold) and $\mathbf{R}_{P+1}$ is the $(P+1)$th-order coefficient of the Taylor series expansion of the residual $\mathbf{R}$ [17, 26]. The idea behind this is to choose the continuation step length $\alpha_{\text{max}}$ in such a way that the residual remains within the specified tolerance $\varepsilon_{\text{tol}}$, with the intent to completely avoid correction iterations. Ideally, the ANM threshold $\varepsilon_{\alpha}$ should be set to the tolerance $\varepsilon_{\text{tol}}$. Since the expansion point already has a non-zero residual in practice, and Eq. (17) is only an approximation of the true radius of convergence (truncating the residual to the $(P+1)$th-order term of its Taylor series), $\|\mathbf{R}(\mathbf{X}(\alpha_{\text{max}}))\| > \varepsilon_{\text{tol}}$ in general, such that correction steps are occasionally needed. To avoid this, the ANM threshold $\varepsilon_{\alpha}$ is commonly set to a smaller value, $\varepsilon_{\alpha} < \varepsilon_{\text{tol}}$. This is further discussed in the result sections of this paper. Once $\alpha_{\text{max}}$ is determined, $\mathbf{X}(\alpha_{\text{max}})$ is evaluated, which then serves as the next expansion point and so forth, cf. Fig. 1c. The ANM thus yields a section-wise power series expansion of the solution path, which together form an accurate continuous representation of the entire solution branch. The power series expansion of the solution path can be useful for bifurcation analysis [27]. Experience has shown that for most problems, an order $P$ of around 20 yields reasonable computational performance [17], which is a compromise between the effort of computing higher order terms and the reward of having an enlarged radius of convergence, and this value was selected throughout this study.

An essential ingredient of the ANM is a procedure to compute the series coefficients $\mathbf{X}_p$. For this purpose,
automatic differentiation can be used. However, this is known to be computationally costly when more than a few hundreds of equations are considered [28]. The computation of the series coefficients is particularly efficient when the algebraic equation system is in quadratic form,

$$R(X) = C + L(X) + Q(X, X), \qquad (18)$$

where C, L, Q are constant, linear and bilinear (quadratic) functions. For a given solution point $X_0$ and tangent $X_1$ determined as in Eq. (13), the higher-order series coefficients can then be determined successively from the linear equations (see e. g. [15])

$$J_{X_0} X_p + \sum_{i=1}^{p-1} Q(X_i, X_{p-i}) = 0 \quad p = 2, \dots, P. \qquad (19)$$

Note that the sought $X_p$ only occur in the first term, while the sum only involves $X_i$ with $i < p$. The equation systems (7)-(8) are generally not in quadratic form. However, if $f_{nl}$ is an analytic function in $q$ (and more generally also in $\dot{q}$), the equation system can usually be recast into quadratic form. To this end, extra unknowns $x_{aux}$ are introduced besides the main unknowns $\dot{q}_H$, along with additional equations $R_{aux} = 0$,

$$x = \begin{bmatrix} \dot{q}_H \\ x_{aux} \end{bmatrix}, \quad R = \begin{bmatrix} \dot{r}_H \\ R_{aux} \end{bmatrix} = 0. \qquad (20)$$

This recast is particularly simple if the nonlinear terms, $f_{nl}$, in the equation of motion are polynomials or rational functions [15]. One can then introduce auxiliary variables $v$ defined by algebraic equations $r_{aux}(q, v) = 0$ in such a way that the differential-algebraic equation system contains only quadratic nonlinear terms in $q$ and $v$. Both $q$ and $v$ are then represented by truncated Fourier series. Application of HB (either the cHB or the AFT variant) leads to a quadratic algebraic equation system with regard to the unknowns $\dot{q}_H, \dot{v}_H$. The described formalism has been extended to some transcendental functions [16, 17], but this is not relevant for the benchmark problems within the present study.

Note that the $X_p$ are the solutions of linear equation systems (19) that share the same coefficient matrix, $J_{X_0}$, for each $p$. Hence, this matrix only needs to be factorized once per expansion point $X_0$. Moreover, if one defines every new auxiliary variable only by the main and lower-index auxiliary variables, it can be shown that the Jacobian block associated to the extra unknowns $x_{aux}$ is triangular (and thus its inverse is easy to determine). A block elimination can then be applied to solve Eq. (19) such that the critical problem dimension reduces to that of the main variables [17]. For the efficient evaluation of C, L and Q in large-dimensional problems, a sparse tensorial formalism can be used. The associated tensor lists can be generated automatically using polarization formulae. The user then only has to define the equation of motion (1), and the algebraic equations $r_{aux}(q, v) = 0$ defining the auxiliary variables $v$ [17].

HB theoretically does not fail if $f_{nl}$ is only differentiable up to a certain order ($C^m$ function with $m < \infty$), i.e. representing a non-smooth function. However, the convergence of the Fourier series is generally slower than in the case of analytic functions [3]. Urabe [4] established sufficient conditions for the convergence of HB to an exact solution. Here, the nonlinear forces, along with their first-order derivatives with respect to $q$ (and more generally also $\dot{q}$), were only required to be continuously differentiable. As one can easily see from the left part of Eq. (5), finite differentiability of $f_{nl}$ with respect to $q$ implies (by the chain rule) finite differentiability of $R(X(\alpha))$ with respect to $\alpha$. Hence, the ANM can only be applied to a certain order. More importantly, cHB can usually not be applied, and the system cannot be brought into quadratic form. To apply cHB and the ANM properly, a regularization is needed, where the non-smooth functions are approximated by analytic ones. This introduces additional errors and comes with the task to properly choose the parameter(s) of the regularization.

## 2.3. Morphological box of methods

In principle, any continuation algorithm (PreCo method vs. ANM) can be combined with any Fourier method (cHB vs. AFT scheme). On the other hand, the ANM only makes sense if one has an efficient
procedure for computing the required Taylor series coefficients up to high order. Here, the aforementioned
recast of the equation of motion (1) into a differential-algebraic equation system with auxiliary variables $v$
is particularly useful. Both the main and the auxiliary variables are then approximated by $H$-order Fourier
series. To evaluate the residual of the quadratically recast system in Eq. (18), one has to determine the
Fourier coefficients of the nonlinear terms. This can, in principle, also be done using the discrete Fourier
transform. For at most quadratic nonlinear terms, the discrete Fourier transform yields the exact value of
the Fourier coefficients for $N \ge 3H + 1$, as shown in Appendix A. In this sense, cHB and the AFT variant
(applied to the quadratically recast problem) are theoretically fully equivalent. However, the numerical
implementation of the nonlinear terms differs (discrete convolutions within cHB vs. FFT/iFFT within the
AFT scheme), leading to different computational effort in general.

In many vibration problems, the generalized coordinates $q$ are dominated by only one or a few harmonics.
In contrast, more harmonics are commonly required to represent the nonlinear terms $f_{nl}$ in the equation
of motion Eq. (1), where $f_{nl}$ depends on the generalized coordinates $q$. It is therefore expected that more
harmonics are required to accurately approximate the nonlinear forces $f_{nl}$, and also the auxiliary variables
$v$. Considering the quadratically recast system, it is therefore interesting to use a higher truncation order
$H_v > H_q$ for the auxiliary variables than for the generalized coordinates. Analogously, considering instead
the initial system of equations of motion (1), it is interesting to apply the AFT scheme with a large number
of samples $N \gg 2H + 1$.

In this work, we focus on two algorithms:

• AFT-PreCo: HB is applied to the initial equation of motion (1) using the AFT scheme (harmonic truncation order $H$, number of time samples $N$); continuation is done with tangent predictor (Eq. (12)), pseudo arc-length constraint (Eq. (14)) and Newton-type correction iterations (Eq. (15))

• cHB-ANM: classical HB is applied to the quadratically recast differential-algebraic equation system (harmonic truncation order $H$ used for both generalized coordinates $q$ and auxiliary variables $v$); continuation is done with $P$-order Taylor series predictor (Eq. (16)), where the step length $\alpha$ is set to $\alpha_{max}$ in accordance with Eq. (17)

It must be emphasized that a *crucial difference* between these two algorithms is whether HB is applied to
the *ordinary differential equation of motion* or a *quadratic recast*, and, thus, how the *nonlinear terms are
represented*. The problem dimension of cHB-ANM is larger due to the quadratic recast. To reach the same
accuracy, it is expected that cHB-ANM needs a higher harmonic truncation order $H$ to properly represent
the auxiliary variables. These two aspects could lead to a larger computational effort for cHB-ANM as
compared with AFT-PreCo. On the other hand, the high-order expansion of the solution path within the
ANM might permit larger step lengths. Moreover, only a single Jacobian matrix has to be factorized for each
continuation step, as compared with typically 3 to 10 for the Newton iterations within the PreCo method.
These two aspects could (over-)compensate the aforementioned computational disadvantages. This will be
analyzed for several benchmark problems in the next section.

The algorithms were implemented in Matlab. The point of departure were the open source tools MAN-
LAB 4.0 (cHB-ANM) [29, 30] and NLvib 1.0 (AFT-PreCo) [3, 31]. Our intent is obviously to compare the
algorithms, not the tools. We therefore placed these implementations into a unified framework. The com-
putation of the initial solution point, pre- and post-processing tasks were excluded from the computational
effort comparisons. Where appropriate, the same linear algebra operations (e. g. matrix factorizations) were
used. Where appropriate, loops were vectorized. Nevertheless, a certain *implementation bias* cannot be
excluded, i. e., when a different person implements the methods, the computational efficiency might differ
to some degree. Thus, the computational effort comparisons must be regarded with care; only substantial
changes (e. g. an order of magnitude) should be considered relevant.

## 2.4. Summary of error sources

It is useful to recall at this point the causes of potential errors made by the described algorithms,
and by what parameters these are affected, as listed in Tab. 1. **The PreCo algorithm just generates a
sequence of solution points, without information about the course of the solution branch between two**
adjacent points. Hence, when the step length is large, important features (e.g. super-harmonic regime, primary resonance peak) along the solution branch might be resolved poorly or overlooked completely. Then, approximation quality can be increased by recovering a continuous representation by interpolation (here quadratic). However, the level of accuracy is then not just driven by HB approximation quality (given by $\epsilon_{tol}$), but also by an additional interpolation error (related to $\alpha$). In contrast, the ANM yields a continuous representation of the solution branch; the accuracy of this representation is limited by the truncation error (controlled by Taylor series order $P$), and the assumed range of utility (controlled by the ANM threshold $\epsilon_{\alpha}$).

<table><caption>Table 1: Error sources in AFT-PreCo and cHB-ANM algorithms</caption><thead><tr><th></th><th>error source</th><th>control parameter</th></tr></thead><tbody><tr><td></td><td>truncation</td><td>harmonic order <i>H</i></td></tr><tr><td></td><td>aliasing (for non-polynomial<br>nonlinearities within the AFT scheme)</td><td>number of samples per period <i>N</i></td></tr><tr><td></td><td>regularization (to approximate<br>non-smooth terms by analytic ones)</td><td>parameter <i>&epsilon;</i><sub>reg</sub></td></tr><tr><td></td><td>numerical tolerance</td><td>threshold <i>&epsilon;</i><sub>tol</sub>, ANM: threshold <i>&epsilon;</i><sub>&alpha;</sub></td></tr><tr><td></td><td>interpolation</td><td>PreCo: step length <i>&alpha;</i></td></tr></tbody></table>

### 3. Results for representative benchmark problems

In the course of the presented study, various benchmark problems were analyzed, with different types of nonlinearities and different underlying linear base structures having different numbers of degrees of freedom. In the following, we present only selected representative results of this comprehensive study. Three different types of nonlinearity are presented: polynomial stiffness, unilateral stiffness and dry friction.

#### 3.1. System with polynomial stiffness

First, we consider the equations of motion [32]

$$q_k'' + 0.1k^2 q_k' + k^4 q_k + \sum_{j=1}^{N_{\text{DOF}}} k^2 j^2 q_k q_j^2 = \hat{f}_c(1) \sin(k \frac{\pi}{2}) \cos(\eta \tau) \quad k = 1, \dots, N_{\text{DOF}}. \quad (21)$$

Herein, $'$ denotes derivative with respect to normalized time $\tau$. For $N_{\text{DOF}} = 1$, Eq. (21) simplifies to the equation of motion of the harmonically driven and damped Duffing oscillator. For $N_{\text{DOF}} > 1$, Eq. (21) can be interpreted as the equations of motion of a pinned-pinned beam: It is derived e. g. in [32], considers the von-Karman beam theory and neglects the longitudinal inertia. The resulting momentum balance condensed to the transverse direction is then projected onto the basis composed of the first $N_{\text{DOF}}$ lowest-frequency odd-numbered (linear) normal modes. This results in the set of nonlinear ordinary differential equations (21), where $q_k$ are the associated modal coordinates of the mass-normalized modes. In summary, the cubic polynomial stiffness terms in Eq. (21) model the geometric stiffening caused by the midplane stretching due to bending under axial constraints. The harmonic forcing with normalized angular frequency $\eta$ is applied as a concentrated load at the center of the beam with magnitude $\hat{f}_c(1)$. Only the odd-numbered modes receive external forcing, because the remaining modes have a node at the center of the symmetric beam. Yet, both even- and odd-numbered modes are considered in the following analyses. For the quadratic recast, one auxiliary variable $v_k$ is introduced per modal coordinate along with the additional equation $r_{\text{aux},k} = v_k - q_k^2 = 0$.

As mentioned before, the sampling procedure within the AFT scheme yields the Fourier coefficients of polynomial nonlinear terms without aliasing errors (exactly) for sufficiently many samples per period, $N$.
In the considered examples, $q(t)$ only responds with odd-numbered harmonics. As shown in Appendix A, $N \ge 4H+1$ is sufficient to completely avoid aliasing errors in the relevant $H$ lowest-order Fourier coefficients of the cubic nonlinear terms. For the results depicted in Fig. 2-6, the number of time samples was set to $N = 4H+1$ accordingly.

### 3.1.1. Effect of PreCo step length $\alpha$ and ANM threshold $\varepsilon_\alpha$

<figure><img src="image_3.png" /><figcaption>Figure 2: Duffing oscillator, effect of nominal continuation step length of the AFT-PreCo method on amplitude-frequency curves; (b), (c), (d) illustration of cyclic jumping between solution branches in super-harmonic regime for large step lengths (do not represent zooms into (a))</figcaption></figure>

In the following, we first focus on the case $N_{\text{DOF}} = 1$ (Duffing oscillator), and the excitation level is set to $\dot{f}_c(1) = 1.5$. The amplitude $A_{\text{RMS}}$ is defined as the root-mean-square value

$$A_{\text{RMS}} = \sqrt{\dot{q}^2} = \sqrt{\dot{q}^2(0) + \frac{1}{2} \sum_{k=1}^{H} \dot{q}_c^2(k) + \dot{q}_s^2(k)}. \qquad (22)$$

Since an exact solution is not available, the reference solution for this benchmark is generated by cHB-ANM with a large number of harmonics (here $H = 200$ chosen). This seems justified since for polynomials ($C^\infty$ function), the truncated Fourier series is known to converge rapidly towards the reference solution with increasing $H$. Besides that the continuous representation of the solution path allows an analytical determination of the resonance peak.

The step length $\alpha$ within the PreCo method is generally a critical parameter and has to be specified by the user, at least to a certain extent. For the considered PreCo method, only a nominal step length $\alpha_{\text{nom}}$ is specified by the user. The lower and upper bounds are set to $\alpha_{\text{nom}}/5$ and $5\alpha_{\text{nom}}$, respectively. This was empirically determined as reasonable tradeoff between adaptability in ranges with strong gradient changes of the solution path and probability that important features of the solution path are not overlooked. If the number of Newton iterations for consecutive continuation steps exceeded 9, the step length is halved, and it is doubled if this number fell below 6.
<figure><img src="image_4.png" /></figure>

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

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

<figure><img src="image_7.png" /><figcaption>Figure 3: Duffing oscillator, effect of nominal step length $\alpha_{nom}$ on AFT-PreCo: (a) number of continuation steps, (b) averaged number of Newton iterations, (c) computation time; effect of ANM threshold $\varepsilon_{\alpha}$: (d) number of continuation steps, (e) averaged number of Newton iterations, (f) computation time</figcaption></figure>

For very large nominal step lengths, the PreCo method might jump between solution branches. As a consequence, important features of the bifurcation diagram might be overlooked, or the algorithm cycles back and forth on the solution path as illustrated in Fig. 2b-d. In general, a very large step length may also lead to divergence or unreasonably slow convergence of the Newton method. For moderately large nominal step lengths, such problems do no longer occur. However, the solution path is only poorly represented in ranges with rapid gradient changes, as indicated in the zooms within Fig. 2a.

In Fig. 3a-c, the effect of the nominal step length $\alpha_{nom}$ (AFT-PreCo) on the number of solution points (equal to the number of continuation steps), average number of Newton iterations and computation time is shown ($H = 20$). For small step lengths, the number of continuation steps $N_{pt}$ and the time $t_{comp}$ for computing the whole branch are approximately proportional to $1/\alpha_{nom}$. At the same time, the average number of Newton iterations, $\bar{N}_{newt}$, increases with the nominal step length. The proper choice of the nominal step length requires some level of experience and knowledge of the considered problem. The value $\alpha_{nom} = 10^{-2}$ yields a reasonable tradeoff between computational effort, robustness and accuracy in this case. For the analyzed benchmarks in this article the nominal step length was set to $\alpha_{nom} = 10^{-2}$, if not explicitly stated otherwise.

The performance of the cHB-ANM algorithm depends on the ANM threshold $\varepsilon_{\alpha}$, as illustrated in Fig. 3d-f. As one can easily deduce from Eq. (17), the larger $\varepsilon_{\alpha}$, the larger the step lengths, leading to fewer continuation steps $N_{pt}$. Throughout this work, the tolerance $\varepsilon_{tol}$ (residual norm) is set to $10^{-10}$. For this benchmark, setting $\varepsilon_{\alpha} = \varepsilon_{tol}$ leads to occasional Newton iterations. This comes from the accumulation of errors over the continuation steps. To avoid Newton iterations, we set $\varepsilon_{\alpha} = \varepsilon_{tol}/10 = 10^{-11}$ in the following if not explicitly stated otherwise. It is interesting to see that the computational effort can be further reduced when the ANM threshold is further increased. The ANM then degenerates to a PreCo algorithm with a high-order predictor and conventional Newton corrector iterations.
<figure><img src="image_8.png" /><figcaption>Figure 4: Duffing oscillator, effect of harmonic truncation order $H$: (a) overview of amplitude-frequency curve, (b) zoom into the super-harmonic resonance regime; $N = 4H + 1$ for AFT-PreCo</figcaption></figure>

### 3.1.2. Effect of harmonic truncation order $H$

In Fig. 4, the amplitude-frequency curves are depicted for different harmonic truncation orders $H$. To better analyze the convergence behavior, error measures are defined. Note that the overall goal of each algorithm is to determine a whole solution branch, not just a single point. It is therefore appropriate to introduce global error measures that are representative for a complete branch. Because of its high technical importance, we consider only the point of maximum amplitude (resonance peak), max $A_{RMS}(\eta)$, and determine the relative errors with respect to amplitude and frequency, $\varepsilon_A$ and $\varepsilon_\omega$, respectively. Various alternative error measures were analyzed as well, including the amplitude $A_{RMS}$ for a fixed frequency, and it was found that considering only the resonance peak yields representative results. To determine the maximum amplitude point, the solution path was expressed as a continuous function. For cHB-ANM, the available Taylor series expansion with respect to the arc length is used. For AFT-PreCo, a quadratic interpolation polynomial is constructed between the three solution points with largest amplitude. The convergence of these error measures with the harmonic truncation order is depicted in Fig. 5.

As expected, the results of both methods converge to the reference (cHB-ANM with $H = 200$) for sufficiently large harmonic truncation order $H$. An exponential convergence rate is theoretically expected for systems with analytic nonlinearities [33]. This would correspond to an error which has a slope in the log-log plots (Fig. 5a-b) that approaches minus infinity in the asymptotic limit case for large $H$. Due to finite numerical precision, however, a round-off plateau is reached in practice. For both methods, this appears for about $H = 10^1$. With cHB-ANM, the level $\varepsilon_A \approx \varepsilon_\omega \approx 10^{-12}$ is reached for an ANM threshold $\varepsilon_\alpha = 10^{-11}$. For a much higher ANM threshold $\varepsilon_\alpha = 10^{-6} \gg \varepsilon_{tol}$, the ANM degenerates to a predictor-corrector method with $P$-order predictor. This raises the error plateau to about $10^{-9}$. To reach a similar order of the interpolation error with AFT-PreCo, the nominal continuation step length $\alpha_{nom}$ must be decreased to a relatively small value. This reflects the advantage of the ANM of yielding a high-order representation of the whole solution path.

For the same $H$, one might think that both methods yield the same truncation error. However, the results generally differ. It is interesting to note that a reasonable accuracy can be achieved with the AFT HB already for $H = 1$ near the primary resonance and for $H = 3$ near the 1 : 3 super-harmonic resonance

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
(cf. Fig. 4). In contrast, the corresponding results for cHB considerably deviate from the reference. This can be attributed to the fact that the nonlinear terms are represented in a different way when AFT HB is applied immediately to the equations of motion, as compared to applying cHB to the quadratically recast system. In AFT HB, the nonlinear terms are represented by the $N$ time samples (recall that $N = 4H + 1$ here, where $H$ is the harmonic truncation order). In cHB, the nonlinear terms are represented by the $2H + 1$ Fourier coefficients (since the nonlinear terms are also approximated by a $H$-order Fourier series). Near the primary resonance, the main variable $q_1$ is dominated by the fundamental harmonic, and consequently $q_1^2$ should be dominated by the zeroth and the second harmonic. Therefore, the defining equation $r_{\text{aux},1} = v_1 - q_1^2 = 0$

<figure><img src="image_9.png" /><figcaption>Figure 5: Duffing oscillator, effect of harmonic order $H$: relative error of (a) amplitude and (b) frequency of the resonance peak, (c) computation time; $N = 4H + 1$ for AFT-PreCo; cHB-ANM\* shall indicate that in this case $\epsilon_\alpha \gg \epsilon_{\text{tol}}$ such that the ANM degenerates to a predictor-corrector method with $P$-order predictor</figcaption></figure>

for the auxiliary variable $v_1$ is not represented well for $H = 1$. In general, $q_1$ only responds with odd harmonics in this benchmark problem, such that the auxiliary variable $v_1$ should only respond with the zeroth and higher even harmonics. Thus, the error for the quadratically recast system generally changes when $H$ is changed, while additional higher even harmonics are not relevant when only the equation of motion is directly considered (as done with AFT-PreCo). For instance, the AFT-PreCo error remains the same when increasing $H$ from 1 to 2, but decreases for $H = 3$ (highlighted by pink diamonds in Fig. 5a-b).

The computational effort is depicted in Fig. 5c. To estimate the time complexity of the analyzed methods, we provide asymptotes with different scaling behavior (e. g. $O(H)$, $O(H \log H)$) in corresponding plots throughout this study. Here, the computation time increases significantly with $H$, for large $H$ slightly super-linearly for both methods. The number of Jacobian factorizations, $N_{\text{fact}}$, is also given for some points in Fig. 5c. Recall that in the case of the ANM this equals the number of continuation steps, i. e., the number of Taylor series expansion points, $N_{\text{fact}} = N_{\text{pt}}$, and the number of continuation steps (equal to the number of generated solution points) times the average number of Newton iterations in the case of the PreCo method, $N_{\text{fact}} = N_{\text{pt}} \bar{N}_{\text{newt}}$. Here, the number of continuation steps for the PreCo method is about 300, which is slightly larger than that of the ANM. The larger step lengths within the ANM are explained by the higher order of the predictor. In accordance, $\bar{N}_{\text{newt}} \approx 2-3$ here. Hence, the ANM requires fewer Jacobian factorizations, here by a factor of about 4 to 9. But apparently this does not lead to substantially shorter computation times. This is explained by the larger problem dimension: The quadratic recast introduces
an auxiliary variable $v_1$ (here $N_{\text{DOF}} = 1$). Thus, we have $N_{\text{eq}} = 2(2H + 1) + 1$ for cHB-ANM vs. only $N_{\text{eq}} = (2H + 1) + 1$ for AFT-PreCo.

As mentioned before, the computational effort should be viewed in the light of *implementation bias*. The efficiency of each algorithm is affected by how well it is implemented (use of parallelization, vectorization of loops within Matlab, etc.) for the given computing architecture. Although we unified the implementation where appropriate, we are convinced that the implementation bias can easily account for an order of magnitude of computational effort. Here, the computation time is always in the same order of magnitude for both methods (for $\alpha_{\text{nom}} = 10^{-2}$ and $\varepsilon_{\alpha} = 10^{-11}$). However, we conclude that the cHB-ANM is computationally more robust for this benchmark problem, since it is less sensitive to the threshold $\varepsilon_{\alpha}$, while AFT-PreCo is quite sensitive to the nominal step length $\alpha_{\text{nom}}$.

### 3.1.3. Effect of number of degrees of freedom $N_{\text{DOF}}$

<figure><img src="image_10.png" /><figcaption>Figure 6: System with polynomial stiffness: (a) amplitude-frequency curves, (b) computation time</figcaption></figure>

Now the effect of the number of modal degrees of freedom, $N_{\text{DOF}} \ge 1$, is analyzed, and the excitation level is set to $\hat{f}_c(1) = 10$. In Fig. 6a, the frequency response is illustrated for $N_{\text{DOF}} = 5$ and $H = 9$. Compared to the Duffing oscillator, additional primary and secondary resonances occur, as expected.

In Fig. 6b, the computation time is depicted vs. $N_{\text{DOF}}$. Analogously to the above discussed case of $N_{\text{DOF}} = 1$, the quadratic recast now introduces one auxiliary variable per modal coordinate. Thus, we have $N_{\text{eq}} = 2N_{\text{DOF}}(2H + 1) + 1$ for cHB-ANM vs. only $N_{\text{eq}} = N_{\text{DOF}}(2H + 1) + 1$ for AFT-PreCo. The problem dimension of cHB-ANM is, hence, approximately twice that of AFT-PreCo. A change of slope can be seen in Fig. 6b for AFT-PreCo at about $N_{\text{DOF}} = 10$. This indicates that the bottleneck of the algorithms changed. For large problem dimensions, the factorization of the Jacobian becomes the bottleneck. The effort for this grows quadratically with $N_{\text{DOF}}$ for both methods. Consequently, AFT-PreCo is about one order of magnitude faster than cHB-ANM for $N_{\text{DOF}} = 10$, but this difference reduces to a factor of about 4 up to $N_{\text{DOF}} = 99$ (1400 s for AFT-PreCo and about 6000 s for cHB-ANM).
### 3.2. System with unilateral stiffness

We now consider the equation of motion

$$q'' + 0.1q' + q + f_u(q) - 0.2 \cos(\eta t) = 0 \quad (23)$$

with the unilateral stiffness term $f_u(q) = 100 \max(q-1, 0)$. If $q < 1$, the spring is in separated state, $f_u = 0$. If $q > 1$, the spring is in contact state, and a force proportional to the spring deformation $q-1$ is built up. $f_u$ is only continuous in $q$, its derivative has jumps, it is therefore non-smooth. As discussed before, cHB cannot be applied since the Fourier coefficients of the nonlinear force cannot be expressed analytically¹. Also, it is not reasonable to apply the ANM directly, since the high-order derivatives required to expand the Taylor series simply do not exist. To apply cHB-ANM, $f_u$ must be regularized to an analytic function. Here, we use

$$f_u(q) \approx f_{u,\varepsilon_{\text{reg}}}(q) = \frac{100(q-1)}{2} \pm \sqrt{\left(\frac{100(q-1)}{2}\right)^2 + \varepsilon_{\text{reg}}}. \quad (24)$$

The smaller the regularization parameter $\varepsilon_{\text{reg}}$, the better the approximation of the unilateral stiffness term. For the quadratic recast, one auxiliary variable $f_{u,\varepsilon_{\text{reg}}}$ is introduced, along with one additional equation $r_{\text{aux}} = [f_{u,\varepsilon_{\text{reg}} - 100(q-1)}] f_{u,\varepsilon_{\text{reg}} - \varepsilon_{\text{reg}}}$.

The reference solution was developed as follows: Without regularization, the differential equation (23) is piecewise linear. Hence, an analytical expression for the general solution in the separated and the contact state can be easily derived. The main challenge is then to determine the contact state transition times (contact↔separation). The condition of a contact state transition is $q-1=0$. Substituting each general solution into this equation yields an implicit equation for the next contact state transition time. This equation is nonlinear and transcendental, and numerical root finding is therefore used to determine each transition time. This procedure is well-known and often called *stitching*. To find periodic solutions, it is required that the coordinate and velocity at the end of an excitation period are the same as those at the beginning; i. e., the problem is cast into a two-point periodic boundary value problem (shooting method). Finally, the periodic solutions are continued using the PreCo method.

#### 3.2.1. Effect of regularization parameter $\varepsilon_{\text{reg}}$

We studied the influence of the regularization parameter on the cHB-ANM solution, $H = 20$, by applying different values down to $\varepsilon_{\text{reg}} = 10^{-10}$. In this benchmark, cHB-ANM was able to continue the whole branch even for very steep regularizations, when setting the error tolerances properly. Moreover, the computational effort showed no significant dependence on $\varepsilon_{\text{reg}}$. For the following analyses, we choose $\varepsilon_{\text{reg}} = 10^{-9}$, which yields a highly accurate approximation of the non-smooth force law. It is important to understand that the quadratic recast of $f_u$ given below Eq. (24) generates an ambiguity: Besides the desired solution with the + sign in Eq. (24), which approximates the unilateral stiffness behavior, we get an artificial solution with the (-) sign in Eq. (24), which gives $f_{u,\varepsilon_{\text{reg}}}(q) \approx 0$. During the continuation procedure, this ambiguity comes to light by bifurcation points, cf. Fig. 7, where the unilateral spring would be activated. The artificial solution yields the quasi-linear behavior indicated as dashed curve in Fig. 7. The cHB-ANM implementation requires the user to manually select one of the emanating branches to be further continued. When the quasi-linear branch is continued further, we encountered additional bifurcation points not depicted in Fig. 7. To avoid these kind of branching problems in subsequent analyses, we apply a predefined rule to automatically choose the correct (physical) solution path. The rule comprises that as soon as cHB-ANM detects a bifurcation point, the emanating branch has to be continued in positive direction with regard to excitation frequency $\eta$. Note that for switching reliably the branch at a bifurcation point we had to reduce the ANM threshold to $\varepsilon_{\text{reg}} = 10^{-13}$.
<figure><img src="image_11.png" /><figcaption>Figure 7: SDOF oscillator with unilateral stiffness: bifurcations and non-physical solutions caused by the selected regularization and quadratic recast</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8: SDOF oscillator with unilateral stiffness, effect of samples per period *N* applied to the non-smooth system: error of AFT HB using (a) *H* = 1 and (b) *H* = 7, (c) computation time for AFT-PreCo</figcaption></figure>

### 3.2.2. Effect of number of samples per period *N*

In contrast to the system with polynomial stiffness, the sampling procedure within the AFT scheme is not exact beyond a certain value *N*: The unilateral stiffness function $f_u(q)$ generates an infinite sequence

<footnote n="1">¹A piecewise analytical integration is possible. However, the transition times are the roots of high-order polynomial equations, which can only be approximated numerically, see e. g. [34].</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 9: SDOF oscillator with unilateral stiffness: amplitude-frequency curves for (a) different harmonic truncation orders *H* and (b) different numbers of time samples *N*; reference and AFT-PreCo consider the non-regularized, cHB-ANM the regularized force law</figcaption></figure>

of (decaying but) nonzero Fourier coefficients $\hat{f}_u(k)$ already for harmonic input $q$. Hence, aliasing errors cannot be avoided within the AFT scheme. In Fig. 8a-b, the influence of $N$ on the relative error of the amplitude $\epsilon_A$ is depicted. For this we solved the AFT HB residual at a fixed excitation frequency close to the resonance peak. This way, interpolation errors inevitable with the PreCo method are avoided. The error decreases with $N$ but does not reach any plateau at the numerical precision level in the depicted range. $N$ also influences the computational effort, as shown in Fig. 8c. For large $N$, the computational effort is dominated by the calculation of the Fourier coefficients of the nonlinear force. The bottleneck of the algorithm are then the forward and inverse fast Fourier transforms, the effort of which increases approximately as $N \log(N)$.

In Fig. 9b, results are shown for different $N$ when using the AFT-PreCo algorithm with $H = 100$. Aliasing errors may lead to divergence and numerical oscillations of the HB approximation. In this case, $N \ge 603$ was sufficient to ensure convergence. Even in the converged results, numerical oscillations occur, caused by the aliasing-induced non-physical energy transfer from high to low frequencies. To better see these oscillations, the nominal continuation step length was decreased to $\alpha_{nom} = 10^{-4}$. This effect decreases for higher sampling rates. To the authors' knowledge, there is no rigorous theory available to a priori determine a useful value for $N$ that ensures convergence and reduces the aliasing error below a specified threshold. Clearly, the optimal choice of $N$ depends on the mechanical system, the considered dynamic regime (affected e. g. by the excitation level). For the following analyses, we considered two different settings for $N$: a fixed, high number $N = 2^{13}$ and the empirical rule $N = N_{\text{emp}}(H) = \min(500 + 25 \cdot H, 2000)$. A lower bound of $N_{\text{emp}}(0) = 500$ ensures convergence and reasonable accuracy for low numbers of $H$, while the upper bound $N_{\text{emp}}(H \ge 60) = 2000$ avoids unnecessary computational effort for very high $H$.

### 3.2.3. Effect of harmonic truncation order $H$

In Fig. 9a, the amplitude-frequency curve is depicted for different harmonic truncation orders $H$. The error measures, $\epsilon_A$ and $\epsilon_\omega$ (relative errors of resonance peak amplitude and frequency), introduced for the system with polynomial stiffness are again considered to illustrate the convergence behavior in Fig. 10a-b. Again, we observe that for the same (small) $H$, cHB-ANM yields larger errors than AFT-PreCo, since the auxiliary variable is only represented by the $2H + 1$ Fourier coefficients within cHB, whereas the nonlinear
<figure><img src="image_14.png" /><figcaption>Figure 10: SDOF oscillator with unilateral stiffness, effect of harmonic truncation order H: relative error of (a) amplitude and (b) frequency of the resonance peak, (c) computation time</figcaption></figure>

forces are represented with more samples $N > 2H + 1$ within the AFT scheme. Using AFT-PreCo with $N_{\text{emp}}(H)$ leads to a more oscillatory decrease of the error than with $N = 2^{13}$. This oscillatory behavior can already be inferred from Fig. 9b. For $N = 2000$, the numerical oscillations have still a relative amplitude of around $10^{-4}$ both in the amplitude and the frequency direction, respectively. This numerical oscillation amplitude is reflected in the error level depicted in Fig. 10a-b.

Compared to the system with polynomial stiffness, the truncation error decreases much more slowly. For about $H = 10^1$, the errors decrease below 1% for this benchmark problem, while $H = 1$ or 2 was sufficient for the system with cubic stiffness. While the error of the cHB-ANM still decreases (in average) beyond $H > 100$, the AFT-PreCo reaches a plateau. As opposed to the case of the system with polynomial stiffness, this plateau is driven by both interpolation and aliasing errors.

As can be seen in Fig. 10c, AFT-PreCo is faster than cHB-ANM when using the empirical rule $N_{\text{emp}}(H)$, but largely slower when setting $N = 2^{13}$ in the depicted range. In both cases, the number of Jacobian factorizations $N_{\text{fact}}$ are larger for AFT-PreCo compared to cHB-ANM (cf. Section 3.1.2). But here the selected number of AFT time samples $N$ decides whether the lower number of Jacobian factorizations $N_{\text{fact}}$ (but also larger problem size) within the cHB-ANM can be compensated or not. Overall, the differences are again not considered substantial (less than one order of magnitude). One should remark the ANM threshold was set to $\epsilon_{\alpha} = 10^{-13}$, as opposed to the other benchmarks, not only to avoid Newton iterations but also to ensure that the physical solution branch is continued, cf. Section 3.2.1.

### 3.3. System with dry friction

We now consider a single degree-of-freedom oscillator with elastic dry friction element, described with the equation of motion

$$q'' + 0.02q' + q + f_{\rho} - \hat{f}_{c}(1) \cos(\eta t) = 0, \qquad (25)$$
where the friction force $f_\rho$ is defined by the differential law

$$d f_{\rho} = \begin{cases} \kappa d q & |f_{\rho} + \kappa d q| < \rho \\ 0 & \text{otherwise} \end{cases}. \qquad (26)$$

The elastic dry friction element can be interpreted as a linear spring with stiffness $\kappa$ connected to a Coulomb slider. If the magnitude of the force in the spring remains smaller than $\rho$, the Coulomb slider is sticking, and the friction element acts as preloaded linear spring. Once the magnitude of the force reaches $\rho$, the Coulomb slider slides in the direction of $q'$ and the friction element produces a constant force $\rho$ in the opposite direction. The current value of $f_\rho$ not only depends on $q$ and $q'$, but also on their time history (hysteresis), hence, $f_\rho$ cannot be expressed as a unique function of $q$ and $q'$. It is, however, still possible to account for the differential law in Eq. (26) within the AFT scheme [35]: One starts with arbitrary initial conditions, e. g., $f_\rho = 0$, and determines the evolution of $f_\rho$ in a time-discrete way until the steady hysteresis cycle is reached (two periods were sufficient here). Thus, only $q$ has to be expanded in a truncated Fourier series within AFT-PreCo.

As in the case of the unilateral stiffness, neither cHB nor the ANM can be directly applied to Eqs. (25)-(26), since the Fourier coefficients of the nonlinear force cannot be expressed in closed form, and the system is non-smooth, respectively. Therefore, a regularization is applied, which is set up in two steps: First, the coordinate of the Coulomb slider, $g$, is introduced, and the signum-function describing the Coulomb law is approximated as $\operatorname{sgn} g' \approx g'/\sqrt{\varepsilon_{\text{reg}} + (g')^2}$,

$$q'' + 0.02q' + q + \kappa (q - g) - \hat{f}_c(1) \cos(\eta t) = 0, \qquad (27)$$

$$\kappa (g - q) + \rho \frac{g'}{\sqrt{\varepsilon_{\text{reg}} + (g')^2}} = 0. \qquad (28)$$

For the quadratic recast, the four auxiliary variables $g$, $u$, $f_\rho$ and $v_4$ are introduced, and Eq. (28) is replaced by the four equations

$$\kappa(g - q) + f_{\rho} = \qquad 0, \qquad (29)$$

$$u - g' = \qquad 0, \qquad (30)$$

$$v_4 f_{\rho} - \rho u = \qquad 0, \qquad (31)$$

$$v_4^2 - u^2 - \varepsilon_{\text{reg}} = \qquad 0. \qquad (32)$$

### 3.3.1. Effect of regularization parameter $\varepsilon_{\text{reg}}$

In Fig. 11, the amplitude-frequency curves are depicted for different values of the regularization parameter $\varepsilon_{\text{reg}}$, and three different excitation levels $\hat{f}_c(1)$. Here and in the following, we set $\kappa = 3$ and $\rho = 1$. As for the system with unilateral stiffness, the shooting method is used as reference, with analytical integration of the piecewise linear system, Eqs. (25)-(26), and direct calculation of the transition times between stick and slip friction (not regularized).

In contrast to the unilateral stiffness benchmark, we observed that as the regularization becomes steeper, i. e. with decreasing regularization parameter $\varepsilon_{\text{reg}}$, the number of harmonics in the cHB-ANM has to be increased to continue the whole branch. This is analyzed in more detail in Section 3.3.3. In Fig. 11, $H = 20$, and cHB-ANM continues the whole branch only if $\varepsilon_{\text{reg}}(\hat{f}_c(1) = 0.5) \ge 3 \cdot 10^{-4}$. (Using $H = 100$ steeper regularizations down to $\varepsilon_{\text{reg}}(\hat{f}_c(1) = 0.5) = 5 \cdot 10^{-5}$ are possible.) For the higher excitation level, we even need $\varepsilon_{\text{reg}}(\hat{f}_c(1) = 1.1) \ge 5.5 \cdot 10^{-3}$. Compared to the unilateral stiffness, the remaining regularization error is not negligible, and particularly severe for higher excitation levels (more sliding friction). For $\varepsilon_{\text{reg}} = 10^{-5}$ (finest considered $\varepsilon_{\text{reg}}$), the regularization error is substantial for the highest excitation level only. We also applied AFT-PreCo to the regularized (but not quadratically recast) system, and did not encounter any convergence problems down to $\varepsilon_{\text{reg}} = 10^{-5}$ for $N \ge 45$ ($H = 20$). We thus conclude that the convergence problems are associated with the representation of the nonlinear terms induced by applying cHB to the quadratically recast differential-algebraic equation system.
<figure><img src="image_15.png" /><figcaption>Figure 11: SDOF oscillator with elastic dry friction element, effect of regularization $\varepsilon_{\text{reg}}$ on amplitude frequency curve for different forcing levels; reference considers the non-regularized force law</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12: SDOF oscillator with dry friction, effect of samples per period $N$ applied to the non-smooth system: error of AFT HB using (a) $H = 1$ and (b) $H = 7$, (c) computation time for AFT-PreCo</figcaption></figure>

### 3.3.2. Effect of number of samples per period $N$

If the non-smooth friction law in Eq. (26) is used, the AFT scheme inevitably leads to aliasing errors, as in the case of the unilateral stiffness. In Fig. 12a-b we analyzed the influence of $N$ on the relative error of the amplitude $\varepsilon_A$ for a fixed frequency $\eta = 1.7$. Here, aliasing effects can be minimized using much fewer
<figure><img src="image_17.png" /><figcaption>Figure 13: SDOF oscillator with elastic dry friction element, effect of harmonic truncation order *H* on amplitude-frequency curves: (a) some successful and (b) some failed computations</figcaption></figure>

time samples *N* compared to the unilateral stiffness benchmark. While *N* ≈ 60 (*H* = 1) is sufficient to decrease the error to ||ε<sub>A</sub> - ε<sub>A</sub>(*N* = 2·10<sup>4</sup>)|| = 10<sup>-3</sup>, about *N* ≈ 400 (*H* = 1) were required to reduce the error to about the same level in the unilateral stiffness benchmark.

Analogous to the benchmark with the unilateral stiffness, we use two settings of *N* for the following analyses: a fixed, high number *N* = 2<sup>11</sup> and the empirical rule *N* = *N*<sub>emp</sub>(*H*) = min(50 + 10·*H*; 500). Compared to the benchmark with the unilateral stiffness, much smaller sampling rates are thus used for the benchmark with dry friction.

3.3.3. Effect of harmonic truncation order *H*

For the results depicted in the following, the intermediate excitation level $\hat{f}_c(1) = 0.5$ is considered, and the regularization parameter is set to $\epsilon_{reg} = 3 \cdot 10^{-4}$ (smallest feasible with cHB-ANM for $H = 20$). In Fig. 13, the amplitude-frequency curve is depicted for different harmonic truncation orders $H$. The convergence behavior is illustrated in Fig. 14a-b, in terms of the error measures introduced earlier.

Remarkably, $H = 1$ with $N_{\text{emp}} = 60$ is sufficient to achieve reasonable accuracy (errors < 1%) with AFT-PreCo applied to the non-smooth system. A plateau is reached for about $H \ge 13$, at an error level $\epsilon_A \approx 5 \cdot 10^{-5}$. Because of aliasing, $N_{\text{emp}}(H)$ leads to a slightly larger error $\epsilon_A$ than for $N = 2^{11}$ in a wide range, and to more pronounced numerical oscillations. Note that as opposed to $\epsilon_A$, $\epsilon_\omega$ stays almost constant. The reason for this is that, by coincidence, the PreCo continuation generates a point very close to the actual resonance peak (for the particular choice of the continuation limits and the specified nominal step length). Errors < 1% are never reached with cHB-ANM. The accuracy is mainly limited by the regularization error. The cHB-ANM error stabilizes beyond about $H = 26$.

It is important to note that the harmonic truncation order has to be set to $H \ge 62$ to ensure that the whole solution branch is computed with cHB-ANM. Fig. 13b shows some cases where the computation of the (physical) solution branch is unsuccessful. Branches of *ghost solutions* appear, which satisfy the condition $||\mathbf{R}|| < \epsilon_{\text{tol}}$, but are far from any reference solution. In other cases, the continuation gets stuck or turns around. Interestingly and in contrast to the benchmark with unilateral stiffness, the ANM did not detect any branching points. In general, the existence of a HB approximation is only guaranteed for sufficiently large
<figure><img src="image_18.png" /><figcaption>Figure 14: SDOF oscillator with elastic dry friction element, effect of harmonic truncation order $H$: relative error of (a) amplitude and (b) frequency of the resonance peak, (c) computation time</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: SDOF oscillator with elastic dry friction element, generalized and auxiliary variables computed with $H = 500$ using regularized force law ($\epsilon_{reg} = 3 \cdot 10^{-4}$) close to resonance peak: (a) one period of vibration and (b) decay of corresponding Fourier coefficients for $k \ge 1$</figcaption></figure>

$H$ [4]. This explains why the cHB can fail to converge in some cases. Apparently, the harmonic truncation order must be set to a comparatively large number, $H \ge 62$, to ensure convergence. This number is cer-
tainly driven by the required accuracy for representing certain auxiliary variables of the quadratically recast problem. To provide further insight into this, the time evolution and the Fourier spectrum of the auxiliary variables are depicted in Fig. 15. In the non-regularized case, the velocity of the Coulomb slider jumps at a stick-to-slip transition. This results in the almost discontinuous behavior of $u(\tau)$ (and consequently $v_4(\tau)$) in the regularized case, associated with a Gibbs-type phenomenon and slow decay of the corresponding Fourier coefficients. In contrast, the friction force and the deformation of the spring $\kappa$ are continuous ($C^0$) functions, even in the non-regularized case, yielding slightly better convergence of the Fourier series. Even better is the Fourier representation of the coordinate $q(\tau)$, which is a $C^2$ function.

We also observed that the ANM algorithm applied Newton correction iterations to reach some solution points, even in those runs where the whole branch was successfully continued. This means that the high-order Taylor series predictor produced unacceptably large residuals (although the ANM threshold $\epsilon_\alpha$ is already set to one tenth of the tolerance $\epsilon_{tol}$).

The non-negligible regularization error and the poor numerical robustness are already strong arguments against applying cHB-ANM to the considered dry frictional system. Yet another argument is the computational effort (see Fig. 14c): AFT-PreCo takes less than 0.5 s to compute the frequency response of the non-regularized problem with $H = 1$, $N_{emp} = 60$. This is about 25 times faster than cHB-ANM applied to the regularized problem with $H = 26$. As opposed to the previous benchmarks, not only the increased problem size (quadratic recast and higher number of harmonics) of cHB-ANM prevents better computational performance compared to AFT-PreCo but also the number of Jacobian factorizations $N_{fact}$: 491 needed for cHB-ANM vs. only 164 for AFT-PreCo when using $H = 100$. We observed, when applying the same regularization according to Eq. (28) to AFT-PreCo the required number of Jacobian factorizations also increases significantly with $H$ (e. g. about $N_{fact}(H = 100) = 714$).

## 4. Concluding guidelines

In this work, we compared the performance of the cHB-ANM and the AFT-PreCo algorithm for harmonically driven systems with either polynomial stiffness, unilateral spring or elastic dry friction nonlinearity. On the basis of the results, we derive the following conclusions.

If the considered differential-algebraic equation system contains non-smooth terms, arising e. g. from contact constraints, the AFT-PreCo algorithm is highly recommended. To apply cHB or the ANM, a regularization would be necessary to approximate the non-smooth terms by analytic ones. Additionally, a very steep regularization can be required to limit the error made by this approximation, as in the case of the benchmark problem with dry friction. A prohibitively high harmonic truncation order might then be needed to ensure convergence of the cHB-ANM algorithm. In contrast, reasonable accuracy can often be achieved with the AFT-PreCo already with a very low truncation order, e. g. $H = 1$ or $H = 3$. This is because the AFT-PreCo algorithm treats only the main variables (usually the generalized coordinates) as unknowns, which can often be represented with only a few Fourier terms, while the nonlinear terms are treated internally only, and represented by (potentially a much larger number of) time samples.

There are two important drawbacks of the AFT-PreCo algorithm: first, the AFT scheme yields inevitable aliasing errors for non-polynomial nonlinear terms. If the number of samples per period is too small, numerical oscillations or even divergence are possible. To the authors' knowledge there is no rigorous theory available to select this number in accordance with a priori estimated aliasing errors, and empirical rules are commonly used. Second, the PreCo continuation step length has to be properly selected by the user. Large step lengths may lead to a poor resolution (leading to inevitable interpolation errors when behavior between discrete solution points is of interest) or even jumps over important features of the solution path. Small step lengths increase the computational effort without providing significant additional information.

If the considered differential equation system contains only analytic terms (which are not designed to approximate non-smooth terms), the cHB-ANM algorithm is recommended. To compute the required coefficients of the high-order Taylor series expansion efficiently, a quadratic recast of the differential equation system is commonly applied (which has to be done by the user). For this recast, auxiliary variables and additional equations need to be introduced, which increases the problem dimension. In some cases, this recast can
introduce mathematical ambiguities (it is not always bijective) leading to artificial bifurcations, which requires some user interaction. Recently developed algorithms can then be used to automatically set up the harmonic balance equations and carry out the continuation. The important advantage of this procedure is that the continuation step length is automatically chosen depending on the range of utility of the truncated Taylor series. This makes the procedure numerically robust (the step length does not have to be selected by the user) and computationally competitive (since correction iterations are avoided). The cHB-ANM algorithm appears useful for smooth nonlinearities, for instance geometrically nonlinear finite element models. For large numbers of degrees of freedom ($N_{DOF}$), the Jacobian factorization becomes the bottleneck of both algorithms, the computational effort then increases approximately quadratically with $N_{DOF}$.

In the future, it would be interesting to extend the cHB-ANM algorithm in such a way that different variables are resolved with individual harmonic truncation orders. This would permit to combine a more precise representation of the nonlinear terms (idea behind AFT HB) with the sophisticated continuation framework of the ANM. Furthermore, the aforementioned conclusions on the performance of the cHB-ANM (high harmonic truncation order needed, mathematical ambiguities) for initially non-smooth problems were drawn for the considered regularization; it remains to be investigated, how different regularizations affect the cHB-ANM performance. Moreover, it would be interesting to extend the comparison of cHB-ANM and AFT-PreCo algorithms to bifurcation analyses.

## Acknowledgement

This work presents results of the Tribomechadynamics Research Camp (TRC), formerly called Nonlinear Dynamics of Coupled Structures and Interfaces summer school, an annual month long international research collaboration for graduate students and postdoctoral researchers. For more information, visit http://tmd.rice.edu.

The author Fabia Tubita is grateful for the funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 721865.

The authors Bruno Cochelin, Louis Guillot and Christophe Vergez would like to thank the Labex MEC (ANR-10-LABX-0092) and 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).

## Appendix A. On the exactness of the AFT scheme for polynomials

The goal of the AFT scheme is to compute the Fourier coefficients of the nonlinear forces up to the considered truncation order $H$. In this appendix, we show that the AFT gives the exact Fourier coefficients in the special case of polynomial nonlinear forces for a sufficiently large number of samples $N$. To this end, we consider a scalar function

$$f(q) = q^P, \qquad (A.1)$$

i.e., a monomial in $q$ of order $P$. The following argument can be extended to arbitrary polynomials in a straight-forward way. In accordance with HB, we use as ansatz for $q$ a Fourier series truncated to order $H$,

$$q = \sum_{k=-H}^{H} \hat{q}(k)e^{ik\tau}, \qquad (A.2)$$

with the normalized time $\tau = \Omega t$. In this appendix, we use the complex-exponential representation of a Fourier series, which is, of course, fully equivalent to the sine-cosine representation e. g. in Eq. (2). To ensure that the sum in Eq. (A.2) yields a real result, we must have $\hat{q}(-k) = \hat{q}(k)$.

One can easily see that substituting the $H$-order Fourier series in Eq. (A.2) into the $P$-order monomial in Eq. (A.1) generates a $PH$-order Fourier series,

$$f = \sum_{k=-PH}^{PH} \hat{f}(k)e^{ik\tau}. \qquad (A.3)$$
For instance, a cubic term ($P = 3$) generates at most Fourier terms with frequency $3\Omega$.
The Fourier coefficients $\hat{f}(k)$ in Eq. (A.3) are defined by the integral

$$ \hat{f}(k) = \frac{1}{2\pi} \int_{0}^{2\pi} f(q(\tau))e^{-ik\tau} d\tau. \qquad (A.4) $$

The discrete Fourier transform (DFT) generally gives an approximation of the continuous integral in Eq. (A.4),

$$ \hat{f}(k) \approx \frac{1}{N} \sum_{n=0}^{N-1} \tilde{f}(n) e^{-ik \frac{2\pi n}{N}}, \qquad (A.5) $$

where $\tilde{f}(n)$ are samples, $\tilde{f}(n) = f(q(2\pi n/N))$, and $N$ is the number of samples. Eq. (A.5) can be identified as applying the trapezoidal rule to approximate the continuous integral in Eq. (A.4) [3].

We now show that if $f$ is a truncated Fourier series, Eq. (A.5) gives the exact value of the continuous integral in Eq. (A.4), and therefore the exact Fourier coefficients, for sufficiently large $N$.

$$ \begin{aligned} \frac{1}{N} \sum_{n=0}^{N-1} \tilde{f}(n) e^{-ik \frac{2\pi}{N} n} &= \frac{1}{N} \sum_{n=0}^{N-1} \left( \sum_{m=-PH}^{PH} \hat{f}(m) e^{im \frac{2\pi}{N} n} \right) e^{-ik \frac{2\pi}{N} n} \\ &= \sum_{m=-PH}^{PH} \hat{f}(m) \sum_{n=0}^{N-1} \frac{1}{N} e^{i2\pi \frac{m-k}{N} n} \\ &= \sum_{m=-PH}^{PH} \hat{f}(m) \sum_{n=0}^{N-1} \frac{1}{N} \left( e^{i2\pi \frac{m-k}{N}} \right)^n \\ &= \sum_{m=-PH}^{PH} \hat{f}(m) \begin{cases} 1 & |m-k| \bmod N = 0 \\ 0 & \text{otherwise} \end{cases}. \end{aligned} \qquad (A.6) $$

Herein, a b denotes the modulo operation which determines the remainder of the division of a by b. Thus, the expression $|m-k| \bmod N = 0$ means that $|m-k|$ is an integer multiple of $N$. To verify the last step of Eq. (A.6), we use the geometric regression lemma:

$$ \begin{aligned} s_N &= \sum_{n=0}^{N-1} \alpha z^n, \\ \Rightarrow s_N - zs_N &= \alpha - \alpha z^N, \\ \Rightarrow s_N &= \alpha \frac{1-z^N}{1-z}. \end{aligned} \qquad (A.7) $$

With $z = e^{i \frac{2\pi}{N} \ell}$ and $\alpha = 1/N$ we obtain

$$ \sum_{n=0}^{N-1} \frac{1}{N} \left( e^{i \frac{2\pi}{N} \ell} \right)^n = \frac{1}{N} \frac{1 - e^{i 2\pi \ell}}{1 - e^{i \frac{2\pi}{N} \ell}} = \begin{cases} 1 & |\ell| \bmod N = 0 \\ 0 & \text{otherwise} \end{cases}. \qquad (A.8) $$

The expression in the middle is actually indeterminate for $|\ell| \bmod N = 0$, but its limit for $|\ell| \bmod N \to 0$ is well-defined and given on the right. As alternative to taking the limit, the sum can be easily calculated noticing that $e^{i \frac{2\pi}{N} \ell} = 1$ if $|\ell|$ is an integer multiple of $N$.

Now, the last line of Eq. (A.6) can be verified by setting $\ell = m - k$. To avoid aliasing, we need to ensure that $N$ is sufficiently large so that $|m-k|$ never reaches $N$ (or any higher integer multiple), and thus
$|m-k| \mod N = 0$ is only achieved for $m=k$. To set up the HB equations, we need only the Fourier coefficients of $f$ up to order $H$. Thus, the largest possible value of $|m-k|$ is obtained when $m = +PH$ and $k = -H$, or when $m = -PH$ and $k = +H$. In both cases, $|m-k| = (P+1)H$. As a consequence, for $N \ge (P+1)H$ the sum on the right hand side of Eq. (A.5) yields exactly $\hat{f}(k)$ for $k = -H, \dots, H$ (no aliasing errors). It is important to emphasize that $N = (P+1)H$ is not sufficient to correctly compute the higher-order Fourier terms of $f$, i. e., $\hat{f}(k)$ with $H < |k| \le PH$. But these are not needed for setting up the HB equations.

In summary, we can say that the AFT scheme for computing the Fourier coefficients of the nonlinear forces is exact if $N \ge 3H + 1$ for polynomial nonlinear terms of up to quadratic order. Analogously, it is exact if $N \ge 4H + 1$ for polynomial nonlinear terms of up to cubic order.

## References

[1] E. Sarrouy, J.-J. Sinou, Non-linear periodic and quasi-periodic vibrations in mechanical systems-on the use of the harmonic balance methods, in: Advances in Vibration Analysis Research, Intech, 2011.

[2] T. Detroux, L. Renson, G. Kerschen, The harmonic balance method for advanced analysis and design of nonlinear mechanical systems, Nonlinear Dynamics, Volume 2: Proceedings of the 32nd IMAC, A Conference and Exposition on Structural Dynamics (2014).

[3] M. Krack, J. Gross, Harmonic Balance for Nonlinear Vibration Problems, Springer, 2019. doi:10.1007/978-3-030-14023-6.

[4] M. Urabe, Galerkin's procedure for nonlinear periodic systems, Archive for Rational Mechanics and Analysis 20 (2) (1965) 120-152.

[5] E. P. Petrov, Stability analysis of multiharmonic nonlinear vibrations for large models of gas turbine engine structures with friction and gaps, Journal of Engineering for Gas Turbines and Power 139 (2) (2016) 022508-022508-10. doi: 10.1115/1.4034353.

[6] L. Salles, L. Blanc, F. Thouverez, A. M. Gouskov, P. Jean, Dual time stepping algorithms with the high order harmonic balance method for contact interfaces with fretting-wear, Paper GT2011-46488, Proc. of GT2011, ASME Turbo Expo 2011: Advancing Clean and Efficient Turbine Technology, June 7-10, Vancouver, Canada (2011).

[7] S. Zucca, C. M. Firrone, Nonlinear dynamics of mechanical systems with friction contacts: Coupled static and dynamic multi-harmonic balance method and multiple solutions, Journal of Sound and Vibration 333 (3) (2014) 916-926. doi: 10.1016/j.jsv.2013.09.032.

[8] D. Süß, K. Willner, Investigation of a jointed friction oscillator using the multiharmonic balance method, Mechanical Systems and Signal Processing 52-53 (2015) 73-87. doi:10.1016/j.ymssp.2014.08.003.

[9] D. Laxalde, F. Thouverez, Complex non-linear modal analysis for mechanical systems: application to turbomachinery bladings with friction interfaces, Journal of Sound and Vibration 322 (4-5) (2009) 1009-1025. doi:10.1016/j.jsv.2008.11.044.

[10] G. Jacquet-Richardet, M. Torkhani, P. Cartraud, F. Thouverez, T. Nouri Baranger, M. Herran, C. Gibert, S. Baguet, P. Almeida, L. Peletan, Rotor to stator contacts in turbomachines. review and application, Mechanical Systems and Signal Processing 40 (2) (2013) 401-420. doi:10.1016/j.ymssp.2013.05.010.

[11] 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 (4-5) (2009) 520-540. doi:10.1016/j.jsv.2009.08.017.

[12] K. C. Hall, J. P. Thomas, W. S. Clark, Computation of unsteady nonlinear flows in cascades using a harmonic balance technique, AIAA Journal 40 (5) (2002) 879-886.

[13] R. J. Gilmore, M. B. Steer, Nonlinear circuit analysis using the method of harmonic balance—a review of the art. part i. introductory concepts, International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering 1 (1) (1991) 22-37.

[14] 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 (1) (1989) 149-154. doi:10.1115/1.3176036.

[15] B. Cochelin, C. Vergez, A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions, Journal of Sound and Vibration 324 (1-2) (2009) 243-262. doi:10.1016/j.jsv.2009.01.054.

[16] S. Karkar, B. Cochelin, C. Vergez, A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities, Journal of Sound and Vibration 332 (4) (2013) 968-977. doi:10.1016/j.jsv.2012.09.033.

[17] L. Guillot, B. Cochelin, C. Vergez, A generic and efficient taylor series-based continuation method using a quadratic recast of smooth nonlinear systems, International Journal for Numerical Methods in Engineering doi:10.1002/nme.6049.

[18] L. Guillot, B. Cochelin, C. Vergez, A taylor series-based continuation method for solutions of dynamical systems, Nonlinear Dynamics doi:10.1007/s11071-019-04989-5.

[19] M. Urabe, A. Reiter, Numerical computation of nonlinear forced oscillations by galerkin's procedure, Journal of Mathematical Analysis and Applications 14 (1) (1966) 107-140.

[20] E. Moussi, S. Bellizzi, B. Cochelin, I. Nistor, Nonlinear normal modes of a two degrees-of-freedom piecewise linear system, Mechanical Systems and Signal Processing (2013) 1-26.

[21] R. W. Klopfenstein, Zeros of nonlinear functions, Journal of the ACM (JACM) 8 (3) (1961) 366-373.
[22] E. Riks, The application of newton's method to the problem of elastic stability, Journal of Applied Mechanics 39 (4) (1972) 1060-1065. doi:10.1115/1.3422829.

[23] R. Seydel, Practical Bifurcation and Stability Analysis: from Equilibrium to Chaos, Springer New York, 1994.

[24] A. Cardona, T. Coune, A. Lerusse, M. Geradin, A multiharmonic method for non-linear vibration analysis, International Journal for Numerical Methods in Engineering 37 (9) (1994) 1593-1608.

[25] M. Jerschl, D. S. K. Willner, Numerical Continuation Methods for the Concept of Non-linear Normal Modes, 2014, pp. 19-26. doi:10.1007/978-3-319-04546-7_3.

[26] B. Cochelin, N. Damil, M. Potier-Ferry, The asymptotic-numerical method: an efficient perturbation technique for non-linear structural mechanics, Revue Européenne des Éléments Finis 3 (2) (1994) 281-297. doi:10.1080/12506559.1994.10511124.

[27] B. Cochelin, M. Medale, Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations, Journal of Computational Physics 236 (1) (2013) 594-607. doi:10.1016/j.jcp.2012.11.016.

[28] I. Charpentier, B. Cochelin, Towards a full higher order ad-based continuation and bifurcation framework, Optimization Methods and Software 33 (4-6) (2018) 945-962. doi:10.1080/10556788.2018.1428604.

[29] Manlab - an interactive path-following and bifurcation analysis software, available at https://manlab.lma.cnrs-mrs.fr.

[30] R. Arquier, B. Cochelin, C. Vergez, Manlab logiciel de continuation interactif (2005).

[31] Nlvib - a matlab tool for nonlinear vibration problems, available via https://www.ila.uni-stuttgart.de/nlvib.

[32] A. H. Nayfeh, D. T. Mook, Nonlinear Oscillations, Vol. 1979, John Wiley & Sons, New York.

[33] E. Tadmor, The exponential accuracy of fourier and chebyshev differencing methods, SIAM J. Numer. Anal. 23 (1) (1986) 1-10. doi:10.1137/0723001.

[34] M. Krack, L. Panning-von Scheidt, J. Wallaschek, A high-order harmonic balance method for systems with distinct states, Journal of Sound and Vibration 332 (21) (2013) 5476-5488. doi:10.1016/j.jsv.2013.04.048.

[35] J. Guillen, C. Pierre, An efficient hybrid frequency-time domain method for the dynamics of large-scale dry-friction damped structural systems, in: C. Glocker, F. Pfeiffer (Eds.), IUTAM Symposium on Unilateral Multibody Contacts, Vol. 72 of Solid mechanics and its applications, Kluwer Academic Publishers, Dordrecht, 1999, pp. 169-178.