Thomas Hélie
email: thomas.helie@ircam.fr

Béatrice Laroche
email: beatrice.laroche@inrae.fr

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

Keywords: beam, partial differential equation, nonlinear systems, Volterra series, modal truncation, convergence, error bounds

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 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 subclass. 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 [START_REF] Géradin | Théorie des vibrations[END_REF][START_REF] Graff | Wave motion in elastic solids[END_REF], 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 [START_REF] Liang | Proper orthogonal decomposition and its applicationspart i: Theory[END_REF], [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF]) coupled with Galerkin projection are very popular (see for instance [START_REF] Kappagantu | Part 1: Dynamical characterization of a frictionally excited beam[END_REF], [START_REF] Kappagantu | Part 2: Proper orthogonal modal modeling of a frictionally excited beam[END_REF] and more recently [START_REF] Hassanabadi | Recursive principal component analysis for model order reduction with application in nonlinear bayesian filtering[END_REF]). Koopman based methods (see for instance [START_REF] Peitz | Koopman operator-based model reduction for switched-system control of pdes[END_REF], [START_REF] Chen | Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses[END_REF]) 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 [START_REF] Volterra | Theory of Functionnals and of Integral and Integro-Differential Equations[END_REF], 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 [START_REF] Gilbert | Functional expansions for the response of nonlinear differential systems[END_REF][START_REF] Fliess | An algebraic approach to nonlinear functional expansions[END_REF][START_REF] Rugh | Nonlinear System Theory, The Volterra/Wiener approach[END_REF][START_REF] Crouch | The observation space and realizations of finite Volterra series[END_REF][START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF]. Volterra series expansions can also be applied to some classes of nonlinear PDEs, using space-time kernels (see e.g. [START_REF] Roze | Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string[END_REF] for Green-Volterra series expansions) or using infinite-dimensional (semigroup) Volterra kernels [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF] (see also [START_REF] Cheng | Volterraseries-based nonlinear system modeling and its engineering applications: A state-of-the-art review[END_REF] 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 [START_REF] Zhu | A new convergence analysis for the volterra series representation of nonlinear systems[END_REF] and references therein (see also [START_REF] Chatterjee | Convergence analysis of Volterra series response of nonlinear systems subjected to harmonic excitation[END_REF][START_REF] Peng | On the convergence of the volterraseries representation of the duffing's oscillators subjected to harmonic excitations[END_REF][START_REF] Li | Analysis of nonlinear oscillators using volterra series in the frequency domain[END_REF] 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 [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF] 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 [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF] 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.

Preamble

In [14,[START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF], 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 by ẋ(t) = A x(t) + B u(t) + K k=2 A k x(t), . . . , x(t)

(multilinear) k , (1)
with initial condition

x(0) = x ini ∈ X, (2)
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 1 Every eigenvalue of A has strictly negative real part.

flow S(t) = exp(At) on [0, T ], whose operator norm is bounded by β exp(αt) for some β > 0 and a negative growth bound2 α < 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

u ∈ U = L ∞ (T, U) and x ∈ X = L ∞ (T, X), (3)
case for which well-posed definitions of solutions (afternamed mild solutions) can be introduced (see details in appendix A, see also [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF][START_REF] Curtain | An introduction to Infinitedimensional Linear Systems Theory[END_REF] and [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF][START_REF] Hélie | Modélisation physique d'instruments de musique et de la voix : systèmes dynamiques, problèmes directs et inverses[END_REF] 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

x(t) = +∞ m=1 x m (t), (4)
where the terms satisfy the sequence of linear problems

ẋ1 (t) = Ax 1 (t) + B u(t) with x 1 (0) = x ini , (5)
ẋm (t) = Ax m (t) + χ m (t) with x m (0) = 0, for m ≥ 2, ( 6
) with χ m (t) = min(m,K) k=2 p∈M k m A k x p1 (t), . . . , x p k (t) , (7)
where the index set

M k m := p ∈ (N * ) k p 1 + • • • + p k = m selects occurrences x pi 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 ini + t 0 S(t -τ ) B u(τ ) dτ, (8)
x m (t) = t 0 S(t -τ ) χ m (τ ) dτ, for m ≥ 2, (9)
that exactly generate a Volterra series expansion (see [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF]Rk1]), whose kernels can be deduced from formula [START_REF] Fliess | An algebraic approach to nonlinear functional expansions[END_REF][START_REF] Géradin | Théorie des vibrations[END_REF].

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 [START_REF] Hélie | Modélisation physique d'instruments de musique et de la voix : systèmes dynamiques, problèmes directs et inverses[END_REF] 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 = ε m xm ).
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 [14,[START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF][START_REF] Cheng | Volterraseries-based nonlinear system modeling and its engineering applications: A state-of-the-art review[END_REF]

): if x 1 X < ρ, then x X < Φ x 1 X , (convergence in norm)(10) x- M m=1 x m X < R M Φ x 1 X , ( error
bound) [START_REF] Graff | Wave motion in elastic solids[END_REF] 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, ζ k ≥ T S(t)A k ML k (X,X) dt, (12)
Step 2: Define the function F , analytic at z = 0,

F (z) = z z -ζ(z) with ζ(z) = K k=2 ζ k z k . ( 13
)
Step 3: Compute the unique root σ defined by

σ > 0 such that σ F (σ) -F (σ) = 0. ( 14
)
Step 4: Compute the convergence radius ρ, given by

ρ = σ/F (σ) > 0. ( 15
)
Step 5: Find the unique function analytic at z = 0

Φ(z) = ∞ m=1 φ m z m such that Φ(z) = z F Φ(z) . ( 16
)
Step 6: Compute the remainder function given by

R M Φ(z) = Φ(z) - M m=1 φ m z m . ( 17
)
Note that the coefficients φ m are given by φ 1 = 1 and, for all m ≥ 2,

φ m = m k=2 ζ k p∈M k m k i=1 φ pi .
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 ≤ γ u U with γ ≥ T S(t)B L(U,X) dt. (18)
So, a sufficient condition for the series convergence is

u U ≤ ρ u with ρ u = ρ/γ. ( 19
)
It is convenient in practice, since computing x 1 is not required. But it is conservative compared to x 1 X < ρ.

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

Nonlinear boundary problem: governing equations

The considered Euler-Bernoulli model assumes the following hypotheses (see [START_REF] Géradin | Théorie des vibrations[END_REF][START_REF] Graff | Wave motion in elastic solids[END_REF]):

(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) ∈ [0, 1] × T → R is given by (T = R + ), ∂ 2 t w+2 (a+b∂ 4 z ) ∂ t w+∂ 4 z w-η 1 0 ∂ z w 2 dz ∂ 2 z w = f, (20)
for all (z, t) ∈ [0, 1] × T, with pinned-type boundary conditions (fixed extremities and no momentum)

w(z, t) = 0 and ∂ 2 z w(z, t) = 0 for all (z, t) ∈ {0, 1} × T, (21)
and zero initial conditions

w(z, 0) = 0 and ∂ t w(z, 0) = 0 for all z ∈ [0, 1]. ( 22
)
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 [START_REF] Nayfeh | Nonlinear Oscillations[END_REF].

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

H := L 2 (0, 1) (square integrable functions), (23)
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 terms3 , 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

∂ 2 t w + 2 aI + bB (∂ t w) + B(w) = f, (24)
for all (z, t) ∈ [0, 1]×T, with zero initial conditions [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], where B is denotes the bi-Laplacian (B(w) = w (4) ), operating on the domain4 of sufficiently regular functions that satisfy the boundary conditions [START_REF] Liang | Proper orthogonal decomposition and its applicationspart i: Theory[END_REF]. This problem admits the state-space representation ẋ = Ax + Bu, for all t ∈ T, with x(0) = 0, [START_REF] Peng | On the convergence of the volterraseries representation of the duffing's oscillators subjected to harmonic excitations[END_REF] where the state x, input u, linear operators A and B are

x(t) = w(•, t) ∂ t w(•, t) , u(t) = f (•, t), (26)
Ax = 0 I -B -2 (a I + b B) x, Bu = 0 u . ( 27
)
This formulation is well-posed in the spaces of bounded signals [START_REF] Chen | Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses[END_REF] for u(t) and x(t) living in, respectively,

U = H (in accordance with H4), (28)
X = H 1 2 × H with x X = x 1 2 H 1 2 + x 2 2 H 1 2 , (29)
and norm x 1

H 1 2 = x (2) 1
H , where the complete functional setting is detailed in appendix A. [START_REF] Chen | Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses[END_REF]. 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 [START_REF] Peng | On the convergence of the volterraseries representation of the duffing's oscillators subjected to harmonic excitations[END_REF][START_REF] Roze | Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string[END_REF][START_REF] Rugh | Nonlinear System Theory, The Volterra/Wiener approach[END_REF] is in the class defined in section 2.1.

Nonlinear problem and Volterra series convergence

The nonlinear problem (η > 0) described by [START_REF] Li | Analysis of nonlinear oscillators using volterra series in the frequency domain[END_REF][START_REF] Liang | Proper orthogonal decomposition and its applicationspart i: Theory[END_REF] can be recast in the state-space formulation (1-2) based on definitions [START_REF] Roze | Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string[END_REF][START_REF] Rugh | Nonlinear System Theory, The Volterra/Wiener approach[END_REF] and by complementing [START_REF] Peng | On the convergence of the volterraseries representation of the duffing's oscillators subjected to harmonic excitations[END_REF] as follows

∂ t x = Ax + A 3 (x, x, x) + Bu, (30)
with multi-linear operator

A 3 (a, b, c) = 0 -η 1 0 a 1 (z) b 1 (z) dz c 1 . (31)
This operator defined from X 3 to X belongs to ML 3 (X). Indeed, its norm is bounded (proof in appendix B.1) as

A 3 ML3 ≤ a 3 := η 3 √ 10 . (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 ρ 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 [START_REF] Kappagantu | Part 2: Proper orthogonal modal modeling of a frictionally excited beam[END_REF]) and to a 3 (defined in (32)), and F is given by

F (z) = 1 1 -ζ 3 z 2 with ζ 3 := γa 3 . ( 33
)
Solving F (σ)σ F (σ) = 0 leads to σ = (3ζ 3 ) -1 2 . Finally, the convergence radius ρ = σ/F (σ) is given by

ρ = 2 3 √ 3ζ 3 , (34)
and ρ u = ρ/γ in remark 3. Expressions of γ, ζ 3 and ρ (eq. 18, 33, 34) w.r.t. physical coefficients are provided in section 4.2. Moreover (steps 5-6), solving Φ(z) = zF Φ(z) , using the Cardano's method and isolating the positive solution on [0, ρ[ which is zero at 0 yields, for all z ∈ [0, ρ[,

Φ(z) = ρ Ψ z ρ with Ψ (Z) = 3 cos π+arccos Z 3 . (35)
Its first Taylor coefficients are given by φ 2p = 0 for all p ≥ 1 and

φ 1 = 1, φ 3 = 2 2 3 3 ρ 2 , φ 5 = 2 4 3 5 ρ 4 , φ 7 = 2 8 3 8 ρ 6 , φ 9 = 2 8 ×5×11 3 12 ρ 8 (etc.), from which R M φ(z) = Φ(z) - M m=1 φ m z m , (36)
is deduced. The gain bound function Φ and its remainder R M Φ are displayed in figure 1.

Φ/ρ M = 1 M = 3 M = 5 M = 7 M = 9 z/ρ R M Φ/ρ a ⃝ b ⃝ Figure 2: a ⃝ Function D ϵ Φ/ρ where D ϵ Φ(z) = Φ z(1 + ϵ) -Φ(z)
for relative errors ϵ of 0.1%, 1% and 10%. b ⃝ functions R M Φ/ρ (M = 1, . . . , 9) and Φ/ρ.

Numerical results

Two sets of parameters are studied.

Configuration 1. We consider a single mode (n = N = 1), a fluid damping only (b = 0) tuned to be close to the critical regime (a = 0.999π 2 ). Parameter

η is chosen such that ρ = 1, which corresponds to η = 13.857. The excitation duration is θ = 3. Four amplitudes are tested, chosen such that ∥x 1 ∥ X /ρ ∈ {0.8; 1; 1.2; 2}. The truncated Volterra expansion (8-9) that solves the nonlinear beam model [START_REF] Zhu | A new convergence analysis for the volterra series representation of nonlinear systems[END_REF](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.

Modal decomposition and projection

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 b/a, almost reached when k 2 n ∼ a/b.

Assuming that λ ± n are always distinct, the associated eigenfunctions of A are

e ± n (z) = 1/λ ± n 1 s n (z) with s n (z) = √ 2 sin(k n z) and k n = nπ (n ≥ 1). ( 38
)
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

H = U := span(s 1 , . . . , s N ) ⊂ H = U, (39)
H 1 2 := span(s 1 , . . . , s N ) ⊂ H 1 2 , (40)
S i = s i 0 if 1 ≤ i ≤ N,
and

S i = 0 s i-N if N + 1 ≤ i ≤ 2N. ( 42
)
This basis is the one used in the sequel. The associated functional spaces U and X are summarized in table 1 (column 2).

5 characterized by the smallest absolute value of the real part of λ ± n .

Convergence bound estimate

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

γ = T sup n≥1 |k 4 n h n (t) 2 + h n (t) 2 | dt, ( 43
)
where the impulse responses h n of second order linear systems are given according to the damping ratio ξ n by

h n (t) =          e -k 2 n ξnt sin(k 2 n √ 1-ξ 2 n t) k 2 n √ 1-ξ 2 n if ξ n < 1, te -k 2 n t if ξ n = 1, e -k 2 n ξnt sinh(k 2 n √ ξ 2 n -1t) k 2 n √ ξ 2 n -1 if ξ n > 1.
Moreover, equations (33-34) become

ζ 3 = η 3 √ 10 γ and ρ = 2 4 √ 10 3 1 ηγ . (44)
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

ρ u = 2 4 √ 10 /(3(γ) 3/2 √ η).

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 A : x ∈ X -→ Ax ∈ X and B : u ∈ U -→ Bu ∈ X, 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 v = [V 1 , . . . , V N ] T ∈ V, s.t. u(t) = N i=1 V i (t)s i , (45)
y ini = [Y (0) 1 , . . . , Y (0) 2N ] T ∈ Y, s.t. x ini = 2N i=1 Y (0) i S i , (46)
y = [Y 1 , . . . , Y 2N ] T ∈ Y, s.t. x = 2N i=1 Y i S i , (47)
where spaces V, Y and Y are defined in table 1 (column 3) and Y is equipped with the norm built from X

y Y = 2N i=1 Y i S i X = y T Qy, (48)
with

Q = diag( S 1 2 X , . . . , S 2N 2
X ) = π 2 diag(1, . . . , N ) 2 0 N ×N 0 N ×N I N . ( 49
)
The coordinates y ∈ Y of the trajectories are governed by

ẏ(t) = A Y y(t) + B Y v(t), y(0) = y ini ∈ Y, (50)
with

A Y = 0 N×N I N -π 4 diag(1, . . . , N) 4 -2 a I N +bπ 4 diag(1, . . . , N) 4 , (51)
B Y = 0 N ×N I N , (52)
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 for all (a, b, c) in X 3 , A 3 (a, b, c) ∈ X. ( 53
)
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 [START_REF] Zhu | A new convergence analysis for the volterra series representation of nonlinear systems[END_REF] 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 problem

˙ x = A x + BΠ H u + A 3 ( x, x, x), (54)
x ini = Π X x ini ∈ X.
Note that x does not identify with x := Π X x for a solution x in general. For 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 section 4.3, the finite-dimensional problem (54) satisfied by x on X can be restated on Y for the coordinates in the basis S as

ẏ(t) = A Y y(t)+B Y v(t)+A Y3 (y, y, y), y(0) = y ini ∈ Y, (55)
where

A Y3 ∈ ML 3 (Y) is given by A Y3 (a, b, c) = -η a T π 2 diag(1, . . . , N, 0 1×N ) 2 b × 0 N ×N 0 N ×N -π 2 diag(1, . . . , N ) 2 0 N ×N c. ( 56
)
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

x 1 = x 1 = Π X (x 1 ), (57)
x m = t 0 e A(t-τ ) p∈M 3 m A 3 x p1 (τ ), x p2 (τ ), x p3 (τ ) dτ, for all m such that 2 ≤ m ≤ M. (58)
Then, we define the pseudo-modal Volterra approximation of x as

X M N = M m=1 x m . ( 59
)
Original spaces Projection sub-spaces Modal coordinates spaces

U = H U = H V = R N (equipped with 2-norm) X = H × H 1 2 X = H × H 1 2 (equipped with (29)) Y = R 2N (equipped with • Y , see (48)) U = L ∞ (T, U) U = L ∞ (T, U) V = L ∞ (T, V) X = L ∞ (T, X) X = L ∞ (T, X) Y = L ∞ (T, Y)
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

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 sections 3.3 and 4.2 to (54) and the Volterra series decomposition (57-58). This yields the estimates

γ N = T N sup n=1 | k 4 n h n (t) 2 + h n (t) 2 |dt, (60)
and

A 3 ML3( X) = A Y3 ML3(Y) ≤ a 3N with a 3N = η π 2 N n=1 1 n 4 1/2 . (61)
Then, we build

ζ 3N = γ N a 3N , ρ N = 2 3 √ 3ζ 3N , (62)
and the gain bound function (see (35))

∀z ∈ [0, ρ N ), Φ N (z) = ρ N Ψ z ρ N = ρ N ρ Φ ρ ρ N z . ( 63
)
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

r N := ρ N ρ = ζ 3 ζ 3N = a 3 a 3N γ γ N > 1, (64)
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 γ/γ 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 ρ N and gain bound function

Φ N , so that y = ∞ m=1 y m is convergent in norm if y 1 Y = < ρ N . Moreover, we have y Y = x X < Φ N x 1 X = Φ N y 1 Y .

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

(A1) u = u+u ⊥ where u = Π U u and u ⊥ = (Id-Π U )u, (A2) x ⊥ 1 X ≤ x 1 X with 1, (A3) x 1 X + x ⊥ 1 X ≤ (1 + ) x 1 X < ρ,
where x 1 and x ⊥ 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 = xx the error on x due to the pseudo-modal truncation on X. Then, e is bounded as

e X ≤ D Φ x 1 X ( 65
)
where D Φ(z)

:= Φ (1 + )z -Φ(z).
The proof is given in C. Bound estimates are given in figure 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 pseudo-modal Volterra approximation

X M N is bounded by x -X M N X ≤ x -x X + x -X M N X ≤ D Φ x 1 X + R M Φ N x 1 X . (66)
Φ/ρ M = 1 M = 3 M = 5 M = 7 M = 9 z/ρ R M Φ/ρ a ⃝ b ⃝ Figure 2: a ⃝ Function D ϵ Φ/ρ where D ϵ Φ(z) = Φ z(1 + ϵ) -Φ(z) for relative errors ϵ of 0.1%, 1% and 10%. b ⃝ functions R M Φ/ρ (M = 1, . . . , 9
) and Φ/ρ.

Numerical results

Two sets of parameters are studied.

Configuration 1. We consider a single mode (n = N = 1), a fluid damping only (b = 0) tuned to be close to the critical regime (a = 0.999π 2 ). Parameter η is chosen such that ρ = 1, which corresponds to η = 13.857. The excitation duration is θ = 3. Four amplitudes are tested, chosen such that ∥x 1 ∥ X /ρ ∈ {0.8; 1; 1.2; 2}.

The nonlinear system is simulated using a standard ODE solver (ode15s, Matlab). These trajectories are compared to those obtained from truncated versions of ( 7) by simulating a sequence of linear systems (lsim, Matlab)

given in proposition 1. Signals w(z = 0.5, t) = X 1 (z = 0.5, t) are displayed in figure 3 for several truncation orders M = 1, 3, 5, 7. For ∥x∥ X = 0.8, a good approximation is obtained as soon as M ≥ 3. For ∥x∥ X = 1, the approximations are significantly more accurate when increasing the order M . This is no longer true for ∥x∥ X = 1.2, for which the convergence seems 24

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 ( 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 section 2.1. It follows that we obtain a convergence bound ρ and a gain bound function Φ 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 ini and input u, we define its pseudo-modal truncation at order N x m as in equation (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

A 3 -A 3 X = sup (a,b,c)∈( X) 3 a X , b X , c X =1 A 3 (a, b, c) -A 3 (a, b, c) X , (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 = xx the error on x due to the pseudo-modal truncation on X. Then, e is bounded as

e X ≤ D Φ x 1 X + A 3 -A 3 X a 3 Φ x 1 X , ( 68
)
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 x -X M N X ≤ D Φ x 1 X + R M Φ N x 1 X + A 3 -A 3 X a 3 Φ x 1 X .
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

A3-A3 X a3 by max k∈{2,••• ,K} A k -Π X A k X a k .

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

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 → X only act on subspaces of X and admits a sandwich decomposition

A k (x 1 , . . . , x k ) = B k A † k Cx 1 , . . . , Cx k , (69)
where

C : X → X † , A † k : (X † ) k → W † k and B k : W † k → X lie respectively in L(X, X † ), ML k (X † , W † k ) and L(W † k , X).
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

x † = Cx, defined in X † = L ∞ (T, X † ). ( 70
)
Indeed, its Volterra series is the sum of x † m = C x m , straightforwardly derived from (8-9) and (69) as

x † 1 (t) = CS(t)x ini + t 0 S † 1 (t -τ )u(τ ) dτ, (71)
x † m (t) = min(m,K) k=2 p∈M k m t 0 S † k (t -τ )A † k x † p1 (τ ), . . . , x † p k (τ ) dτ, for m ≥ 2, ( 72
)
where we have defined, for 1 ≤ k ≤ K and

B 1 = B, S † k (t) = C S(t) B k . (73)
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 [START_REF] Hassanabadi | Recursive principal component analysis for model order reduction with application in nonlinear bayesian filtering[END_REF]. Thus, bounds ζ k are replaced in step 1 by the new estimates, for 2

≤ k ≤ K, ζ † k ≥ T S † k (t)A † k ML k (X † ,X † ) dt, ( 74
) some overestimated values of which are γ † k a † k with a † k ≥ A † k ML k (X † ,W † k ) , γ † k ≥ T S † k (t) L(W † k ,X † ) dt. ( 75
)
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:

x † X † < Φ † x † 1 X † (76) x † - M m=1 x † m X † < R M Φ † x † 1 X † , for M ≥ 1, (77)
x-

M m=1 x m X < r R M Φ † C x 1 X † , for M ≥ 1, (78) x X < x 1 X + r R 1 Φ † Cx 1 X † , (79)
where

r = max 2≤k≤K (ζ k /ζ † k ).
The proof follows exactly the same steps as in [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF].

Remark 5 (Supplement to remark 3) An alternative convergence radius on the input is

ρ † u = ρ † /γ † 1 with γ † 1 = T S † 1 L(U,X † ) 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)ξ

+ 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, |U | < ρ † / |H(iω)|.

Application to the beam model

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

w =: x † = Cx with C : w v ∈ X -→ w ∈ H 1 2 , (80)
from which the nonlinear term in (20) is built using

A † 3 (w a , w b , w c ) = -η 1 0 w a (z) w b (z) dz w c , (81)
and contributes to the state equation through (as input u)

B 3 = B : H → X. (82)
Operator A † 3 belongs to ML 3 (X † := H Then, the nonlinear systems (30) governing the beam state x and (54) governing its pseudo-modal truncation x can both benefit from the refined convergence results presented in section 6.1. The estimates replacing γ N -→ N →+∞ γ are given by (see Appendix B.2)

γ † N := T N sup n=1 k 2 n |h n (t)| dt,
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

ρ † N = 2 3 3a 3N γ † N > ρ N > ρ and ρ † = 2 3 3a 3 γ † > ρ, (83)
and, in theorem 1, factors r N (for x) and r (for x) are given by

r N = ζ 3N ζ † 3N = γ N γ † N > 1 and r = ζ 3 ζ † 3 = γ γ † 3 > 1.
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

x † -C X M N X † ≤ E M N and x -X M N X ≤ r E M N ,( 84
)
with E M N = D Φ † x † 1 X † + R M Φ † N x † 1 X † .
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

ρ † 1 = 2 3 3 2 |α 3 | h 1 . ( 85
)
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 ( [START_REF] Li | Analysis of nonlinear oscillators using Volterra series in the frequency domain Part I : convergence limits[END_REF], [START_REF] Li | Analysis of nonlinear oscillators using volterra series in the frequency domain[END_REF]), established for a Duffing oscillator, with impulse response h. In the harmonic regime, we obtain a bound

ρ † 1 /|H(iω)| = 2/ 3 3 2 |α 3 | h 1 |H(iω)| .
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.

7 Numerical experiments

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 timedomain 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 [START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF]) for the more general class (1-2)), in which the finite dimensional linear systems (W in grey) are simulated using a standard linear 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 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.

equations, Volterra series expansions, convergence domain Email addresses: thomas.helie@ircam.fr (Thomas Hélie), blaroche@jouy.inra.fr (Béatrice Laroche)

1. Introduction v y 1 y 1 y 3 y 3 y 5 y 5 y 7 y trunc x trunc 2M 3 2M 3 2M 3 M 3 M 3 M 3 M 3 M 3 C W Q W W W a b c M 3 (

Preprint submitted to Elsevier

June 20, 2019 WB R 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 Y3 =: M 3 is symmetric with respect to the two first variables (A Y3 (a, b, c) = A Y3 (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 W B respectively denote the simulation of the semigroup generated by matrix A Y and the one with input matrix B Y . The space-time function xtrunc (z, t) = R(z) y trunc (t) is built by aggregating the activated eigenfunctions as in formula (47).

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

η = 4π 2 27γ † 1 (config. 1) and η = 4 √ 10 9γ † (config. 2). ( 86
)
Moreover, in order to excite the beam by a relevant rich spectral content, we choose inputs u : t → z → f (z, t) in (30) with step time-shapes

f (z, t) = f tot g(z) 1 [0,Te] (t), (87)
where the total force f tot is spatially distributed according to g > 0 (where 1 0 g(z)dz = 1). As mentioned in section 7.1, we simulate the pseudomodal Volterra terms xm 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π 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 response x † 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). Table 3: Summary of the main physical parameters of the configurations 1 and 2 used for numerical experiments.

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, centered at z c = 1/7 described by g(z) = G cos(π(zz 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,

g n = cos(nπδ/2) 1-(nδ) 2 e n (z c ) if nδ = 1,

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) = [ x trunc z = 0.5, t ] 1 ), for truncation orders M = 1, 3, 5, 7. Results are displayed in figure 5.

For x † 1 X † = 0.8, a good approximation is obtained as soon as M ≥ 3. As guaranteed by [START_REF] Gilbert | Functional expansions for the response of nonlinear differential systems[END_REF], 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 figure 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. In figure 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 [START_REF] Li | Analysis of nonlinear oscillators using volterra series in the frequency domain[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF] and [START_REF] Zhu | A new convergence analysis for the volterra series representation of nonlinear systems[END_REF]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 ( x †

1 X † = 3 > ρ † 2 > ρ † 16 > ρ † = 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 ( 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 x † x † 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 [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. 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 = 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.

1 X † = 0.8 x † 1 X † = 1 0 1 2 3 4 5
1 X † = 1.2 x † 1 X † = 2
Finally, figure 11 illustrates the effect of the fluid damping coefficient decrease: the parameters of config- 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 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 section 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 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:

• 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 oper- ators from X × • • • × X k to X, equipped with norm A k = sup (x 1 ,...,x k )∈X k x 1 =•••= x k =1 A k (x 1 , . . . , x k ) X ,
• U = L ∞ (T, U) and X = L ∞ (T, 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 [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF][START_REF] Curtain | An introduction to Infinitedimensional Linear Systems Theory[END_REF].

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

x(t) = S(t)x ini + t 0 S(t -τ )Bu(τ ) dτ.
A similar definition is given for nonlinear systems.

Definition 2 (Mild solution of a nonlinear system) Let u ∈ L ∞ loc (T, U). Then, x is said to be a mild solution of (1-2) iff x ∈ C 0 (T, X) and satisfies, for all t ∈ T,

x(t) = S(t)x ini + t 0 S(t -τ )× Bu(τ ) + P x(τ ) + Q x(τ ), u(τ ) dτ. ( 89
)
As a standard result from [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF], 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 problem8 (24), the bi-Laplacian B is defined as the unbounded operator on H = L 2 (0, 1) with domain9 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 :

D(B) = w ∈ H 4 (0, 1) s.t. w(0) = w(1) = 0, w (0) = w (1) = 0 ,
H 1 2 := D(K) with . H 1 2 = K . H .
Then, according to [START_REF] Jacob | Analyticity and Riesz basis property of semigroups associated to damped vibrations[END_REF], the linearized problem (24) admits the state-space representation [START_REF] Peng | On the convergence of the volterraseries representation of the duffing's oscillators subjected to harmonic excitations[END_REF][START_REF] Roze | Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string[END_REF][START_REF] Rugh | Nonlinear System Theory, The Volterra/Wiener approach[END_REF], which is well-posed for the functional setting U and X introduced in (28-29) and defining operator A on domain D(A) as

A : D(A) -→ X with D(A) = {(u, v) ∈ H 1 2 × H 1 2 s.t. 2(a + b B)v + Bu ∈ H}.
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 section 2.1. 10 This justifies the notation introduced in (29) where

x 1 H 1 2 = K x 1 H = x (2) 1 H .

B Estimates for the nonlinear beam model

B.1 Proof of equation (32)

For all X 1 ∈ H

1 2 and z ∈ [0, 1], X 1 (z) = X 1 (0) z + z 0 (z - ζ)X 1 (ζ) dζ which is zero at z = 1. Hence, X 1 (z) = 1 0 K(z, ζ)X 1 (ζ) dζ, where K(z, ζ) = ζ(z -1) if 0 ≤ ζ ≤ z and K(z, ζ) = z(ζ -1) if z ≤ ζ ≤ 1. Using the Cauchy-Schwartz inequal- ity, it comes X 1 (z) ≤ µ(z) 1 0 X 1 (ζ) 2 dζ 1 2
, where µ(z)

= 1 0 K(z, ζ) 2 dζ 1 2 = z(1-z) √ 3 . Then, 1 0 X 1 (z) 2 dz 1 2 ≤ ν X 1 2 H 1 2
, where ν = On each modal subspace, the linearized system behaves like second order system. Indeed, for all x ∈ H

1 0 µ(z) 2 dz 1 2 = 1 3 √ 10 . Finally, for all (X, Y, Z) ∈ X 3 , A 3 (X, Y, Z) X ≤ η 1 0 X 1 (z) Y 1 (z) dz 1 0 Z 1 (z) dz ≤ η 1 0 X 1 (z) Y 1 (z) dz 1 0 Z 1 (z) 2 dz 1 2 ≤ η 1 0 X 1 (z) 2 dz 1 2 Y 1 H 1 2 Z 1 H 1 2 ≤ ν η X 1 H 1 2 Y 1 H 1 2 Z 1 H 1 2 ≤ ν η X X Y X Z X ,
S(t)Bx = ∞ n=1 κ - n S(t)E - n (90)
where κ - n =< Bx, E - n > X =< x, s n > H , and

E + n (z) = (e + n -e - n )λ + n λ - n k 2 n (λ - n -λ + n ) = 1 k 2 n 1 0 s n (z) E - n (z) = (λ + n e + n -λ - n e - n ) λ + n -λ - n = 0 1 s n (z) (91)
From (91) we obtain Now we notice that for all (x, y, z) in X, A 3 (x, y, z) = Bg(x, y, z) where g is a third order multilinear operator from X 3 to H, such that g ML(X,H) = a 3 . We therefore obtain that for all (a, b, c) in X We assume (A1-A3) in section 5.

S
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 equation (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 multilinearity of A 3 yield the following expansion (omitting variable τ for sake of legibility)

F p = A 3 x p 1 ,
≤ (β m -φ m ) x 1 m X = ψ m x 1 m X ,
so that C m is satisfied. An immediate consequence is that

e X ≤ ∞ m=1 ψ m x 1 m X = ∞ m=1 β m x 1 m X - ∞ m=1 φ m x 1 m X ≤ Φ ((1 + ) x 1 X ) -Φ( x 1 X ),
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 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 m X , x m X ≤ φ m x 1 m X ≤ φ m ((1 + ) x 1 X ) m ≤ β m x 1 m X ,
where the sequence (β m ) m∈N * was defined in appendix C. Now, for all m ≥ 1, denote the error terms e m := x mx 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 e m (t) = 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 (τ ) and G p (τ ) = (A 3 -A 3 ) x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) Setting

G p X ≤ A 3 -A 3 X x 1 m X φ p 1 φ p 2 φ p 3
Then, following the same steps as in appendix C, we introduce ψ m = β mφ m , and obtain by induction that for all m ≥ 1, e m X ≤ (ψ m + A

3 -A 3 X a 3 φ m ) x 1 m X . An immediate consequence is that e X ≤ ∞ m=1 (ψ m + A 3 -A 3 X a 3 φ m ) x 1 m X ≤ ∞ m=1 β m x 1 m X -(1 - A 3 -A 3 X a 3 ) ∞ m=1 φ m x 1 m X ≤ Φ ((1 + ) x 1 X ) -(1 - A 3 -A 3 X a 3 )Φ( x 1 X ),
which concludes the proof.

Fig. 1 :

 1 Fig. 1: Normalized gain bound function Φ/ρ and normalized remainder functions R M Φ/ρ.

Fig. 2 :

 2 Fig. 2: Function D Φ/ρ (= D Ψ ) for relative errors of 0.1%, 1% and 10% (see (35) and (65)).

1 2 , W † 3 :† 3 on H 1 2 1 2

 1311 = H) and has the same norm value as A 3 . Similarly, its restriction A , which generates A 3 = C A † 3 (B•, B•, B•), belongs to ML 3 ( H , H) and has the same norm value as A 3 .

Figure 1 :

 1 Figure 1: Zoom for the case x X = 1 in figure ??.

Configuration 1 (

 1 non oscillating) 2 (oscillating+hammer) Parameters: a, b a = 0.999π 2 and b = 0 a = 0.1 and b = 10 -6 Impulse response bound γ † 1 ≈ 0.1013 γ † ≈ 8.7749 η (see 86) η ≈ 14.434 η ≈ 0.1602 Input (simulation on [0, 2T e ]): N = 1 and T e = 3 N ∈ {2, 16} and T e = 40 f (z, t) = f tot 1 [0,T e ] (t) N n=1 g n e n (z) g 1 = 1 g n : see (88) and fig. 4a

π 4 e

 4 n (z c ) otherwise. (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 ρ † 2 ≈ 1.0498 and (b) N = 16 modes for which ρ † 16 ≈ 1 + 2.64 × 10 -5 ≈ ρ † . Three amplitudes are tested. They are chosen such that x † 1 X † = A with A ∈ {1; 3; 5} (still corresponding to f tot = A/ x † 1 X † where computations yield f tot ≈ 8.2002×A for N = 2 and f tot ≈ 7.9533×A for N = 16). The deflection wave is observed at the center of the distributed force (z = z c = 1/7).

Configuration 2 (

 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.

Figure 4 :Configuration 1 ( 14 Fig. 4 :

 41144 Figure 4: Configuration 2. (a) First coe cients g n . (b) Reference unitary spatial distribution g(z) (-) and its modal reconstructions with N = 2 modes (• • • ) and N = 16 modes ( ).

Fig. 5 :

 5 Fig. 5: Configuration 1. Signals of the deflection wave at center.

Fig. 6 :

 6 Fig. 6: Configuration 1, x † 1 X † = 1. Zoom of figure 5.

Figure 7 :

 7 Figure 7: (Configuration 2, ky 1 k X = 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.

17

 17

Fig. 7 :

 7 Fig. 7: (Configuration 2, x † 1 X † = 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.

Figure 8 :Fig. 8 :

 88 Figure 8: (Configuration 2, ky 1 k X = 3) To be compared to fig. 7.18

Fig. 9 :

 9 Fig. 9: (Configuration 2 x † 1 X † = 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.

Figure 9 : 19 Fig. 10 :

 91910 Figure 9: (Configuration 2, ky 1 k X = 5) To be compared to fig. 7 and 8.19Fig. 10: (Configuration 2, x † 1 X † = 5) To be compared to fig. 7 and 8.

Fig. 1 :Fig. 11 :

 111 Fig. 1: (Configuration 2 with smaller damping a = 0.05, 0.02, 0.01) 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.

t 0 S 0 S

 00 (tτ )Bg(a, b, c)(τ )dτ X ≤ t (t)Bx L(H,X) a 3 dτ a X b X c X ≤ γa 3 a X b X c X (94)This proves that ζ 3 ≤ γa 3 .In the same way, for operator C, we obtain from (92) thatCS(t)Bx = ∞ n=1 κ - n e λ + n te λ - n t

3 m 3 m β p 1 β p 2 β p 3 . 0 S(t -τ ) p∈M 3 mA 3 x 0 S(t -τ ) p∈M 3 mF

 33303303 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 p∈M φ p 1 φ p 2 φ p 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=1x m are normally convergent and such that, for all m ≥ 1,x m X ≤ φ m x 1 m X , x m X ≤ φ m x 1 m X ≤ φ m ((1 + ) x 1 X ) m ≤ β m x 1 m X ,where the sequence (β m ) m∈N * defined by β m = (1 + ) m φ m is such thatβ 1 = 1 + and for all m ≥ 2, β m = γa 3 p∈MNow, for all m ≥ 1, denote the error terms e m := x mx 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 e m (t) = t p 1 (τ ), x p 2 (τ ), x p 3 (τ )-A 3 x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) dτ = t p (τ )dτ,

t 0 S(t -τ ) p∈M 3 mA 3 0 S(t -τ ) p∈M 3 mF

 03303 x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) -A 3 x p 1 (τ ), x p 2 (τ ), x p 3 (τ ) dτ = t

A 3 -

 3 A 3 X = sup (a,b,c)∈( X) 3 a X , b X , c X =1 A 3 (a, b, c) -A 3 (a, b, c) X ,we have

Table 2 :

 2 Coefficient a 3 /a 3N -1 for 1 ≤ N ≤ 8 (lower rounded values in per cent).

	N	1	2	3	4	5	6	7	8
	10 2 ( a 3 /a 3N -1)	1.99 0.463 0.173 0.0826 0.0455 0.0277 0.0181 0.0124

  (t)Bx =

	∞ n=1	κ -n	∞ n=1 k 2 n e λ + κ -n n t -e λ -λ + n e λ + n t λ + n -λ -n n t λ + n -λ -n E + n + e + n -λ + n e λ + λ -n e λ -n t λ + n n -λ -n t -λ -e -n n e λ -n t = λ + n n -λ -E -n	(92)
	Setting h n (t) = e λ + n t -e λ -n t λ + n -λ -n sions	and using the standard expres-
	h n (t) =	        	e -k 2 n ξ n t sin(k 2 n k 2 n √ √ 1-ξ 2 1-ξ 2 n t) n te -k 2 n t e -k 2 n ξ n t sinh(k 2 n √ ξ 2 n -1t) k 2 n √ ξ 2 n -1	if ξ n < 1 if ξ n = 1 if ξ n > 1
	we obtain that
				∞
	S(t)Bx X =	(κ -
				n=1

n ) 2 (k 4 n h n (t) 2 + h n (t) 2 ) ≤ sup n |k 4 n h n (t) 2 + h n (t) 2 | x H so that we can set γ = T sup n |k 4 n h n (t) 2 + h n (t) 2 |dt (93)

  T sup n |k 2 n h n (t) 2 + h n (t) 2 |dt T sup n (k 2 n |h n (t)|) dt

		and
						∞
		CS(t)Bx	L(H,H	1 2 )	=	(κ -n ) 2 k 4 n h n (t) 2
						n=1
						≤ sup n	k 2 n |h n (t)| x H
		It follows that CS(t)Bx we can set	L(H,H	1 2 )	≤ sup n k 2 n |h n (t)|, so that
		γ † =	T	sup n	k 2 n |h n (t)| dt
		r = C Proof of proposition 1 (truncation error
		estimate)
	λ + n -λ -n	s n ,

  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 .Thomas Hélie, Béatrice LarocheThen, 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 thate m X ≤ γa 3 p 1 φ p 2 ψ p 3 + φ p 1 ψ p 2 φ p 3 + ψ p 1 φ p 2 φ p 3 + φ p 1 ψ p 2 ψ p 3 + ψ p 1 φ p 2 ψ p 3 + ψ p 1 ψ p 2 φ p 3 + ψ p 1 ψ p 2 ψ p 3 ≤ γa 3 x 1 β p 1 β p 2 β p 3φ p 1 φ p 2 φ p 3

	x 1	(p 1 +p 2 +p 3 ) X
	p∈M 3 m
	m
	X
	p∈M 3 m

φ

The maximal real part of the eigenvalues of A is the smallest optimal bound α.

few oscillatory signals, sweeps or sequences of impulses...

see appendix A.3 for the detailed construction of the functional setting.

This case involves some so-called secular modes due to internal and undamped resonances.

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.

For recall, ∂ 2 t w +2 aI +bB (∂ t w)+B(w) = f on [0, 1]×T.

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

Conflict of Interest

The authors declare that they have no conflict of interest.


Sahana Bhattacharjee
email: sahana.bhattacharjee@hotmail.com

Subrata Chakraborty
email: subratastats@dibru.ac.in

INFLATED LOG LINDLEY DISTRIBUTION

Keywords: Inflated beta Distribution, Inflated Kumaraswamy Distribution, Inflated Unit Lindley Distribution, Zero-one inflated data MSC 2010: 60E05, 62E15

This paper introduces the Zero and Zero-One Inflated versions of the Log Lindley distribution. Important distributional properties are investigated. Special case of both the proposed distributions are shown to be members of the exponential family. Moment methods and maximum likelihood estimation is discussed in detail with derivation of the exact expression for the information matrix.

Introduction

Many researchers in the field of applied statistics often come across data on standardized rates, proportions or fractions, which assume values in the range (0, 1) . For instance, data on proportion of successful (or unsuccessful candidates in an examination), fraction of an hour required to accomplish a task by a group of workers etc. However, it is not rare that data on proportions or fractions contain a number of zeroes and/or ones, i.e., the data arising in the range [0, 1) , (0, 1] or [0, 1] . In such cases, the commonly known distributions on (0, 1) such as the Beta distribution [START_REF] Johnson | Continuous univariate distributions[END_REF], Kumaraswamy distribution [START_REF] Kumaraswamy | A generalized probability density function for double-bounded random processes[END_REF], Unit Lindley distribution [START_REF] Mazucheli | On the one parameter unit-Lindley distribution and its associated regression model for proportion data[END_REF] fail to model the data in question as all of them allow modeling of values lying strictly between 0 and 1. This calls upon the requirement to develop models which is able to capture the 1 probability mass concentrated at the points 0, 1 or both, as the case may be. One of the most popular approaches is to mix a continuous distribution in the range (0, 1) with a degenerate distribution whose probability mass is concentrated at either 0 or 1 (in case of data observed in the range [0, 1) or (0, 1]) and that with the Bernoulli distribution which assigns non-negative probability to both 0 and 1. These models are members of the inflated family of distributions. The case when a significantly high number of zeroes or ones in the data arises, it is specifically referred to as zero-inflated or oneinflated data. Zero inflated data rise very commonly in nature and can be modelled suitably using a number of continuous distributions (See [START_REF] Burch | Zero-inflated beta distribution applied to word frequency and lexical dispersion in corpus linguistics[END_REF], [START_REF] Hashimoto | Zero-spiked regression models generated by gamma random variables with application in the resin oil production[END_REF], [START_REF] Liu | Statistical analysis of zero-inflated nonnegative continuous data: A review[END_REF] and [START_REF] Tomazella | Zeroadjusted reparameterized Birnbaum-Saunders regression model[END_REF]).The data where a high number of both zeroes and ones arise is referred to as zero-one inflated data. The first work on development of above-mentioned models was reported in [START_REF] Ospina | Inflated beta distributions[END_REF], where the authors developed the zero-and/or inflated versions of the Beta distribution and showed its utility through a real-life data set on percentage of qualified nurses in 645 Brazilian municipal districts. [START_REF] Cribari-Neto | Inflated Kumaraswamy distributions[END_REF] introduced the inflated Kumaraswamy distribution and showed its advantage over the inflated Beta distribution through modeling of a data set on proportion of inhabitants in each of the 5,566 Brazilian municipalities that lived in homes with at least one bathroom and piped water in 2010. The inflated Unit Lindley distribution was recently developed by [START_REF] Chakraborty | Modeling proportion of success in high school leaving certificate examination-a comparative study of inflated unit Lindley and inflated beta distribution[END_REF], who studied their various statistical properties and it proved better in modeling data on proportion of successful students in high school leaving certificate examination of Manipur state of India as compared to the inflated Beta distribution. The Log Lindley distribution was proposed by [START_REF] Goméz-Deniz | The Log-Lindley distribution as an alternative to the beta regression model with applications in insurance[END_REF] as an alternative to the classical beta probability model having support in the unit interval. One of the main advantages of this distribution is that the cdf is obtainable in a closed form, unlike the Beta distribution and also, the expression for pdf and cdf does not involve any special function. It also presented an alternative to the Beta regression model which models response variable distributed in the range (0, 1) and allows more flexibility of the covariates. So far, no work has been attempted on developing the inflated version of the Log Lindley distribution, which may prove to be a viable alternative of the inflated Beta distribution. Also, only a handful of zero-or/and-one inflated distributions are available in the literature. Hence, this paper attempts to develop the zero-or/and-one inflated Log Lindley distribution. Section 2 introduces the zero-or-one inflated and zero-and-one inflated Log Lindley distribution and explores some of the distributional properties of both the distributions. Section 3 deals with the parameter estimation of both the zero-or-one and zero-and-one inflated Log Lindley distribution. Final section provides concluding remarks.

f (y; σ, λ) = σ 2 1 + λσ (λ -log y) y σ-1 ; 0 < y < 1 (1)
For data arising from some real life phenomena, a number of zeroes and/or ones may be present, which is not then suitable to be modeled by any distribution having support on (0, 1). Hence, it becomes necessary to incorporate a discrete component into the continuous data generating process so that the zeroes and/or ones are observed with a positive probability. Consequently, we need to look for a probability law which will merge both the continuous and discrete data generating processes. One way of achieving this is by mixing two distributions -the continuous Log Lindley distribution on (0, 1) and the degenerate distribution with the entire probability mass concentrated at the point c, where either c = 0 or c = 1. We then say that the data is inflated at one of two the end points of the interval (0, 1). We name this mixture distribution as the inflated Log Lindley (ILL) distribution, whose cumulative distrbution function (c.d.f.) is given by

ILL c (y; α, σ, λ) = αI [c,1] (y) + (1 -α) F (y; σ, λ) (2)
where I A is the indicator function which assumes the value 0 if y ∈ A and 1 if y / ∈ A, 0 < α < 1 is the mixture parameter and F (; σ, λ) is the c.d.f. of the Log Lindley distribution with parameters σ and λ. The random variable Y follows the Log Lindley distribution with parameters σ and λ with probability (1 -α) and follows the degenerate distribution at c with probability α. Let Y be a r.v. with c.d.f. given by equation [START_REF] Chakraborty | Modeling proportion of success in high school leaving certificate examination-a comparative study of inflated unit Lindley and inflated beta distribution[END_REF]. Then the p.d.f. of Y is given by

iLL c (y; α, θ) = α, if y = c (1 -α) f (y; σ, λ) , if y ∈ (0, 1) (3)
where f (; σ, λ) is the Log Lindley density given in equation ( 1) and α ∈ (0, 1) is the probability mass at c, which shows the probability of observing 0 (when c = 0) or 1 (when c = 1). We denote Y ∼ ILL c (α, σ, λ).

Definition 2.1 Let Y ∼ ILL c (α, σ, λ).
1. If c = 0, the distribution in ( 3) is called the zero-inflated Log Lindley distribution (LLZI) and we specifically write Y ∼ LLZI (α, σ, λ), where α = P (Y = 0).

2. If c = 1, the distribution in ( 3) is called the one-inflated unit Lindley distribution (LLOI) and we specifically write Y ∼ LLOI (α, σ, λ), where α = P (Y = 1).

Figure [START_REF] Burch | Zero-inflated beta distribution applied to word frequency and lexical dispersion in corpus linguistics[END_REF] displays the LLZI and LLOI densities inflated at the points c = 0 and c = 1 for different combinations of the values of σ and λ, with the value of the mixing parameter α kept fixed at 0.5. It is evident from the figure that the p.d.f. of both the LLZI and LLZO distribution behave differently for different for different values of σ, λ and c. σ and λ seem to control the shape of the probability curve, whose skewness decreases with an increase in the value of σ. Further, the peakedness of the curve also seem to decrease with an increase in the value of λ. In each sub plots in figure [START_REF] Burch | Zero-inflated beta distribution applied to word frequency and lexical dispersion in corpus linguistics[END_REF], the vertical bar with a circle above represents α = 0.5 (P r (Y = 0) or P r (Y = 1)). The functional shape of the both the zero inflated and one inflated distribution seem to vary with the combination of σ and λ values.

The zero-and-one inflated Log Lindley distribution

The zero-or-one inflated Log Lindley distribution proposed in the previous section is not suitable for moedling data that may not only contain values in the range [0, 1) or (0, 1] but also in [0, 1]. We refer to this case as double inflation. For modeling double inflated data, a mixture of Log Lindley distribution and Bernoulli distribution, which assigns non-zero probabilities to both the end points 0 and 1, is to be considered. The c.d.f. of this mixture distribution, referred to as zero-and-one inflated Log Lindley distribution (LLZOI), is given by

LLZOI (y; α, p, σ, λ) = αBer (y; p) + (1 -α) F (y; σ, λ) (4)
where y ∈ [0, 1], Ber (; p) denotes the c.d.f. of a Bernoulli r.v. with parameter p and F (; σ, λ) is the c.d.f. of the LogLindley distribution with parameters σ and λ. Further, α is the mixing parameter which lies between 0 and 1. It follows that the p.d.f. corresponding to the c.d.f. in (4) is

llzoi (y; α, p, σ, λ) =      αp, if y = 1 α (1 -p) , if y = 0 (1 -α) f (y; σ, λ) , if y ∈ (0, 1) (5)
We then say that the r.v. Y having the p.d.f. given in [START_REF] Johnson | Continuous univariate distributions[END_REF]

Properties

Result 1: The r th raw moment of the zero-or-one inflated Log Lindley distribution defined in (3) is given by

E (Y r ) = αc + (1 -α) µ r ; r = 1, 2, ... where µ r = σ 2 1+λσ 1+λ(σ+r)
(σ+r) 2 is the r th raw moment of the Log Lindley distribution.

In particular the mean and variance of Y are

E (Y ) = αc + (1 -α) σ 2 1 + λσ 1 + λ (σ + 1) (σ + 1) 2 (6)
V ar

(Y ) = αc (1 -αc) + σ 2 1 + λσ (1 -α) 1 + λ (σ + 2) (σ + 2) 2 - 1 + λ (σ + 1) (σ + 1) 2 1 + (1 -α) σ 2 1 + λσ 1 + λ (σ + 1) (σ + 1) 2 (7)
Result 2: The r th raw moment of Y ∼ ZOILL is

E (Y r ) = αp + (1 -α) µ r ; r = 1, 2, ... (8)
In particular, the mean and variance of Y are

E (Y ) = αp + (1 -α) σ 2 1 + λσ 1 + λ (σ + 1) (σ + 1) 2 (9)
V ar

(Y ) = αp (1 -αp) + σ 2 1 + λσ (1 -α) 1 + λ (σ + 2) (σ + 2) 2 - 1 + λ (σ + 1) (σ + 1) 2 1 + (1 -α) σ 2 1 + λσ 1 + λ (σ + 1) (σ + 1) 2 (10)
Result 3: A special case of zero-or-one inflated Log Lindley distribution defined in (3) (specifically, for fixed value of λ) belongs to the two-parameter exponential family.

Proof: The p.d.f. in (3) can be re-written as

f (y; α, σ, λ) = α Ic(y) (1 -α) 1-Ic(y) σ 2 1 + λσ (λ -log y) y σ-1 1-Ic(y) (11)
where

I c (y) = 1, if y = c (0 or 1) 0, if y ∈ (0, 1) Denoting T 1 (y) = I c (y) T 2 (y) = log y {1 -I c (y)} η 1 = log α (1 + λσ) σ 2 (1 -α) η 2 = σ
We can again write the p.d.f. in expression [START_REF] Tomazella | Zeroadjusted reparameterized Birnbaum-Saunders regression model[END_REF] as

f (y; η 1 , η 2 ) = exp I c (y) log α 1 -α -log σ 2 1 + λσ + {1 -I c (y)} σ log y + log (1 -α) + log σ 2 1 + λσ λ -log y y 1-Ic(y)
Now, putting

T (y) = (T 1 (y), T 2 (y)) , η = (η 1 , η 2 ) , B(η) = log σ 2 (1 -α) 1 + λσ , h(y) = λ -log y y 1-Ic(y)
the p.d.f. can finally be expressed as

f (y; η 1 , η 2 ) = exp η T (y) + B (η) h (y)
Here, the function B(η) is a real valued function of η 1 and η 2 , h(y) is a positive real valued function. The transformation from (α, σ) to (η 1 , η 2 ) is obviously one-one from (0, 1) × R + to R × R + . Hence the p.d.f. in (3) belongs to a two parameter exponential family distribution when λ is held fixed.

Result 4: A special case of zero-and-one inflated Log Lindley distribution defined in (5) (specifically, for fixed value of λ) belongs to the threeparameter exponential family.

Proof: The p.d.f. in ( 5) can be re-written as

f (y; p, α, σ, λ) = α I {0,1} (y) (1 -α) 1-I {0,1} (y) p yI {0,1} (y) (1 -p) I {0,1} (y)-yI {0,1} (y) (1 -α) σ 2 1 + λσ (λ -log y) y σ-1 1-I {0,1} (y) (12)
where

I {0,1} (y) = 1, if y = 1 0, if y = 0 Denoting T 1 (y) = I {0,1} (y)
T 2 (y) = 1 -I {0,1} (y)
T 3 (y) = yI {0,1} (y)

η 1 = log α (1 -p) (1 + σλ) σ 2 (1 -α) 2 η 2 = σ η 3 = log p 1 -p
We can re-write the p.d.f. in expression (12) as

f (y; η 1 , η 2 , η 3 ) = exp I {0,1} (y) log α (1 -p) (1 + λσ) (1 -α) 2 σ 2 + 1 -I {0,1} (y) σ log y +yI {0,1} (y) log p 1 -p + log (1 -α) 2 σ 2 (1 + λσ) λ -log y y 1-I {0,1} (y)
Now, putting

T (y) = (T 1 (y), T 2 (y), T 3 (y)) , η = (η 1 , η 2 , η 3 ) , B * (η) = log σ 2 (1 -α) 2 1 + λσ , h(y) = λ -log y y 1-I {0,1} (y)
the p.d.f. can finally be expressed as

f (y; η 1 , η 2 , η 3 ) = exp η T (y) + B * (η) h (y)
Here, the functions B * (η) is a real valued function of (η 1 , η 2 , η 3 ), h(y) is a positive real valued function. Further, the transformation from (α, p, σ) to (η 1 , η 2 , η 3 ) is one-one from (0, 1) × (0, 1) × R + to R × R + × R. Hence the p.d.f. in (5) belongs to a three parameter exponential family distribution when λ is held fixed.

Estimation

In this section, the estimation of parameters through the maximum likelihood (ML) method and moment estimation method of both the zero-or-one and zero-and-one inflated Log Lindley distribution are considered and the construction of the Fisher Information matrix for ML estimators are discussed.

MLE: Zero-or-one inflated Log Lindley distribution

The likelihood function for ν = (α, σ, λ) of the random sample y = (y 1 , y 2 , . . . , y n ) from iLL c is given by

L (ν, y) = n i=1 iuL c (y i ; α, σ, λ) = n i=1 α I [c] (y i ) (1 -α) 1-I [c] (y i ) f (y; σ, λ) 1-I [c] (y i ) = L 1 (α, y) × L 2 (σ, λ, y)
It can be observed that the likelihood function has been factorized into two terms L 1 and L 2 where L 1 depends on α and L 2 only on σ and λ. Now,

L 1 (α, y) = n i=1 α I [c] (y i ) (1 -α) 1-I [c] (y i ) = α n i=1 I [c] (y i ) (1 -α) n-n i=1 I [c] (y i )
and

L 2 (σ, λ, y) = n i=1 y i ∈(0,1) f (y i ; σ, λ) 1-I [c] (y i )
The corresponding log-likelihood function of the Zero-or-one inflated Log Lindley distribution is given by: l (ν) = log L (ν, y) = l 1 (α; y) + l 2 (σ, λ; y) where

l 1 (α; y) = log L 1 (α, y) = log α n i=1 I [c] (y i ) + log (1 -α) n - n i=1 I [c] (y i ) and l 2 (σ, λ; y) = log L 2 (σ, λ, y) = n - n i=1 I [c] (y i )     2n log σ -n log (1 + λσ) + n i=1 y i ∈(0,1) log (λ -log y i ) + (σ -1) n i=1 y i ∈(0,1) log y i    
The score function is then obtained by differentiating the log-likelihood function and is denoted by

U (ν) = [U α (α) , U σ (σ) , U λ (λ)],
where

U α (α) = ∂l 1 (α; y) ∂α = 1 α n i=1 I [c] (y i ) - 1 1 -α n - n i=1 I [c] (y i ) U σ (σ) = ∂l 2 (σ, λ; y) ∂σ = n - n i=1 I [c] (y i )     2n σ - nλ 1 + λσ + n i=1 y i ∈(0,1) log y i     U λ (λ) = ∂l 2 (σ, λ; y) ∂λ = n - n i=1 I [c] (y i )     -nσ 1 + λσ + n i=1 y i ∈(0,1) 1 λ -log y i    
The m.l.e's of the parameters are obtained by solving the equations

U α (α) = 0 =⇒ α = 1 n n i=1 I [c] (y i ) U σ (σ) = 0 =⇒ -nλσ -2n σ (1 + λσ) = n i=1 log y i U λ (λ) = 0 =⇒ n i=1 1 λ -log y i = nσ 1 + λσ
Since the second and third maximum likelihood equations are non-linear in σ and λ, so, numerical techniques are to be employed to get a solution to these equations. The Fisher Information matrix for the zero-or-one inflated Log Lindley distribution is

K (ν) =   k αα k αλ k ασ k λα k λλ k λσ k σα k σλ k λλ  
where

k αα = -E ∂U α (α) ∂α = n α (1 -α) k αλ = -E ∂U λ (λ) ∂α = 0 = k λα k ασ = -E ∂U α (α) ∂θ = 0 = k σα k λλ = -E ∂U λ (λ) ∂λ = n 2 σ 2 e λσ (1 -α) 1 + λσ Γ (0, λσ) - (1 -α) n 2 σ 2 (1 + λσ) 2 - n (n -1) α 2 (λ -log c) 4 - n (n -1) ασ 2 e λσ (1 -α) (1 + λσ) (λ -log c) 2 Γ (0, λσ) + nα (n -1) (λ -log c) 2 k λσ = -E ∂U σ (σ) ∂λ = n 2 (1 -α) (1 + λσ) 2 k σλ = -E ∂U λ (λ) ∂σ = -n 2 σ 2 (1 -α) (1 + λσ) 2 k σσ = -E ∂U σ (σ) ∂σ = (1 -α) 2n 2 σ 2 + n 2 λ 2 (1 + λσ) 2
Suppose ν = α, σ, λ denote the m.l.e. of ν = (α, σ, λ). In large samples, ν is asymptotically normally distributed, i.e., ν D -→ N 3 ν, k (ν) -1 where k (ν) is the Fisher Information Matrix. This result can be used to approximate the confidence intervals for the parameters α, σ and λ. Let δ ∈ (0, 0.5). Then (1 -δ)×100% asymptotic confidence intervals for α, σ and λ are given respectively by

α ± Z 1-δ 2 SE (α), σ ± Z 1-δ 2 SE (σ) and λ ± Z 1-δ 2 SE λ
where SE (.) denotes standard error and Z 1-δ 2 is the standard normal quantile at δ level of significance.

MLE: Zero-and-one inflated Log Lindley distribution

An alternate expression of the p.d.f. of the zero-and-one inflated Log Lindley distribution is llzoi (y; α, p, θ) = αp y (1 -p) 1-y I {0,1} (y)

[(1 -α) f (y; σ, λ)] 1-I {0,1} (y)

= α I {0,1} (y) (1 -α) 1-I {0,1} (y) p y (1 -p) 1-y I {0,1} (y)
[f (y; σ, λ)] 1-I {0,1} (y)

where I {0,1} (y) is an indicator function defined as

I {0,1} (y) = 1, if y ∈ {0, 1} 0, if y / ∈ {0
, 1} It can be observed from the above expression that the first term depends only on α, the second term depends only on p and the third term depends only on σ and λ. The likelihood function for ν = (α, p, σ, λ) on the random sample (y 1 , y 2 , . . . , y n ) drawn from the llzoi (y; α, p, σ, λ) distribution is given by

L (ν, y) = n i=1 llzoi (y i ; α, p, θ) = L 1 (α; y) × L 2 (p; y) × L 3 (σ, λ; y)
where

L 1 (α; y) = n i=1 α I {0,1} (y i ) (1 -α) 1-I {0,1} (y i ) = α n i=1 I {0,1} (y i ) (1 -α) n-I {0,1} (y i ) L 2 (p; y) = n i=1 p y i (1 -p) 1-y i I {0,1} (y i ) = p n i=1 y i I {0,1} (y i ) (1 -p) n i=1 (1-y i )I {0,1} (y i ) L 3 (σ, λ; y) = n i=1 y i ∈(0,1) f (y i ; σ, λ) = σ 2n (1 + λσ) n n i=1 (λ -log y i ) n i=1 y σ-1 i
The corresponding log-likelihood function is given by l (ν) = log L (ν, y) = l 1 (α; y) + l 2 (p; y) + l 3 (σ, λ; y) where

l 1 (α; y) = log α n i=1 I {0,1} (y i ) + log (1 -α) n - n i=1 I {0,1} (y i ) l 2 (p; y) = log p n i=1 y i I {0,1} (y i ) + log (1 -p) n i=1 I {0,1} (y i ) - n i=1 y i I {0,1} (y i ) l 3 (θ; y) = 2n log σ -n log     1 + λσ + n i=1 y i ∈(0,1) log (λ -log y i ) + (σ -1) n i=1 y i ∈(0,1) log y i     Again, the score function is U (ν) = [U α (α) , U p (p) , U σ (σ) , U λ (λ)] where U α (α) = ∂l 1 (α; y) ∂α = 1 α n i=1 I {0,1} (y i ) - 1 1 -α n -I {0,1} (y i ) U p (p) = 1 p n i=1 y i I {0,1} (y i ) - 1 1 -p I {0,1} (y i ) -I [1] (y i ) U σ (σ) = 2n σ - nλ 1 + λσ + n i=1 y i ∈(0,1) log y i U λ (λ) = ∂l 2 (σ, λ; y) ∂λ = -nσ 1 + λσ + n i=1 y i ∈(0,1) 1 λ -log y i
The maximum likelihood estimates of the parameters are obtained by solving the following equations:

U α (α) = 0 =⇒ α = 1 n n i=1 y i ∈(0,1) I {0,1} (y i ) U p (p) = 0 =⇒ p = p = n i=1 I {1} (y i ) n i=1 I {0,1} (y i ) U σ (σ) = 0 =⇒ 2n σ - nλ 1 + λσ + n i=1 y i ∈(0,1) log y i = 0 U λ (λ) = 0 =⇒ -nσ 1 + λσ + n i=1 y i ∈(0,1) 1 λ -log y i = 0
Since the last two m.l. equations are non-linear in σ and λ, they are solved using numerical computational techniques. Now, the Fisher Information matrix for the Zero-and-One inflated Log Lindley probability law is

K (ν) =     k αα k αp k αλ κ ασ k pα k pp k pλ k pσ k θα k θp k λλ k λσ k σα k σp k σλ k σσ    
where

k αα = n α (1 -α) k αp = k αλ = k ασ = k pα = k pλ = k pσ = k λα = k λp = k σα = k σp = 0 k pp = nα p (1 -p) k λλ = n 2 σ 2 e λσ (1 -α) 1 + λσ Γ (0, λσ) - (1 -α) n 2 σ 2 (1 + λσ) 2 - n (n -1) α 2 (λ -log c) 4 k λσ = n 2 (1 -α) (1 + λσ) 2 k σλ = -n 2 σ 2 (1 -α) (1 + λσ) 2 k σσ = (1 -α) 2n 2 σ 2 + n 2 λ 2 (1 + λσ) 2 Suppose ν = α,

Method of Moments

The method of moments consists in estimating the parameters by equating the expressions for E (Y r ) to the corresponding sample moments and then proceeding to solve the resulting equations for those parameters. Let us denote the first two sample moments by m 1 and m 2 respectively.

MoM: Zero-or-one inflated Log Lindley distribution

The two equations to be solved for obtaining the moment estimates of σ and λ are: and following the calculations in [START_REF] Goméz-Deniz | The Log-Lindley distribution as an alternative to the beta regression model with applications in insurance[END_REF], the moment estimates of σ and λ can be similarly obtained as in the previous section.

E (Y ) =

  follows the zero-and-one inflated Log Lindley (ZOILL) distribution and we denote it by Y ∼ [0, 1] . Here, αp = P r(Y = 1) and α (1 -p) = P r(Y = 0).

Figure 1 : 5 Figure 2 5 Figure ( 2 )

 15252 Figure 1: p.d.f. plots for the LLZI and LLOI distributions for different values of σ and λ; α = 0.5

2 SE 2 SE

 22 p, σ, λ . In large samples, ν is asymptotically normally distributed, i.e., ν D -→ N 4 ν, k (ν)-1 where k (ν) is the Fisher Information Matrix. Let δ ∈ (0, 0.5). Then (1 -δ) × 100% asymptotic confidence intervals for α, p, σ and λ are given respectively byα ± Z 1-δ ( α), p ± Z 1-δ (p), σ ± Z 1-δ 2 SE θ and λ ± Z 1-δ 2 SE λwhere the symbols have their usual meaning.

m 1 =⇒ 2 - 1 - 3 . 5 1 =⇒

 121351 c] (y i ), we get, following the calculations in[START_REF] Goméz-Deniz | The Log-Lindley distribution as an alternative to the beta regression model with applications in insurance[END_REF],λ = ȳ2 (σ + 2) 2 -ȳ1 (σ + 1) 2 ȳ1 (σ + 1) 2 (σ + 2) -ȳ2 (σ + 2) 2 (σ + 1ȳ1 ) ± 2 (ȳ 1 -2ȳ 2 + ȳ1 ȳ2 )2ȳ MoM: Zero-and-one inflated Log Lindley distributionThe two equations to be solved for obtaining the moment estimates of σ and λ are:E (Y ) = m c] (y i )and p by n i=1 I {1} (y i ) n i=1 I {0,1} (y i )

Concluding remarks and discussion

Inflated Log-Lindley distributions is introduced and some of their distributional properties are investigated. Maximum likelihood and moment method of parameter estimation is outlined. Further results, simulation to check the efficacy of estimation and real life examples will be added in the next version.