# Input/output reduced model of a damped nonlinear beam based on Volterra series and modal decomposition with convergence results

Thomas Hélie, Béatrice Laroche

► To cite this version:

Thomas Hélie, Béatrice Laroche. Input/output reduced model of a damped nonlinear beam based on Volterra series and modal decomposition with convergence results. Nonlinear Dynamics, 2021, 105 (1), pp.515-540. 10.1007/s11071-021-06529-6. hal-03306198

HAL Id: hal-03306198

https://hal.science/hal-03306198

Submitted on 1 Oct 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Input/output reduced model of a damped nonlinear beam based on Volterra series and modal decomposition with convergence results

Thomas Hélie · Béatrice Laroche

Received: date / Accepted: date

**Abstract** This paper addresses the model reduction and the simulation of a damped Euler-Bernoulli-von Kármán pinned beam excited by a distributed force. This nonlinear problem is formulated as a PDE and reformulated as a well-posed state-space system. The model order reduction and simulation are derived by combining two approaches: a Volterra series expansion and truncation; a pseudo-modal truncation defined from the eigenbasis of the linearized problem. The interest of this approach lies in the large class of input wave-shapes that can be considered and in the simplicity of the simulation structure. This structure only involves cascades of finite dimensional decoupled linear systems and multilinear functions. Closed-form bounds depend-ing on the model coefficients and the truncation orders are provided for the Volterra convergence domain and the approximation error. These theoretical results are generalized to a large class of nonlinear models and re-finement of bounds are also proposed for a large sub-class. Numerical experiments confirm that the beam model is well approximated by the very first Volterra terms inside the convergence domain.

**Keywords** beam · partial differential equation · nonlinear systems · Volterra series · modal truncation · convergence · error bounds

## 1 Introduction

This paper addresses the model reduction and simulation of a beam with external excitation for fast simulation purpose. The model under consideration, proposed in [9,11], represents a damped nonlinear pinned beam. Under the simplifying assumptions of Euler-Bernoulli and Von Karman's kinematics coupled with viscous and structural damping, the resulting model is expressed as a nonlinear partial differential equation. Our objective is to design and simulate an accurate reduced causal model when the excitation is not known in advance and may have any waveshape: the excitation is then considered as the input of a beam system.

Model reduction for nonlinear systems, in particu-lar nonlinear PDEs, is an important problem in engi-neering sciences, that has been and is currently still thoroughly studied through a wide variety of methods. Among them, POD based approaches (including SVD, PCA and KLD decomposition [21],[1]) coupled with Galerkin projection are very popular (see for instance [17], [18] and more recently [12]). Koopman based meth-ods (see for instance [24], [3]) are also becoming increas-ingly used. These methods are data driven in the sense that they build approximate models based on values collected from numerous real or numerical experiments: to generate these data, the system is fed by a selected set of inputs or forcing signals and a set of initial con-ditions.

<footnote>T Hélie
CNRS, STMS laboratory (UMR 9912, IRCAM-CNRS-Sorbonne Université), Paris, France
Tel.: +33-1-44784824
Fax: +33-1-44781540
E-mail: thomas.helie@ircam.fr</footnote>

<footnote>B. Laroche
Université Paris-Saclay, INRAE, MaIAGE, 78350 Jouy-en-Josas, France
Université Paris-Saclay, Inria, Inria Saclay-Île-de-France, 91120, Palaiseau, France
Tel.: +33-1-34652843
E-mail: beatrice.laroche@inrae.fr</footnote>
In this work, we propose a different approach, based on the exact knowledge of the systems equations, and exploit two standard tools.

The first one is the Volterra series expansion. First introduced by Vito Volterra [29], these series have been widely used for signal processing, electronics, mechanics, acoustics (etc), for model order reduction and real-time simulation purposes, since they transform a non-linear dynamical system into a series of linear systems cascaded with multilinear interconnecting functions. There exists a vast literature concerning Volterra series for ODEs. They were studied for instance in [10, 8, 27, 6, 28]. Volterra series expansions can also be applied to some classes of nonlinear PDEs, using space-time kernels (see e.g. [26] for Green-Volterra series expansions) or using infinite-dimensional (semigroup) Volterra kernels [15] (see also [4] for a review). One important issue raised by this approach in practical applications is to find a bound on the system's input magnitude for which the series representation is convergent, and on the series truncation error. This is still an active domain of research as evidenced in the recent paper [30] and references therein (see also [2, 25, 20] for frequency domain criteria). This paper relies on computable convergence results valid for the class of bounded input signals, first established for linear-analytic ODE systems in [14], and improved in [15] where tighter convergence bounds and criterion are provided for both ODEs and PDEs.

The second standard tool that we have at hand for dynamical model reduction is the modal decomposition of linear PDEs, followed by truncation and projection on a finite dimension modal subspace spanned by the first modes. In this context, the external excitation is projected on this subspace, and the resulting system's trajectories correspond to a time varying excitation of the finite modal basis. The resulting reduced model becomes a set of finite dimension ODE systems, whose simulation is easy. This ability to approximate infinite dimensional linear systems by finite dimension linear ODE systems constitutes the main advantage of this approach.

In this paper, as our model is both nonlinear and infinite dimensional, we propose a method for combining both tools and exploit their respective advantages in order to obtain an approximated model composed of finite dimensional linear ODE systems and nonlinear interconnections functions.

The paper is organized as follows: the notations and mathematical results borrowed from [15] are briefly recalled in section 2, and the beam model is introduced in section 3, together with its state space representation. The linearized version of the model and its modal

decomposition are studied in 4. In this section, we also introduce the notion of pseudo-modal truncation that we use to formulate an approximated finite dimensional nonlinear model in the modal sub-spaces. Section 5 presents the reduced model and provides a convergence bound as well as an approximation error estimate, depending on the input signal. These results also are generalized for a wide class of nonlinear models. Then, in section 6, we develop new theoretical results that refine convergence bounds for specific classes of nonlinearities, and apply them to the beam model. Finally, numerical simulations showing the application of our method to two different parameter configurations are provided in section 7.

## 2 Preamble

In [14,15], we have proved that the trajectories of a class of input/output nonlinear systems can be decomposed into Volterra series expansions for any bounded input in a computable convergence domain. A truncation error bound can be also computed. This section reminds these results, based on which is addressed the case of a damped nonlinear pinned beam excited by a distributed force of any waveshape.

### 2.1 Class of nonlinear systems

We consider systems excited on positive times $T = \mathbb{R}_+$ by inputs $u: T \to \mathbb{U}$, with state-space representation of state $x: T \to \mathbb{X}$, governed by

$$ \dot{x}(t) = A x(t) + B u(t) + \sum_{k=2}^{K} A_k (x(t), \dots, x(t)), \quad (1) $$

with initial condition

$$ x(0) = x_{\text{ini}} \in \mathbb{X}, \quad (2) $$

where $A: \mathbb{X} \to \mathbb{X}$ and $B: \mathbb{U} \to \mathbb{X}$ are linear operators, $A_k: \mathbb{X}^k \to \mathbb{X}$ are multilinear operators (accounting for nonlinearities of homogeneous degree $k$ w.r.t. $x$) and where $A$ is assumed to generate stable dynamics for the linearized system (hypothesis H0).

Finite dimensional case ($\mathbb{X} = \mathbb{R}^n, \mathbb{U} = \mathbb{R}^m$). The state $x$ and the input $u$ can be vectors, in which case $A$ and $B$ can be represented by $n \times n$ and $n \times m$ matrices and the $A_k$'s are multilinear functions. Hypothesis (H0) means that $A$ is a Hurwitz stable matrix<fnref n="1" />: $A$ generates the

<footnote n="1">¹ Every eigenvalue of $A$ has strictly negative real part.</footnote>
flow $S(t) = \exp(At)$ on $[0, T]$, whose operator norm is bounded by $\beta \exp(\alpha t)$ for some $\beta > 0$ and a negative growth bound$^2 \alpha < 0$.

*Infinite dimensional case (X, U: Banach spaces).* The state $x$ and the input $u$ can also be space-dependent functions, in which case $A, B$ (resp., $A_k$) are linear (resp., multilinear) space operators in the spaces setting $\mathcal{L}(X)$, $\mathcal{L}(U, X)$ (resp., $\mathcal{ML}_k(X)$) detailed in appendix A. Technically, hypothesis (H0) means that $A$ generates a flow (called a strongly continuous semigroup) $S$ with negative growth bound $\alpha < 0$ and with $\beta > 0$ such that for all $t \in \mathbb{T}$, $\|S(t)\|_{\mathcal{L}(X,X)} \le \beta \exp(\alpha t)$.

*Systems under consideration.* The class of problems (1-2) with (H0) is examined in the framework of bounded signals, namely

$$u \in \mathcal{U} = L^{\infty}(\mathbb{T}, \mathbb{U}) \text{ and } x \in \mathcal{X} = L^{\infty}(\mathbb{T}, \mathbb{X}), \quad (3)$$

case for which well-posed definitions of solutions (after-named mild solutions) can be introduced (see details in appendix A, see also [23,7] and [15,13] for more general classes of systems).

## 2.2 Volterra series expansion

The trajectory $x$ can be decomposed as a sum of contributions with homogeneous of order $m$ w.r.t. to input $u$ and initial condition $x_{\text{ini}}$, given by

$$x(t) = \sum_{m=1}^{+\infty} x_m(t), \quad (4)$$

where the terms satisfy the sequence of linear problems

$$\begin{align}
\dot{x}_1(t) &= Ax_1(t) + B u(t) \quad \text{with } x_1(0) = x_{\text{ini}}, \tag{5} \\
\dot{x}_m(t) &= Ax_m(t) + \chi_m(t) \quad \text{with } x_m(0) = 0, \quad \text{for } m \ge 2, \tag{6}
\end{align}$$

$$\text{with } \chi_m(t) = \sum_{k=2}^{\min(m,K)} \sum_{p \in \mathbb{M}_m^k} A_k(x_{p_1}(t), \dots, x_{p_k}(t)), \quad (7)$$

where the index set $\mathbb{M}_m^k := \{p \in (\mathbb{N}^*)^k \mid p_1 + \dots + p_k = m\}$ selects occurrences $x_{p_i}$ such that their combination through the multi-linear operator $A_k$ is of order $m$.

Solutions are given by a sequence of convolutions with $S$

$$x_1(t) = S(t)x_{\text{ini}} + \int_0^t S(t-\tau) B u(\tau) \, d\tau, \quad (8)$$

$$x_m(t) = \int_0^t S(t-\tau) \chi_m(\tau) \, d\tau, \quad \text{for } m \ge 2, \quad (9)$$

that exactly generate a Volterra series expansion (see [15, Rk1]), whose kernels can be deduced from formula (8-9).

