M Di

A Pirrotta

M Zingales

Itô calculus extended to systems driven by alphastable lévy white noises (a novel clip on the tails of lévy motion)

Keywords: α-Stable Processes, Einstein-Smoluchowsky Equation, Itô Calculus, Truncated Lévy Motion

The paper deals with probabilistic characterization of the response of non-linear systems under α-stable Lévy white noise input. It is shown that, by properly selecting a clip in the probability density function of the input, the moments of the increments of Lévy motion process remain all of the same order ( ) dt , like the increments of the Compound Poisson process. It follows that the Itô calculus extended to Poissonian input, may also be used for α-stable Lévy white noise input processes. It is also shown that, when the clip on the tails of the probability of the increments of the Lévy motion approaches to infinity, the Einstein-Smoluchowsky equation is restored. Once these concepts are outlined extension to single oscillator is readily obtained. A discussion on the proper way to perform Monte Carlo simulation is also exploited.

INTRODUCTION

Normal white noises are very popular stochastic processes and they have been used to model several types of physical phenomena. Such processes may be defined as formal time derivative of Wiener processes. In this setting the powerful machinery of Itô stochastic differential calculus may be used to yield the response probabilistic characterization of systems driven by normal However, several real phenomena observed in physics, seismology, electrical engineering, economics and in some other research fields show evident non-Gaussianity either in heavy tails distribution or in the impulsive nature of the recorded samples. According to this new need, the Itô stochastic differential calculus has been extended to Poissonian white noises too, providing the equation governing the evolution of the probability density, known as Kolmogorov-Feller equation [START_REF] Iwanckievicz | Dynamic Response of non-linear systems to Poisson distributed random pulses[END_REF][START_REF] Pirrotta | Non-linear systems under parametric white noise input: digital simulation and response[END_REF][START_REF] Pirrotta | Multiplicative cases from additive cases: Extension of Kolmogorov-Feller equation to parametric Poisson white noise processes[END_REF][START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF]. The need for non-Gaussian models, to describe the fluctuations exhibited by non-Gaussian phenomena, has also raised the interest in the so-called α -stable Lévy processes [START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF][START_REF] Grigoriu | Linear and NonLinear Systems with non-Gaussian White Noise Input[END_REF][START_REF] Metzler | Boundary Value Problems for Fractional Diffusion Equations[END_REF][START_REF] Grigoriu | Equivalent Linearization for Systems driven by Lévy White Noise[END_REF].

This kind of stochastic process is characterized by the knowledge of four parameters, which are, respectively, the stability index α , the scale factor σ , the skewness β and the shift μ The main challenge in the analysis of dynamical systems in presence of α -stable Lévy white noise is related to the divergence of statistical moments of the α -stable random variable ( )

L t α , namely ( ) p E L t α ⎡ ⎤= ∞ ⎣ ⎦ if p α ≥ .
Then unless the case of normal white noise [START_REF] Pirrotta | Non-linear systems under parametric white noise input: digital simulation and response[END_REF] α = all the moments starting from the variance would diverge because of the heavy tails of the distribution of the random variable ( )

L t α .
In order to overcome this drawback some truncation of the PDF of α -stable random variable ( ) L t α have recently been performed [START_REF] Sokolov | Fractional diffusion equation for a power-lawtruncated Lévy process[END_REF] in which a power law with exponent (5-α ) has been used to truncate the PDF of α -stable distribution. In this context a modified fractional FPK equation has been obtained and the PDF of the response converges towards a Gaussian density in the central part. Tails of solution of the response of nonlinear systems are discussed in [11].

A c c e p t e d m a n u s c r i p t 3

In this paper, the Itô calculus will be extended to non-linear systems under α-stable Lévy white noise. In order to aim at this, a new form of truncation is introduced, different from that proposed in [START_REF] Sokolov | Fractional diffusion equation for a power-lawtruncated Lévy process[END_REF]. We, in fact, exclude the tails on the PDF for value greater than ( )

INTRODUCTORY REMARKS ON ITÔ DIFFERENTIAL CALCULUS

In this section some well known concepts of Itô calculus are briefly summarized for sake of clarity as well as for introducing appropriate symbology.

Let the non-linear system under normal white noise be given in the form:

( ) ( ) ( ) 2 , , Z f Z t g Z t W t = + (1)
where ( )

2
W t is the normal white noise characterized by the correlation function

( ) ( ) 2 1 2 2 2 1 , W R t t q t t δ = - being ( ) δ ⋅ a Dirac delta function. Moreover 2
q is the strength of the white noise,

( )

, f Z t and ( ) , g Z t are deterministic non-linear functions of the random process ( ) Z t and of time t. Eq. ( 1) is a non-linear one excited by a parametric excitation.

In the following we will confine ourselves to the case in which ( ) ( )

, g Z t g t =
, that is the function g does not depend on the response ( ) Z t , Eq. ( 1) is referred as a non-linear one (for the presence of the non-linear term ( ) , f Z t ) driven by an external excitation.

A c c e p t e d m a n u s c r i p t 4

From the definition of ( ) 2 W t , as formal time derivative of Wiener process ( ) B t having independent increments, it follows that Eq.( 1) may be rewritten in Itô form as follows:

( ) ( ) ( )

, dZ f Z t dt g t dB t = + (2)
Once the original equation is transformed into an Itô type stochastic differential equation, the Itô differential rule may be used for the probabilistic characterization of the response process ( )

Z t : ( ) ( ) 2 2 2 1 , 2 d Z t dt dZ dZ t Z Z ψ ψ ψ ψ ∂ ∂ ∂ = + + ∂ ∂ ∂ (3)
where ( )

, Z t
ψ is any non-linear function of Z and t , continuous and differentiable on t and twice differentiable on Z . The third term in the right-hand side of eq.( 3) is essential since ( ) dB t is of order of magnitude ( ) 1 2 dt and then term ( ) 2 dZ is of the same order of the first term. In order to derive the equation governing the evolution of the characteristic function 3), taking mathematical expectation and accounting for the non-anticipative property of Itô calculus

( ) ( ) ( ) , exp Z t E i Z t φ θ θ ⎡ ⎤ = ⎣ ⎦ we put ( ) ( ) , exp Z t i Z ψ θ = in eq.(
( ) ( ) ( ) ( ) ( ) , , k k E Z t dB t E Z t E dB t ψ ψ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = ⎣ ⎦ ⎣ ⎦ ⎣ ⎦
the differential equation governing the evolution of the CF is readily written as:

( ) ( ) ( ) ( ) ( ) 2 2 2 , exp , , 2
Z Z t q i E i Z f Z t t g t t φ θ θ θ θ φ θ ∂ È = - Î ∂ (4)
An inverse Fourier transform of eq.( 4) yields the so-called Fokker-Planck-Kolmogorov (FPK) equation governing the evolution of the PDF of the response:

( ) ( ) ( ) ( ) ( ) 2 2 2 2 , , , , 2
Z Z Z p z t p z t q p z t f z t g t t z z ∂ ∂ ∂ È = - + Î ∂ ∂ ∂ (5)
Itô rule extended to the case of Poissonian white noise is reported in Appendix A.

THE CASE OF α-STABLE WHITE NOISE EXTERNAL EXCITATION

Let us now suppose that ( )

2 W t
( ) ( ) ( ) ( ) 1 ,0,0 L t L s S t s α α α α ∼ - - (6a)
ii)

