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

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 and a pseudo-modal truncation defined from the eigenbasis of the linearized problem. The interest of this approach lies in the large class of input waveshapes 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 depending 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 refinement 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.


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 particular nonlinear PDEs, is an important problem in engineering 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 [1,21]) coupled with Galerkin projection are very popular (see, for instance, [17,18] and more recently [12]). Koopman-based methods (see, for instance, [3,24]) are also becoming increasingly 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 conditions.
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 realtime simulation purposes, since they transform a nonlinear 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 [6,8,10,27,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,20,25] 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-dimensional 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 Sect. 2, and the beam model is introduced in Sect. 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 Sect. 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 Sect. 7.

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.

Class of nonlinear systems
We consider systems excited on positive times T = R + by inputs u : T → U, with state-space representation of state x : T → X, governed bẏ with initial condition where A : X → X and B : U → X are linear operators, A k : X k → 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 (X = R n , U = R m ). The state x and the input u can be vectors, in which case A and B can be represented by n × n and n × m matrices and the A k 's are multilinear functions. Hypothesis (H0) means that A is a Hurwitz stable matrix 1 : A generates the flow S(t) = exp(At) on [0, T ], whose operator norm is bounded by β exp(αt) for some β > 0 and a negative growth bound 2 α < 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 L(X), L(U, X) (resp., 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 α < 0 and with β > 0 such that for all t ∈ T, S(t) L(X,X) ≤ β exp(αt). Systems under consideration. The class of problems (1-2) with (H0) is examined in the framework of bounded signals, namely case for which well-posed definitions of solutions (afternamed mild solutions) can be introduced (see details in "Appendix A", see also [7,23] and [13,15] for more general classes of systems).

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 ini , given by 1 Every eigenvalue of A has strictly negative real part. 2 The maximal real part of the eigenvalues of A is the smallest optimal bound α.
where the terms satisfy the sequence of linear problemṡ with where the index set M k m := p ∈ (N * ) k p 1 +· · ·+ 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 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) Eq. (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 ε (u = εũ, x ini = ε x ini ), by sorting terms along the powers of this marker ε (x m = ε mx 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ẏ = Ay + v) and instantaneous nonlinearities (multilinear operators A k with k ≤ m fed by subsets of yet-simulated signals x j with j ≤ m − 1) [15, Rk2 and Fig.1] (see also Fig. 3).

Convergence domain and truncation error bound
The series (4-9) converges towards a mild solution of (1-2) in the space X of bounded trajectories, if the solution (8) of the linearized problem is bounded as follows (see [5,14,15]): (11) where the convergence radius ρ, gain bound function Φ and error bound function R M Φ are characteristic quantities of the system. In practice, they can be derived as follows: Step 1 Compute the upper bounds ζ k , for 2 ≤ k ≤ K , Step 2 Define the function F, analytic at z = 0, Step 3 Compute the unique root σ defined by Step 4 Compute the convergence radius ρ, given by Step 5 Find the unique function analytic at z = 0 Step 6 Compute the remainder function given by Note that the coefficients φ m are given by φ 1 = 1 and, for all 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 So, a sufficient condition for the series convergence is It is convenient in practice, since computing x 1 is not required. But it is conservative compared to x 1 X < ρ.

Physical model: nonlinear boundary problem, state-space representation and well-posedness
We consider an Euler-Bernoulli model of a damped nonlinear pinned beam, initially at rest. First, the governing equations of the nonlinear boundary problem are presented in Sect. 3.1. Then, the problem is reformulated as a well-posed state-space representation (1)(2). The linearized problem is first examined in Sect. 3.2, and then, the full nonlinear problem in Sect. 3.3.
A dimensionless model governing the deflection waves w : (z, t) ∈ [0, 1] × T → R is given by (T = R + ), for all (z, t) ∈ [0, 1] × T, with pinned-type boundary conditions (fixed extremities and no momentum) and zero initial conditions Coefficients a > 0 and b > 0 are fluid and structural damping parameters and η > 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 (·, t) belongs to the Lebesgue space for each fixed time t. Moreover, we assume that (H4) these integrals are bounded in time, meaning that functions t → f (·, t) belong to L ∞ (T, H). Note that, compared to structured forcing terms 3 , the range largeness of such excitation signals gives the following method one of its main interests.

Linearized problem and state-space representation
The governing equations of the linearized boundary problem ((20-21) with η = 0) can be reformulated as for all (z, t) ∈ [0, 1] × T, with zero initial conditions (22), where B is denotes the bi-Laplacian (B(w) = w (4) ), operating on the domain 4 of sufficiently regular functions that satisfy the boundary conditions (21). This problem admits the state-space representatioṅ where the state x, input u, linear operators A and B are This formulation is well-posed in the spaces of bounded signals (3) for u(t) and x(t) living in, respectively, and norm x 1 H H , where the complete functional setting is detailed in "Appendix A.3". In this setting, A belongs to L(X) and generates a strongly continuous semigroup S with negative growth bound, B belongs to L(U, X) and B = 1, proving that system (25)(26)(27) is in the class defined in Sect. 2.1. 3 few oscillatory signals, sweeps or sequences of impulses. 4 see "Appendix A.3" for the detailed construction of the functional setting.

Nonlinear problem and Volterra series convergence
The nonlinear problem (η > 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: with multi-linear operator This operator defined from X 3 to X belongs to ML 3 (X). Indeed, its norm is bounded (proof in "Appendix B.1") as Hence, the state-space representation (30) belongs to the class of the well-posed nonlinear problems introduced in Sect. 2.1. According to Sect. 2.3, there exists a positive ρ such that the Volterra series expansion (4-9) is convergent in norm if x 1 X < ρ. Such a convergence radius ρ can be derived using steps 1-4. Computations in "Appendix B.2" provide a bound estimate ζ 3 related to γ (defined in (18)) and to a 3 (defined in (32)), and F is given by Solving F(σ )−σ F (σ ) = 0 leads to σ = (3ζ 3 ) − 1 2 . Finally, the convergence radius ρ = σ/F(σ ) is given by and ρ u = ρ/γ in Remark 3. Expressions of γ , ζ 3   Its first Taylor coefficients are given by φ 2 p = 0 for all p ≥ 1 and is deduced. The gain bound function Φ and its remainder R M Φ are displayed in Fig. 1.
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.

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 λ ± n of the characteristic polynomials P n (λ) = λ 2 + 2(a + bk 4 n )λ + k 4 n with k n = nπ, (37) for all n ≥ 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 ξ 1 := (a+bπ 4 )/π 2 (damping ratio) is smaller than 1, (H6) higher modal numbers n ≥ 2 yield faster decays 5 than n = 1, which is true if 2b(a + bπ 4 ) ≤ 1. Note that when b > 0 the corresponding second order systems can be oscillating at low n (if ξ n := (a + bk 4 n )/k 2 n is smaller than 1) and eventually become overdamped. A lower bound for the damping ratios is given by 2 Assuming that λ ± n are always distinct, the associated eigenfunctions of A are and k n = nπ (n ≥ 1).
Functions s n form a orthonormal basis of H and functions e ± n form a basis of X. For N > 0 (modal truncation order), we introduce the approximation subspaces of H = U (N -dimensional), H 1 2 (N -dimensional) and X (2N -dimensional) generated by the first related basis functions, as follows: and . X , respectively. Using previous notations, we introduce the orthogonal basis {S 1 , · · · , S 2N } of X defined by This basis is the one used in the sequel. The associated functional spaces U and X are summarized in Table 1 (column 2).

Convergence bound estimate
A bound γ satisfying (18) can be derived as (see "Appendix B.2" for details) where the impulse responses h n of second-order linear systems are given according to the damping ratio ξ n by Moreover, Eqs. (33-34) become The convergence radius denominator is proportional to the square root of the nonlinear stiffness η and to that of γ 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 √ b/a constant or increasing will also help the convergence, and that no damping (γ → +∞) prevents convergence 6 . Note also that, from Remark 3, a possibly conservative bound on the system input is

Projection of the linearized problem
We denote Π X : X → X the orthogonal projection on X. By construction of the modal basis, if x ∈ X is the solution of (25) with input u ∈ U and initial condition x ini ∈ X, then x = Π X (x) ∈ X is the solution of the same problem with input u = Π H (u) ∈ U and initial condition x ini = Π X (x ini ) ∈ X. This corresponds to a finite-dimensional problem in which operators A and B can be replaced by their restrictions respectively. This problem can be restated for the coordinates in the basis {S 1 , · · · , S 2N }, by introducing input v, initial condition y ini and state y where spaces V, Y and Y are defined in Table 1 (column 3) and Y is equipped with the norm built from X The coordinates y ∈ Y of the trajectories are governed bẏ with Note that by definition of the modal basis, when N → +∞, the truncated system trajectory given by (47) converges in norm towards the solution of (25) in X .

Projection of the nonlinear terms
Considering operator A 3 defined in (31), we notice that in the modal subspace X 3 we have Notations of the spaces involved in: (1) the original state-space problem, (2) the projected problem, (3) the representation of (2) with coordinates Original spaces Projection sub-spaces Modal coordinates spaces This means that the restriction of A 3 on X defines a multilinear operator on this space, that we denote A 3 .
Consider the trajectory x of the nonlinear beam model (30) with initial condition x ini and input u, respectively, in X and U. Then, we define the pseudomodal truncation of x on X as the solution x of probleṁ Note that x does not identify with x := Π X x for a solution x in general.
For the sake of conciseness, in the sequel x will be referred to as the pseudo-modal truncation of x omitting the input u and initial condition x ini . As in Sect. 4.3, the finite-dimensional problem (54) satisfied by x on X can be restated on Y for the coordinates in the basis S aṡ where Remark 4 (Duffing oscillators) If the N first damping ratios ξ n are less than 1, then the pseudo-modal system (55-56) exactly corresponds to N coupled damped Duffing oscillators.

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 pseudomodal approximation, with modal truncation order N ≥ 1 (defined in (54)) and Volterra series truncation order M ≥ 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.

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 ≥ 2 its pseudo-modal truncation x governed by system (54) and introduce the first M terms of the Volterra expansion of x as Then, we define the pseudo-modal Volterra approximation of x as We first examine the convergence of this series, and second, the error between x and the approximation X M N .

Convergence
We apply the same steps as in Sects. 3.3 and 4.2 to (54) and the Volterra series decomposition (57-58). This yields the estimates and Then, we build and the gain bound function (see (35)) Note that γ N , a 3N and ζ 3N are all strictly lower than, respectively, γ , a 3 and ζ 3 and tend towards them as N → ∞. Then, the modal truncation makes the convergence radius ρ N increase according to the dilatation factor through which the gain bound functions are related as Φ N (z) = r N Φ(z/r N ). This factor combines a first dilatation factor √ a 3 /a 3N > 1 depending on N only (see values in Table 2) and a second one 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 ρ N and gain bound func-

Truncation error estimates
Based on these results, we can now address the derivation of an estimate of the error bound of the pseudomodal Volterra approximation. Consider now an input u = f ∈ U and an initial condition x ini ∈ X such that for relative errors of 0.1%, 1% and 10% (see (35) and (65)) x 1 X < ρ. Assume that u and the linear response x 1 can be accurately described by their N -order modal decomposition u and x 1 in the sense that 1 denote the responses of the linearized system excited by u and u ⊥ , respectively.
As x 1 and x 1 are in the convergence domain and after (64), the series expansions x = ∞ m=1 x m for the system excited by u and initial condition x ini , and x = ∞ m=1 x m for the pseudo-modal truncation excited by u and initial condition x ini are both convergent in norm. The error on the approximated trajectory x is bounded as stated in the following proposition.

Proposition 1 (pseudo-modal truncation error bound)
Assume (A1-A3). Denote e = x − x the error on x due to the pseudo-modal truncation on X. Then, e is bounded as The proof is given in Appendix C. Bound estimates are given in Fig. 2 for several values of . 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

Generalization of pseudo-modal Volterra approximation
The pseudo-modal Volterra approximation (57-58) can be generalized to models for which A 3 does not fulfil (53), or to more general models (1) such that A k ( X k ) ⊆ 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 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 X is not stable, that is, A 3 ( X) is not a subset of 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 Sect. 2.1. It follows that we obtain a convergence bound ρ and a gain bound function Φ given by Eqs. For a solution of this new beam model x, with initial condition x ini and input u, we define its pseudo-modal truncation at order N x m as in Eq. (54), where A 3 is now defined as A 3 = Π X A 3 , the projection of A 3 on the modal subspace X. Note that thanks to this definition, x belongs to X.
The pseudo-modal Volterra approximation of x is now X M N = M m=1 x m , the sum of the first M terms of the Volterra expansion of x defined as in (57-58). As in previous section, the Volterra series expansion of x on X has a convergence radius ρ greater than ρ since, as Π X is an orthogonal projection, A 3 ≤ A 3 ≤ a 3 , and a gain bound function Φ. Then, defining 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 − x the error on x due to the pseudo-modal truncation on X. Then, e is bounded as where D Φ(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 X M N is bounded by 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 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 X stable for the linear part of the system, then Proposition 2 and Corollary 2 hold replacing

Refinements on convergence bounds
This section presents new results to improve convergence and error bound estimates. We refine the theoretical results in Sect. 2 for specific classes of nonlinearities (Sect. 6.1) and apply them to the beam (Sect. 6.2).

New theoretical results for a class of nonlinearities
The convergence bounds in Sect. 2.3 can be improved if the multilinear operators A k : X k → X only act on subspaces of X and admit a sandwich decomposition where C : Typically, X † and W † k refer to spaces (R n † ≤n or Banach) smaller than X, and A † k to a reduced version of A k . For clarity, we stamp all the labels of the smaller or reduced objects with the dagger symbol † .
The idea is to examine the convergence of (4) on x through that on Indeed, its Volterra series is the sum of x † m = C x m , straightforwardly derived from (8-9) and (69) as where we have defined, for 1 ≤ k ≤ K and B 1 = B, Then, we can compute new convergence estimates that benefit from the focused action of operators S † k , A † k : compared to S, A k , these operators are circumscribed to smaller spaces and expected to yield lower bound estimates than ζ k in (12). Thus, bounds ζ k are replaced in step 1 by the new estimates, for 2 ≤ k ≤ K , some overestimated values of which are γ † k a † k with Using ζ † k instead of ζ k in steps 2 to 6, we successively define F † , ζ † , σ † , ρ † and Φ † and are ready to state the following theorem.
Theorem 1 (Refined Volterra convergence and error bounds) Assume that x † 1 X † < ρ † . Then, the Volterra series of x † and x converge in norm in X † and X , respectively. Moreover, the following inequalities hold: The proof follows exactly the same steps as in [15].

Remark 5 (Supplement to Remark 3) An alternative convergence radius on the input is
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)ξ + K k=2 α k ξ k = u, for which we can choose x † = ξ and A † k (ξ 1 , . . . , ξ k ) = α k ξ 1 . . . ξ 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 ω (t) = U e iωt (harmonic regime), then ξ 1 (t) = H (iω)u ω (t) leads to a convergence bound for amplitude U that depends on the pulsation, namely,

Application to the beam model
Operator A 3 defined in (31) admits a sandwich decomposition (69). It involves the deflection wave from which the nonlinear term in (20) is built using and contributes to the state equation through (as input u) with limit γ † , which are strictly lower than γ N and γ , respectively, and from which we build ζ † 3N = a 3N γ † N and ζ † = a 3 γ † . Finally, the improved convergence radii are then and, in Theorem 1, factors r N (forx) and r (for x) are given by We may observe that the convergence criterion on x † 1 is actually weaker than the one on x 1 . Indeed, if x 1 X < ρ, then x † 1 X † < C L(X,X † ) x 1 X = x 1 X < ρ < ρ † , and the convergence radius dilatation is ρ † /ρ = √ r > 1. Moreover, following Theorem 1 and assuming (A1-A3) as in Proposition 1 for u, x † 1 , ρ † instead of u, x 1 , ρ, Corollary 1 is improved as follows.
Corollary 3 (Improved pseudo-modal Volterra approximation error) The error due to the pseudo-modal Volterra approximation X M N is bounded by In the specific situation where N = 1 (the beam is modelled as a single Duffing oscillator), setting h(t) = k 2 1 h 1 (t), we obtain γ † 1 = h 1 and the convergence bound ρ † 1 is expressed as Remark 6 (Single Duffing oscillator in harmonic regime) Following Remark 5, the resulting conservative bound on the systems input ρ † 1u = ρ † 1 / 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 As |H (iω)| < h 1 , this last bound is close, but slightly less than the bound obtained by these authors, whose expression is 2/ (3|H (iω)|) 3 2 √ |α 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. 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 Fig. 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 v y 1 y 1 y 3 y 3  (47) solver (lsim, Matlab) and the evaluation of instantaneous operators is exact. Moreover, to serve as a reference, the trajectory of the nonlinear problem (54) on 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 Fig. 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.

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 ρ † 1 (see (83)). Configuration 2 corresponds to the beam in oscillating regime, examined with respect to ρ † (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 η is chosen to provide a unitary convergence radius (ρ † 1 for configuration 1 and ρ † for the configuration 2). This yields where the total force f tot is spatially distributed according to g > 0 (where 1 0 g(z)dz = 1).  As mentioned in Sect. 7.1, we simulate the pseudomodal Volterra termsx m governed by (57-58) to build the approximation (59), through the realization in Fig. 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π 2 ).
The excitation (87) is spatially distributed on the first mode only (g(z) = s 1 (z) = √ 2 sin(π z) defined in (38)) and the step duration is T e = 3. We obtain γ † 1 = 0.1013 and (86) yields η = 14.434. Four amplitudes are tested. They are chosen such that x † 1 X † = A with A ∈ {0.8; 1; 1.2; 2}. They correspond to excitation forces f tot = A/ x † 1 X † where x † 1 X † is computed for the linear responsex † 1 to unit force. This yields f tot ≈ 9.8696 × 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 γ † = 8.774986 and η ≈ 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π f 0 /( f s ms 1 ) where (s 1 , s 1 ) are the poles associated with the first mode (T ≈ 2.4536 × 10 −3 for f 0 = 185Hz, f s = 48000Hz).
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 δ = 1/8, centred at z c = 1/7 described by g(z) = G cos(π(z−z c )/δ) on [z c −δ/2, z c −δ/2] and zero outside, with G = π/(2δ). The coefficients of the modal decomposition of g are, for all n ≥ 1,
For x † 1 X † = 0.8, a good approximation is obtained as soon as M ≥ 3. As guaranteed by (10), the convergence is numerically observed.
For x † 1 X † = 1 = ρ < ρ 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 Fig. 6). This is no longer true for x † 1 X † = 1.2, for which the convergence seems to be lost. For x † 1 X † = 2, the divergence is so fast (on the range 1 ≤ t ≤ 3) that the best approximation is the linear approximation. More precisely, nonlinear contributions y m≥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: x † 1 X † = 1, 3, 5 in figs. 7, 8 and 10, respectively.
In Fig. 7 ( x † 1 X † = 1 < ρ † 16 < ρ † 2 ), the convergence is guaranteed. The signals appear to be well approximated as soon as M ≥ 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 Fig. 8 , 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 (Fig. 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 Fig. 10 ( x † 1 X † = 5), the convergence is lost. As in the case x † 1 X † = 2 for configuration 1, the series approximation is efficient at the very beginning of the 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 artefacts. 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 polyno-mial combinations of fixed-frequency oscillating signals.
In the case of a conservative problem, this kind of approximation produces signal envelopes that increase as t p−1 for orders m = 2 p +1. As a consequence, there is no possible convergence over T = R and ρ is zero. In the case of the damped beam, the damping makes these envelopes behaves as t p−1 exp e(λ n )α 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 ρ 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, Fig. 11 illustrates the effect of the fluid damping coefficient decrease: the parameters of configuration 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 x † 1 X † = 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 Sect. 4.1, H6, the lower bound for damping ratio is 2 √ b/a, reached for natural frequencies close to √ a/b, whose linear response may significantly excite the nonlinear part of the system.

Conclusion and perspectives
We proposed and analysed a method, the pseudomodal 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 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.

Conflict of Interest
The authors declare that they have no conflict of interest.

A Well-posedness and definitions
A.1 General setting The class of systems (1-2) is considered for the general following setting: -T denotes the time set R + , -U and X are Banach spaces on the field R, -L(U, X) and L(X) are the sets of bounded linear operators from U to X, and from X to X, respectively, -ML k (X) (k ≥ 2) is the set of bounded multilinear operators from X × · · · × X k to X, equipped with norm

A.2 Mild solutions
For the linearized version of systems (1-2), mild solutions are defined in the following way [7,23].
Definition 1 (Mild solution of a linear system) Let u ∈ L ∞ loc (T, U). The mild solution of the linearized system is the function x ∈ C 0 (T, X) defined for all t ∈ T as A similar definition is given for nonlinear systems.

Definition 2 (Mild solution of a nonlinear system) Let
Then, x is said to be a mild solution of (1-2) iff x ∈ C 0 (T, X) and satisfies, for all t ∈ T, As a standard result from [23], for a fixed u ∈ C 0 (T, U), there exist t max ∈ (0, T ] and a unique function x ∈ C 0 ([0, t max ), 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 L ∞ loc (T, U), as stated in Definition 2.

A.3 Functional setting for the linearized beam problem
In the linear problem 8 (24), the bi-Laplacian B is defined as the unbounded operator on H = L 2 (0, 1) with domain 9 D(B) = w ∈ H 4 (0, 1) s.t. w(0) = w(1) = 0, w (0) = w (1) = 0 , such that B(w) = w (4) for all w ∈ D(B). This operator is closed, densely defined, self-adjoint and positive on H. Hence, we can introduce its uniquely defined positive square root K, with domain D(K) = w ∈ H 2 (0, 1) s.t. w(0) = w(1) = 0 , such that K(w) = −Δw. The domain D(K) equipped with the K-norm defines a Hilbert space 10 : Then, according to [16], the linearized problem (24) admits the state-space representation (25)(26)(27), which is well-posed for the functional setting U and X introduced in (28)(29)  In this setting, A generates a C 0 contraction semigroup on X ([16, (A1-A2), p. 6]), so that there exists a negative growth bound α < 0. More precisely (see [16, corollary 5.2]), it is a Riesz spectral operator on X which generates an analytic semigroup S, provided that − 1 b is not in the point spectrum of A. In addition, operator B belongs L(U, X) and B = 1. Therefore, this linear problem is well-posed and is in the class of systems defined in Sect. 2.1.

C Proof of Proposition 1 (truncation error estimate)
We assume (A1-A3) in Sect. 5.3 and set Φ(z) = ∞ m=1 φ m z m where the sequence (φ m ) m∈N * is defined by (see step 6 in page 4) φ 1 = 1 and, for all m ≥ 2, φ m = γ a 3 with M 3 m = p ∈ (N * ) 3 p 1 + p 2 + p 3 = m . Then, the trajectories x = ∞ m=1 x m and x = ∞ m=1 x m are normally convergent and such that, for all m ≥ 1, where the sequence (β m ) m∈N * defined by β m = (1 + ) m φ m is such that Now, for all m ≥ 1, denote the error terms e m := x m − x m . The first terms are such that e 1 X = x 1 X (from (A2)), e 2 X = 0. Moreover, for m ≥ 3, it follows from (9) that where we set F p (τ ) = A 3 x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) − A 3 x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) . It should be noted that (95) crucially depends on Eq. (53), from which in X, A 3 is identical to A 3 . Now, replacing x p i by x p i + e p i and exploiting the multi-linearity of A 3 yield the following expansion (omitting variable τ for sake of legibility) F p = A 3 x p 1 , x p 2 , e p 3 + A 3 x p 1 , e p 2 , x p 3 + A 3 e p 1 , x p 2 , x p 3 + A 3 x p 1 , e p 2 , e p 3 + A 3 e p 1 , x p 2 , e p 3 + A 3 e p 1 , e p 2 , x p 3 + A 3 e p 1 , e p 2 , e p 3 .
Then, introducing ψ m = β m − φ m for all m ≥ 1, we prove by induction that (claim C m ) e m X ≤ ψ m x 1 m X : (m = 1): the claim (C 1 ) is true for m = 1 by construction; (m ≥ 2):assume that C p holds for all p ≤ m − 1, then using expressions above, it follows that so that C m is satisfied. An immediate consequence is that which concludes the proof.

D Proof of generalized result
We assume (A1-A3) in Sect. 5.3 and consider a beam model with a third order nonlinearity A 3 for which X is not invariant. We assume that a bound a 3 of A 3 X is available, and define Φ(z) = ∞ m=1 φ m z m as in "Appendix C".
We denote A 3 = Π X A 3 . Since Π X is an orthogonal projection, A 3 X ≤ A 3 X ≤ a 3 , and therefore, as in "Appendix C", the trajectories x = ∞ m=1 x m and x = ∞ m=1 x m are convergent and satisfy, for all m ≥ 1, x m X ≤ φ m x 1