**Remark 1 (Link with the regular perturbation method)** Equation (8) isolates the classical solution of the linearized problem (1-2 without the $A_k$'s). Solution (4-9) formally results from the regular perturbation method [13] applied to the input and the initial condition, both marked by a common scalar $\varepsilon$ ($u = \varepsilon \tilde{u}, x_{\text{ini}} = \varepsilon \tilde{x}_{\text{ini}}$), by sorting terms along the powers of this marker $\varepsilon$ ($x_m = \varepsilon^m \tilde{x}_m$).

**Remark 2 (Simulation issues)** The truncated series is well-adapted to time-domain simulation, since the signals $x_m$ can be processed by combining linear filters (simulation of several occurrences of $\dot{y} = Ay + v$) and instantaneous nonlinearities (multilinear operators $A_k$ with $k \le m$ fed by subsets of yet-simulated signals $x_j$ with $j \le m-1$) [15, Rk2 and Fig.1] (see also fig. 3).

## 2.3 Convergence domain and truncation error bound

The series (4-9) converges towards a mild solution of (1-2) in the space $\mathcal{X}$ of bounded trajectories, if the solution (8) of the linearized problem is bounded as follows (see [14,15,5]):
if $\|x_1\|_{\mathcal{X}} < \rho$, then

$$\begin{align}
& \|x\|_{\mathcal{X}} < \Phi(\|x_1\|_{\mathcal{X}}), \quad (\text{convergence in norm}) \tag{10} \\
& \|x - \sum_{m=1}^{M} x_m\|_{\mathcal{X}} < R_M \Phi(\|x_1\|_{\mathcal{X}}), \quad (\text{error bound}) \tag{11}
\end{align}$$

where the convergence radius $\rho$, gain bound function $\Phi$ and error bound function $R_M\Phi$ are characteristic quantities of the system.

In practice, they can be derived as follows:

**Step 1:** Compute the upper bounds $\zeta_k$, for $2 \le k \le K$,

$$\zeta_k \ge \int_{\mathbb{T}} \|S(t) A_k\|_{\mathcal{ML}_k(X, X)} \, dt, \quad (12)$$

**Step 2:** Define the function $F$, analytic at $z = 0$,

$$F(z) = \frac{z}{z - \zeta(z)} \quad \text{with} \quad \zeta(z) = \sum_{k=2}^{K} \zeta_k z^k. \quad (13)$$

<footnote>$^2$ The maximal real part of the eigenvalues of $A$ is the smallest optimal bound $\alpha$.</footnote>
Step 3: Compute the unique root $\sigma$ defined by

$$ \sigma > 0 \text{ such that } \sigma F'(\sigma) - F(\sigma) = 0. \quad (14) $$

Step 4: Compute the convergence radius $\rho$, given by

$$ \rho = \sigma/F(\sigma) > 0. \quad (15) $$

Step 5: Find the unique function analytic at $z = 0$

$$ \Phi(z) = \sum_{m=1}^{\infty} \phi_m z^m \text{ such that } \Phi(z) = z F(\Phi(z)). \quad (16) $$

Step 6: Compute the remainder function given by

$$ R_M \Phi(z) = \Phi(z) - \sum_{m=1}^{M} \phi_m z^m. \quad (17) $$

Note that the coefficients $\phi_m$ are given by $\phi_1 = 1$ and, for all $m \ge 2$, $\phi_m = \sum_{k=2}^m \zeta_k \sum_{p \in M_m^k} \prod_{i=1}^k \phi_{p_i}$.

Remark 3 (Convergence radius on the input) Following [14,19], a convergence condition can be proposed on the system input for zero initial conditions. From (8), $\|x_1\|_X \le \gamma \|u\|_U$ with

$$ \gamma \ge \int_{\mathbb{T}} \|S(t)B\|_{L(U,X)} dt. \quad (18) $$

So, a sufficient condition for the series convergence is

$$ \|u\|_U \le \rho_u \text{ with } \rho_u = \rho/\gamma. \quad (19) $$

It is convenient in practice, since computing $x_1$ is not required. But it is conservative compared to $\|x_1\|_X < \rho$.

### 3 Physical model: nonlinear boundary problem, state-space representation and well-posedness

We consider an Euler-Bernoulli model of a damped non-linear pinned beam, initially at rest. First, the governing equations of the nonlinear boundary problem are presented in § 3.1. Then, the problem is reformulated as a well-posed state-space representation (1-2): the linearized problem is first examined in (§ 3.2), and then, the full nonlinear problem in (§ 3.3).

#### 3.1 Nonlinear boundary problem: governing equations

The considered Euler-Bernoulli model assumes the following hypotheses (see [9,11]):

(H1) Euler-Bernoulli kinematics (any cross-section before deformation remains straight after deformation);

(H2) Von Karman's coupling (between the axial and the bending movements) introducing a nonlinearity;

(H3) viscous and structural damping phenomena.

A dimensionless model governing the deflection waves $w: (z,t) \in [0,1] \times \mathbb{T} \to \mathbb{R}$ is given by $(\mathbb{T} = \mathbb{R}_+)$,

$$ \partial_t^2 w + 2(a+b\partial_z^4)\partial_t w + \partial_z^4 w - \eta \left( \int_0^1 (\partial_z w)^2 dz \right) \partial_z^2 w = f, \quad (20) $$

for all $(z,t) \in [0,1] \times \mathbb{T}$, with pinned-type boundary conditions (fixed extremities and no momentum)

$$ w(z,t) = 0 \text{ and } \partial_z^2 w(z,t) = 0 \text{ for all } (z,t) \in \{0,1\} \times \mathbb{T}, \quad (21) $$

and zero initial conditions

$$ w(z,0) = 0 \text{ and } \partial_t w(z,0) = 0 \text{ for all } z \in [0,1]. \quad (22) $$

Coefficients $a > 0$ and $b > 0$ are fluid and structural damping parameters and $\eta > 0$ is the nonlinear coupling coefficient between the bending momentum and the displacement under the von Kármán assumption [22].

The excitation is assumed to be a time dependent distributed force $f(z,t)$, whose space profile is square integrable, that is, $f(\cdot, t)$ belongs to the Lebesgue space

$$ \mathbb{H} := L^2(0,1) \text{ (square integrable functions),} \quad (23) $$

for each fixed time $t$. Moreover, we assume that

(H4) these integrals are bounded in time,

meaning that functions $t \mapsto f(\cdot, t)$ belong to $L^\infty(\mathbb{T}, \mathbb{H})$.

Note that, compared to structured forcing terms<fnref n="3" />, the range largeness of such excitation signals gives the following method one of its main interests.

#### 3.2 Linearized problem and state-space representation

The governing equations of the linearized boundary problem ((20-21) with $\eta = 0$) can be reformulated as

$$ \partial_t^2 w + 2[aI + b\mathcal{B}](\partial_t w) + \mathcal{B}(w) = f, \quad (24) $$

for all $(z,t) \in [0,1] \times \mathbb{T}$, with zero initial conditions (22), where $\mathcal{B}$ is denotes the bi-Laplacian $(\mathcal{B}(w) = w^{(4)})$, operating on the domain<fnref n="4" /> of sufficiently regular functions that satisfy the boundary conditions (21).

This problem admits the state-space representation

$$ \dot{x} = Ax + Bu, \text{ for all } t \in \mathbb{T}, \text{ with } x(0) = 0, \quad (25) $$

<footnote n="3">³ few oscillatory signals, sweeps or sequences of impulses...</footnote>

<footnote n="4">⁴ see appendix A.3 for the detailed construction of the functional setting.</footnote>
where the state $x$, input $u$, linear operators $A$ and $B$ are

$$x(t) = \begin{bmatrix} w(\cdot, t) \\ \partial_t w(\cdot, t) \end{bmatrix}, \qquad u(t) = f(\cdot, t), \qquad (26)$$

$$Ax = \begin{bmatrix} 0 & I \\ -B & -2(aI + bB) \end{bmatrix} x, \quad Bu = \begin{bmatrix} 0 \\ u \end{bmatrix}. \qquad (27)$$

This formulation is well-posed in the spaces of bounded signals (3) for $u(t)$ and $x(t)$ living in, respectively,

$$U = \mathbb{H} \text{ (in accordance with H4)}, \qquad (28)$$

$$X = \mathbb{H}^{\frac{1}{2}} \times \mathbb{H} \text{ with } \|x\|_{X} = \left( \|x_1\|_{\mathbb{H}^{\frac{1}{2}}}^2 + \|x_2\|_{\mathbb{H}}^2 \right)^{\frac{1}{2}}, \qquad (29)$$

and norm $\|x_1\|_{\mathbb{H}^{\frac{1}{2}}} = \|x_1^{(2)}\|_{\mathbb{H}}$, where the complete functional setting is detailed in appendix A.3. In this setting,

A belongs to $\mathcal{L}(X)$ and generates a strongly continuous semigroup $S$ with negative growth bound,

B belongs to $\mathcal{L}(U, X)$ and $\|B\| = 1$,

proving that system (25-27) is in the class defined in section 2.1.

### 3.3 Nonlinear problem and Volterra series convergence

The nonlinear problem ($\eta > 0$) described by (20-21) can be recast in the state-space formulation (1-2) based on definitions (26-27) and by complementing (25) as follows

$$\partial_t x = Ax + A_3(x, x, x) + Bu, \qquad (30)$$

with multi-linear operator

$$A_3(a, b, c) = \begin{bmatrix} 0 \\ -\eta \left( \int_0^1 a'_1(z) b'_1(z) \, dz \right) c''_1 \end{bmatrix}. \qquad (31)$$

This operator defined from $X^3$ to $X$ belongs to $\mathcal{ML}_3(X)$. Indeed, its norm is bounded (proof in appendix B.1) as

$$\|A_3\|_{\mathcal{ML}_3} \le a_3 := \frac{\eta}{3\sqrt{10}}. \qquad (32)$$

Hence, the state-space representation (30) belongs to the class of the well-posed nonlinear problems introduced in section 2.1.

According to section 2.3, there exists a positive $\rho$ such that the Volterra series expansion (4-9) is convergent in norm if $\|x_1\|_{X} < \rho$. Such a convergence radius $\rho$ can be derived using steps 1-4. Computations in appendix B.2 provide a bound estimate $\zeta_3$ related to $\gamma$ (defined in (18)) and to $a_3$ (defined in (32)), and $F$ is given by

$$F(z) = \frac{1}{1 - \zeta_3 z^2} \quad \text{with} \quad \zeta_3 := \gamma a_3. \qquad (33)$$

Solving $F(\sigma) - \sigma F'(\sigma) = 0$ leads to $\sigma = (3\zeta_3)^{-\frac{1}{2}}$. Finally, the convergence radius $\rho = \sigma/F(\sigma)$ is given by

$$\rho = \frac{2}{3\sqrt{3\zeta_3}}, \qquad (34)$$

and $\rho_u = \rho/\gamma$ in remark 3. Expressions of $\gamma$, $\zeta_3$ and $\rho$ (eq. 18, 33, 34) w.r.t. physical coefficients are provided in section 4.2.

Moreover (steps 5-6), solving $\Phi(z) = zF(\Phi(z))$, using the Cardano's method and isolating the positive solution on $[0, \rho]$ which is zero at 0 yields, for all $z \in [0, \rho]$,

$$\Phi(z) = \rho \Psi\left(\frac{z}{\rho}\right) \quad \text{with} \quad \Psi(Z) = 3 \cos\left(\frac{\pi + \arccos Z}{3}\right). \qquad (35)$$

Its first Taylor coefficients are given by $\phi_{2p} = 0$ for all $p \ge 1$ and $\phi_1 = 1$, $\phi_3 = \frac{2^2}{3^3\rho^2}$, $\phi_5 = \frac{2^4}{3^5\rho^4}$, $\phi_7 = \frac{2^8}{3^8\rho^6}$, $\phi_9 = \frac{2^8 \times 5 \times 11}{3^{12}\rho^8}$ (etc.), from which

$$R_M \phi(z) = \Phi(z) - \sum_{m=1}^{M} \phi_m z^m, \qquad (36)$$

is deduced. The gain bound function $\Phi$ and its remainder $R_M \Phi$ are displayed in figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Normalized gain bound function $\Phi/\rho$ and normalized remainder functions $R_M \Phi/\rho$.</figcaption></figure>

The truncated Volterra expansion (8-9) that solves the nonlinear beam model (30-31) provides a finite cascade of linear interconnected systems, with a controlled approximation error. However, these linear systems are infinite dimensional PDEs, so that their real-time simulation may still be an issue. This is why in next section we introduce a standard dimension reduction technique, the modal decomposition, and investigate how it can be combined with the Volterra expansion to yield a controlled approximation of the beam model by a finite cascade of finite dimension linear systems.
# 4 Modal decomposition and projection

## 4.1 Spectral decomposition of A and projection spaces

Operator A is a spectral operator (see appendix A.3). Its point spectrum is composed of the roots $\lambda_n^{\pm}$ of the characteristic polynomials

$$P_n(\lambda) = \lambda^2 + 2(a + bk_n^4)\lambda + k_n^4 \quad \text{with } k_n = n\pi, \quad (37)$$

for all $n \ge 1$. In the following, the fluid and structural damping coefficients $a > 0$ and $b > 0$ are supposed to be such that

(H5) the first pair of roots ($n = 1$) corresponds to a (second order) damped oscillator, which is true if $\xi_1 := (a+b\pi^4)/\pi^2$ (damping ratio) is smaller than 1,

(H6) higher modal numbers $n \ge 2$ yield faster decays$^5$ than $n = 1$, which is true if $2b(a + b\pi^4) \le 1$. Note that when $b > 0$ the corresponding second order systems can be oscillating at low $n$ (if $\xi_n := (a + bk_n^4)/k_n^2$ is smaller than 1) and eventually become overdamped. A lower bound for the damping ratios is given by $2\sqrt{b/a}$, almost reached when $k_n^2 \sim \sqrt{a/b}$.

Assuming that $\lambda_n^{\pm}$ are always distinct, the associated eigenfunctions of A are

$$e_n^{\pm}(z) = \begin{bmatrix} 1/\lambda_n^{\pm} \\ 1 \end{bmatrix} s_n(z) \quad \text{with} \quad s_n(z) = \sqrt{2} \sin(k_n z) \\ \text{and} \quad k_n = n\pi \quad (n \ge 1). \quad (38)$$

Functions $s_n$ form an orthonormal basis of $\mathbb{H}$ and functions $e_n^{\pm}$ form a basis of $\mathbb{X}$. For $N > 0$ (modal truncation order), we introduce the approximation subspaces of $\mathbb{H} = \mathbb{U}$ (N-dimensional), $\mathbb{H}^{1/2}$ (N-dimensional) and $\mathbb{X}$ (2N-dimensional) generated by the first related basis functions, as follows

$$\hat{\mathbb{H}} = \hat{\mathbb{U}} := \operatorname{span}(s_1, \dots, s_N) \quad \subset \mathbb{H} = \mathbb{U}, \quad (39)$$

$$\hat{\mathbb{H}}^{1/2} := \operatorname{span}(s_1, \dots, s_N) \quad \subset \mathbb{H}^{1/2}, \quad (40)$$

$$\hat{\mathbb{X}} := \operatorname{span}(e_1^+, \dots, e_N^+, e_1^-, \dots, e_N^-) \subset \mathbb{X}, \quad (41)$$

equipped with norms $\|\cdot\|_{\mathbb{H}}$, $\|\cdot\|_{\mathbb{H}^{1/2}}$ and $\|\cdot\|_{\mathbb{X}}$, respectively. Using previous notations, we introduce the orthogonal basis $\{S_1, \dots, S_{2N}\}$ of $\hat{\mathbb{X}}$ defined by

$$S_i = \begin{pmatrix} s_i \\ 0 \end{pmatrix} \quad \text{if } 1 \le i \le N,$$

$$\text{and } S_i = \begin{pmatrix} 0 \\ s_{i-N} \end{pmatrix} \quad \text{if } N + 1 \le i \le 2N. \quad (42)$$

This basis is the one used in the sequel. The associated functional spaces $\hat{\mathcal{U}}$ and $\hat{\mathcal{X}}$ are summarized in table 1 (column 2).

## 4.2 Convergence bound estimate

A bound $\gamma$ satisfying (18) can be derived as (see Appendix B.2 for details)

$$\gamma = \int_{\mathbb{T}} \sup_{n \ge 1} \sqrt{|k_n^4 h_n(t)^2 + h_n'(t)^2|} dt, \quad (43)$$

where the impulse responses $h_n$ of second order linear systems are given according to the damping ratio $\xi_n$ by

$$h_n(t) = \begin{cases} e^{-k_n^2 \xi_n t} \frac{\sin(k_n^2 \sqrt{1-\xi_n^2} t)}{k_n^2 \sqrt{1-\xi_n^2}} & \text{if } \xi_n < 1, \\ te^{-k_n^2 t} & \text{if } \xi_n = 1, \\ e^{-k_n^2 \xi_n t} \frac{\sinh(k_n^2 \sqrt{\xi_n^2-1} t)}{k_n^2 \sqrt{\xi_n^2-1}} & \text{if } \xi_n > 1. \end{cases}$$

Moreover, equations (33-34) become

$$\zeta_3 = \frac{\eta}{3\sqrt{10}} \gamma \quad \text{and} \quad \rho = \frac{2\sqrt{10}}{3} \sqrt{\frac{1}{\eta\gamma}}. \quad (44)$$

The convergence radius denominator is proportional to the square root of the nonlinear stiffness $\eta$ and to that of $\gamma$ which depends on the damping physical parameters $a$ and $b$ in a complex way. It follows that decreasing the nonlinear coefficient helps the convergence. It can be shown that increasing $a$ and $b$ while keeping the lower bound $2\sqrt{b/a}$ constant or increasing will also help the convergence, and that no damping ($\gamma \to +\infty$) prevents convergence$^6$. Note also that, from remark 3, a possibly conservative bound on the system input is $\rho_u = 2\sqrt{10}/(3(\gamma)^{3/2}\sqrt{\eta})$.

## 4.3 Projection of the linearized problem

We denote $\Pi_{\hat{\mathbb{X}}} : \mathbb{X} \to \hat{\mathbb{X}}$ the orthogonal projection on $\hat{\mathbb{X}}$. By construction of the modal basis, if $x \in \mathcal{X}$ is the solution of (25) with input $u \in \mathcal{U}$ and initial condition $x_{\text{ini}} \in \mathbb{X}$, then $\hat{x} = \Pi_{\mathbb{X}}(x) \in \hat{\mathcal{X}}$ is the solution of the same problem with input $\hat{u} = \Pi_{\mathbb{H}}(u) \in \hat{\mathcal{U}}$ and initial condition $\hat{x}_{\text{ini}} = \Pi_{\mathbb{X}}(x_{\text{ini}}) \in \hat{\mathbb{X}}$. This corresponds to a finite dimensional problem in which operators $A$ and $B$ can be replaced by their restrictions $\hat{A} : x \in \hat{\mathbb{X}} \mapsto Ax \in \hat{\mathbb{X}}$ and $\hat{B} : u \in \hat{\mathbb{U}} \mapsto Bu \in \hat{\mathbb{X}}$, respectively. This problem can be restated for the coordinates in the basis $\{S_1, \dots, S_{2N}\}$, by introducing input $v$, initial condition

<footnote>$^5$ characterized by the smallest absolute value of the real part of $\lambda_n^{\pm}$.</footnote>

<footnote>$^6$ This case involves some so-called secular modes due to internal and undamped resonances.</footnote>
$y_{ini}$ and state $y$

$$v = [V_1, \dots, V_N]^T \in \mathcal{V}, \text{ s.t. } \hat{u}(t) = \sum_{i=1}^{N} V_i(t)s_i, \qquad (45)$$

$$y_{ini} = [Y_1^{(0)}, \dots, Y_{2N}^{(0)}]^T \in \mathbb{Y}, \text{ s.t. } \hat{x}_{ini} = \sum_{i=1}^{2N} Y_i^{(0)} S_i, \qquad (46)$$

$$y = [Y_1, \dots, Y_{2N}]^T \in \mathcal{Y}, \text{ s.t. } \hat{x} = \sum_{i=1}^{2N} Y_i S_i, \qquad (47)$$

where spaces $\mathcal{V}$, $\mathbb{Y}$ and $\mathcal{Y}$ are defined in table 1 (column 3) and $\mathbb{Y}$ is equipped with the norm built from $\mathbb{X}$

$$\|y\|_{\mathbb{Y}} = \left\| \sum_{i=1}^{2N} Y_i S_i \right\|_{\mathbb{X}} = \sqrt{y^T Q y}, \qquad (48)$$

$$\begin{aligned} & \text{with } Q = \operatorname{diag}(\|S_1\|_{\mathbb{X}}^2, \dots, \|S_{2N}\|_{\mathbb{X}}^2) \\ & = \begin{bmatrix} \pi^2 \operatorname{diag}(1, \dots, N)^2 & 0_{N \times N} \\ 0_{N \times N} & I_N \end{bmatrix}. \end{aligned} \qquad (49)$$

The coordinates $y \in \mathcal{Y}$ of the trajectories are governed by

$$\dot{y}(t) = A_{\mathbb{Y}}y(t) + B_{\mathbb{Y}}v(t), \quad y(0) = y_{ini} \in \mathbb{Y}, \qquad (50)$$

with

$$A_{\mathbb{Y}} = \begin{bmatrix} 0_{N \times N} & I_N \\ -\pi^4 \operatorname{diag}(1, \dots, N)^4 & -2(a I_N + b \pi^4 \operatorname{diag}(1, \dots, N)^4) \end{bmatrix}, \qquad (51)$$

$$B_{\mathbb{Y}} = \begin{bmatrix} 0_{N \times N} \\ I_N \end{bmatrix}, \qquad (52)$$

Note that by definition of the modal basis, when $N \to +\infty$, the truncated system trajectory given by (47) converges in norm towards the solution of (25) in $\mathcal{X}$.

## 4.4 Projection of the nonlinear terms

Considering operator $A_3$ defined in (31), we notice that in the modal subspace $\hat{\mathbb{X}}^3$ we have

$$\text{for all } (a, b, c) \text{ in } \hat{\mathbb{X}}^3, \quad A_3(a, b, c) \in \hat{\mathbb{X}}. \qquad (53)$$

This means that the restriction of $A_3$ on $\hat{\mathbb{X}}$ defines a multilinear operator on this space, that we denote $\hat{A}_3$.

Consider the trajectory $x$ of the nonlinear beam model (30) with initial condition $x_{ini}$ and input $u$ respectively in $\mathbb{X}$ and $\mathbb{U}$. Then, we define the pseudo-modal truncation of $x$ on $\hat{\mathbb{X}}$ as the solution $\tilde{x}$ of problem

$$\dot{\tilde{x}} = \hat{A}\tilde{x} + \hat{B}\Pi_{\mathbb{H}}u + \hat{A}_3(\tilde{x}, \tilde{x}, \tilde{x}), \qquad (54)$$

$$\hat{x}_{ini} = \Pi_{\mathbb{X}} x_{ini} \in \hat{\mathbb{X}}.$$

Note that $\tilde{x}$ does not identify with $\hat{x} := \Pi_{\mathbb{X}}x$ for a solution $x$ in general.

For sake of conciseness, in the sequel $\tilde{x}$ will be referred to as the pseudo-modal truncation of $x$ omitting the input $u$ and initial condition $x_{ini}$. As in section 4.3, the finite-dimensional problem (54) satisfied by $\tilde{x}$ on $\hat{\mathbb{X}}$ can be restated on $\mathbb{Y}$ for the coordinates in the basis $\mathcal{S}$ as

$$\dot{y}(t) = A_{\mathbb{Y}}y(t) + B_{\mathbb{Y}}v(t) + A_{\mathbb{Y}_3}(y, y, y), \quad y(0) = y_{ini} \in \mathbb{Y}, \qquad (55)$$

where $A_{\mathbb{Y}_3} \in \mathcal{ML}_3(\mathbb{Y})$ is given by

$$A_{\mathbb{Y}_3}(a, b, c) = -\eta \left( a^T \pi^2 \operatorname{diag}(1, \dots, N, 0_{1 \times N})^2 b \right) \times \begin{bmatrix} 0_{N \times N} & 0_{N \times N} \\ -\pi^2 \operatorname{diag}(1, \dots, N)^2 & 0_{N \times N} \end{bmatrix} c. \qquad (56)$$

Remark 4 (Duffing oscillators) If the $N$ first damping ratios $\xi_n$ are less than 1, then the pseudo-modal system (55-56) exactly corresponds to $N$ coupled damped Duffing oscillators.

# 5 Reduced beam model combining Volterra expansion and pseudo-modal truncation

This section presents the approximation of the beam model by the truncated Volterra series of its pseudo-modal approximation, with modal truncation order $N \ge 1$ (defined in (54)) and Volterra series truncation order $M \ge 1$. Truncation error estimates are given, that take into account both the modal and the Volterra truncations. Moreover, hints to generalize the approach to other nonlinear models are also provided.

## 5.1 Pseudo-modal Volterra approximation

For a solution of the beam model $x$, with initial condition $x_{ini}$ and input $u$, we consider for order $N \ge 2$ its pseudo-modal truncation $\tilde{x}$ governed by system (54), and introduce the first $M$ terms of the Volterra expansion of $\tilde{x}$ as

$$\tilde{x}_1 = \hat{x}_1 = \Pi_{\mathbb{X}}(x_1), \qquad (57)$$

$$\tilde{x}_m = \int_0^t e^{\hat{A}(t-\tau)} \sum_{p \in \mathbb{M}_m^3} \hat{A}_3(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau)) d\tau, \qquad \text{for all } m \text{ such that } 2 \le m \le M. \qquad (58)$$

Then, we define the pseudo-modal Volterra approximation of $x$ as

$$\tilde{X}_N^M = \sum_{m=1}^{M} \tilde{x}_m. \qquad (59)$$
<table><caption>Table 1: Notations of the spaces involved in: (1) the original state-space problem, (2) the projected problem, (3) the representation of (2) with coordinates</caption><thead><tr><th>Original spaces</th><th>Projection sub-spaces</th><th>Modal coordinates spaces</th></tr></thead><tbody><tr><td>U = H</td><td>Û = Ĥ</td><td>V = R<sup>N</sup> (equipped with 2-norm)</td></tr><tr><td>X = H × H<sup>1/2</sup></td><td>X̂ = Ĥ × Ĥ<sup>1/2</sup> (equipped with (29))</td><td>Y = R<sup>2N</sup> (equipped with ||·||<sub>Y</sub>, see (48))</td></tr><tr><td>U = L<sup>∞</sup>(T, U)</td><td>Û = L<sup>∞</sup>(T, Û)</td><td>V = L<sup>∞</sup>(T, V)</td></tr><tr><td>X = L<sup>∞</sup>(T, X)</td><td>X̂ = L<sup>∞</sup>(T, X̂)</td><td>Y = L<sup>∞</sup>(T, Y)</td></tr></tbody></table>

We first examine the convergence of this series, and second, the error between $x$ and the approximation $\tilde{X}_N^M$.

## 5.2 Convergence

We apply the same steps as in sections 3.3 and 4.2 to (54) and the Volterra series decomposition (57-58). This yields the estimates

$$ \gamma_N = \int_{\mathbb{T}} \sup_{n=1}^{N} |\sqrt{k_n^4 h_n(t)^2 + h'_n(t)^2}| dt, \quad (60) $$

$$ \text{and } \|\hat{A}_3\|_{\mathcal{ML}_3(\hat{X})} = \|A_{\mathbb{Y}3}\|_{\mathcal{ML}_3(\mathbb{Y})} \le a_{3N} \text{ with} $$

$$ a_{3N} = \frac{\eta}{\pi^2} \left( \sum_{n=1}^{N} \frac{1}{n^4} \right)^{1/2} . \quad (61) $$

Then, we build

$$ \zeta_{3N} = \gamma_N a_{3N}, \quad \rho_N = \frac{2}{3\sqrt{3\zeta_{3N}}}, \quad (62) $$

and the gain bound function (see (35))

$$ \forall z \in [0, \rho_N), \quad \hat{\Phi}_N(z) = \rho_N \Psi\left(\frac{z}{\rho_N}\right) = \frac{\rho_N}{\rho} \Phi\left(\frac{\rho}{\rho_N} z\right). \quad (63) $$

Note that $\gamma_N$, $a_{3N}$ and $\zeta_{3N}$ are all strictly lower than, respectively, $\gamma$, $a_3$ and $\zeta_3$ and tend towards them as $N \to \infty$. Then, the modal truncation makes the convergence radius $\rho_N$ increase according to the dilatation factor

$$ r_N := \frac{\rho_N}{\rho} = \sqrt{\frac{\zeta_3}{\zeta_{3N}}} = \sqrt{\frac{a_3}{a_{3N}}} \sqrt{\frac{\gamma}{\gamma_N}} > 1, \quad (64) $$

through which the gain bound functions are related as $\hat{\Phi}_N(z) = r_N \Phi(z/r_N)$. This factor combines a first dilatation factor $\sqrt{a_3/a_{3N}} > 1$ depending on $N$ only (see values in table 2) and a second one $\sqrt{\gamma/\gamma_N} > 1$ depending on $a$ and $b$ in a complex way but that can be even more effective in practice.

Note also that, obviously, the Volterra series expansion (4-9) of the pseudo-modal coordinate system has the same convergence radius $\rho_N$ and gain bound function $\hat{\Phi}_N$, so that $y = \sum_{m=1}^{\infty} y_m$ is convergent in norm if $\|y_1\|_Y = \| < \rho_N$. Moreover, we have $\|y\|_Y = \|\tilde{x}\|_{\hat{X}} < \hat{\Phi}_N(\|\tilde{x}_1\|_{\hat{X}}) = \hat{\Phi}_N(\|y_1\|_Y)$.

## 5.3 Truncation error estimates

Based on these results, we can now address the derivation of an estimate of the error bound of the pseudo-modal Volterra approximation. Consider now an input $u = f \in U$ and an initial condition $x_{\text{ini}} \in X$ such that $\|x_1\|_X < \rho$. Assume that $u$ and the linear response $x_1$ can be accurately described by their $N$-order modal decomposition $\hat{u}$ and $\hat{x}_1$ in the sense that

$$ \begin{align*} \text{(A1)} \quad & u = \hat{u} + u^{\perp} \text{ where } \hat{u} = \Pi_{\mathcal{U}} u \text{ and } u^{\perp} = (\text{Id} - \Pi_{\mathcal{U}})u, \\ \text{(A2)} \quad & \|x_1^{\perp}\|_{\mathcal{X}} \le \epsilon \|\hat{x}_1\|_{\mathcal{X}} \text{ with } \epsilon \ll 1, \\ \text{(A3)} \quad & \|\hat{x}_1\|_{\mathcal{X}} + \|x_1^{\perp}\|_{\mathcal{X}} \le (1 + \epsilon) \|\hat{x}_1\|_{\mathcal{X}} < \rho, \end{align*} $$

where $\hat{x}_1$ and $x_1^{\perp}$ denote the responses of the linearized system excited by $\hat{u}$ and $u^{\perp}$, respectively.

As $x_1$ and $\hat{x}_1$ are in the convergence domain and after (64), the series expansions $x = \sum_{m=1}^{\infty} x_m$ for the system excited by $u$ and initial condition $x_{\text{ini}}$, and $\tilde{x} = \sum_{m=1}^{\infty} \tilde{x}_m$ for the pseudo-modal truncation excited by $\hat{u}$ and initial condition $\hat{x}_{\text{ini}}$ are both convergent in norm. The error on the approximated trajectory $\tilde{x}$ is bounded as stated in the following proposition.

**Proposition 1 (pseudo-modal truncation error bound)** Assume (A1-A3). Denote $e = x - \tilde{x}$ the error on $x$ due to the pseudo-modal truncation on $\tilde{X}$. Then, $e$ is bounded as

$$ \|e\|_{\mathcal{X}} \le D_{\epsilon} \Phi(\|\hat{x}_1\|_{\mathcal{X}}) \quad (65) $$

$$ \text{where } D_{\epsilon}\Phi(z) := \Phi((1 + \epsilon)z) - \Phi(z). $$

The proof is given in C. Bound estimates are given in figure 2 for several values of $\epsilon$.

Moreover, as a consequence of (64), we have the following corollary.

**Corollary 1 (Pseudo-modal Volterra approximation error)** With the above hypotheses, the error due to the pseudo-modal Volterra approximation $\tilde{X}_N^M$ is bounded by

$$ \begin{align} \|x - \tilde{X}_N^M\|_{\mathcal{X}} & \le \|x - \tilde{x}\|_{\mathcal{X}} + \|\tilde{x} - \tilde{X}_N^M\|_{\mathcal{X}} \nonumber \\ & \le D_{\epsilon} \Phi(\|\hat{x}_1\|_{\mathcal{X}}) + R_M \hat{\Phi}_N(\|\hat{x}_1\|_{\mathcal{X}}). \tag{66} \end{align} $$
<table><caption>Table 2: Coefficient $\sqrt{a_3/a_{3N}} - 1$ for $1 \le N \le 8$ (lower rounded values in per cent).</caption><thead><tr><th>N</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>10<sup>2</sup>(√a<sub>3</sub>/a<sub>3N</sub> - 1)</td><td>1.99</td><td>0.463</td><td>0.173</td><td>0.0826</td><td>0.0455</td><td>0.0277</td><td>0.0181</td><td>0.0124</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Function $D_\epsilon\Phi/\rho (=D_\epsilon\Psi)$ for relative errors $\epsilon$ of 0.1%, 1% and 10% (see (35) and (65)).</figcaption></figure>

## 5.4 Generalization of pseudo-modal Volterra approximation

The pseudo-modal Volterra approximation (57-58) can be generalized to models for which $A_3$ does not fulfill (53), or to more general models (1) such that $A_k(\hat{\mathbf{X}}^k) \subseteq \hat{\mathbf{X}}$ does not hold.

Indeed, we considered so far a beam model for which the definition of the pseudo-modal decomposition (54) and the definition of the pseudo-modal Volterra approximation (57-58) crucially rely on (53), which is a particular feature of the model: operator $A_3$ (the nonlinear part of the model) defines a multilinear operator in the modal subspace $\hat{\mathbf{X}}$.

We now consider a beam model with the same linear part as above, but with a multilinear operator $A_3$ such that the order $N$ modal subspace $\hat{\mathbf{X}}$ is not stable, that is, $A_3(\hat{\mathbf{X}})$ is not a subset of $\hat{\mathbf{X}}$ anymore. We assume that a bound $a_3$ of $\|A_3\|$ is available so that the new model is in the class of well-posed problems defined in section 2.1. It follows that we obtain a convergence bound $\rho$ and a gain bound function $\Phi$ given by equations (34-35) computed with the new value of $a_3$.

For a solution of this new beam model $x$, with initial condition $x_{\text{ini}}$ and input $u$, we define its pseudo-modal truncation at order $N$ $\tilde{x}_m$ as in equation (54), where $\hat{A}_3$ is now defined as $\hat{A}_3 = \Pi_{\mathbf{X}}A_3$, the projection of $A_3$ on the modal subspace $\hat{\mathbf{X}}$. Note that thanks to this definition, $\tilde{x}$ belongs to $\hat{\mathbf{X}}$.

The pseudo-modal Volterra approximation of $x$ is now $\tilde{X}_N^M = \sum_{m=1}^M \tilde{x}_m$, the sum of the first $M$ terms of

the Volterra expansion of $\tilde{x}$ defined as in (57-58). As in previous section, the Volterra series expansion of $\tilde{x}$ on $\hat{\mathbf{X}}$ has a convergence radius $\hat{\rho}$ greater than $\rho$ since, as $\Pi_{\mathbf{X}}$ is an orthogonal projection, $\|\hat{A}_3\| \le \|A_3\| \le a_3$, and a gain bound function $\hat{\Phi}$. Then, defining

$$ \|A_3 - \hat{A}_3\|_{\hat{\mathbf{X}}} = \sup_{\substack{(a,b,c) \in (\hat{\mathbf{X}})^3 \\ \|a\|_{\hat{\mathbf{X}}}, \|b\|_{\hat{\mathbf{X}}}, \|c\|_{\hat{\mathbf{X}}}=1}} \|A_3(a,b,c) - \hat{A}_3(a,b,c)\|_{\mathbf{X}}, \quad (67) $$

modified versions of proposition 1 and corollary 1 are stated below, whose proof is given in appendix D.

**Proposition 2 (pseudo-modal truncation error bound)** Assume (A1-A3). Denote $e = x - \tilde{x}$ the error on $x$ due to the pseudo-modal truncation on $\hat{\mathbf{X}}$. Then, $e$ is bounded as

$$ \|e\|_{X} \le D_{\epsilon}\Phi\left(\|\hat{x}_{1}\|_{X}\right) + \frac{\|A_{3} - \hat{A}_{3}\|_{\hat{\mathbf{X}}}}{a_{3}}\Phi\left(\|\hat{x}_{1}\|_{X}\right), \quad (68) $$

where $D_\epsilon\Phi(z)$ is defined in proposition 1.

**Corollary 2 (Pseudo-modal Volterra approximation error)** With the above hypotheses, the error due to the pseudo-modal Volterra truncation $\tilde{X}_N^M$ is bounded by

$$ \|x - \tilde{X}_N^M\|_{X} \le D_{\epsilon}\Phi\left(\|\hat{x}_1\|_X\right) + R_M\hat{\Phi}_N\left(\|\hat{x}_1\|_X\right) + \frac{\|A_3 - \hat{A}_3\|_{\hat{\mathbf{X}}}}{a_3}\Phi\left(\|\hat{x}_1\|_X\right). $$

Note that for the beam model under consideration in the rest of the paper, because of (53), (67) is zero and proposition 2 and corollary 2 boil down to proposition 1 and corollary 1. Note also that (67) is easy to compute in the finite dimension space $\hat{\mathbf{X}}$ and that its value tends to zero when the order $N$ of the modal truncation goes to infinity.

Finally, these results can be further extended to any system in the class of well-posed problems (1-2) defined in Hilbert spaces: if one can find a subspace $\hat{\mathbf{X}}$ stable for the linear part of the system, then proposition 2 and corollary 2 hold replacing $\frac{\|A_3 - \hat{A}_3\|_{\hat{\mathbf{X}}}}{a_3}$ by

$$ \max_{k \in \{2, \dots, K\}} \frac{\|A_k - \Pi_{\mathbf{X}} A_k\|_{\hat{\mathbf{X}}}}{a_k}. $$

## 6 Refinements on convergence bounds

This section presents new results to improve convergence and error bound estimates. We refine the theoretical results of section 2 for specific classes of nonlinearities (§ 6.1) and apply them to the beam (§ 6.2).
## 6.1 New theoretical results for a class of nonlinearities

The convergence bounds in section 2.3 can be improved if the multilinear operators $A_k: X^k \to X$ only act on subspaces of $X$ and admits a sandwich decomposition

$$A_k(x_1, \dots, x_k) = B_k A_k^\dagger(Cx_1, \dots, Cx_k), \quad (69)$$

where $C: X \to X^\dagger$, $A_k^\dagger: (X^\dagger)^k \to W_k^\dagger$ and $B_k: W_k^\dagger \to X$ lie respectively in $\mathcal{L}(X, X^\dagger)$, $\mathcal{ML}_k(X^\dagger, W_k^\dagger)$ and $\mathcal{L}(W_k^\dagger, X)$. Typically, $X^\dagger$ and $W_k^\dagger$ refer to spaces $(\mathbb{R}^{n^\dagger} \le n$ or Banach) smaller than $X$, and $A_k^\dagger$ to a reduced version of $A_k$. For clarity, we stamp all the labels of the smaller or reduced objects with the dagger symbol $^\dagger$.

The idea is to examine the convergence of (4) on $x$ through that on

$$x^\dagger = Cx, \text{ defined in } \mathcal{X}^\dagger = L^\infty(\mathbb{T}, X^\dagger). \quad (70)$$

Indeed, its Volterra series is the sum of $x_m^\dagger = C x_m$, straightforwardly derived from (8-9) and (69) as

$$x_1^\dagger(t) = CS(t)x_{\text{ini}} + \int_0^t S_1^\dagger(t-\tau)u(\tau) \, d\tau, \quad (71)$$

$$x_m^\dagger(t) = \sum_{k=2}^{\min(m,K)} \sum_{p \in \mathbb{M}_m^k} \int_0^t S_k^\dagger(t-\tau) A_k^\dagger(x_{p_1}^\dagger(\tau), \dots, x_{p_k}^\dagger(\tau)) \, d\tau, \\ \text{for } m \ge 2, \quad (72)$$

where we have defined, for $1 \le k \le K$ and $B_1 = B$,

$$S_k^\dagger(t) = C S(t) B_k. \quad (73)$$

Then, we can compute new convergence estimates that benefit from the focused action of operators $S_k^\dagger$, $A_k^\dagger$: compared to $S$, $A_k$, these operators are circumscribed to smaller spaces and expected to yield lower bound estimates than $\zeta_k$ in (12). Thus, bounds $\zeta_k$ are replaced in step 1 by the new estimates, for $2 \le k \le K$,

$$\zeta_k^\dagger \ge \int_{\mathbb{T}} \|S_k^\dagger(t) A_k^\dagger\|_{\mathcal{ML}_k(X^\dagger, X^\dagger)} dt, \quad (74)$$

some overestimated values of which are $\gamma_k^\dagger a_k^\dagger$ with

$$a_k^\dagger \ge \|A_k^\dagger\|_{\mathcal{ML}_k(X^\dagger, W_k^\dagger)}, \quad \gamma_k^\dagger \ge \int_{\mathbb{T}} \|S_k^\dagger(t)\|_{\mathcal{L}(W_k^\dagger, X^\dagger)} dt. \quad (75)$$

Using $\zeta_k^\dagger$ instead of $\zeta_k$ in steps 2 to 6, we successively define $F^\dagger$, $\zeta^\dagger$, $\sigma^\dagger$, $\rho^\dagger$ and $\Phi^\dagger$ and are ready to state the following theorem.

**Theorem 1 (Refined Volterra convergence and error bounds)** Assume that $\|x_1^\dagger\|_{X^\dagger} < \rho^\dagger$. Then, the Volterra series of $x^\dagger$ and $x$ converge in norm in $X^\dagger$

and $X$, respectively. Moreover, the following inequalities hold:

$$\|x^\dagger\|_{X^\dagger} < \Phi^\dagger(\|x_1^\dagger\|_{X^\dagger}) \quad (76)$$

$$\|x^\dagger - \sum_{m=1}^{M} x_m^\dagger\|_{X^\dagger} < R_M \Phi^\dagger(\|x_1^\dagger\|_{X^\dagger}), \text{ for } M \ge 1, \quad (77)$$

$$\|x - \sum_{m=1}^{M} x_m\|_{X} < r R_M \Phi^{\dagger}(\|Cx_1\|_{X^{\dagger}}), \text{ for } M \ge 1, (78)$$

$$\|x\|_{X} < \|x_1\|_{X} + r R_1 \Phi^{\dagger}(\|Cx_1\|_{X^{\dagger}}), \quad (79)$$

where $r = \max_{2 \le k \le K}(\zeta_k / \zeta_k^\dagger)$.

The proof follows exactly the same steps as in [15].

**Remark 5 (Supplement to remark 3)** An alternative convergence radius on the input is $\rho_u^\dagger = \rho^\dagger / \gamma_1^\dagger$ with $\gamma_1^\dagger = \int_{\mathbb{T}} \|S_1^\dagger\|_{\mathcal{L}(U, X^\dagger)} dt$. Moreover, this radius can be adapted and improved for specific inputs. As a simple example, consider a scalar ODE of the form $L(d/dt)\xi + \sum_{k=2}^K \alpha_k \xi^k = u$, for which we can choose $x^\dagger = \xi$ and $A_k^\dagger(\xi_1, \dots, \xi_k) = \alpha_k \xi_1 \dots \xi_k$. Denote $H$ the Laplace transfer function and $h$ the impulse response of its linear part. If the system is excited by the specific input $u_\omega(t) = U e^{i\omega t}$ (harmonic regime), then $\xi_1(t) = H(i\omega)u_\omega(t)$ leads to a convergence bound for amplitude $U$ that depends on the pulsation, namely, $|U| < \rho^\dagger / |H(i\omega)|$.

## 6.2 Application to the beam model

Operator $A_3$ defined in (31) admits a sandwich decomposition (69). It involves the deflection wave

$$w := x^\dagger = Cx \text{ with } C: \begin{bmatrix} w \\ v \end{bmatrix} \in X \mapsto w \in \mathbb{H}^{\frac{1}{2}}, \quad (80)$$

from which the nonlinear term in (20) is built using

$$A_3^\dagger(w_a, w_b, w_c) = -\eta \left( \int_0^1 w_a'(z) w_b'(z) dz \right) w_c'', \quad (81)$$

and contributes to the state equation through (as input $u$)

$$B_3 = B: \mathbb{H} \to \mathbb{X}. \quad (82)$$

Operator $A_3^\dagger$ belongs to $\mathcal{ML}_3(X^\dagger := \mathbb{H}^{\frac{1}{2}}, W_3^\dagger := \mathbb{H})$ and has the same norm value as $A_3$. Similarly, its restriction $\tilde{A}_3^\dagger$ on $\tilde{\mathbb{H}}^{\frac{1}{2}}$, which generates $\tilde{A}_3 = C \tilde{A}_3^\dagger(B, B, B)$, belongs to $\mathcal{ML}_3(\tilde{\mathbb{H}}^{\frac{1}{2}}, \tilde{\mathbb{H}})$ and has the same norm value as $\tilde{A}_3$.

Then, the nonlinear systems (30) governing the beam state $x$ and (54) governing its pseudo-modal truncation $\tilde{x}$ can both benefit from the refined convergence
results presented in section 6.1. The estimates replacing $\gamma_N \xrightarrow[N \to +\infty]{} \gamma$ are given by (see Appendix B.2)

$$ \gamma_N^\dagger := \int_T \sup_{n=1}^N (k_n^2 |h_n(t)|) \, dt, $$

with limit $\gamma^\dagger$, which are strictly lower than $\gamma_N$ and $\gamma$ respectively, and from which we build $\zeta_{3N}^\dagger = a_{3N}\gamma_N^\dagger$ and $\zeta^\dagger = a_3\gamma^\dagger$. Finally, the improved convergence radii are then

$$ \rho_N^\dagger = \frac{2}{3\sqrt{3a_{3N}\gamma_N^\dagger}} > \rho_N > \rho \quad \text{and} \quad \rho^\dagger = \frac{2}{3\sqrt{3a_3\gamma^\dagger}} > \rho, \qquad (83) $$

and, in theorem 1, factors $r_N$ (for $\tilde{x}$) and $r$ (for $x$) are given by

$$ r_N = \frac{\zeta_{3N}}{\zeta_{3N}^\dagger} = \frac{\gamma_N}{\gamma_N^\dagger} > 1 \quad \text{and} \quad r = \frac{\zeta_3}{\zeta_3^\dagger} = \frac{\gamma}{\gamma_3^\dagger} > 1. $$

We may observe that the convergence criterion on $x_1^\dagger$ is actually weaker than the one on $x_1$. Indeed, if $\|x_1\|_X < \rho$, then $\|x_1^\dagger\|_{X^\dagger} < \|C\|_{L(X,X^\dagger)} \|x_1\|_X = \|x_1\|_X < \rho < \rho^\dagger$, and the convergence radius dilatation is $\rho^\dagger/\rho = \sqrt{r} > 1$.

Moreover, following theorem 1 and assuming (A1-A3) as in proposition 1 for $u, x_1^\dagger, \rho^\dagger$ instead of $u, x_1, \rho$, corollary 1 is improved as follows.

**Corollary 3 (Improved pseudo-modal Volterra approximation error)** The error due to the pseudo-modal Volterra approximation $\tilde{X}_N^M$ is bounded by

$$ \begin{gathered} \|x^\dagger - C\tilde{X}_N^M\|_{X^\dagger} \le E_N^M \quad \text{and} \quad \|x - \tilde{X}_N^M\|_X \le r E_N^M, (84) \\ \text{with } E_N^M = D_\epsilon \Phi^\dagger (\|\hat{x}_1^\dagger\|_{X^\dagger}) + R_M \hat{\Phi}_N^\dagger (\|\hat{x}_1^\dagger\|_{X^\dagger}). \end{gathered} $$

In the specific situation where $N=1$ (the beam is modelled as a single Duffing oscillator), setting $h(t) = k_1^2 h_1(t)$, we obtain $\gamma_1^\dagger = \|h\|_1$ and the convergence bound $\rho_1^\dagger$ is expressed as

$$ \rho_1^\dagger = \frac{2}{3^{\frac{3}{2}} \sqrt{|\alpha_3| \|h\|_1}}. \qquad (85) $$

**Remark 6 (Single Duffing oscillator in harmonic regime)** Following remark 5, the resulting conservative bound on the systems input $\rho_{1u}^\dagger = \rho_1^\dagger/\|h\|_1$ deduced from (85) allows to recover the one given in ([19],[20]), established for a Duffing oscillator, with impulse response $h$. In the harmonic regime, we obtain a bound

$$ \rho_1^\dagger/|H(i\omega)| = 2/(3^{\frac{3}{2}}\sqrt{|\alpha_3|\|h\|_1}|H(i\omega)|). $$

As $|H(i\omega)| < \|h\|_1$, this last bound is close, but slightly less than the bound obtained by these authors, whose

expression is $2/((3|H(i\omega)|)^{3/2}\sqrt{|\alpha_3|})$. This illustrates that for specific inputs and specific systems, the bounds we provide can be improved, however our bound is relevant for inputs with rich frequency content.

# 7 Numerical experiments

## 7.1 Numerical method

The trajectories of the Euler-Bernoulli beam are simulated, based on the truncated Volterra series of the pseudo-modal solution described by (57-58). Numerical experiments are processed inside the convergence domain, where the error estimate (84) is available (see figures 1-2 and (63)), and outside the convergence domain to investigate the behaviour of the approximated solution in this case.

The truncated series expansion is convenient for time-domain simulation as it can be achieved by combining linear systems and instantaneous nonlinearities. This is described by the block-diagram in figure 3 (see remark 2 and figure 1 in [15]) for the more general class (1-2)), in which the finite dimensional linear systems (W in grey) are simulated using a standard linear solver (lslm, Matlab) and the evaluation of instantaneous operators is exact. Moreover, to serve as a reference, the trajectory of the nonlinear problem (54) on $\tilde{X}$ is computed using a standard ODE solver on a refined time grid (here, ode15s, Matlab).

Note that inside the guaranteed convergence domain (constraint characterized in this paper), there can be some benefits to using a Volterra series instead ODE solvers. It yields explicit computation (no need for iterative solver) which, through W in figure 3, reproduces the exact spectral values at low amplitude, guarantees stability (also in the nonlinear regimes) and, through the expansion with respect to homogeneous orders, does not introduce any cumulative error in the contributions due to the (pre-computed) kernels. It also provides a convenient form for application in control. Moreover, when necessary in such applications (or also in telecommunications, audio, etc.), it also allows to reject the aliasing due to nonlinearities by sandwiching the multilinear function $M_k$ between an oversampler (of factor $k$) and a Shannon anti-aliasing filter$^7$. This strategy extracts the correct signal in the frequency baseband, whereas applying it on the ODE (or its field) modifies the baseband behaviour.

<footnote>$^7$ Thanks to filter commutations and using an oversampler (starting in the block-chain at input $v$ or at the linear contribution $y_1$), this can be done by simply combining each W with the anti-aliasing filter in the grey blocks.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3: Input-Output simulation based the Volterra series expansion (57-58) truncated at order $M = 7$. This structure exploits that even order terms are zero and that $A_{Y_3} =: M_3$ is symmetric with respect to the two first variables ($A_{Y_3}(a, b, c) = A_{Y_3}(b, a, c)$). Triplets of numbers pqr above blocks $M_3$ account for the terms $M_3(y_p, y_q, y_r)$. $W$ and $WB$ respectively denote the simulation of the semigroup generated by matrix $A_Y$ and the one with input matrix $B_Y$. The space-time function $\tilde{x}_{\text{trunc}}(z, t) = R(z)y_{\text{trunc}}(t)$ is built by aggregating the activated eigenfunctions as in formula (47).</figcaption></figure>

## 7.2 Parameters and input signal

Two configuration sets of experiments are examined. Configuration 1 corresponds to a damped Duffing oscillator near the critical regime: only the first mode of the beam is excited and we examine the nonlinear system with respect to the convergence radius $\rho_1^\dagger$ (see (83)). Configuration 2 corresponds to the beam in oscillating regime, examined with respect to $\rho^\dagger$ (see also (83)) and an increasing number $N$ of modes, on which the force excitation is decomposed. The finite dimensional system corresponds to coupled damped Duffing oscillators.

For presentation purposes (without loss of generality), systems are all built such that, given the linearized model (here, given the damping coefficients $a$ and $b$), the nonlinear coefficient $\eta$ is chosen to provide a unitary convergence radius ($\rho_1^\dagger$ for configuration 1 and $\rho^\dagger$ for the configuration 2). This yields

$$ \eta = \frac{4\pi^2}{27\gamma_1^\dagger} \quad (\text{config. 1}) \quad \text{and} \quad \eta = \frac{4\sqrt{10}}{9\gamma^\dagger} \quad (\text{config. 2}). \qquad (86) $$

Moreover, in order to excite the beam by a relevant rich spectral content, we choose inputs $u: t \mapsto (z \mapsto f(z, t))$ in (30) with step time-shapes

$$ f(z, t) = f_{\text{tot}} g(z) 1_{[0, T_e]}(t), \qquad (87) $$

where the total force $f_{\text{tot}}$ is spatially distributed according to $g > 0$ (where $\int_0^1 g(z) dz = 1$).

As mentioned in section 7.1, we simulate the pseudomodal Volterra terms $\tilde{x}_m$ governed by (57-58) to build

the approximation (59), through the realisation in figure 3 that uses the coordinate description $y$ in (55-56).

The two configurations are now detailed below, and all the parameters (damping and nonlinear coefficients) as well as the excitation (shape and duration) are summarized in table 3.

**Configuration 1 (non oscillating single mode)** We consider a fluid damping only ($b=0$) tuned to be close to the critical regime ($a = 0.999\pi^2$). The excitation (87) is spatially-distributed on the first mode only ($g(z) = s_1(z) = \sqrt{2}\sin(\pi z)$ defined in (38)) and the step duration is $T_e = 3$. We obtain $\gamma_1^\dagger = 0.1013$ and (86) yields $\eta = 14.434$. Four amplitudes are tested. They are chosen such that $\|\tilde{x}_1^\dagger\|_{X^\dagger} = A$ with $A \in \{0.8; 1; 1.2; 2\}$. They correspond to excitation forces $f_{\text{tot}} = A/\|\tilde{x}_1^{\dagger*}\|_{X^\dagger}$ where $\|\tilde{x}_1^{\dagger*}\|_{X^\dagger}$ is computed for the linear response $\tilde{x}_1^{\dagger*}$ to unit force. This yields $f_{\text{tot}} \approx 9.8696 \times A$. The simulation time-step is chosen as $T = 10^{-3}$.

**Configuration 2** The second configuration corresponds to an oscillating beam, which can be used for sound synthesis of e.g. vibraphones, xylophones, marimbas, etc. Here, the damping parameters $a=0.1$ and $b=10^{-6}$ (for which we obtain $\gamma^\dagger = 8.774986$ and $\eta \approx 0.1602$) are chosen so that it sounds like a “wooden beam”. Note that to listen to a result with a first mode at frequency $f_0$ on a sound card with sampling frequency $f_s$, the simulated trajectories must be sampled at $T = 2\pi f_0 / (f_s \Im s_1)$ where $(s_1, \bar{s}_1)$ are the poles associated with the first mode ($T \approx 2.4536 \times 10^{-3}$ for $f_0 = 185$Hz, $f_s = 48000$Hz).
<table><caption>Table 3: Summary of the main physical parameters of the configurations 1 and 2 used for numerical experiments.</caption><thead><tr><th>Configuration</th><th>1 (non oscillating)</th><th>2 (oscillating+hammer)</th></tr></thead><tbody><tr><td>Parameters: a, b</td><td>a = 0.999π² and b = 0</td><td>a = 0.1 and b = 10<sup>-6</sup></td></tr><tr><td>Impulse response bound</td><td>γ<sub>1</sub> ≈ 0.1013</td><td>γ<sup>†</sup> ≈ 8.7749</td></tr><tr><td>η (see 86)</td><td>η ≈ 14.434</td><td>η ≈ 0.1602</td></tr><tr><td>Input (simulation on [0, 2T<sub>e</sub>]):</td><td>N = 1 and T<sub>e</sub> = 3</td><td>N ∈ {2, 16} and T<sub>e</sub> = 40</td></tr><tr><td>f(z, t) = f<sub>tot</sub> 1<sub>[0, T<sub>e</sub>]</sub>(t) Σ<sub>n=1</sub><sup>N</sup> g<sub>n</sub>e<sub>n</sub>(z)</td><td>g<sub>1</sub> = 1</td><td>g<sub>n</sub>: see (88) and fig. 4a</td></tr></tbody></table>

The reference excitation (87) is chosen with a step duration of $T_e = 3$. It is spatially distributed according to a cosinusoidal activation of width $\delta = 1/8$, centered at $z_c = 1/7$ described by $g(z) = G \cos(\pi(z - z_c)/\delta)$ on $[z_c - \delta/2, z_c - \delta/2]$ and zero outside, with $G = \pi/(2\delta)$. The coefficients of the modal decomposition of $g$ are, for all $n \ge 1$,

$$g_n = \begin{cases} \frac{\cos(n\pi\delta/2)}{1-(n\delta)^2} e_n(z_c) & \text{if } n\delta \ne 1, \\ \frac{\pi}{4} e_n(z_c) & \text{otherwise.} \end{cases} \quad (88)$$

These coefficients $g_n$, function $g$ and its approximation with $N = 2$ and $N = 16$ modes are summarized in figure 4.

Two sets of numerical experiments are performed: (a) $N = 2$ modes, for which the improved convergence radius on the modal subspace is $\rho_2^\dagger \approx 1.0498$ and (b) $N = 16$ modes for which $\rho_{16}^\dagger \approx 1 + 2.64 \times 10^{-5} \approx \rho^\dagger$.

Three amplitudes are tested. They are chosen such that $\|\tilde{x}_1^\dagger\|_{X^\dagger} = A$ with $A \in \{1; 3; 5\}$ (still corresponding to $f_{tot} = A/\|\tilde{x}_1^\dagger\|_{X^\dagger}$ where computations yield $f_{tot} \approx 8.2002 \times A$ for $N = 2$ and $f_{tot} \approx 7.9533 \times A$ for $N = 16$). The deflection wave is observed at the center of the distributed force ($z = z_c = 1/7$).

## 7.3 Numerical results

### Configuration 1 (non oscillating beam, single mode)

Deflection signals are simulated according to fig. 3 and observed at the beam center ($w(z = 0.5, t) = [\tilde{x}_{\text{trunc}}(z = 0.5, t)]_1$), for truncation orders $M = 1, 3, 5, 7$. Results are displayed in figure 5.

For $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 0.8$, a good approximation is obtained as soon as $M \ge 3$. As guaranteed by (10), the convergence is numerically observed.

For $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 1 = \rho < \rho_N$ ($N = 1$ mode), the approximations are significantly more accurate when increasing the truncation order $M$ of the Volterra series (see also the zoom in figure 6).

This is no longer true for $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 1.2$, for which the convergence seems to be lost. For $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 2$, the divergence is so fast (on the range $1 \le t \le 3$) that the best approximation is the linear approximation. More precisely, nonlinear contributions $y_{m \ge 3}$ build unrealistic

high amplitudes signals after the very beginning of the trajectory. High order contributions $x_m$ have increasing amplitudes with $m$.

Thus, for this single mode system, the guaranteed bound is close to the exact convergence radius.

### Configuration 2 (multiple modes and damped oscillating waves)

Three figures represent the deflection waves (at $z = z_c = 1/7$) for $N = 2$ compared to $N = 16$ for three amplitudes: $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 1, 3, 5$ in figs. 7, 8 and 10, respectively.

In figure 7 ($\|\tilde{x}_1^\dagger\|_{X^\dagger} = 1 < \rho_{16}^\dagger < \rho_2^\dagger$), the convergence is guaranteed. The signals appear to be well approximated as soon as $M \ge 3$ (slight errors on the linear approximations are visible in the zoomed figures on periods [20,22] and [30,32]). As expected, signals built with $N=16$ modes (right column) are sharper (with a richer spectral content) than with $N=2$ modes (left column).

In figure 8 ($\|\tilde{x}_1^\dagger\|_{X^\dagger} = 3 > \rho_2^\dagger > \rho_{16}^\dagger > \rho^\dagger = 1$), the convergence is not guaranteed for any input waveform. However, for the tested excitation, the nonlinear approximations at order $M = 7$ all appear to be numerically accurate: the signal shape, the amplitude and the synchronization with the reference solution are correct. This suggests that the generic convergence bound may happen to be conservative when specific excitations are applied to the beam model. For the linear approximation, the shape and the synchronization are lost (accounting for an accumulated delay due to frequency shifts as expected in Duffing oscillators). Approximations at orders 3 and 5 are not sufficient but illustrate how the summation of contributions operate to reproduce the complete (hardening) nonlinear effect. Moreover, in the frequency domain, approximations at order 1 and 7 (figure 9) make clearly appear the frequency structure: the eigenfrequencies (imaginary part of the eigenvalues) are activated at order 1 and 7; the approximation at order 7 makes clearly appear additional super-harmonics and inter-modulation (additive and difference combination of these frequencies), generated by the nonlinear contributions.

In figure 10 ($\|\tilde{x}_1^\dagger\|_{X^\dagger} = 5$), the convergence is lost. As in the case $\|\tilde{x}_1^\dagger\|_{X^\dagger} = 2$ for configuration 1, the series approximation is efficient at the very beginning of the
<figure><img src="image_5.png" /><figcaption>Fig. 4: Configuration 2. (a) First coefficients $g_n$. (b) Reference unitary spatial distribution $g(z)$ (-) and its modal reconstructions with $N = 2$ modes (··) and $N = 16$ modes (---).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Configuration 1. Signals of the deflection wave at center.</figcaption></figure>

simulation (over [0,12] at order 7) but the nonlinear contributions eventually grow (the higher is order $m$, the worse is $y_m$) and produce large artifacts. This corresponds to the so-called phenomenon of *secular modes*. This phenomenon is well-known for the Duffing oscillator [22]. It is due to the nonlinear effect type: in this system, the nonlinearity is mainly responsible for a frequency modulation of the input signal for large magnitudes. The Volterra series expansion is a regular perturbation method, which attempts to represent and approximate this frequency modulation with polynomial combinations of fixed-frequency oscillating signals. In the case of a conservative problem, this kind of approx-

imation produces signal envelopes that increase as $t^{p-1}$ for orders $m = 2p + 1$. As a consequence, there is no possible convergence over $T = \mathbb{R}$ and $\rho$ is zero. In the case of the damped beam, the damping makes these envelopes behaves as $t^{p-1} \exp(\Re(e^{\lambda_n} \alpha_n t))$ for each mode $n$, so that signals $x_m$ are all bounded and the convergence radius is nonzero. From these results and our numerical simulations, the convergence radius $\rho$ appears to be a guaranteed bound on the linear response, and therefore on the input signal, for which high order contributions $x_m$ are not secular over $T$.

Finally, figure 11 illustrates the effect of the fluid damping coefficient decrease: the parameters of config-
<figure><img src="image_7.png" /><figcaption>Fig. 6: Configuration 1, $||\tilde{x}_1^\dagger||_{\mathcal{X}^\dagger} = 1$. Zoom of figure 5.</figcaption></figure>

uration 2 are kept except that $a$ is progressively decreased to 0.05, 0.02 and 0.01. In practice, it corroborates the previous observations for $||\tilde{x}_1^\dagger||_{\mathcal{X}^\dagger} = 1, 3, 5$ and also makes appear that the smaller is the damping, the larger are the difference between the linear approximation and the nonlinear response. This may be explained by the presence of modes with slightly higher, but still low damping and lower natural frequency than with higher values of $a$ as above, since according to section 4.1, H6, the lower bound for damping ratio is $2\sqrt{b/a}$, reached for natural frequencies close to $\sqrt{a/b}$, whose linear response may significantly excite the nonlinear part of the system.

## 8 Conclusion and perspectives

We proposed and analyzed a method, the pseudo-modal Volterra approximation, that combines Volterra series expansion and modal decomposition to provide a reduced order model representing the vibrations of a damped nonlinear beam. This reduced model has a simple structure (finite dimension linear systems with nonlinear interconnection), well suited for real time simulation.

The convergence analysis of the resulting expansion provides a computable domain, available for bounded excitation signals with any waveshape. Our theoretical results guarantee that in this domain the series expansion provides accurate approximations of the solution: this is well-adapted to account for distortions and timbres modification of sounds and vibrations resulting from the external input. Moreover, we showed that our theoretical results on pseudo-modal Volterra approximation can be extended to a large class of systems.

However, as Volterra series approximation is a regular perturbation method, it cannot accommodate situations where the system exhibits secular modes, in which case, other nonlinear approximation techniques such as POD can be best suited. A perspective is then to propose extensions of this work to efficiently represent such modulations, based on other perturbation methods and

approaches that are adapted to any input signal waveform.

## A Well-posedness and definitions

### A.1 General setting

The class of systems (1-2) is considered for the general following setting:

* $\mathbb{T}$ denotes the time set $\mathbb{R}_+$,

* $\mathbb{U}$ and $\mathbb{X}$ are Banach spaces on the field $\mathbb{R}$,

* $\mathcal{L}(\mathbb{U}, \mathbb{X})$ and $\mathcal{L}(\mathbb{X})$ are the sets of bounded linear operators from $\mathbb{U}$ to $\mathbb{X}$, and from $\mathbb{X}$ to $\mathbb{X}$, respectively,

* $\mathcal{ML}_k(\mathbb{X})$ ($k \ge 2$) is the set of bounded multilinear operators from $\mathbb{X} \times \cdots \times \mathbb{X}$ to $\mathbb{X}$, equipped with norm

$$ \|A_k\| = \sup_{\substack{(x_1, \dots, x_k) \in \mathbb{X}^k \\ \|x_1\| = \dots = \|x_k\| = 1}} \|A_k(x_1, \dots, x_k)\|_{\mathbb{X}}, $$

* $\mathcal{U} = L^\infty(\mathbb{T}, \mathbb{U})$ and $\mathcal{X} = L^\infty(\mathbb{T}, \mathbb{X})$ are standard Lebesgue spaces (in which trajectories $u$ and $x$ live).

### A.2 Mild solutions

For the linearized version of systems (1-2), mild solutions are defined in the following way [23, 7].

**Definition 1 (Mild solution of a linear system)** Let $u \in \mathcal{L}_{\text{loc}}^{\infty}(\mathbb{T}, \mathbb{U})$. The mild solution of the linearized system is the function $x \in C_0(\mathbb{T}, \mathbb{X})$ defined for all $t \in \mathbb{T}$ as

$$ x(t) = S(t)x_{\text{ini}} + \int_{0}^{t} S(t - \tau)Bu(\tau) \, d\tau. $$

A similar definition is given for nonlinear systems.

**Definition 2 (Mild solution of a nonlinear system)** Let $u \in \mathcal{L}_{\text{loc}}^{\infty}(\mathbb{T}, \mathbb{U})$. Then, $x$ is said to be a mild solution of (1-2) iff $x \in C_0(\mathbb{T}, \mathbb{X})$ and satisfies, for all $t \in \mathbb{T}$,

$$ x(t) = S(t)x_{\text{ini}} + \int_{0}^{t} S(t - \tau) \times (Bu(\tau) + P(x(\tau)) + Q(x(\tau), u(\tau))) \, d\tau. \quad (89) $$

As a standard result from [23], for a fixed $u \in C_0(\mathbb{T}, \mathbb{U})$, there exist $t_{\max} \in (0, T]$ and a unique function $x \in C_0([0, t_{\max}), \mathbb{X})$ such that $x$ is a mild solution in the sense of (89). For the class of systems (1-2), it can be easily shown that the local existence and uniqueness of mild solutions still holds when the input $u$ is taken in $\mathcal{L}_{\text{loc}}^{\infty}(\mathbb{T}, \mathbb{U})$, as stated in definition 2.

### A.3 Functional setting for the linearized beam problem

In the linear problem<fnref n="8" /> (24), the bi-Laplacian $\mathcal{B}$ is defined as the unbounded operator on $\mathbb{H} = L^2(0, 1)$ with domain<fnref n="9" />

$$ D(\mathcal{B}) = \{ w \in H^4(0, 1) \text{ s.t. } w(0) = w(1) = 0, \\ w''(0) = w''(1) = 0 \}, $$

<footnote n="8">⁸ For recall, $\partial_t^2 w + 2[aI + b\mathcal{B}](\partial_t w) + \mathcal{B}(w) = f$ on $[0, 1] \times \mathbb{T}$.</footnote>

<footnote n="9">⁹ $H^k(0, 1)$ is the standard Sobolev space of functions in $L^2(0, 1)$ whose first $k$ derivatives are in $L^2(0, 1)$.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 7: (Configuration 2, $||\tilde{x}_1^\dagger||_{X^\dagger} = 1$) Deflection waves at $z = 1/7$ for $N = 2$ modes (top and left figures) compared to $N = 16$ modes (right and bottom figures). Note that Zoomed parts (starting at 0, 10, 20, 30, 40, 70 respectively, with duration 2) have the same legend as the plain figures.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8: (Configuration 2, $||\tilde{x}_1^\dagger||_{X^\dagger} = 3$) To be compared to fig. 7.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9: (Configuration 2 $||\tilde{x}_1^\dagger||_{X^\dagger} = 3$, N=16 modes) Spectrograms of $w$: (top) linear approximation ; (bottom) order 7. Zooms in the low-frequency range (right column) detail the spectral differences.</figcaption></figure>

such that $\mathcal{B}(w) = w^{(4)}$ for all $w \in D(\mathcal{B})$. This operator is closed, densely defined, self-adjoint and positive on $\mathbb{H}$. Hence, we can introduce its uniquely defined positive square root $K$, with domain

$$D(\mathcal{K}) = \{w \in H^2(0, 1) \text{ s.t. } w(0) = w(1) = 0\},$$

such that $\mathcal{K}(w) = -\Delta w$. The domain $D(\mathcal{K})$ equipped with the $K$-norm defines a Hilbert space<fnref n="10" />:

$$\mathbb{H}^{\frac{1}{2}} := D(\mathcal{K}) \text{ with } \|\cdot\|_{\mathbb{H}^{\frac{1}{2}}} = \|\mathcal{K} \cdot\|_{\mathbb{H}}.$$

Then, according to [16], the linearized problem (24) admits the state-space representation (25-27), which is well-posed for the functional setting $\mathbb{U}$ and $\mathbb{X}$ introduced in (28-29) and defining operator $A$ on domain $D(A)$ as

$$A: D(A) \to \mathbb{X} \text{ with } D(A) = \{(u, v) \in \mathbb{H}^{\frac{1}{2}} \times \mathbb{H}^{\frac{1}{2}} \text{ s.t.} \\ 2(a + b \mathcal{B})v + \mathcal{B}u \in \mathbb{H}\}.$$

In this setting, $A$ generates a $C_0$ contraction semigroup on $\mathbb{X}$ ([16, (A1-A2), p. 6]), so that there exists a negative growth bound $\alpha < 0$. More precisely (see [16, corollary 5.2]), it is a Riesz spectral operator on $\mathbb{X}$ which generates an analytic semigroup $S$, provided that $-\frac{1}{b}$ is not in the point spectrum of $A$. In addition, operator $B$ belongs $\mathcal{L}(\mathbb{U}, \mathbb{X})$ and $\|B\| = 1$. Therefore, this linear problem is well-posed and is in the class of systems defined in section 2.1.

# B Estimates for the nonlinear beam model

## B.1 Proof of equation (32)

For all $X_1 \in \mathbb{H}^{\frac{1}{2}}$ and $z \in [0, 1]$, $X_1(z) = X_1'(0)z + \int_0^z (z - \zeta)X_1''(\zeta) d\zeta$ which is zero at $z = 1$. Hence,

$$X_1(z) = \int_0^1 K(z, \zeta) X_1''(\zeta) d\zeta,$$

where $K(z, \zeta) = \zeta(z - 1)$ if $0 \le \zeta \le z$ and $K(z, \zeta) = z(\zeta - 1)$ if $z \le \zeta \le 1$. Using the Cauchy-Schwartz inequality, it comes $|X_1(z)| \le \mu(z)(\int_0^1 X_1''(\zeta)^2 d\zeta)^{\frac{1}{2}}$, where $\mu(z) = (\int_0^1 K(z, \zeta)^2 d\zeta)^{\frac{1}{2}} = \frac{z(1-z)}{\sqrt{3}}$. Then, $(\int_0^1 |X_1(z)|^2 dz)^{\frac{1}{2}} \le \nu \|X_1\|_{\mathbb{H}^{\frac{1}{2}}}^2$, where $\nu = (\int_0^1 \mu(z)^2 dz)^{\frac{1}{2}} = \frac{1}{3\sqrt{10}}$. Finally, for all $(X, Y, Z) \in \mathbb{X}^3$,

$$\begin{aligned} \|A_3(X, Y, Z)\|_{\mathbb{X}} &\le \eta \left| \int_0^1 X_1'(z) Y_1'(z) dz \right| \int_0^1 |Z_1''(z)| dz \\ &\le \eta \left| \int_0^1 X_1(z) Y_1''(z) dz \right| \left( \int_0^1 |Z_1''(z)|^2 dz \right)^{\frac{1}{2}} \\ &\le \eta \left( \int_0^1 |X_1(z)|^2 dz \right)^{\frac{1}{2}} \|Y_1\|_{\mathbb{H}^{\frac{1}{2}}} \|Z_1\|_{\mathbb{H}^{\frac{1}{2}}} \\ &\le \nu \eta \|X_1\|_{\mathbb{H}^{\frac{1}{2}}} \|Y_1\|_{\mathbb{H}^{\frac{1}{2}}} \|Z_1\|_{\mathbb{H}^{\frac{1}{2}}} \\ &\le \nu \eta \|X\|_{\mathbb{X}} \|Y\|_{\mathbb{X}} \|Z\|_{\mathbb{X}}, \end{aligned}$$

proves the result.

<footnote n="10">¹⁰ This justifies the notation introduced in (29) where $||x_1||_{\mathbb{H}^{\frac{1}{2}}} = ||\mathcal{K} x_1||_{\mathbb{H}} = ||x_1^{(2)}||_{\mathbb{H}}$.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 10: (Configuration 2, $||\hat{x}_1^\dagger||_{\mathcal{X}^\dagger} = 5$) To be compared to fig. 7 and 8.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11: (Configuration 2 (N=16 modes) with smaller damping coefficient $a = 0.05, 0.02, 0.01$) To be compared to fig. 7, 8 and 10.</figcaption></figure>
## B.2 Proof of equations (33,43)

Functions $e_n^±$ form an orthogonal basis of $\mathcal{X}$, from which we define an orthonormal basis

$$E_n^+(z) = \frac{1}{k_n^2} \begin{bmatrix} 1 \\ 0 \end{bmatrix} s_n(z) \quad \text{and} \quad E_n^-(z) = \begin{bmatrix} 0 \\ 1 \end{bmatrix} s_n(z)$$

On each modal subspace, the linearized system behaves like second order system. Indeed, for all $x \in \mathbb{H}$

$$S(t)Bx = \sum_{n=1}^{\infty} \kappa_n^- S(t)E_n^- \quad (90)$$

where $\kappa_n^- = \langle Bx, E_n^- \rangle_{\mathcal{X}} = \langle x, s_n \rangle_{\mathbb{H}}$, and

$$\begin{aligned} E_n^+(z) &= \frac{(e_n^+ - e_n^-)\lambda_n^+ \lambda_n^-}{k_n^2 (\lambda_n^- - \lambda_n^+)} = \frac{1}{k_n^2} \begin{bmatrix} 1 \\ 0 \end{bmatrix} s_n(z) \\ E_n^-(z) &= \frac{(\lambda_n^+ e_n^+ - \lambda_n^- e_n^-)}{\lambda_n^+ - \lambda_n^-} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} s_n(z) \end{aligned} \quad (91)$$

