Stochastic Stability of Damped Mathieu Oscillator Parametrically Excited by a Gaussian Noise

Abstract

This paper analyzes the stochastic stability of a damped Mathieu oscillator subjected to a parametric excitation of the form of a stationary Gaussian process, which may be both white and coloured. By applying deterministic and stochastic averaging, two Itô’s differential equations are retrieved. Reference is made to stochastic stability in moments. The differential equations ruling the response statistical moment evolution are written by means of Itô’s differential rule. A necessary and sufficient condition of stability in the moments of order r is that the matrix 𝐀_𝑟 of the coefficients of the ODE system ruling them has negative real eigenvalues and complex eigenvalues with negative real parts. Because of the linearity of the system the stability of the first two moments is the strongest condition of stability. In the case of the first moments (averages), critical values of the parameters are expressed analytically, while for the second moments the search for the critical values is made numerically. Some graphs are presented for representative cases.

1. Introduction

The so-called Mathieu equation has attracted the attention of the scholars in the past, and many papers and books are available on this subject such as [1, 2]. The interest in this equation stems from two facts: first, its dynamics is very rich; second, it describes the vibrations of important mechanical systems such as a prismatic bar stretched or compressed by a sinusoidal axial force, a pendulum whose support is subjected to a sinusoidal motion, and an elliptic membrane [3, 4]. In general, attention is focused on the stability of the motion.

When the excitation acting on the Mathieu oscillator is a stochastic process, its response is a stochastic process too. Thus, the analyst is concerned with the problem of the statistical characterization of the response, and in the case of a system prone to instability the problem of the stochastic stability must be solved. Nevertheless, studies on the Mathieu systems with stochastic excitation are not numerous. Probably, the first authors who addressed the problem were Stratonovich and Romanovski [5]. Analyses of the Mathieu-type equations with external random excitation were presented by Dimentberg [6] and by Cai and Lin [7].

A Mathieu system with random parametrical excitation was considered by Rong et al. [8] and by Xie [9] (see the references cited herein too). The first authors use the method of multiple scale to determine the equations of the amplitude and the phase of the response. The check for almost-sure stability is made by computing the largest Lyapunov exponent. On the contrary, in ([9, Chapter 7]) the method of stochastic averaging is used, and the stability in moments is analyzed. It recalled that there are different definitions of stochastic stability, for which the reader is referred to [10, 11]. According to the selected stability definition, different stability bounds are found. Thus, it is not surprising that the results of [8, 9] are different. However, the conclusions that are reached are opposite: in [8] it is claimed that the stochastic parametric excitation is always detrimental with respect to the stability, while in [9] it is affirmed that there are cases in which the stochastic excitation stabilizes the system.

In this paper, the Mathieu damped oscillator parametrically excited by a Gaussian noise is considered again. A suitable coordinate transformation is made, and deterministic and stochastic averaging are applied in such a way that the second-order motion equation is replaced by two first-order equations, are interpreted in Itô’s sense. As the Fokker-Planck-Kolmogorov equation for the response probability density function (PDF) is not analytically solvable, the stochastic stability in moments is considered as in [9, 12]: the differential equations governing the response moment evolution are written by means of the rules of Itô’s stochastic differential calculus [13–16]. A necessary and sufficient condition of stability in the moments of order 𝑟 is that the matrix 𝐀_𝑟 of the coefficients of the ODE system governing them has negative real eigenvalues and complex eigenvalues with negative real parts. Because of the linearity of the system the stability of the first two moments is the strongest condition of stability so that only the first two moments are considered ([11, Chapter 7]). Analytical stability bounds are established for the first moments, while the stability of the second moments is searched numerically.

2. Statement of the Problem

Consider the Mathieu-type stochastic differential equation𝑋̈(𝑡)+2𝜀𝜁₀𝜔₀𝑋̇(𝑡)+𝜔²₀1+𝜀𝛽sinΩ𝑡+√𝜀𝐹(𝑡)𝑋(𝑡)=0, (2.1) where 𝜁₀ is the ratio of critical damping, 𝜔₀ is the undamped pulsation of the oscillator, 𝜀 is a small parameter, β and Ω are the amplitude and the pulsation of the sinusoidal term, respectively, and F(t) is a stochastic stationary Gaussian process with zero mean. Because of the parametric excitation F(t) the response X(t) is a stochastic process too.

The pulsations Ω and 𝜔₀ can be linked by writing𝜔²₀=Ω²4+𝜀𝛿, (2.2) where δ is a detuning parameter. By inserting (2.2) in (2.1), this becomes𝑋̈(𝑡)+Ω²4𝑋(𝑡)=−𝜀2𝜁₀𝜔₀𝑋̇(𝑡)+𝛿+𝜔²₀𝛽sinΩ𝑡𝑋(𝑡)−√𝜀𝜔²₀𝑋(𝑡)𝐹(𝑡). (2.3)