The CF of an increment of ( ) α dL t , takes the form:

( ) exp dL dt α α φ θ θ È = - Î ˚ (6b)
where ( )

dL
( ) L t α .
In the case of Lévy white noise the equation governing the evolution of the PDF, for ( ) ( ) , g Z t g t = (external excitation), is the so-called Einstein-Smoluchowsky (ES) equation [START_REF] Chechkin | Stationary States of Non-linear Oscillators Driven by Lévy Noise[END_REF] involving Riesz-Weil fractional derivative in the diffusion term, that is:

( ) ( ) ( ) ( ) ( ) ( ) , , , , Z Z
Z p z t p z t f z t g t p z t t z z α α α ∂ ∂ ∂ = - + ∂ ∂ ∂ (7)
in which z α ∂ ∂ is the symmetric functional space derivative [START_REF]Fractional Calculus in Physics[END_REF][START_REF] Samko | Marichev Fractional Integrals and Derivatives: Theory and Applications[END_REF] which is defined for a "sufficiently well-behaved" function through its Fourier transform [ ]

i ¡ : ( ) ( ) , , Z Z p z t t z α α α θ φ θ È ∂ ¡ =- Í ∂ Í Î ˚ (8)
or in terms of the Riemann-Liouville derivatives as:

( ) ( ) ( ) ( ) , 1 , , 2 cos 2 Z Z Z p z t D p z t D p z t z α α α α πα + - ∂ È = - + Î ∂ (9)
where 0 α > and if 0 1 α < < then Riemann-Liouville derivatives reads: are guaranteed since the area under the PDF is one. The derivatives in eqs. [START_REF] Sokolov | Fractional diffusion equation for a power-lawtruncated Lévy process[END_REF]11) are characterized in the Fourier transform space as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , 1 , 1 , 1 , 1 z Z Z Z Z z p t D p z t d z z p t D p z t d z z α α α α ξ ξ α ξ ξ ξ α ξ Γ Γ + -• • - ∂ = -∂ - ∂ = - -∂ - Ú Ú ( 10
( ) ( ) ( ) ( ) 1 1 , , n z n n Z Z n D p z t p z t d n z α α ξ ξ ξ α Γ ∓ -- ± -• ± ∂ = -∂ Ú (11)
( ) ( ) ( ) , , Z Z D p z t i t α α θ φ θ ∓ ± È ¡ = Î ˚ (12)
with:

( ) ( )

exp sgn 2 i i α α απ θ θ θ ∓ È - = Í Î ˚ (13) and
( )

sgn i is the well-known signum function. Combining eqs.(9, 11) with eq.( 12), eq.( 8) is recovered. Fourier transform of eq.( 7) yields the CF equation in the form:

( ) ( ) ( ) , , , i Z Z Z t i E f Z t e t t α θ φ θ θ θ φ θ ∂ È = - Î ∂ (14)
in the following termed as spectral ES. Eq. ( 14) has also been obtained in [START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF][START_REF] Ditlevsen | Invalidity of the Spectral Fokker-Planck Equation for Cauchy Noise Driven Langevin Equation[END_REF] using semimartingales.

In the next section it will be presented an alternative way for finding ES equation or the spectral counterpart as well as for properly performing Monte Carlo simulation.

ITÔ CALCULUS FOR EXTERNAL LÉVY WHITE NOISE

Let us consider the dynamic system reported in eq.( 1) (with function ( )

1 g t = )
enforced by Lévy white noise and let us cast the differential equation of motion in incremental form as:

( ) ( ) ( ) ( ) ( ) 1 , k k k k k Z t t Z t f Z t t t L t α α Δ Δ Δ + = + + (15)
in which t Δ is a small time increment and

( ) k L t
α is a realization of an α-stable random variable, so that ( )

1 k t L t α α Δ is the area of the impulse in the interval k k t t t Δ ∏ + .
( ) ( ) ,0,0 k L t S α α σ ∼ by: ( ) { } Pr ob k L t C α α α ν ν - > = ; ν AE • (16)
where ( )

1 0 sin C x x dx α α - ∞ - ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ ∫
is a real number dependent on the stability index α , (for instance, [ ]

0 for 2, 2 for 1, 1 (2) =1 for 0 8 C C C α α α α π α Γ α = = = = = =
). It follows that the probability of occurrence of values of ( )

k L t α larger than ( ) 1 t α - Δ is C t α Δ , that is ( ) ( ) { } 1 Pr ob ; 0 k L t t C t t α α α - > Δ Δ Δ → (17)
Based on this consideration in [START_REF] Di Paola | Stochastic response of linear and non-linear systems to αstable Lévy white noises[END_REF], realization of ( )

k L t α greater than 1 t α - Δ may be neglected.
Looking at eq.( 17), we may say that, introducing a real value n α , ( 1 n α ≥ ∈R ), in such a way:

( ) 1 Pr ob ; 0 k t t L t C t n n α α α α α Δ Δ Δ - Ï Ȩ Ô Ô > AE Ì Á Ë Ô Ô Ó ˛ (18)
the probability remains of order t Δ also for ( )

1 t n α α - Δ
. This means that neglecting outcomes of ( )

k L t α larger than ( ) 1 t n α α - Δ
, we introduce an approximation of order ( ) Δ and in the limit as 0 t Δ → becomes an infinitesimal of higher order thus it may be really neglected. It is worth stressing that, at the limit for α → ∞ n the real PDF is restored.

t
Previous considerations lead us to represent α-stable increments of the Lévy process ( )

dL t α in the form: ( ) ( ) ( ) 1 1 0 lim n t t dL t dt L t L t n α α α α α α α →∞ Δ → ⎛ ⎞ Δ = = ⎜ ⎟ ⎝ ⎠ (19) It is worth noting, that if 2 α = ( ) ( ) ( ) L t B t α →
, then C α =0 and the area of each impulse that is ( )

1 2 k t B t Δ remains finite setting n α =1 in eq.(19).
In order to better clarify the role played by eq.( 19) we may write

( ) ( ) ( ) j j j j j k k L E dL t E L t dt dt x p x dx α α α α α • -• È ˘È = = Î ˚Î ˚Ú (20)
Notice that in eq.( 20) it is present the product of an infinitesimal quantity j dt α and a divergent integral term. The question is: what is the order of magnitude of ( )

j
( )

1 1 0 lim t n j j j k L j t t n E dL t t x p x dx K n dt α α α α α α α α Δ Δ Δ Δ - - AE - È ˘= = Î ˚Ú (21)
where ( ) j K n α are real finite numbers and then from eq.( 21) we may state that the increments of the Lèvy motion are of order dt . This remarkable result allows us to assess that in this perspective, once n α has been selected, the increments of the Lèvy motion behave exactly as the increments of the compound Poisson process [START_REF] Iwanckievicz | Dynamic Response of non-linear systems to Poisson distributed random pulses[END_REF][START_REF] Pirrotta | Non-linear systems under parametric white noise input: digital simulation and response[END_REF][START_REF] Pirrotta | Multiplicative cases from additive cases: Extension of Kolmogorov-Feller equation to parametric Poisson white noise processes[END_REF][START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF] and then Itô rule for Poissonian white noise may be used (see eq. (A5) in the Appendix A.).

As an example we may assume that ( )

L
( ) ( ) 1 2 1 2 2 2 2 0 1 1 lim exp 2 2 t j j j j j t t E dB t t x x dx t K π - - Δ Δ → -Δ ⎛ ⎞ ⎡ ⎤= Δ - =Δ ⎜ ⎟ ⎣ ⎦ ⎝ ⎠ ∫ (23)
where 2 1 0,

j K + = and
( )

2 2 2 2 1 !! j j K q j = -
, (being ( ) ( ) ( )

2 1 !! 2 1 2 3 3 1 j j j - = -⋅ - ⋅ ; with (- 1)!!=1
), that is dB(t) is of order 1 2 dt . It follows that the Itô rules according eq.( 3) contains terms up to the order ( )

2
dZ since the latter remains of order dt . Conversely, when the system is driven by Lèvy α -stable white noise, for fixed n α , all the increments ( )

j
( ) ( ) ( ) ( ) ( ) ( ) 1 exp , exp exp ! j j j i d i Z i f Z t i Z dt i Z dL t j α θ θ θ θ θ • = È ˘= + Î ˚Â (24)
By inserting eq.( 21) into eq.( 24), performing mathematical expectation, dividing by dt and using the non-anticipative property of Itô calculus, we get:

( ) ( ) ( ) ( ) ( ) ( ) 1 , , e x p , ! j Z Z j j i t i E f Z t i Z t K n t j α θ φ θ θ θ φ θ ∞ = ∂ ⎡ ⎤ = + ⎣ ⎦ ∂ ∑ (25)
Summation in eq.( 25) for 1 when α → ∞ n coincides with eq.( 14), and its Fourier transform fully restores the ES equation.

α = leads to ( ) ( ) ( ) ( ) 2 2 1 1 1 2 ! 2! 2 1 j j j j j j i i n K n j j j α α θ θ π - ∞ ∞ = = = - ∑ ∑ (
( ) ( ) ( ) ( ) 2 1 1 1 2 ! ! 2 1 j j j j j j i i n K n j j j α α θ θ π - ∞ ∞ = = = - ∑ ∑ (
On the other hand once these results have been obtained we may also assert that for every value of α the following relation

( ) ( ) ( ) ( ) 1 e x p 1 1 ! j j dL j i E dL t dt dt j α α α α θ θ φ θ θ ∞ = ⎡ ⎤= - -= -- ⎣ ⎦ ∑ (28)
holds true.

At this stage some further observations will be exploited about the validity of eq.( 28) as claimed in [START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF][START_REF] Ditlevsen | Invalidity of the Spectral Fokker-Planck Equation for Cauchy Noise Driven Langevin Equation[END_REF]. It may be observed that formal analytical derivations led to eq.( 28) but along way we introduced two fundamental assumptions:

• Setting the value of n α , moments of increments

j E dL α ⎡ ⎤ ⎣
⎦ exist for every value of j and they are of order dt (eq.21).

• We consider that the series in eq.( 28) converges to the CF of the increment ( ) , for assessing that the error in clipping the tails of the PDF of ( )

dL
k L t α is of order 2 t Δ .
Obviously if we select n α and we compute ( ) j K n α and insert these values in eq.( 24) we get ( )

, E Z t ψ ⎡ ⎤ ⎣ ⎦ .
The same result is obtained by performing MCS excluding realizations of ( )

k L t α larger than ( ) 1 t n α α Δ -
. In this framework some numerical applications will be reported in the following sections.

NUMERICAL APPLICATIONS

In this section applications to linear and cubic non-linear system under αstable Lèvy white noise ( ) 1 α = will be presented. The problem involving cubic non-linear oscillator has been selected since the solution in terms of spectral ES equation has been extensively studied (see i.e.

[4], [START_REF] Chechkin | Stationary States of Non-linear Oscillators Driven by Lévy Noise[END_REF], [START_REF] Ditlevsen | Invalidity of the Spectral Fokker-Planck Equation for Cauchy Noise Driven Langevin Equation[END_REF]) and the existent solutions may be used as benchmark to highlight the role played by the truncation coefficient n α . The analyzed system is ( )

3 Z a Z b Z W t α = + + (29)

Linear system under Cauchy white noise

The linear system, which is the hitherto studied case, is obtained setting 0 b = in eq.( 29). In this case ( ) Z t is also an αstable process having the same stability index as the Lèvy white noise, but with different scale (or amplitude). For this case the spectral ES may be written as:

( ) ( ) ( ) , , , Z Z Z t t a t t φ θ φ θ θ θφ θ θ ∂ ∂ = - ∂ ∂ ( 30 a
( ) exp Z a φ θ θ ⎡ ⎤ = ⎣ ⎦ (31)
The Ornstein-Uhelembeck CF described in eq.( 31) has been reported in fig. = Δ = Δ have been considered for numerical analysis, but the obtained results match the reported CF for the limit value 1 t Δ . This aspect is much more important for numerical simulation of non-linear dynamical systems as it will be shown in the following section.

Non-linear system under Cauchy white noise

Let us now consider that the non-linear coefficient 0 b ≠ in eq.( 29). The associated spectral ES equation is [START_REF] Grigoriu | Stochastic Calculus, Application in Science and Engineering[END_REF], [START_REF] Chechkin | Stationary States of Non-linear Oscillators Driven by Lévy Noise[END_REF], [START_REF] Ditlevsen | Invalidity of the Spectral Fokker-Planck Equation for Cauchy Noise Driven Langevin Equation[END_REF].

( ) ( ) ( ) ( ) ( ) ( ) ( ) 3 3 0 , , , , , , 1 ; 0
; lim , 0 Z Z Z Z Z Z Z t t t a b t t φ θ φ θ θ = →∞ ∂ ∂ ∂ = + - ∂ ∂ ∂ ∂ = = = ∂ (32)
The boundary conditions in eq. (32 c) has been used in [START_REF] Chechkin | Stationary States of Non-linear Oscillators Driven by Lévy Noise[END_REF] and it is valid only if

( ) Z φ θ is differentiable for 0 θ = .
For the steady-state case the boundary value problem described in eq.( 32a, b) may be solved for 0 and 0 (

θ θ > < assuming that ( ) [ ]
1 2 2 1 1 2 1 2 1 2 2 1 1 2 1 2 exp exp ; 0 exp exp ; 0 Z Z λ λ φ θ λ θ λθ θ λ λ λ λ λ λ φ θ λ θ λθ θ λ λ λ λ ⎛ ⎞ = - ≥ ⎜ ⎟ - - ⎝ ⎠ ⎛ ⎞ = - - + - < ⎜ ⎟ - - ⎝ ⎠ ) ( ) ( ) ( ) ( ) ( )
In order to assess the validity of eqs.(34 a,b) MCS has been performed estimating the CF with 747.000.000 deviates and for different values of n α . The time step of numerical integration has been set to In more detail we notice that for values of the truncation parameter too small, say 1 n α = , the estimated CF does not coincide with the exact expression in eq.(34 a,b). However as soon as the value of n α reaches 1 t Δ the estimated CF perfectly match the exact expression and also for higher values of the truncation parameter, exceeding the limit value 1 n t α = Δ , results do not change. Moreover, it must be stressed that Monte Carlo simulation of non-linear systems fails due to numerical overflows as soon as the truncation parameter exceeds the value 20 n t α = Δ and then a clipping control parameter 1 n t α = Δ is the right choice, being the committed error of order 2 t Δ . Moreover at the limit when 0 t Δ → and the n α → ∞ the ES equation is fully restored. Question raised in [START_REF] Ditlevsen | Invalidity of the Spectral Fokker-Planck Equation for Cauchy Noise Driven Langevin Equation[END_REF] about the validity of the spectral ES equation remains clarified in sense that the spectral ES and hence the ES equation are fully valid and the Itô calculus may still be applied.

k k k L t L t L t α α α > > > ( ) ( ) ( ) 1000 (stars) ; 50000 (triangles) ; 1000000 (squares)
K W t E R K W t W t E R t t K W t W t W t E R t t t t λ λ δ λ δ δ ⎡ ⎤= ⎣ ⎦ ⎡ ⎤ ⎡ ⎤= - ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ = - - ⎣ ⎦ ⎣ ⎦ (A2)
White noise process may be considered as the formal derivative of the Compound Poisson process ( )

C t with increments characterized by the expressions:

( ) r r E dC t E R dt λ ⎡ ⎤ ⎡ ⎤ = ⎣ ⎦ ⎣ ⎦ (A3)
Eq.(A1) may be reported in incremental form as:

( ) ( ) ( )

, dZ t f Z t dt dC t = + (A4)
and since ( )

r
dC t is of order dt r ∀ the Taylor expansion of ( ) , Z t ψ will be rewritten in its complete form (omitting arguments) as:

( ) ( ) ( ) 1 1
,

j j j j j j j j d dt dZ t Z dt f Z t dt dC t Z Z ψ ψ ψ ψ ψ ψ ∞ = ∞ = ∂ ∂ = + = ∂ ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∑ ∑ (
( ) ( ) ( ) ( ) ( ) j Z j Z j 1 , t i i E f Z,t exp i Z ,t E R t j ! φ ϑ ϑ ϑ ϑ λ φ ϑ ∞ = ∂ ⎡ ⎤ ⎡ ⎤ = + ⎣ ⎦ ⎣ ⎦ ∂ ∑ (A6)
( ) ( ) j j R j 1 i 1 E R j! ϑ φ ϑ ∞ = ⎡ ⎤ = + ⎣ ⎦ ∑ (A7)
Then eq.(A6) may also be rewritten

( ) ( ) ( ) ( ) ( ) ( ) Z Z R , t i E f Z,t exp i Z ,t 1 t φ ϑ ϑ ϑ λ φ ϑ φ ϑ ∂ ⎡ ⎤ = + - ⎣ ⎦ ∂ (A8)
An inverse Fourier transform of eq. (A8) yields the Kolmogorov-Feller equation ruling the evolution of the PDF

( ) ( ) ( ) ( ) ( ) ( ) ( ) Z Z Z Z R p z,t f z,t p z p z p z p d t z λ λ ξ ξ ξ ∞ -∞ ∂ ∂ = - - + - ∂ ∂ ∫ (A9)
A c c e p t e d m a n u s c r i p t

1 1 2 2 2 2 2 8 1 1 2 2 2 2 8 a b a b Y t Y t Y t Y t Y t W t a b a b Y Y t Y t Y t Y t α ε ε ⎧ + - ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ = - - + + + + ⎜ ⎟ ⎜ ⎟ ⎪ ⎜ ⎟ ⎪ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎨ + + ⎛ ⎞ ⎛ ⎞ ⎪ = -+ - + + + ⎜ ⎟ ⎜ ⎟ ⎪ ⎝ ⎠ ⎝ ⎠ ⎩ (B4)
Eqs.(B4 a,b) may be cast in an equivalent vector form as:

( ) ( ) ( ) , d t t dt dL t α = + Z g Z v (B5)
where ( ) ( ) ( )

1 2 T t Y t Y t = ⎡ ⎤ ⎣ ⎦ Z , ( ) ( ) ( ) 1 2
, , ,

T t g t g t = ⎡ ⎤ ⎣ ⎦ g Z Z Z , [ ]
1 0 T = v
with the non-linear functions: ( )

( ) ( ) ( ) ( ) 3 1 1 2 1 2 3 2 1 2 1 2 1 1 , 2 2 2 2 8 1 1 , 2 2 2 2 8 a b a b g t Y Y Y Y a b a b g t Y Y Y Y ε ε + - ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ = - - + + + ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ + + ⎛ ⎞ ⎛ ⎞ = -+ - + + + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ Z Z (
( ) [ ] [ ] ( ) 1 exp , exp exp ! T T T T j j j j T j d i i t i i dL i j α ∞ = ⎡ ⎤ ⎡ ⎤ = ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ + ⎣ ⎦ T j j j j j d t i E t i dt i E d L t j α φ φ ∞ = ⎡ ⎤ ⎡ ⎤ = ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ + ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ T j j j j j j j n j j i i E dL lim K n dt dt j j α α α α θ θ ∞ ∞ →∞ = = ⎡ ⎤ = = - ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ∑ ∑ v θ (

A c c e p t e d m a n u s c r i p t 2 white

 2 noises, in terms of moments or of conditional probability density function (PDF), that is the Fokker-Planck-Kolmogorov (FPK) equation. The Fourier transform, of the FPK equation, rules the evolution of the characteristic function (CF) of the response, often termed as spectral FPK equation.

2 α<

 2 [START_REF] Samorodnitsky | Stable non-Gaussian Random Processes: Stochastic Models with Infinite Variance[END_REF]. The choice for different values of the stability index α ( 0 ≤ ) leads to formidable variety of stochastic processes including the Gaussian white noise obtained for 2 α = . Linear and non-linear systems driven by external α -stable Lévy white noise processes (formal time derivative of the Lévy motion processes labeled as treated, in the past, either in terms of PDF Einstein-Smoluchowsky (ES) equation or in terms of CF (see e.g. [9]]). Moreover closed-form expressions of the probability density function of dynamical systems driven by external α -stable Lévy noises have been obtained for limited values of the stability index.

  OE . In the specific case 1 α = the Riemann-Liouville derivative coalesces with the Hilbert transform of the function

A c c e p t e d m a n u s c r i p t 7 .

 7 Peculiarity of α-stable random variables is long tail PDF. This feature will be reflected on experiencing eventual very high value This leads to numerical overflow in the response evaluation. To avoid the latter problem, it should be selected the highest value of a large real value ν of an α-stable random variable

A c c e p t e d m a n u s c r i p t 8 2 t

 2

  α = for j odd and for j even.

A c c e p t e d m a n u s c r i p t 9

 9 Moreover

(

  26) And in the limit when α → ∞ n converges to θ as shown in fig.(1a) in which the sum in the series (26) has been reported for different values of n α . Summation in eq.(Lèvy distribution) is given as

A c c e p t e d m a n u s c r i p t 10 limit

 10 α have been tested for different values of α and always they give that at the may conclude that eq.(25)

A c c e p t e d m a n u s c r i p t 11

 11

  eq.(30 b) representing the attendant boundary condition. Steady-state solution of the differential equation described in eqs.(30) reads:

  (2) with continuous line. Such exact solution has been contrasted with results obtained by MCS for different values of n α . The values of the time interval has been selected to value eq.(31) with the decaying factor 1 a = -, has been depicted in fig.(2a) with continuous line. The CF obtained via MCS have been contrasted in fig. (2a) for several values of the truncation parameter n α set equal to 1 n α = (star)

A c c e p t e d m a n u s c r i p t 13 disregarded

 13 λ obtained by the solution of the algebraic equation: latter positive real root must be since it does not fulfill the boundary condition

  the results of the estimation have been reported in fig.(3a) contrasted with solution in eq.(34). Values of the truncation parameter n α has been set to 1 n α = (star), 50 n α = (triangle) and for the limit value 1 1000 n t α = Δ = (square). The selected values of n α

according to sec. 4 .

 4 Close observation of fig.(3a) shows that in the limiting case 1 n t α = Δ the estimated CF totally overlaps the exact CF in eqs.(34 a,b) leading to conclude that the assumption reported in [9] about the differentiability of the CF at 0 θ = is correct. This consideration is still more evident in fig. (3b) where the behavior of the CF in the close vicinity of the origin has been investigated.

A c c e p t e d m a n u s c r i p t 14 6 .A c c e p t e d m a n u s c r i p t 15 order 2 dt

 6152 CONCLUSIONSIn this paper the Itô stochastic differential rule holding for normal and Poissonian white noise has been extended to analysis of Lévy white noise in presence of external excitation. Extension of the differential rule to analysis of Lévy white noise has led to the formulation of the Einstein-Smoluchowsky fractional differential equation in terms of characteristic function. This remarkable result relies on the assumption to neglect increments of the Lévy white noise larger than ( context the Einstein-Smoluchowsky equation is restored. This latter consideration allows to perform Monte-Carlo simulation of dynamical systems under stable αexternal excitation with an opportune clip on the increments of the stable α -Lévy flights. Some numerical analyses, contrasting the proposed formulation with a Monte-Carlo simulation obtained neglecting values of the realization of the Lévy motion to assess the validity of the concepts here exploited. Analyses have been conducted for linear and non-linear dynamical systems excited by a particular class of Lévy flights with stability index 1 α = , namely Cauchy white noises, which allows for an exact solution in terms of spectral ES. In this context a simple Langevin equation has been investigated to yield the stationary characteristic function of the response contrasted with the proposed Monte-Carlo simulation method. The numerical analyses conducted have shown that the proposed clip on tails of the probability density function of the increments of the Lévy flights yields estimates which are in good agreement with the estimate via the Monte-Carlo simulation, contrasted by solving the spectral ES. Remarkably it has been shown that the use of stochastic differential calculus yields the same differential equation for the characteristic function obtained via Fourier transform of the fractional Einstein-Smoluchowsky differential equation for the pdf of the dynamic response. This latter consideration may raise some comments about the range of validity of the Einstein-Smoluchowsky equation or of its spectral counterpart in the analysis of non-linear dynamical systems under stable α -Lévy flights. As in fact we may conclude that the ES equation may be considered as the limit case when the clip on the tails reaches infinity and the error made is of . In this perspective question raised in [14] about the validity of the spectral ES equation remains fully clarified. That is ES equation and own spectral counterpart are valid as soon as differentials of order greater than 2 dt are neglected as usual done.

A c c e p t e d m a n u s c r i p t 18

 18 Since the Taylor expansion of CF of the random variable R is given as

B6) A c c e p t e d m a n u s c r i p t 20 The

 20 Itô rule for Poissonian white noise is written as: vector, the exponent into brackets is the Kronecker power and ( ) ,t ψ Z is any real-valued function continuosly differentiable in t and times ∞differentiable in Z . Particularization of eq.(B7) with

A c c e p t e d m a n u s c r i p t 21 Figure Caption Figure 1a :

 21Caption1a Figure Caption

Figure 1b :

 1b Figure 1b: Asymptotic trend of the series ( ) , f n α θ

Figure 2a :

 2a Figure 2a: Characteristic function of the linear system with 1 a =contrasted with Monte Carlo estimates with different values of the truncation parameter n α .

Figure 2b :

 2b Figure 2b: Behavior of the Characteristic function of linear system in the neighborhood of the origin contrasted with different values of the truncation parameter n α .

Figure 3a :

 3a Figure 3a: Characteristic function of the non-linear system with 1, 0.1 a b = -=contrasted with Monte Carlo estimates with different values of the truncation parameter n α .

Figure 3b :A c c e p t e d m a n u s c r i p t 22 Figure 1a :Figure 1b :Figure 2a :

 3b221a1b2a Figure 3b: Behavior of the Characteristic function of non-linear system in the neighborhood of the origin contrasted with different values of the truncation parameter n α .

Figure 2b :u s c r i p t 24 Figure 3a :

 2b243a Figure 2b: Behavior of the Characteristic function of linear system in the neighborhood of the origin contrasted with different values of the truncation parameter n α

Figure 3b :

 3b Figure 3b: Behavior of the Characteristic function of non-linear system in the neighborhood of the origin contrasted with different values of the truncation parameter n α .

A c c e p t e d m a n u s c r i p t

	i)	It has independent, stationary increments following the	α -	stable	distribution, that is
			is substituted by	( ) W t α	, being	( ) W t α	the α -stable Lévy white
	noise process. In analogy to the definition of the previous white noise, an α -stable Lévy white
	noise may be defined as the formal derivative of a corresponding α -stable Lévy motion	( ) L t α
	(α is the stability index of the	α -	stable	process) having the following properties:
						5

  t

		α	represents the increment of α -stable Lévy motion	( ) L t α	and	( ) k L t α	, at a selected
	time ( ) k t , is a symmetric α-stable ( ) S S α	random variable.
	iii)	For	2 α = ,	( ) L t α	reverts to ( ) B t and then the non-normal	α -	stable	process reverts to
		normal white noise process, or in other words the Wiener process is a particular case of
		the random process

  .4. It may be observed that for small values of n α the Monte-Carlo estimate of the CF does not coincide with CF in eq.(31). A different scenario is observed increasing the value of n α until its upper bound 1 t Δ and it may be noticed that the CF obtained via numerical simulation coincides with the exact reported in eq.(31). A more precise description of the trend of the CF, for different n α , is furnished by observation of fig.(2b).In more detail, it may be observed that, a smooth CF at the origin appears, assuming 1 n α = . On the other hand for values of the coefficient n α approaching the limit 1 t Δ the peak of the exact CF at the origin is captured. Several other

	A c c e p t e d clipping values of 2 or 3 n t n t α α	m a n u s c r i p t
	The selected values of		n α	and	Δ	t	led us to exclude values of
	( ) k L t α	>	1000 (stars) ;	( ) k L t α	>	50000 (triangles);	( ) k L t α	>	1000000 (squares)	according to

sec

  B10) yielding the spectral equation, ruling the evolution of the CF as: Fourier transform of eq.(B11) yields the Einstein-Smoluchowsky equation extended for the single non-linear oscillator reported in eq.(B1) in the form:

	( ) , t t ∂ Z φ ∂ θ	=	i	T θ		E	⎡ ⎣	( ) , exp t g Z	⎡ ⎣	i	θ	T	Z	⎤ ⎦	⎤ ⎦	-	α θ φ Z	( ) , t θ	(B11)
	p ∂	Z	( ) , t t ∂ Z	= -∇	Z	T	⎡ ⎣	( ) ( ) , , t p t Z g Z Z	⎤ ⎦	+	∂	1 Z α ∂	α	p	Z	( ) , t Z	(B12)
	Eq.(B12) represents the extension of the ES equation for the single degree of freedom non-linear
	system.

Acknoledgements:

The author are very grateful to the Italian bureau of research and technology trough grant: PRIN 2004 coordinated by Prof. A. Materazzi. This financial support has been strongly appreciated.

Appendix A: Itô Rule for Poissonian White Noise

Let the equation of motion be given as:

( ) ( ) ( )

Appendix B: Extension to a Single Oscillator System

In this appendix the analysis of a non-linear system with two degrees of freedom has been reported. Let the equation of motion of the system represented in the form:

Introducing the Lagrangian parameters

the second-order equation of motion in eq.(B1) may be represented by a system of first-order differential equations in the form:

( ) ( ) ( ) ( )

System of differential equations in eqs.(B2 a,b) may be cast in an alternative form by means of a linear transform which reads:

which may be substituted into system in eqs.(B2 a,b) yielding the equivalent non-linear system:


B S Lazarov

J S Jensen

Low-frequency band gaps in chains with attached non-linear oscillators

Keywords: Non-linear wave propagation, Local resonators, Low frequency band gaps

The aim of this article is to investigate the wave propagation in one-dimensional chains with attached non-linear local oscillators by using analytical and numerical models. The focus is on the influence of non-linearities on the filtering properties of the chain in the low frequency range. Periodic systems with alternating properties exhibit interesting dynamic characteristics that enable them to act as filters. Waves can propagate along them within specific bands of frequencies called pass bands, and attenuate within bands of frequencies called stop bands or band gaps. Stop bands in structures with periodic or random inclusions are located mainly in the high frequency range, as the wave length has to be comparable with the distance between the alternating parts. Band gaps may also exist in structures with locally attached oscillators. In the linear case the gap is located around the resonant frequency of the oscillators, and thus a stop band can be created in the lower frequency range. In the case with non-linear oscillators the results show that the position of the band gap can be shifted, and the shift depends on the amplitude and the degree of non-linear behaviour.

Introduction

The filtering properties of periodic structures with alternating characteristics have been investigated by many authors by using analytical, numerical and experimental methods. Such systems possess interesting filtering properties. Waves can propagate unattenuated along these structures within specific bands of frequencies called propagation or pass bands, and attenuate within bands of frequencies called stop bands, attenuation zones or band gaps. Band gaps in linear systems can also be created by introducing random inclusions or geometric disturbances. External excitation in such structures results in localised response surrounding the external input. The majority of the texts consider linear systems [START_REF] Brillouin | Wave propagation in periodic structures[END_REF][START_REF] Mead | Wave-propagation and natural modes in periodic systems .1. monocoupled systems[END_REF][START_REF] Mead | Wave-propagation and natural modes in periodic systems .2. multicoupled systems, with and without damping[END_REF]. Localisation phenomena can also be observed in perfectly periodic non-linear structures [START_REF] Vakakis | Forced localization in a periodic chain of nonlinear oscillators[END_REF][START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF][START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF]14], where spring-mass chains are studied and the non-linear behaviour is introduced either in the spring between the two neighbour masses, or by adding non-linear springs between the ground and the masses. The applications of these filtering phenomena are mainly in the high frequency range, as the distance between the inclusions has to be comparable with the wave length.

In the beginning of the twentieth century Frahm discovered the vibration absorber as a very efficient way to reduce the vibration amplitude of machinery

A c c e p t e d m a n u s c r i p t

and structures by adding a spring with a small mass to the main oscillatory body [START_REF] Den | Mechanical vibrations 4th Edition[END_REF]. The additional spring-mass system is tuned to be in resonance with the applied load. When the natural frequency of the attached absorber is equal to the excitation frequency, the main structure does not oscillate at all, as the attached absorber provides force equal and with opposite sign to the applied one. The idea can be exploited further by attaching multiple absorbers on a wave guide. Waves are attenuated in a frequency band located around the resonant frequency of the local oscillators, and thus stop bands can be created in the lower frequency range, which is often more important in mechanical applications. The effect has been studied experimentally and analytically in [START_REF] Yu | Low frequency torsional vibration gaps in the shaft with locally resonant structures[END_REF][12][13].

If the attached oscillators are non-linear, the response displays a dependency between the amplitude and the frequency. Very little is known about the filtering properties of the systems in this case. Periodic spring-mass system with attached non-linear pendulums are investigated in [15]. The attached pendulums are considered to be stiffer than the main chain, and they do not introduce band gaps in the lower frequency range. The aim of this paper is to investigate the behaviour of one-dimensional infinite spring-mass chain with locally attached oscillators with linear or non-linear behaviour. The oscillators are considered to be relatively soft compared to the main chain, and to create band gaps in the lower frequency range. The non-linearities in the attached oscillators are considered to be cubic. First the mechanical system together with the equations of motion is presented. The appearance of band gaps is shown in the linear case, and in the case with non-linear attached oscillators the method of harmonic balance is utilised to obtain a system of equations for the wave amplitude, as well as an approximate expression for the wave

Mechanical system and equations of motion

The considered spring-mass system with attached local oscillators is shown in Figure 1. Using the dynamic equilibrium condition for mass j, the equations describing its motion together with the equation of motion of the attached oscillator can be written as

(m + M) d 2 u j dτ 2 + 2ku j -ku j-1 -ku j+1 + M d 2 q j dτ 2 = 0 (1)
M d 2 q j dτ 2 + c qj + k l f (q j ) + M d 2 u j dτ 2 = 0 ( 2
)
where k is the spring stiffness between two masses, m is the mass at position j = 0, 1, . . . , ∞, M is the mass of the attached local oscillator, k l f (q j ) is its restoring force and c is a viscous damping coefficient. u j is the displacement of mass j, with positive direction shown in Figure 1, and q j is the relative displacement of the local oscillator. Both u j and q j are functions of the time and their positions j. The function f (q j ) is assumed to be in the form

f (q j ) = q j + γ j q 3 j ( 3
)
where γ j is a parameter controlling the degree of non-linearity. By introducing normalised time t = ω 0 τ , where ω 0 = k/m, the coefficients β = M/m and α = k l /k, the equations ( 1) and (2) can be written as üj + 2u j -u j-1 -u j+1 -βκ 2 f (q j ) -2βζκ qj = 0 (4) qj + 2ζκ qj + κ 2 f (q j ) + üj = 0 (5) where κ = α/β and () denotes the second derivative with respect to the normalised time t.

3 Band gaps in spring-mass chains with attached oscillators

Linear undamped oscillators

Using the idea for the vibration absorber, an efficient filter for waves propagating along a spring-mass chain can be created by attaching multiple absorbers to the chain as shown in Figure 1. The equations of motion are given by ( 4) and [START_REF] Jensen | Phononic band gaps and vibrations in one-and two-dimensional mass-spring structures[END_REF], where the non-linear parameter γ and the damping parameter ζ are set to zero. In order to study the transmission properties of the chain, the solution is sought in the form [START_REF] Brillouin | Wave propagation in periodic structures[END_REF] u j = Be jµu+iωt [START_REF] Mead | Wave-propagation and natural modes in periodic systems .1. monocoupled systems[END_REF] where µ u is the so-called wave propagation constant of harmonic wave with frequency ω, which is equal to the wave number multiplied by the distance between the masses in the main chain, and j is the mass index. The displacements of the neighbour masses j + 1 and j -1 become

A c c e p t e d m a n u s c r i p t

u j+1 = e µu u j , u j-1 = e -µu u j [START_REF] Mead | Wave-propagation and natural modes in periodic systems .2. multicoupled systems, with and without damping[END_REF] By inserting ( 6) and ( 7) into the equations of motion ( 4) and ( 5), and requiring the resulting equations to be valid at all time instances, the following relation between µ u and ω is obtained

cosh(µ u ) = 1 + 1 2 ω 4 -(1 + β) κ 2 ω 2 κ 2 -ω 2 (8)
Eq. ( 8) is also known as the dispersion relation.

Non-linear local oscillators

A solution of the non-linear system of equations ( 4) and ( 5) is not known in a closed form. An approximate solution can be obtained by using the method of harmonic balance, the method of multiple scales or the method of averaging, e.g. [START_REF] Thomsen | Vibrations and stability[END_REF][START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Non-linear oscillators with third order non-linearities, as the one described by [START_REF] Jensen | Phononic band gaps and vibrations in one-and two-dimensional mass-spring structures[END_REF], are well studied in the literature. Based on these solutions the method of harmonic balance is applied here, in order to study the wave propagation along the chain with attached damped non-linear oscillators.

The time response of q j (t) is assumed to be periodic in the form of complex Fourier series

q j (t) = k ε k-1 2 A k,j e ikωt + ε k-1 2
Āk,j e -ikωt k = 1, 3, . . .

where ε is a dimensionless bookkeeping parameter showing the order of the amplitude of the motion. By substituting (9) into ( 5) and integrating twice with respect to time, the following expression for the time response of mass j can be obtained

u j (t) = A 1,j -1 + κ 2 ω 2 + 3ε γ j κ 2 ω 2 A 1,j Ā1,j + 2iζκ ω e iωt + +ε -A 3,j + 2 3 iζκ ω A 3,j + 1 9 κ 2 ω 2 A 3,j + 1 9 γ j κ 2 ω 2 A 3 1,j e i3ωt + c.c. + O(ε 2 ) (10)

A c c e p t e d m a n u s c r i p t

The expression in the above equation is obtained by truncating the time series and preserving only the first two terms. By substituting [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF] and ( 9) into ( 4)

and equating to zero the coefficients in front of e ikωt , a system of algebraic equations for the amplitudes A k,j can be obtained

2 κ 2 ω 2 -1 + ω 2 -β 1 κ 2 + 2iζ κ ω 2 -β 1 ω 2 A 1,j +3γ j κ 2 ω 2 2 -β 1 ω 2 A 2 1,j Ā1,j + 1 - κ 2 ω 2 -2iζ κ ω (A 1,j+1 + A 1,j-1 ) -3γ j+1 κ 2 ω 2 A 2 1,j+1 Ā1,j+1 -3γ j-1 κ 2 ω 2 A 2 1,j-1 Ā1,j-1 = 0 (11) 2 κ 2 ω 2 3 -1 + ω 2 3 -β 1 κ 2 + 2iζ κ ω 3 2 -β 1 ω 2 3 A 3,j + 2 κ 2 ω 2 3 -β 1 κ 2 γ j A 3 1,j + 1 - κ 2 ω 2 3 -2iζ κ ω 3 (A 3,j+1 + A 3,j-1 ) -γ j+1 κ 2 ω 2 3 A 3 1,j+1 -γ j-1 κ 2 ω 2 3 A 3 1,j-1 = 0 (12)
where β 1 = β +1 and ω 3 = 3ω. By specifying the amplitudes for two neighbour masses and respectively for the attached oscillators, a solvable system of equations for the response amplitudes of the attached oscillators can be obtained.

The system is non-linear and there can be multiple solutions satisfying it in certain cases. The amplitude of the wave travelling along the main chain can be obtained by inserting the amplitude of the attached oscillators into [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF].

The wave travelling along the spring-mass chain can also be investigated by using the wave number multiplied by the distance between the masses in the main chain µ u . The displacements of the masses neighbour to mass j can be expressed as

u j±1 = k B k,j e ±µu k e ikωt + c.c (13)

A c c e p t e d m a n u s c r i p t

where B k,j are the coefficients in front of e ikωt from equation [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF]. It can be clearly seen that the solution in the non-linear case consists of a zero order wave with frequency ω and additional high order waves. As the nonlinearities are assumed to be small and the focus is on the filtering properties of the chain around the linear natural frequency of the attached oscillators, the contributions to the solution of the harmonics of order equal to and higher than k = 3 are neglected. Substituting (13) and ( 9) into (4), equating the coefficients in front of e iωt to zero and solving the resulting equation with respect to µ u , the dispersion relation for the non-linear case is obtained in the

form cosh(µ u ) = 1 - 1 2 ω 2 - 1 2
ω 2 β κ 2 + 2iζκω + 3κ 2 γ j A 1,j Ā1,j κ 2 + 2iζκω + 3κ 2 γ j A 1,j Ā1,j -ω 2 (14)
The dispersion relation for the linear case can be obtained from ( 14 of the stop band can be increased by increasing the attached mass.

For the non-linear spring-mass chain, µ u depends on the amplitude of the attached oscillator at position j. Plots of the dispersion relation for different

A c c e p t e d m a n u s c r i p t

values of the non-linear parameter γA 1,j Ā1,j are shown in Figure 5. As can be seen, the position of the maximal value of the attenuation rate varies with the value of the nonlinear parameter. The maximum is shifted above the linear natural frequency of the attached oscillators, and the shift increases for larger values of the amplitude A 1,j . For a negative value of the non-linear parameter γ, the shift will occur in the opposite direction, below the linear natural frequency. As the wave propagates, the amplitude A 1,j decreases, due to the reflections from the attached oscillators, as well as due to the damping in the system, and the position of the maximal value of Re (µ u ) moves toward ω/κ = 1. The amplitudes vary along the chain, and µ u becomes a function of the spatial coordinate j. A fixed wave propagation constant cannot be defined in the non-linear case.

Transmission properties based on analytical calculations

For given values of the amplitudes B 1,j and A 1,j of mass j and the corresponding attached oscillator, the value of µ u can be calculated using (14).

The amplitude of the mass at position j + 1 is calculated as B 1,j+1 = e µu B 1,j .

The amplitude A 1,j+1 of the attached oscillator can be calculated by solving the following equation

B 1,j+1 = A 1,j+1 -1 + κ 2 ω 2 + 3 γ j+1 κ 2 ω 2 A 1,j+1 Ā1,j+1 + 2iζκ ω (15)
which is obtained from ( 10) and ( 13) by equating the coefficients of e iωt .

Eq. ( 15

9γ 2 r 4 R 6 A + 6γ j+1 r 2 r 2 -1 R 4 A + r 2 -1 2 + 4ζ 2 r 2 R 2 A -R 2 B = 0 (16)
where r = κ/ω. Eq. ( 16) can be represented as a cubic equation with respect to R 2 A . It has either one, or three real roots for R 2 A . In the second case, two of the solutions are stable and one of them is unstable, which is well known property of the Duffing oscillator [START_REF] Thomsen | Vibrations and stability[END_REF]. There are two values of µ u corresponding to the two stable solutions of ( 16), and a unique value for B 1,j+2 cannot be obtained. The two solutions for µ u,j+1 can be ordered by the magnitude of their real part. The one with larger absolute real value, µ l u,j+1 , produces an amplitude with a smaller absolute value, and the other one, µ u u,j+1 , produces an amplitude with a larger absolute value. Continuing in this way, an estimate for the upper and the lower bounds of the absolute value of the wave amplitude can be obtained.

B u 1,j+k = n=1...k e Re(µ u u,j+n ) B u 1,j B l 1,j+k = n=1...k e Re(µ l u,j+n ) B l 1,j (17)
In the case where only a single real solution for R 2 A exists, µ l u,j+k is equal to

µ u u,j+k and B u 1,j+k = B l 1,j+k .

Numerical simulations

The theoretical analysis derived in the previous section is validated by using numerical simulations for finite spring-mass chain, as shown in Figure 1. Ab- The system is integrated numerically, until steady state is reached. The RMS (root-mean-squared) value of the response amplitude

[B j ] = 1 √ 2 1 T t+T t u 2 j (t) dt (18)
is calculated at the right end of the chain and is compared with the results obtained by using the analytical prediction derived in the previous section. The RMS value is calculated for a finite time interval T = 20λ, where λ = 2π/ω and ω is the frequency of the generated wave. The contribution of the higher order harmonics to the RMS value of the response amplitude is assumed to be small.

Influence of the nonlinearities

The transmission properties of a chain with 2000 attached oscillators calculated by using numerical simulations, and the procedure in section 3.

Chains with variable non-linear coefficients

If the wave frequency is located in the band gap, the amplitude is decaying fast, and the band gap is shifted back to the one obtained for linear attached oscillators. The shift depends on the value of γ j A 1,j Ā1,j . As the amplitude decreases along the chain, γ j can be changed, in order to keep the shift of the band gap at a desirable frequency. For very small amplitudes the value of γ needs to be very large and thus, there will always be a wave with finite amplitude propagating after the part of the chain with attached oscillators.

Plots of the transmitted amplitude for chains with 400 attached oscillators are shown in Figure 8. A chain with variable non-linearities is chosen with γ equal to 0.12 exp((i -1)/100), and for comparison the results for chains with

A c c e p t e d m a n u s c r i p t

  propagation properties of the chain. The analytical results are compared with the ones obtained by numerical simulations, some of which are presented in[16].

Fig. 1 .

 1 Fig. 1. Spring mass system with attached local oscillators.

  Plots of the real and the imaginary part of µ u for different frequencies are shown in Figure 2. Around the natural frequency of the attached oscillators there exists a frequency band where waves are attenuated. Outside the band gap, two neighbour masses j and j +1 oscillate with phase difference Im(µ u ). Inside the band gap, below the natural frequency κ of the attached oscillators, two neighbour masses oscillate with an opposite phase Im(µ u ) = π and above κ they oscillate with phase difference Im(µ u ) = 0. The real part Re(µ u ) gives the attenuation rate of the wave amplitude. At ω = κ the real part of µ u is infinity and decreases to zero outside the band gap. It should be pointed out that the system studied in this section has no energy dissipation mechanism.

Fig. 2 .

 2 Fig. 2. Dispersion relation for κ = 1.0 and β = 0.1.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Dispersion relation for β = 0.1, γ = 0 and different damping ratios.

A c c e p t e d m a n u s c r i p t written as

  ) resembles the equation for the amplitude of the stationary response of externally excited Duffing oscillator. By using polar representation for the amplitudes A 1,j+1 = R A (cos(ψ) + i sin(ψ)) and R B = |B 1,j+1 |, (15) can be

A c c e p t e d m a n u s c r i p t

  sorbing boundary conditions are applied at each end of the chain by adding dash-pot with a damping constant c = 1. The damping constant corresponds to the exact absorbing boundary condition for a continuous rod with distributed mass ρA = 1 and stiffness EA = 1, where ρ is the mass density, A is the section area and E is the elastic modulus. A wave travelling in the right direction is generated by applying a harmonic force with frequency ω and unit amplitude at the left end of the chain. The first and the last 10 masses in the chain are without attached oscillators. The wave travels undisturbed, until it reaches the part with the attached oscillators, where a part of it is reflected back, and a part of it propagates until it reaches the right end of the chain.

Fig. 6 .

 6 Fig. 6. Comparison of the theoretical prediction (dotted curves) of the transmitted amplitude with the one obtained by numerical integration for different values of γ = 0.0, 0.1, 0.25, 0.5, 1.0, 1.5, β = 0.1, ζ = 0.01. The RMS value of the amplitude of the first mass is 0.25. γ. The linear normalised frequency κ of the attached oscillators is 0.05. For small values of the non-linear coefficient the estimated upper and lower values of the absolute value of the amplitude bound well the estimated RMS value of the amplitude from the numerical simulations. For large values of the nonlinearity the results start to deviate. The frequency shift in the band gap can

Fig. 7 .A c c e p t e d m a n u s c r i p t 4 . 2

 742 Fig. 7. Comparison of the theoretical prediction (dotted curve) of the transmitted amplitude with the one obtained by numerical integration for different numbers of the attached oscillators n l = 250, 500, 1000, 1500, 2000, 2500 and γ = 0.3.

Fig. 8 .A c c e p t e d m a n u s c r i p t 5 Conclusions

 85 Fig. 8. Comparison of the transmitted amplitude for variable (solid line) and fixed (circular dots γ = 0.3) γ for chain with 400 attached oscillators. The two dashed lines show the transmitted amplitude for fixed γ (0.12, 0.88) equal to the lowest and the highest value used for the chain with varying non-linearities. The second plot is obtained by using the analytical calculations with variable γ and γ = 0.3. constant γ = 0.12 corresponding to the minimal value of γ in the chain with variable non-linearities, γ = 0.88 (maximal value) and γ = 0.3 (intermediate value), are shown in the figure. The input RMS value of the amplitude is equal to 0.25, the damping coefficient is ζ = 0.02 and i is the number of the attached oscillators. The theoretical prediction is calculated by using the procedure in 3.3. In all cases with γ = 0, a shift in the stop band frequency can be clearly observed. The shift of the chain with variable γ is between the shift for the cases with γ = 0.12 and γ = 0.88, and close to the case with the intermediate value γ = 0.3. The theoretical prediction indicates better filtering properties in the case with variable γ and this result is also supported by the numerical simulations.

Acknowledgement

This work was supported by grant 274-05-0498 from the Danish Research Council for Technology and Production Sciences. The authors wish to thank Professor Jon Juel Thomsen for his suggestions and valuable discussions.