From (91) we obtain

$$\begin{aligned} S(t)Bx = \sum_{n=1}^{\infty} \kappa_n^- \left( \frac{\lambda_n^+ e_n^+ t}{\lambda_n^+ - \lambda_n^-} e_n^+ + \frac{\lambda_n^- e_n^- t}{\lambda_n^- - \lambda_n^+} e_n^- \right) &= \sum_{n=1}^{\infty} \kappa_n^- \left( k_n^2 \frac{e_n^+ t - e_n^- t}{\lambda_n^+ - \lambda_n^-} E_n^+ + \frac{\lambda_n^+ e_n^+ t - \lambda_n^- e_n^- t}{\lambda_n^- - \lambda_n^+} E_n^- \right) \end{aligned} \quad (92)$$

Setting $h_n(t) = \frac{e^{\lambda_n^+ t} - e^{\lambda_n^- t}}{\lambda_n^+ - \lambda_n^-}$ and using the standard expressions

$$h_n(t) = \begin{cases} e^{-k_n^2 \xi_n t} \frac{\sin(k_n^2 \sqrt{1-\xi_n^2} t)}{k_n^2 \sqrt{1-\xi_n^2}} & \text{if } \xi_n < 1 \\ te^{-k_n^2 t} & \text{if } \xi_n = 1 \\ e^{-k_n^2 \xi_n t} \frac{\sinh(k_n^2 \sqrt{\xi_n^2-1} t)}{k_n^2 \sqrt{\xi_n^2-1}} & \text{if } \xi_n > 1 \end{cases}$$