The problem of the statistical characterization of (2.3) will be tackled by using the stochastic averaging method [6, 11, 17–20]. However, in the case of a harmonic term like that in (2.3), the classic coordinate transformation 𝑋(𝑡)=𝐴(𝑡)cos𝜙(𝑡), where 𝐴(𝑡) and 𝜙(𝑡) are stochastic processes, does not work well as the joint probability density function (PDF) of the response is not separable ([11, chapter 7]). A suitable coordinate transformation was proposed in [6, 17, 18], which reads as𝑋(𝑡)=𝑋₁(𝑡)cosΩ𝑡2+𝑋₂(𝑡)sinΩ𝑡2. (2.4) According to the principles of deterministic averaging for weakly nonlinear systems [21] it is assumed that 𝑋̇(𝑡)=−12Ω𝑋₁(𝑡)sinΩ𝑡2+12Ω𝑋₂(𝑡)cosΩ𝑡2𝑋̈(𝑡)=−12Ω𝑋̇₁sinΩ𝑡2−14Ω²𝑋₁cosΩ𝑡2+12Ω𝑋̇₂cosΩ𝑡2−14Ω²𝑋₂sinΩ𝑡2. (2.5)

The expressions in (2.5) are not exact but only approximate. The exact expression of the first derivative is −(1/2)Ω𝑋₁sin(Ω𝑡/2)+𝑋̇₁cos(Ω𝑡/2)+(1/2)Ω𝑋₂cos(Ω𝑡/2)+𝑋̇₂sin(Ω𝑡/2). According to the method of deterministic averaging [21] the second and fourth terms are neglected so that the derivative retains the same form as it would have if 𝑋₁ and 𝑋₂ were constant.

By inserting (2.4), (2.5) into (2.3), after some algebra we obtain a pair of first-order stochastic differential equations:𝑋̇₁(𝑡)=−𝜀2𝜁₀𝜔₀𝑋₁sin²Ω𝑡2−12𝑋₂−𝛿sinΩ𝑡Ω𝑋₁sinΩ𝑡−2𝛿Ω𝑋₂sin²Ω𝑡2−𝜔₀𝛽2𝑋₁sin²Ω𝑡−𝜔₀𝛽𝑋₂sin²Ω𝑡2+√sinΩ𝑡𝜀𝜔₀12𝑋₁sinΩ𝑡+𝑋₂sin²Ω𝑡2𝐹(𝑡),𝑋̇₂(𝑡)=𝜀2𝜁₀𝜔₀12𝑋₁sinΩ𝑡−𝑋₂cos²Ω𝑡2−𝛿Ω𝑋₂sinΩ𝑡−2𝛿Ω𝑋₁cos²Ω𝑡2−𝜔₀𝛽2𝑋₂sin²Ω𝑡−𝜔₀𝛽𝑋₁cos²Ω𝑡2−√sinΩ𝑡𝜀𝜔₀𝑋₁cos²Ω𝑡2+𝑋₂2sinΩ𝑡𝐹(𝑡). (2.6)(2.7)

The method of stochastic averaging is applied to (2.6), (2.7). In the field of stochastic dynamics this method of analysis was proposed first by Stratonovich [17, 18] (see also [11, 19, 20]) as an extension to stochastic systems of the deterministic method by Bogoliubov and Mitropolsky [21] (see Mettler [22], too). Then, it was rigorously demonstrated by Khasminskii [23]. It involves two phases that can be performed in either order: in the former, the deterministic terms that do not contain the forcing functions are averaged. In the second phase, the terms containing the forcing functions are reduced to Gaussian stationary white noises.

The deterministic averaging of the terms in brace brackets in (2.6), (2.7) requires the evaluation of integrals like Θ⁻¹⋅∫^Θ₀(•)d𝑡 [6, 17, 18], where Θ=2𝜋/Ω. It is obtained that𝑋̇₁(𝑡)=−𝜀𝜁₀𝜔₀𝑋₁−𝛿Ω𝑋₂−𝜔₀𝛽4𝑋₁+√𝜀𝜔₀12𝑋₁sinΩ𝑡+𝑋₂sin²Ω𝑡2𝐹(𝑡),𝑋̇₂(𝑡)=−𝜀𝜁₀𝜔₀𝑋₂+𝛿Ω𝑋₁+𝜔₀𝛽4𝑋₂−√𝜀𝜔₀𝑋₁cos²Ω𝑡2+12𝑋₂sinΩ𝑡𝐹(𝑡). (2.8)

