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

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 Ar 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 twomoments 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.


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

Statement of the Problem
Consider the Mathieu-type stochastic differential equation where ζ 0 is the ratio of critical damping, ω 0 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 ω 0 can be linked by writing where δ is a detuning parameter.By inserting 2.2 in 2.1 , this becomes 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 X t A t cos φ t , where A t and φ t 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 According to the principles of deterministic averaging for weakly nonlinear systems 21 it is assumed that

2.5
The expressions in 2.5 are not exact but only approximate.The exact expression of the first derivative is − 1/2 ΩX 1 sin Ωt/2 Ẋ1 cos Ωt/2 1/2 ΩX 2 cos Ωt/2 Ẋ2 sin Ωt/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 X 1 and X 2 were constant.
By inserting 2.4 , 2.5 into 2.3 , after some algebra we obtain a pair of first-order stochastic differential equations: sin Ωt F t .

4 Mathematical Problems in Engineering
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

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 Ωt F t is replaceable by a stationary Gaussian white noise W 1 t having the autocorrelation function K 0 δ τ .The constant K 0 is given by the following integral: , where R FF τ and S FF Ω are the auto-correlation function and the power spectral density PSD of F t , respectively.However, the forcing terms in 2.8 contain the contributions F t sin 2 Ωt/2 , F t cos 2 Ωt/2 too.These contributions are replaced by a stationary Gaussian white noise W 2 t , whose intensity is computed by adapting to the present case the derivation of Stratonovich 17, Volume 1, Section 7.2 .We have By expanding the product of the trigonometric functions in the integral 2.9 , we obtain: sin 2 Ωt/2 sin 2 Ω/2 t τ 1/4{1 − cos Ω t τ − cos Ωt 0.5 cos Ωτ 0.5 cos Ω 2t τ }.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 The expansion of ∞ −∞ R FF τ cos 2 Ωt/2 cos 2 Ω/2 t τ dτ leads to the same result.In much the same way, it results that E W 1 t W 2 t τ 0, that is, they are uncorrelated.
Performing these operations, 2.8 simplify into

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 ± √ KW t 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 socalled Wong-Zakai-Stratonovich corrective terms must be added to the drift terms 25, 26 .They are computed according to the following formula: where in the present case n = m 2 and g kj , g ij are elements of the diffusion matrix Hence, the two It ô stochastic differential equations that govern the problem are where S S FF 0 S FF Ω and B 1 t and B 2 t are two standard Wiener processes Brownian motion , for which the formal relationship dB i /dt W i t i 1, 2 holds.
From a theoretical point of view the probabilistic characterization of the random vector X {X 1 , X 2 } t should be obtained by solving the so-called Fokker-Planck-Kolmogorov FPK equation in the joint PDF p X of the state variables 11, 17, 18, 24 where the summation rule with respect to a repeated index has been adopted, and , with the matrix G 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.

Moment Equation Approach
In 2.14 the excitation is multiplicative only, and there is not an external excitation.Thus, zero solution X 1 X 2 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 E X r 1 1 X r 2 2 r 1 r 2 r 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 as In 3.1 Ψ is a nonanticipating function of the state variables x i and the summation rule is adopted.It is recalled that 3.1 retains the terms of order dt only, as dB is of order √ dt.In order to write the moment equations, the appropriate nonanticipating function to be chosen is ψ X i 1 X j 2 , where i, j are positive integers or zero with the constraint i j r, with r being the order of the moments to be computed.
By choosing ψ X i 1 X j 2 with i j 1, the ODEs for the first order moments are where the dots mean derivatives with respect to time and S S FF 0 S FF Ω .Analogously, the ODEs ruling the second moments are found to be

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 where m r t is a vector collecting all the moments of order r of the system states and A r is a matrix of the coefficients a ij as shown in 3.3 .The solution to 3.4 is where m 0 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 t grows, the response moments decay to zero whenever the matrix A r 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 A r have negative real parts.In Mathematical Problems in Engineering this way, the stochastic problem is led to the classic deterministic problem of studying the eigenvalues of a matrix.Since the matrix A r 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 A r , its characteristic equation is formed: where I r is a unit matrix having the same order as A r .Equation 3.6 must specialized for the order r 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 where from 2.2 δ ε −1 ω 2 0 − Ω 2 /4 .The roots of 3.7 are 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 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

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 ω 0 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, Solving 3.11 with respect to S, we obtain the critical value of this quantity: 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 S S FF 0 S FF Ω , 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 A 2 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%.