we obtain that

$$\begin{aligned} \|S(t)Bx\|_{\mathcal{X}} &= \sqrt{\sum_{n=1}^{\infty} (\kappa_n^-)^2 (k_n^4 h_n(t)^2 + h_n'(t)^2)} \\ &\le \sqrt{\sup_n |k_n^4 h_n(t)^2 + h_n'(t)^2|} \|x\|_{\mathbb{H}} \end{aligned}$$

so that we can set

$$\gamma = \int_{\mathbb{T}} \sqrt{\sup_n |k_n^4 h_n(t)^2 + h_n'(t)^2|} dt \quad (93)$$

Now we notice that for all $(x, y, z)$ in $\mathcal{X}$, $A_3(x, y, z) = Bg(x, y, z)$ where $g$ is a third order multilinear operator from $\mathcal{X}^3$ to $\mathbb{H}$, such that $\|g\|_{\mathcal{ML}(\mathcal{X}, \mathbb{H})} = a_3$. We therefore obtain that for all $(a, b, c)$ in $\mathcal{X}$

$$\begin{aligned} \| \int_0^t S(t-\tau) Bg(a, b, c)(\tau) d\tau \|_{\mathcal{X}} &\le \int_0^t \| S(t) Bx \|_{\mathcal{L}(\mathbb{H}, \mathcal{X})} a_3 d\tau \|a\|_{\mathcal{X}} \|b\|_{\mathcal{X}} \|c\|_{\mathcal{X}} \\ &\le \gamma a_3 \|a\|_{\mathcal{X}} \|b\|_{\mathcal{X}} \|c\|_{\mathcal{X}} \end{aligned} \quad (94)$$

This proves that $\zeta_3 \le \gamma a_3$.

In the same way, for operator $C$, we obtain from (92) that

$$CS(t)Bx = \sum_{n=1}^{\infty} \kappa_n^- \frac{e^{\lambda_n^+ t} - e^{\lambda_n^- t}}{\lambda_n^+ - \lambda_n^-} s_n,$$

and

$$\begin{aligned} \|CS(t)Bx\|_{\mathcal{L}(\mathbb{H}, \mathbb{H}^{\frac{1}{2}})} &= \sqrt{\sum_{n=1}^{\infty} (\kappa_n^-)^2 k_n^4 h_n(t)^2} \\ &\le \sup_n (k_n^2 |h_n(t)|) \|x\|_{\mathbb{H}} \end{aligned}$$

It follows that $\|CS(t)Bx\|_{\mathcal{L}(\mathbb{H}, \mathbb{H}^{\frac{1}{2}})} \le \sup_n k_n^2 |h_n(t)|$, so that we can set

$$\begin{aligned} \gamma^\dagger &= \int_{\mathbb{T}} \sup_n (k_n^2 |h_n(t)|) dt \\ r &= \frac{\int_{\mathbb{T}} \sqrt{\sup_n |k_n^2 h_n(t)^2 + h_n'(t)^2|} dt}{\int_{\mathbb{T}} \sup_n (k_n^2 |h_n(t)|) dt} \end{aligned}$$

## C Proof of proposition 1 (truncation error estimate)

We assume (A1-A3) in section 5.3 and set $\Phi(z) = \sum_{m=1}^{\infty} \phi_m z^m$ where the sequence $(\phi_m)_{m \in \mathbb{N}^*}$ is defined by (see step 6 in page 4)

$$\phi_1 = 1 \text{ and, for all } m \ge 2, \quad \phi_m = \gamma a_3 \sum_{p \in \mathbb{M}_m^3} \phi_{p_1} \phi_{p_2} \phi_{p_3},$$

with $\mathbb{M}_m^3 = \{p \in (\mathbb{N}^*)^3 \mid p_1 + p_2 + p_3 = m\}$. Then, the trajectories $\tilde{x} = \sum_{m=1}^{\infty} \tilde{x}_m$ and $x = \sum_{m=1}^{\infty} x_m$ are normally convergent and such that, for all $m \ge 1$,

$$\begin{aligned} \|\tilde{x}_m\|_{\mathcal{X}} &\le \phi_m \|\hat{x}_1\|_{\mathcal{X}}^m, \\ \|x_m\|_{\mathcal{X}} &\le \phi_m \|x_1\|_{\mathcal{X}}^m \le \phi_m ((1+\epsilon)\|\hat{x}_1\|_{\mathcal{X}})^m \le \beta_m \|\hat{x}_1\|_{\mathcal{X}}^m, \end{aligned}$$

where the sequence $(\beta_m)_{m \in \mathbb{N}^*}$ defined by $\beta_m = (1+\epsilon)^m \phi_m$ is such that

$$\beta_1 = 1 + \epsilon \quad \text{and for all } m \ge 2, \quad \beta_m = \gamma a_3 \sum_{p \in \mathbb{M}_m^3} \beta_{p_1} \beta_{p_2} \beta_{p_3}.$$

Now, for all $m \ge 1$, denote the error terms $e_m := x_m - \tilde{x}_m$. The first terms are such that $\|e_1\|_{\mathcal{X}} = \epsilon \|\hat{x}_1\|_{\mathcal{X}}$ (from (A2)), $\|e_2\|_{\mathcal{X}} = 0$. Moreover, for $m \ge 3$, it follows from (9) that

$$\begin{aligned} e_m(t) &= \int_0^t S(t-\tau) \sum_{p \in \mathbb{M}_m^3} (A_3(x_{p_1}(\tau), x_{p_2}(\tau), x_{p_3}(\tau)) \\ &\quad - \hat{A}_3(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau))) d\tau \\ &= \int_0^t S(t-\tau) \sum_{p \in \mathbb{M}_m^3} F_p(\tau) d\tau, \end{aligned} \quad (95)$$

where we set

$F_p(\tau) = A_3(x_{p_1}(\tau), x_{p_2}(\tau), x_{p_3}(\tau)) - A_3(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau))$. It should be noted that (95) crucially depends on equation (53), from which in $\tilde{X}$, $\hat{A}_3$ is identical to $A_3$. Now, replacing $x_{p_i}$ by $\tilde{x}_{p_i} + e_{p_i}$ and exploiting the multi-linearity of $A_3$ yield the following expansion (omitting variable $\tau$ for sake of legibility)

$$\begin{aligned} F_p = A_3(\tilde{x}_{p_1}, \tilde{x}_{p_2}, e_{p_3}) &+ A_3(\tilde{x}_{p_1}, e_{p_2}, \tilde{x}_{p_3}) + A_3(e_{p_1}, \tilde{x}_{p_2}, \tilde{x}_{p_3}) \\ &+ A_3(\tilde{x}_{p_1}, e_{p_2}, e_{p_3}) + A_3(e_{p_1}, \tilde{x}_{p_2}, e_{p_3}) + A_3(e_{p_1}, e_{p_2}, \tilde{x}_{p_3}) \\ &+ A_3(e_{p_1}, e_{p_2}, e_{p_3}). \end{aligned}$$
Then, introducing $\psi_m = \beta_m - \phi_m$ for all $m \ge 1$, we prove by induction that (claim $C_m$) $\|e_m\|_{\mathcal{X}} \le \psi_m \|\hat{x}_1\|_{\mathcal{X}}^m$:

(m = 1): the claim ($C_1$) is true for $m=1$ by construction; (m ≥ 2): assume that $C_p$ holds for all $p \le m-1$, then using expressions above, it follows that

$$
\begin{aligned}
\|e_m\|_{\mathcal{X}} &\le \gamma a_3 \sum_{p \in \mathbb{M}_m^3} \|\hat{x}_1\|_{\mathcal{X}}^{(p_1+p_2+p_3)} (\phi_{p_1} \phi_{p_2} \psi_{p_3} \\
&\qquad + \phi_{p_1} \psi_{p_2} \phi_{p_3} + \psi_{p_1} \phi_{p_2} \phi_{p_3} + \phi_{p_1} \psi_{p_2} \psi_{p_3} \\
&\qquad + \psi_{p_1} \phi_{p_2} \psi_{p_3} + \psi_{p_1} \psi_{p_2} \phi_{p_3} + \psi_{p_1} \psi_{p_2} \psi_{p_3}) \\
&\le \gamma a_3 \|\hat{x}_1\|_{\mathcal{X}}^m \sum_{p \in \mathbb{M}_m^3} (\beta_{p_1} \beta_{p_2} \beta_{p_3} - \phi_{p_1} \phi_{p_2} \phi_{p_3}) \\
&\le (\beta_m - \phi_m) \|\hat{x}_1\|_{\mathcal{X}}^m = \psi_m \|\hat{x}_1\|_{\mathcal{X}}^m,
\end{aligned}
$$

so that $C_m$ is satisfied.

An immediate consequence is that

$$
\begin{aligned}
\|e\|_{\mathcal{X}} &\le \sum_{m=1}^{\infty} \psi_m \|\hat{x}_1\|_{\mathcal{X}}^m = \sum_{m=1}^{\infty} \beta_m \|\hat{x}_1\|_{\mathcal{X}}^m - \sum_{m=1}^{\infty} \phi_m \|\hat{x}_1\|_{\mathcal{X}}^m \\
&\le \Phi((1+\epsilon)\|\hat{x}_1\|_{\mathcal{X}}) - \Phi(\|\hat{x}_1\|_{\mathcal{X}}),
\end{aligned}
$$

which concludes the proof.

## D Proof of generalized result

We assume (A1-A3) in section 5.3 and consider a beam model with a third order nonlinearity $A_3$ for which $\tilde{X}$ is not invariant. We assume that a bound $a_3$ of $\|A_3\|_{\tilde{X}}$ is available, and define $\Phi(z) = \sum_{m=1}^{\infty} \phi_m z^m$ as in appendix C.

We denote $\tilde{A}_3 = \Pi_{\tilde{X}} A_3$. Since $\Pi_{\tilde{X}}$ is an orthogonal projection, $\|\tilde{A}_3\|_{\tilde{X}} \le \|A_3\|_{\tilde{X}} \le a_3$, and therefore, as in appendix C, the trajectories $\tilde{x} = \sum_{m=1}^{\infty} \tilde{x}_m$ and $x = \sum_{m=1}^{\infty} x_m$ are convergent and satisfy, for all $m \ge 1$,

$$
\begin{aligned}
\|\tilde{x}_m\|_{\mathcal{X}} &\le \phi_m \|\hat{x}_1\|_{\mathcal{X}}^m, \\
\|x_m\|_{\mathcal{X}} &\le \phi_m \|x_1\|_{\mathcal{X}}^m \le \phi_m ((1+\epsilon)\|\hat{x}_1\|_{\mathcal{X}})^m \le \beta_m \|\hat{x}_1\|_{\mathcal{X}}^m,
\end{aligned}
$$

where the sequence $(\beta_m)_{m \in \mathbb{N}^*}$ was defined in appendix C. Now, for all $m \ge 1$, denote the error terms $e_m := x_m - \tilde{x}_m$. The first terms are such that $\|e_1\|_{\mathcal{X}} = \epsilon \|\hat{x}_1\|_{\mathcal{X}}$ (from (A2)), $\|e_2\|_{\mathcal{X}} = 0$. Moreover, for $m \ge 3$, it follows from (9) that

$$
\begin{aligned}
e_m(t) &= \int_0^t S(t-\tau) \sum_{p \in \mathbb{M}_m^3} (A_3(x_{p_1}(\tau), x_{p_2}(\tau), x_{p_3}(\tau)) \\
&\qquad - \hat{A}_3(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau))) d\tau \\
&= \int_0^t S(t-\tau) \sum_{p \in \mathbb{M}_m^3} F_p(\tau)d\tau + \int_0^t S(t-\tau) \sum_{p \in \mathbb{M}_m^3} G_p(\tau)d\tau
\end{aligned}
\tag{96}
$$

where we set

$$
\begin{aligned}
F_p(\tau) &= A_3(x_{p_1}(\tau), x_{p_2}(\tau), x_{p_3}(\tau)) - A_3(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau)) \text{ and } 10. \\
G_p(\tau) &= (A_3 - \hat{A}_3)(\tilde{x}_{p_1}(\tau), \tilde{x}_{p_2}(\tau), \tilde{x}_{p_3}(\tau))
\end{aligned}
$$