Now, the forcing terms must be worked out. In the classic stochastic averaging method it is required that the process F(t) is broadbanded: in order to remove this restriction, a different way will be followed. Stratonovich [17, Volume 1, Chapter 7] suggested that, given a stochastic process F(t), which may be even narrowbanded, the product sinΩ𝑡𝐹(𝑡) is replaceable by a stationary Gaussian white noise 𝑊₁(t) having the autocorrelation function 𝐾₀𝛿(𝜏). The constant 𝐾₀ is given by the following integral: 𝐾₀=1/2∫^+∞_−∞𝑅_𝐹𝐹(𝜏)cosΩ𝜏d𝜏=𝜋𝑆_𝐹𝐹(Ω), where 𝑅_𝐹𝐹(𝜏) and 𝑆_𝐹𝐹(Ω) are the auto-correlation function and the power spectral density (PSD) of 𝐹(𝑡), respectively. However, the forcing terms in (2.8) contain the contributions 𝐹(𝑡)sin²(Ω𝑡/2),𝐹(𝑡)cos²(Ω𝑡/2) too. These contributions are replaced by a stationary Gaussian white noise 𝑊₂(t), whose intensity is computed by adapting to the present case the derivation of Stratonovich [17, Volume 1, Section 7.2]. We have𝐸𝑊₂(𝑡)𝑊₂(𝑡+𝜏)=𝛿(𝜏)^+∞_−∞𝑅_𝐹𝐹(𝜏)sin²Ω𝑡2sin²Ω2(𝑡+𝜏)d𝜏. (2.9) By expanding the product of the trigonometric functions in the integral (2.9), we obtain: sin²(Ω𝑡/2)sin²[Ω/2(𝑡+𝜏)]=1/4{1−cos[Ω(𝑡+𝜏)]−cosΩ𝑡+0.5cosΩ𝜏+0.5cos[Ω(2𝑡+𝜏)]}. In the last expression, the second and the third addenda, when averaged, give rise to zero, while the fifth term is a faster oscillatory one that according to Stratonovich can be neglected. Thus, the integral (2.9) is equal to^+∞_−∞𝑅_𝐹𝐹(𝜏1)4+18𝜋cosΩ𝜏d𝜏=2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω). (2.10) The expansion of ∫^+∞_−∞𝑅_𝐹𝐹(𝜏)cos²(Ω𝑡/2)cos²[Ω/2(𝑡+𝜏)]d𝜏 leads to the same result. In much the same way, it results that 𝐸[𝑊₁(𝑡)𝑊₂(𝑡+𝜏)]=0, that is, they are uncorrelated.

Performing these operations, (2.8) simplify into𝑋̇₁(𝑡)=−𝜀𝜔₀𝜁₀−𝛽4𝑋₁+𝜀𝛿Ω𝑋₂+𝑋₁𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)𝑊₁(𝑡)+𝑋₂𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)𝑊₂(𝑡),𝑋̇₂(𝑡)=−𝜀𝜔₀𝜁₀+𝛽4𝑋₂−𝜀𝛿Ω𝑋₁+𝑋₂𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)𝑊₁(𝑡)+𝑋₁𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)𝑊₂(𝑡). (2.11) For the sake of simplicity, the white noises in the previous equations have unit intensities. More-over, the forcing terms in the second equation have the sign plus as the white noise processes ±√𝐾𝑊(𝑡) have the same probabilistic characteristics ([24, chapter 6, section 6.3.2]).

In order to transform (2.11) into two Itô-type stochastic differential equations, the so-called Wong-Zakai-Stratonovich corrective terms must be added to the drift terms [25, 26]. They are computed according to the following formula:𝑚_𝑐𝑖=12^𝑛₁𝑘^𝑚_1𝑗𝑔_𝑘𝑗𝜕𝑔_𝑖𝑗𝜕𝑧_𝑘, (2.12) where in the present case n = m = 2 and 𝑔_𝑘𝑗, 𝑔_𝑖𝑗 are elements of the diffusion matrix𝑋𝐆=¹𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)𝑋₂𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω𝑋₂)𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)𝑋₁𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω). (2.13) Hence, the two Itô stochastic differential equations that govern the problem ared𝑋₁=−𝜀𝜔₀𝜁₀−𝛽4𝜋+𝜀4𝜔²₀𝑆𝑋₁𝛿d𝑡+𝜀Ω𝑋₂d𝑡+𝑋₁𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)d𝐵₁(𝑡)+𝑋₂𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)d𝐵₂(𝑡),d𝑋₂=−𝜀𝜔₀𝜁₀+𝛽4𝜋+𝜀4𝜔²₀𝑆𝑋₂𝛿d𝑡−𝜀Ω𝑋₁d𝑡+𝑋₂𝜔₀𝜀𝜋4𝑆_𝐹𝐹(Ω)d𝐵₁(𝑡)+𝑋₁𝜔₀𝜀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)d𝐵₂(𝑡), (2.14) where 𝑆=𝑆_𝐹𝐹(0)+𝑆_𝐹𝐹(Ω) and 𝐵₁(t) and 𝐵₂(t) are two standard Wiener processes (Brownian motion), for which the formal relationship d𝐵_𝑖/d𝑡=𝑊_𝑖(𝑡)(𝑖=1,2) holds.