Stability Analyses
The present section is devoted to the applications of the theory previously explained.Three cases are considered for the parametric excitation F t : 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 σ 2 of the white noise.The power spectral density PSD of a stationary stochastic process X t is defined as where R XX τ is the autocorrelation function of the process and i √ −1.Since the autocorrelation function of a stationary Gaussian white noise W t can be expressed as R WW τ 2πσ 2 δ τ δ τ is Dirac's delta , its PSD is σ 2 , constant on the whole real axis.Thus, for a white noise the first criterion of stability for the first-order moments 3.10 gives The critical value of σ 2 increases linearly with the system ratio of critical damping ζ 0 and is inversely proportional to system pulsation ω 0 , 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 F t , that is, where W t has unit strength.The PSD of F t is The critical value σ 2 cr is found by solving 3.10 with respect to S and taking the expression of this quantity into account, that is, The result of 4.5 is plotted in Figure 1 as a function of Ω for α 10, ω 0 2π, and ζ 0 0.05.It is found that σ 2 cr is an increasing function of Ω.The third type of excitation is obtained by passing the white noise through a secondorder linear filter, that is, In this case the PSD of F t is The results of 4.8 are plotted in Figure 2 for ω 0 2π, ζ 0 0.05, and ζ f 0.05; the top plot refers to ω f ω 0 and the bottom plot to ω f 2ω 0 .Both plots show resonance with a marked minimum when Ω equates ω f .As regards the second condition of stability, S cr from 3.12 is plotted in Figure 3 for ω 0 2π, ζ 0 0.05.Since this condition is valid only when β is larger than the right-hand side of 3.9 , the plots are drawn for β κδ/ ω 0 Ω with κ 4.01, 4.10, 4.5, and 5.0.Equation 3.12 gives a result with physical meaning only when S 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ω 0 4π.
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 σ 2 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 secondorder 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 S FF 0 and S FF Ω .The following analyses are of two types: 1 search of the critical value of σ 2 keeping β, the amplitude of the sinusoidal term in 2.1 , constant; 2 search of the critical value of β keeping σ 2 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 σ 2 equates σ 2 cr or β equates β cr .
For white noise excitation the plots of σ 2 cr and β cr as a function of Ω are in Figures 5  and 6, respectively, where ω 0 2π and ζ 0 0.05.Both σ 2 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 S FF 0 and S FF Ω , which are equal and do   not vary with Ω in the case of white noise, whose PSD is constant on the whole frequency axis.σ 2 cr increases as β decreases Figure 5 , and β cr has a similar trend Figure 6 .
Figure 7 shows E X 2 1 as a function of time for β 0.1, σ 2 0.0012 top plot , and σ 2 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 E X 2  1 0 E X 2 2 0 0.1.In the former case σ 2 is lesser than σ 2 cr , and E X 2  1 starts from the prescribed value 0.1, then it decays to zero and does not show any oscillation.In the latter case σ 2 is larger than σ 2 cr , and E X 2 1 grows without limit.The absence of oscillations is due to the fact that in applying the stochastic averaging the oscillatory terms are cancelled.
As regards the excitation given by the first-order filter of 4.3 , the plots of σ 2 cr and β cr as a function of Ω are shown in Figures 8 and 9, respectively, where ω 0 2π, ζ 0 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 σ 2 cr is marked and reaches a 70% as Ω passes from 0.5π to 6.5π.The curve of β cr refers to σ 2 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 σ 2 , Ω or β, Ω .Then, the curves are traced by means of the routine Spline of MAPLE.
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, ω 0 2π, ζ 0 0.05, β 0.10 in Figure 10 for σ 2 cr , σ 2 0.01 in Figure 11 for β cr , ζ f 0.10, black line ω f 2π, blue line ω f 4π.In both plots the curve deriving from ω f 4π is the higher since the excitation strength is smaller, as already explained.In both plots there is a marked valley for Ω ω f : since S FF Ω assumes Mathematical Problems in Engineering its largest value see Figure 4 , there is a kind of stochastic resonance 30 , even if the system remains stable.

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 r are uncoupled from the ODEs for moments of different order.Because of this feature the ODEs for the moments of order r 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 r is that all the eigenvalues of the matrix A r 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 and 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 char-acteristic equation of the matrix A 2 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 σ 2 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 ω f of the PSD of the excitation.In any case, the critical values σ 2 cr and β cr are 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.

Figure 4 :
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, ζ f 0.10, and ω f 2π or 4π.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 ω f 2π, black line ω f 4π.