Setting

$$
\|A_3 - \hat{A}_3\|_{\tilde{X}} = \sup_{\substack{(a,b,c) \in (\tilde{X})^3 \\ \|a\|_{\tilde{X}}, \|b\|_{\tilde{X}}, \|c\|_{\tilde{X}} = 1}} \|A_3(a,b,c) - \hat{A}_3(a,b,c)\|_{\mathcal{X}},
$$

we have

$$
\|G_p\|_{\mathcal{X}} \le \|A_3 - \hat{A}_3\|_{\tilde{X}} \|\hat{x}_1\|_{\mathcal{X}}^m \phi_{p_1} \phi_{p_2} \phi_{p_3}
$$

Then, following the same steps as in appendix C, we introduce $\psi_m = \beta_m - \phi_m$, and obtain by induction that for all $m \ge 1$, $\|e_m\|_{\mathcal{X}} \le (\psi_m + \frac{\|A_3 - \hat{A}_3\|_{\tilde{X}}}{a_3} \phi_m) \|\hat{x}_1\|_{\mathcal{X}}^m$. An immediate consequence is that

$$
\begin{aligned}
\|e\|_{\mathcal{X}} &\le \sum_{m=1}^{\infty} (\psi_m + \frac{\|A_3 - \hat{A}_3\|_{\tilde{X}}}{a_3} \phi_m) \|\hat{x}_1\|_{\mathcal{X}}^m \\
&\le \sum_{m=1}^{\infty} \beta_m \|\hat{x}_1\|_{\mathcal{X}}^m - \left(1 - \frac{\|A_3 - \hat{A}_3\|_{\tilde{X}}}{a_3}\right) \sum_{m=1}^{\infty} \phi_m \|\hat{x}_1\|_{\mathcal{X}}^m \\
&\le \Phi\left((1+\epsilon)\|\hat{x}_1\|_{\mathcal{X}}\right) - \left(1 - \frac{\|A_3 - \hat{A}_3\|_{\tilde{X}}}{a_3}\right) \Phi(\|\hat{x}_1\|_{\mathcal{X}}),
\end{aligned}
$$