From a theoretical point of view the probabilistic characterization of the random vector X = {𝑋₁,𝑋₂}^𝑡 should be obtained by solving the so-called Fokker-Planck-Kolmogorov (FPK) equation in the joint PDF 𝑝_𝐗 of the state variables [11, 17, 18, 24]𝜕𝑝_𝐗𝜕𝜕𝑡=−𝜕𝑥_𝑖𝑚_𝑖𝑝_𝐗+12𝜕²𝜕𝑥_𝑖𝜕𝑥_𝑗𝑏_𝑖𝑗𝑝_𝐗, (2.15) where the summation rule with respect to a repeated index has been adopted, and 𝑚₁=[−𝜀𝜔₀(𝜁₀−𝛽/4)+𝜀(𝜋/4)𝜔²₀𝑆]𝑋₁+(𝜀𝛿/Ω)𝑋₂,𝑚₂=[−𝜀𝜔₀(𝜁₀+𝛽/4)+𝜀(𝜋/4)𝜔²₀𝑆]𝑋₂−(𝜀𝛿/Ω)𝑋₁,𝑏_𝑖𝑗=⌊𝐆𝐆^𝑡⌋_𝑖𝑗, with the matrix 𝐆 being defined in (2.13), while the apex “t” denotes transpose. Unfortunately, as in many other cases the FPK equation (2.15) does not have an analytical solution because the excitation is a parametric one. Thus, another method has to be chosen in order to characterize the probabilistic characteristics of the response. The moment equation approach is used herein.

3. Moment Equation Approach

In (2.14) the excitation is multiplicative only, and there is not an external excitation. Thus, zero solution 𝑋₁ = 𝑋₂ = 0 for all t satisfies them, even if the multiplicative excitation is a stochastic process. We are concerned with the stability of zero solution. There exist different definitions of stochastic stability, for which the reader is referred to Chap. 6 of [11]. Here, the stability of zero solution is analyzed in the moments. In fact, the stability of the first two moments is the strongest for linear autonomous systems under multiplicative Gaussian excitation [11, 27]. Thus, the stability of the first two moments only will be considered here. As the response of (2.14) is not Gaussian, in order to characterize it statistically, from a theoretical point of view the knowledge of the infinite hierarchy of the moments would be necessary. Nevertheless, the system of (2.14) is linearly parametric: it has been shown that the equations for the statistical moments of such a type of systems are a close set [28], and they can be solved in succession. Thus, the convergence of the moments of first and second order is a necessary condition for being stable the moments of higher orders.

In the field of dynamic stability analyses, a system is stable when it comes back to the initial configuration after being subjected to a small perturbation. Zero statistical moments 𝐸⌊𝑋^𝑟1₁𝑋^𝑟2₂⌋(𝑟₁+𝑟₂=𝑟) correspond to a zero solution. When the zero solution is perturbed, the response moments are no longer zero: if the stochastic Mathieu oscillator is stable, once the perturbation is removed, they decay to zero, otherwise they grow without limits. The first step is to write the ordinary differential equations (ODEs) ruling the time evolution of the response moments. Since (2.14) is an Itô system, use is made of Itô’s differential rule [13–15], which reads asd𝜓=𝜕𝜓+𝜕𝑡𝜕𝜓𝜕𝑥_𝑖d𝑥_𝑖+12𝜕²𝜓𝜕𝑥_𝑖𝜕𝑥_𝑗d𝑥_𝑖d𝑥_𝑗. (3.1) In (3.1) Ψ is a nonanticipating function of the state variables 𝑥_𝑖 and the summation rule is adopted. It is recalled that (3.1) retains the terms of order dt only, as dB is of order√d𝑡. In order to write the moment equations, the appropriate nonanticipating function to be chosen is 𝜓=𝑋^𝑖₁𝑋^𝑗₂, where i, j are positive integers or zero with the constraint 𝑖+𝑗=𝑟, with 𝑟 being the order of the moments to be computed.

By choosing 𝜓=𝑋^𝑖₁𝑋^𝑗₂ with i + j = 1, the ODEs for the first order moments are𝐸̇𝑋₁=−𝜀𝜔₀𝜁₀−𝛽4−𝜔²₀𝜋4𝑆𝐸𝑋₁+𝜀𝛿Ω𝐸𝑋₂,𝐸̇𝑋₂=−𝜀𝜔₀𝜁₀+𝛽4−𝜔²₀𝜋4𝑆𝐸𝑋₂−𝜀𝛿Ω𝐸𝑋₁, (3.2) where the dots mean derivatives with respect to time and 𝑆=𝑆_𝐹𝐹(0)+𝑆_𝐹𝐹(Ω). Analogously, the ODEs ruling the second moments are found to be𝐸̇𝑋²₁=−2𝜀𝜔₀𝜁₀−𝛽4+𝜀𝜔²₀𝜋2𝑆_𝐹𝐹(0)+𝜀𝜔²₀32𝜋𝑆_𝐹𝐹(Ω)𝐸𝑋²₁+2𝜀𝛿Ω𝐸𝑋₁𝑋₂+𝜀𝜔²₀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)𝐸𝑋²₂=𝑎₁₁𝐸𝑋²₁+𝑎₁₂𝐸𝑋₁𝑋₂+𝑎₂₂𝐸𝑋²₂,𝐸̇𝑋₁𝑋₂=−𝜀𝛿Ω𝐸𝑋²₁+−2𝜀𝜔₀𝜁₀+𝜀𝜔²₀𝜋2𝑆_𝐹𝐹(0)+𝑆_𝐹𝐹(Ω)𝐸𝑋₁𝑋₂+𝜀𝛿Ω𝐸𝑋²₂=𝑎₂₁𝐸𝑋²₁+𝑎₂₂𝐸𝑋₁𝑋₂+𝑎₂₃𝐸𝑋²₂,𝐸̇𝑋²₂=𝜀𝜔²₀𝜋2𝑆_𝐹𝐹(0)+𝜋4𝑆_𝐹𝐹(Ω)𝐸𝑋²₁−2𝜀𝛿Ω𝐸𝑋₁𝑋₂+−2𝜀𝜔₀𝜁₀+𝛽4+𝜀𝜔²₀𝜋2𝑆_𝐹𝐹(0)+34𝜋𝑆_𝐹𝐹(Ω)𝐸𝑋²₂=𝑎₁₁𝐸𝑋²₁+𝑎₁₂𝐸𝑋₁𝑋₂+𝑎₂₂𝐸𝑋²₂. (3.3)