which concludes the proof.

## Conflict of Interest

The authors declare that they have no conflict of interest.

## References

1. Antoulas, A.C.: Approximation of large-scale dynamical systems. SIAM (2005)

2. Chatterjee, A., Vyas, N.S.: Convergence analysis of Volterra series response of nonlinear systems subjected to harmonic excitation. Journal of Sound and Vibration **236**(2), 339 – 358 (2000). DOI https://doi.org/10.1006/jsvi.2000.2967

3. Chen, K.K., Tu, J.H., Rowley, C.W.: Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses. Journal of nonlinear science **22**(6), 887–915 (2012)

4. Cheng, C., Peng, Z., Zhang, W., Meng, G.: Volterra-series-based nonlinear system modeling and its engineering applications: A state-of-the-art review. Mechanical Systems and Signal Processing **87**, 340–364 (2017). DOI https://doi.org/10.1016/j.ymssp.2016.10.029

5. Cheng, C., Peng, Z., Zhang, W., Meng, G.: Volterra-series-based nonlinear system modeling and its engineering applications: A state-of-the-art review. Mechanical Systems and Signal Processing **87**, 340–364 (2017)

6. Crouch, P.E., Collingwood, P.C.: The observation space and realizations of finite Volterra series. SIAM journal on control and optimization **25**(2), 316–333 (1987)