By inspecting (3.3), it is noted that the forcing terms are absent, which due to the fact that in (2.14) the excitation is purely parametric. Thus, (3.3) can be written in compact matrix form as𝐦̇_𝑟(𝑡)=𝐀_𝑟𝐦_𝑟, (3.4) where 𝐦_𝑟(t) is a vector collecting all the moments of order 𝑟 of the system states and 𝐀_𝑟 is a matrix of the coefficients 𝑎_𝑖𝑗 as shown in (3.3). The solution to (3.4) is𝐦_𝑟(𝑡)=𝐦₀exp𝐀_𝑟𝐭, (3.5) where 𝐦₀ is a vector whose entries are the initial conditions for the moments. These constitute the perturbation to zero solution.

It is well known (e.g., see [4]) that, as 𝑡 grows, the response moments decay to zero whenever the matrix 𝐀_𝑟 has negative real eigenvalues and complex eigenvalues with negative real parts. Otherwise, the moments increase without limits. Thus, the condition of stability in moments is that the eigenvalues of the matrix 𝐀_𝑟 have negative real parts. In this way, the stochastic problem is led to the classic deterministic problem of studying the eigenvalues of a matrix. Since the matrix 𝐀_𝑟 depends on the system parameters (2.1), there exist critical values of these quantities for which the real part of almost an eigenvalue is zero. Increasing them further, a real part becomes positive, and the moments grow to infinity.

In order to study the eigenvalues of 𝐀_𝑟, its characteristic equation is formed:det𝜆𝐈_𝑟−𝐀_𝑟=0, (3.6) where 𝐈_𝑟 is a unit matrix having the same order as 𝐀_𝑟. Equation (3.6) must specialized for the order 𝑟 of the moments. As previously advanced, due to the linearity of the system and to the Gaussianity of the input, the moment stability analysis is limited to the first two moments.

Coming back to the first-order moments, from (3.2) the characteristic equation looks like𝜔det𝜆+𝜀₀𝜁₀−𝛽4−𝜔²₀𝜋4𝑆−𝜀𝛿Ω𝜀𝛿Ω𝜆+𝜀𝜔₀𝜁₀+𝛽4−𝜔²₀𝜋4𝑆=0, (3.7) where from (2.2) 𝛿=𝜀⁻¹(𝜔²₀−Ω²/4). The roots of (3.7) are𝜆_1,2=𝜀−𝜁₀𝜔₀+14𝜋𝜔²₀1𝑆±4𝜔²₀𝛽²𝛿−16²Ω². (3.8) In examining the eigenvalues given by (3.8), it is necessary to distinguish between the case of real eigenvalues and that of complex conjugate ones. The eigenvalues are real numbers when𝛿𝛽>4𝜔₀Ω, (3.9) where for the sake of simplicity β is assumed to be positive. Otherwise, they are complex conjugate numbers. Clearly, in the latter case the stability condition for the first moments is𝜁₀>𝜋4𝜔₀𝑆. (3.10) Equation (3.10) requires that the oscillator is damped and the amount of damping depends on both the oscillator frequency and the intensity of the exciting noise. This requirement may be rather restrictive when 𝜔₀ is not small.

Now, let us consider the case in which (3.9) is satisfied. From (3.8) it is seen that there is passage to instability when the eigenvalue with the sign plus before the square root becomes zero, that is,−𝜁₀𝜔₀+14𝜋𝜔²₀1𝑆+4𝜔₀𝜔²₀𝛽²𝛿−16²Ω²=0. (3.11) Solving (3.11) with respect to 𝑆, we obtain the critical value of this quantity:𝑆_cr=4𝜁₀𝜋𝜔₀−1𝜋𝜔₀𝛽²𝛿−16²𝜔²₀Ω², (3.12) where β must be larger than the value in the right-hand side of (3.9). By comparing (3.12) with (3.10) it is found that the critical value given by the former is always smaller than that given by the latter. Keeping into account that 𝑆=𝑆_𝐹𝐹(0)+𝑆_𝐹𝐹(Ω), the critical intensity of the excitation depends on the form of power spectral intensity of this: in Section 4 the cases of white and coloured noises will be considered.

As regards the second moment stability, since the moment equations are three, (3.3), the characteristic equation associated with the matrix 𝐀₂ is of the third order. The roots of such a type of equation have analytical expressions. However, they are rather cumbersome as many parameters enter them. Thus, it has been preferred to proceed numerically. By using a computer algebra software a parameter is varied till one out of the three eigenvalues becomes zero: the corresponding value of the parameter is the critical one. Routh-Hurwitz criteria [4, 29] are not used. This is why in another study devoted to the stability of elastic columns with memory-dependent damping [12] it has been found that these criteria may overevaluate the critical values by a 30%.

4. Stability Analyses

The present section is devoted to the applications of the theory previously explained. Three cases are considered for the parametric excitation 𝐹(𝑡): (1) stationary Gaussian white noise; (2) coloured Gaussian process obtained by passing a stationary Gaussian white noise through a first-order linear filter; (3) coloured Gaussian process obtained by passing a stationary Gaussian white noise through a second-order linear filter. This choice makes the comparisons easier as the three stochastic processes have a common parameter that is the intensity σ² of the white noise.

The power spectral density (PSD) of a stationary stochastic process 𝑋(𝑡) is defined as𝑆_𝑋𝑋=12𝜋^+∞_−∞𝑅_𝑋𝑋(𝜏)⋅exp(−𝑖𝑤𝜏)d𝜏, (4.1) where 𝑅_𝑋𝑋(𝜏) is the autocorrelation function of the process and 𝑖=√−1. Since the autocorrelation function of a stationary Gaussian white noise 𝑊(𝑡) can be expressed as 𝑅_𝑊𝑊(𝜏)=2𝜋𝜎²𝛿(𝜏) (𝛿(𝜏) is Dirac’s delta), its PSD is 𝜎², constant on the whole real axis. Thus, for a white noise the first criterion of stability for the first-order moments (3.10) gives𝜎²_cr=2𝜁₀𝜋𝜔₀. (4.2) The critical value of σ² increases linearly with the system ratio of critical damping ζ₀ and is inversely proportional to system pulsation 𝜔₀, but it does not depend on Ω.

Now, the white noise is passed through a first-order linear filter (the Langevin equation) to have the excitation 𝐹(𝑡), that is,𝐹̇(𝑡)+𝛼𝐹(𝑡)=√2𝜋𝜎𝑊(𝑡), (4.3) where 𝑊(𝑡) has unit strength. The PSD of 𝐹(𝑡) is𝑆_𝐹𝐹(𝜔)=2𝜎²𝛼𝛼²+𝜔². (4.4) The critical value 𝜎²_cr is found by solving (3.10) with respect to 𝑆 and taking the expression of this quantity into account, that is,𝑆_cr=2𝜎²_cr𝛼+2𝜎²_cr𝛼𝛼²+Ω²=4𝜁₀𝜋𝜔₀. (4.5) The result of (4.5) is plotted in Figure 1 as a function of Ω for 𝛼=10, 𝜔₀=2𝜋, and 𝜁₀=0.05. It is found that 𝜎²_cris an increasing function of Ω.

Figure 1 

First moment stability: plot 𝜎2crof as a function of Ω from (4.5), first-order filter, 𝛼=10, 𝜔0=2π, 𝜁0=0.05.

The third type of excitation is obtained by passing the white noise through a second-order linear filter, that is,𝐹̈(𝑡)+2𝜁_𝑓𝜔_𝑓𝐹̇(𝑡)+𝜔²_𝑓𝐹(𝑡)=√2𝜋𝜎𝑊(𝑡). (4.6) In this case the PSD of 𝐹(𝑡) is𝑆_𝐹𝐹(𝜔)=2𝜎²𝜔²_𝑓−𝜔²²+4𝜁²_𝑓𝜔²_𝑓𝜔². (4.7)𝜎²_crfor the first condition of stability is computed from𝑆_cr=2𝜎²_cr𝜔⁴_𝑓+2𝜎²_cr𝜔²_𝑓−Ω²²+4𝜁²_𝑓𝜔²_𝑓Ω²=4𝜁₀𝜋𝜔₀. (4.8) The results of (4.8) are plotted in Figure 2 for 𝜔₀=2𝜋, 𝜁₀=0.05, and 𝜁_𝑓=0.05; the top plot refers to 𝜔_𝑓=𝜔₀ and the bottom plot to 𝜔_𝑓=2𝜔₀. Both plots show resonance with a marked minimum when Ω equates 𝜔_𝑓.

(a)

(b)

Figure 2 

First moment stability: plot of 𝜎2cras a function of Ω from (4.8), second-order filter: 𝜔0=2π, 𝜁0=0.05, and 𝜁𝑓=0.10; top plot 𝜔𝑓=𝜔0, bottom plot 𝜔𝑓=2𝜔0.

As regards the second condition of stability, 𝑆_cr from (3.12) is plotted in Figure 3 for 𝜔₀=2𝜋, 𝜁₀=0.05. Since this condition is valid only when 𝛽 is larger than the right-hand side of (3.9), the plots are drawn for 𝛽=𝜅𝛿/(𝜔₀Ω) with 𝜅=4.01,4.10,4.5,and5.0. Equation (3.12) gives a result with physical meaning only when 𝑆_cr is positive. From the plots—the bottom is a detail—it is seen that the interval of validity is small, and it narrows as κ increases. All the curves have a peak for Ω=2𝜔₀=4𝜋.