7. Curtain, R.F., Zwart, H.J.: An introduction to Infinite-dimensional Linear Systems Theory. Springer (1995)

8. Fliess, M., Lamnabhi, M., Lamnabhi-Lagarrigue, F.: An algebraic approach to nonlinear functional expansions. IEEE Trans. Circuits Syst. **30**(8), 554–570 (1983)

9. Géradin, M., Rixen, D.: Théorie des vibrations. Masson (1996)

10. Gilbert, E.G.: Functional expansions for the response of nonlinear differential systems. IEEE Trans. Autom. Control **22**, 909–921 (1977)

11. Graff, K.F.: Wave motion in elastic solids. Dover, New York (1991)

12. Hassanabadi, M.E., Heidarpour, A., Azam, S.E., Arashpour, M.: Recursive principal component analysis for model order reduction with application in nonlinear bayesian filtering. Computer Methods in Applied Mechanics and Engineering **371**, 113334 (2020)
13. Hélie, T.: Modélisation physique d'instruments de musique et de la voix : systèmes dynamiques, problèmes directs et inverses. Tech. rep., Université Pierre et Marie Curie-Paris 6 (2013). Habilitation à Diriger des Recherches

14. Hélie, T., Laroche, B.: Computation of convergence bounds for volterra series of linear analytic single input systems. IEEE Transactions on Automatic Control **56-9**, 2062–2072 (2011). URL http://articles.ircam.fr/textes/Helie10b/