(a)

(b)

Figure 3 

First moment stability: plot of 𝑆cr as function of Ω according to (3.12) for 𝜔0=2π, 𝜁0=0.05, and𝛽 =𝜅𝛿/(𝜔0Ω) with 𝜅=4.01 (red line), 4.10 (black), 4.5 (blue), and 5.0 (violet).

Before presenting the results of the stability analyses for the second-order moments, it is necessary to compare the three PSD functions used in the analyses: they are depicted in Figure 4.The parameter 𝜎²is common to the three PDFs, and it constitutes the strength of the white noise. The plots clearly show as the white noise is by far the most severe excitation that may excite a dynamic system, if the colored excitations are obtained by passing it through a linear filter. As a consequence of this only observation, it cannot produce stabilizing effects on a dynamic system. The Wiener process yielded by the first-order filter of (4.3) (the Langevin equation) is the second in order of severity. As regards the processes generated by the second-order filter of (4.6), they have less strength as the abscissae of the peaks move farther from the origin. These remarks allow to explain the results of the analyses keeping in mind that the excitation affects the moments equations through 𝑆_𝐹𝐹(0) and 𝑆_𝐹𝐹(Ω). The following anal-yses are of two types: (1) search of the critical value of 𝜎² keeping β, the amplitude of the sinusoidal term in (2.1), constant; (2) search of the critical value of β keeping 𝜎² constant. For the sake of simplicity in all analyses the parameter 𝜀 is worth one. It is recalled that there is passage to instability when an eigenvalue becomes zero, which happens when 𝜎²equates 𝜎²_cr or β equates 𝛽_cr.

(a)

(b)

(c)

Figure 4 

(a) Comparison among the PSDs of the white noise (green line) and those given by (4.4) (black line) and by (4.7) (red and blue lines) for𝜎2=0.01,𝛼=10,𝜁𝑓=0.10,and𝜔𝑓=2𝜋or4𝜋.(b) Comparison between the PSDs deriving from (4.4) and (4.7): from (4.4) red line, from (4.7) blue and green lines. (c) Comparison between the PSDs deriving from (4.7): blue line 𝜔𝑓=2𝜋, black line𝜔𝑓=4𝜋.

For white noise excitation the plots of 𝜎²_cr and 𝛽_cras a function of Ω are in Figures 5 and 6, respectively, where 𝜔₀ = 2π and ζ₀ = 0.05. Both 𝜎²_cr and 𝛽_cr are constant with respect to Ω, being the stable regions below the straight lines. This result is not surprising: in fact, the excitation enters the moment equations through 𝑆_𝐹𝐹(0) and 𝑆_𝐹𝐹(Ω), which are equal and do not vary with Ω in the case of white noise, whose PSD is constant on the whole frequency axis. 𝜎²_cr increases as 𝛽decreases (Figure 5), and 𝛽_cr has a similar trend (Figure 6).

Figure 5 

Second moment stability, white noise excitation: plot of 𝜎2cr as a function of Ω for 𝜀=1, 𝜔0=2π, 𝜁0=0.05; 𝛽=0.1 green line, and 𝛽=0.15 black line.

Figure 6 

Second moment stability, white noise excitation: plot of 𝛽cras a function of Ω for 𝜀=1, 𝜔0=2π, 𝜁0=0.05; 𝜎2=0.0005green line, and 𝜎2=0.001 black line.

Figure 7 shows 𝐸[𝑋²₁] as a function of time for 𝛽=0.1, 𝜎²=0.0012 (top plot), and 𝜎²=0.00130 (bottom plot). The moment equations (3.3) are numerically integrated by means of a fourth-order Runge-Kutta method. The initial perturbation is given by means of the initial conditions 𝐸[𝑋²₁(0)]=𝐸[𝑋²₂(0)]=0.1. In the former case 𝜎² is lesser than 𝜎²_cr, and 𝐸⌊𝑋²₁⌋ starts from the prescribed value 0.1, then it decays to zero and does not show any oscillation. In the latter case𝜎² is larger than 𝜎²_cr, and 𝐸⌊𝑋²₁⌋ grows without limit. The absence of oscillations is due to the fact that in applying the stochastic averaging the oscillatory terms are cancelled.

(a)

(b)

Figure 7 

White noise excitation: 𝐸[𝑋21]=𝑧1 with respect to time, 𝜔0=2π, 𝜁0=0.05, 𝛽=0.1, 𝜀=1; (a) 𝜎2=0.00123<𝜎2cr; (b) 𝜎2=0.00130>𝜎2cr.

As regards the excitation given by the first-order filter of (4.3), the plots of 𝜎²_cr and 𝛽_cras a function of Ω are shown in Figures 8 and 9, respectively, where 𝜔₀=2𝜋, 𝜁₀=0.05, and 𝛼=10. Both quantities increase as Ω increases. In fact, the PSD (4.4) is a monotonically decreasing func-tion of Ω, which causes the excitation to diminish with Ω. The curve in the plot of Figure 8 is obtained by keeping β equal to 0.1. The increase of 𝜎²_cr is marked and reaches a 70% as Ω passes from 0.5π to 6.5π. The curve of 𝛽_cr refers to 𝜎²=0.0005: in the same interval the increase is less marked and amounts to about 28%. The analyses are made in a discrete series of values of Ω obtaining a set of points in the plane (σ², Ω) or (β, Ω). Then, the curves are traced by means of the routine Spline of MAPLE.

Figure 8 

Second moment stability, excitation given by the first-order filter (4.3): plot of 𝜎2cr as a function of Ω for 𝛽=0.1, 𝛼=10, 𝜀=1, 𝜔0=2𝜋,𝜁0=0.05.

Figure 9 

Second moment stability, excitation given by the first-order filter (4.3): plot of 𝛽cr as a function of Ω for 𝜎2=0.005, 𝛼=10, 𝜀=1, 𝜔0=2π, 𝜁0=0.05.

The plots for the excitation given by the second-order filter (4.6) are shown in Figures 10 and 11 and are obtained for this data set: 𝜀=1, 𝜔₀= 2π, 𝜁₀=0.05, 𝛽=0.10 in Figure 10 for 𝜎²_cr, 𝜎²=0.01 in Figure 11 for 𝛽_cr, 𝜁_𝑓=0.10, black line 𝜔_𝑓=2𝜋, blue line 𝜔_𝑓=4𝜋. In both plots the curve deriving from 𝜔_𝑓=4𝜋is the higher since the excitation strength is smaller, as already explained. In both plots there is a marked valley for Ω=𝜔_𝑓: since 𝑆_𝐹𝐹(Ω) assumes its largest value (see Figure 4), there is a kind of stochastic resonance [30], even if the system remains stable.

Figure 10 

Excitation from the second order filter (4.6): plots of 𝜎2cr as a function of Ω: 𝜀=1, 𝜔0=2π, 𝜁0=0.05, 𝛽=0.10, 𝜁𝑓=0.10, black line 𝜔𝑓=2π, blue line 𝜔𝑓=4π.

Figure 11 

Excitation from the second-order filter (4.6): plots of 𝛽cr as a function of Ω: 𝜀=1, 𝜔0=2π, 𝜁0=0.05, 𝜎2=0.01, 𝜁𝑓=0.10, black line 𝜔𝑓=2π, blue line 𝜔𝑓=4π.

5. Summary and Conclusions

In this paper the issue of stability of the stochastic Mathieu oscillator is addressed. In order to solve the problem, a suitable coordinate transformation and stochastic averaging are applied to the original dynamic system. In this way, two first-order stochastic differential equations are obtained. Then, they are transformed into two Itô-type stochastic equations. For such a type of stochastic differential equations Itô’s differential rule is applicable and allows one to derive the ODEs ruling the time evolution of the response statistical moments.

In the dynamic system considered here the excitation is merely linearly parametric. Hence, the moments ODEs are linear and homogeneous, and the ODEs for the moments of order 𝑟 are uncoupled from the ODEs for moments of different order. Because of this feature the ODEs for the moments of order 𝑟 admit zero solution, which means that the system is at rest, and the parametric excitation does not disturb this state. The system is stable in moments when these return to zero after the application of an external perturbation; otherwise, it is unstable and the moments diverge. The necessary and sufficient condition for the stability of the moments of order 𝑟 is that all the eigenvalues of the matrix 𝐀_𝑟 of the coefficients of the ODEs written for those moments have negative real part. The passage to the instability requires that an eigenvalue becomes zero. In this way, the stochastic problem is reduced to the deterministic problem of the study of the eigenvalues of a matrix.

The procedure developed in this paper is applied to three types of Gaussian stationary stochastic excitations: (1) white noise; (2) white noise filtered through a first-order filter; (3) white noise filtered through a second-order filter. In order to make it possible to compare the results, in cases (2.2) and (2.3) the source white noise is the same as that of case (1), and it is by far the more severe excitation (Figure 4).

As regards the stability of the first-order moments, it is substantially led by the system damping. Several parameters affect the stability of the second-order moments: since the characteristic equation of the matrix 𝐀₂ is of the third order, so that its roots have complicated expressions, it has been chosen to proceed numerically. The analyses were aimed at finding the critical values of the white noise strength σ² and of the amplitude β of the sinusoidal term in the motion equation (2.1). Differently from the deterministic Mathieu oscillator the pulsation Ω of this term has little or no effects on the bounds of stability, which is due to the presence of damping. Vice versa, in the case of excitation deriving from a second-order filter there is a stochastic resonance when Ω equates the peak frequency 𝜔_𝑓 of the PSD of the excitation. In any case, the critical values 𝜎²_crand 𝛽_crare by far smaller for the white noise excitation since as observed it is the strongest excitation. Thus, the statement of [9, 31] that a white noise excitation may cause a stabilizing effect is not found true. The results of these papers agree with those of [8]: a direct comparison of the results is not possible as the system studied in that reference is a little different and a quite different method of analysis is used there.