15. Hélie, T., Laroche, B.: Computable convergence bounds of series expansions for infinite dimensional linear-analytic systems and application. Automatica **50-9**, 2334–2340 (2014)

16. Jacob, B., Trunk, C., Winklmeier, M.: Analyticity and Riesz basis property of semigroups associated to damped vibrations. J Evol Eq **8**(2), 263–281 (2008)

17. Kappagantu, R., Feeny, B.: Part 1: Dynamical characterization of a frictionally excited beam. Nonlinear Dynamics **22**(4), 317–333 (2000)

18. Kappagantu, R., Feeny, B.: Part 2: Proper orthogonal modal modeling of a frictionally excited beam. Nonlinear Dynamics **23**(1), 1–11 (2000)

19. Li, L., Billings, S.: Analysis of nonlinear oscillators using Volterra series in the frequency domain Part I : convergence limits. Research report, acse research report no.988. automatic control and systems engineering, University of Sheffield (2009)

20. Li, L., Billings, S.: Analysis of nonlinear oscillators using volterra series in the frequency domain. Journal of Sound and Vibration **330**(2), 337 – 355 (2011). DOI https://doi.org/10.1016/j.jsv.2010.08.016

21. Liang, Y., Lee, H., Lim, S., Lin, W., Lee, K., Wu, C.: Proper orthogonal decomposition and its applications–part i: Theory. Journal of Sound and vibration **252**(3), 527–544 (2002)

22. Nayfeh, A., Mook, D.: Nonlinear Oscillations. Wiley Classics Library. Wiley (2008). URL https://books.google.fr/books?id=sj3ebg7jRaoC

23. Pazy, A.: Semigroups of Linear Operators and Applications to Partial Differential Equations, *Applied mathematical sciences*, vol. 44. Springer (1983)

24. Peitz, S., Klus, S.: Koopman operator-based model reduction for switched-system control of pdes. Automatica **106**, 184–191 (2019)

25. Peng, Z., Lang, Z.: On the convergence of the volterra-series representation of the duffing's oscillators subjected to harmonic excitations. Journal of Sound and Vibration **305**(1), 322 – 332 (2007). DOI https://doi.org/10.1016/j.jsv.2007.03.062

26. Roze, D., Hélie, T.: Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string. Journal of Sound and Vibration **333-7**, 2073–2086 (2014)

27. Rugh, W.J.: Nonlinear System Theory, The Volterra/Wiener approach. The Johns Hopkins University Press (1981)

28. Schetzen, M.: The Volterra and Wiener theories of non-linear systems. Wiley-Interscience (1989)

29. Volterra, V.: Theory of Functionnals and of Integral and Integro-Differential Equations. Dover Publications (1959)

30. Zhu, Y.P., Lang, Z.Q.: A new convergence analysis for the volterra series representation of nonlinear systems. Automatica **111**, 108599 (2020)