# Stochastic Stability Analysis of Coupled Viscoelastic Systems with Nonviscously Damping Driven by White Noise

Authors: Di Liu; Yanru Wu; Xiufeng Xie
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Nonviscously damped structural system has been raised in many engineering fields, in which the damping forces depend on the past time history of velocities via convolution integrals over some kernel functions. This paper investigates stochastic stability of coupled viscoelastic system with nonviscously damping driven by white noise through moment Lyapunov exponents and Lyapunov exponents. Using the coordinate transformation, the coupled Itô stochastic differential equations of the norm of the response and angles process are obtained. Then the problem of the moment Lyapunov exponent is transformed to the eigenvalue problem, and then the second-perturbation method is used to derive the moment Lyapunov exponent of coupled stochastic system. Lyapunov exponent also can be obtained according to the relationship between moment Lyapunov exponent and Lyapunov exponent. Finally, the effects of various physical quantities of stochastic coupled system on the stochastic stability are discussed in detail. These results are validated by the direct Monte Carlo simulation technique.


---

## Body

## 1. Introduction

Viscoelastic materials are widely used in aerospace, construction, textile, and other industries, because they have a series of excellent properties, including light weight, high strength, wide source, and good shock absorption. These materials’ stress depends on the past time history of stain; that is, the stress will increase correspondingly with the increase of time and the bucking will be more likely to occur. Therefore, the research of dynamical behavior of viscoelastic system has received a lot of interests in recent years [1–4].To better understand the dynamical behavior of the viscoelastic system, many methods had been put forward. McTavish et al. [5] presented GHM method to analyze the linear multiple degree of freedom viscoelastic damping systems. Li et al. [6] studied the decay rate of energy functional of nonlinear viscoelastic full Marguerre-von Kármán shallow shell system. Later, Li [7] also showed the existence and uniqueness of weak solution of this system by applying the Galerkin finite element method. Nieto and Ahmad [8] used the generalized quasi-linearization method to obtain the explicit approximate solutions of an initial and terminal value problem for the forced Duffing equation with nonviscous damping. Recently, Li and Du [9] studied general energy decay of a degenerate viscoelastic Petrovsky-type plate equation with boundary feedback through a priori estimate and analysis of Lyapunov-like functional.The dynamic stability as an important field of nonlinear dynamics has attached increasing attention in recent years. Guo et al. [10] derived a new simple sufficient condition of the global asymptotic stability of the integrodifferential systems with delay through constructing suitable Lyapunov functionals combined with the analytical technique. And then Meng et al. [11, 12] obtained a sufficient condition of uniform stability for nonlinear integrodifferential system. However, when we consider the dynamical behavior of the system in a real environment, the stochastic excitation cannot be ruled out, such as wind gusts, earthquakes, and ocean waves [13–15]. Stochastic stability of the system under the stochastic excitation has attracted more and more attention [16–20]. The Lyapunov exponent specially has been researched by many scholars as an important indication for judging the stability of stochastic system. For example, Potapov [21] derived the sufficient condition of the almost-sure asymptotic stability of elastic and viscoelastic systems by the Lyapunov’s direct method. He [22] also studied the stability of elastic and viscoelastic systems under the non-Gaussian excitation. For stochastic stability of high dimensional coupled system, Pavlović et al. [23] proposed the direct Lyapunov method to investigate the almost-sure stochastic stability of a viscoelastic double-beam system under parametric excitation. Then they used the same method to analyze the instability of coupled nanobeam systems subject to compressive axial loading [24]. These works only indicate that the system is almost-surely stable.According to the theory of large deviation, which was first proposed by Baxendale and Stroock [25], the almost-sure stability of stochastic system does not mean the moment stability. Therefore, it is necessary to study the moment Lyapunov exponent of a stochastic system. The moment Lyapunov exponent can judge not only almost-sure stability but also the moment stability [26]. Nevertheless, analytic solutions of the moment Lyapunov exponent are very difficult to derive. To overcome this, many scholars developed different approximation methods to study the moment Lyapunov exponents of the stochastic system. For example, Sri Namachchivaya and Van Roessel [27] studied the moment Lyapunov exponents of two-degrees-of-freedom coupled elastic oscillators under real noise excitation by combining an asymptotic approximation for the moment Lyapunov exponents. Kozić and his associates [28–30] developed the first-order perturbation approach to obtain weak noise expansion of moment Lyapunov exponents and Lyapunov exponents for a stochastically coupled double-beam system and Timoshenko beam system, respectively. Subsequently, Stojanović and Petković [31] used a perturbation method to study the moment Lyapunov exponents and the Lyapunov exponents of the three elastically connected Euler beams. More recently, Deng [32] applied the averaging method to establish the moment Lyapunov exponent of coupled gyroscopic stochastic system under bounded noise excitation. But it is rare to see studies of the moment Lyapunov exponent of a viscoelastic coupled system with nonviscous damping. Deng [33] investigated stochastic stability of two-degrees-of-freedom coupled viscoelastic system under white noise through moment Lyapunov exponent, but its damping term is viscous.The purpose of this paper is to study stochastic stability of a viscoelastic coupled system with nonviscous damping subject to Gaussian white noise excitation through moment Lyapunov exponents and Lyapunov exponents, in which the nonviscous damping term is expressed by Boltzmanns superposition integral with a hereditary relaxation kernel. The paper is organized as follows. Section2 derives the governing equations of motion of a stochastic viscoelastic system with a nonviscous damping. And then solving the problem of the moment Lyapunov exponent is changed into the eigenvalue and eigenfunction problem through a suitable transformation. In Section 3, the zeroth-order and the first-order and second-order perturbation solutions of the moment Lyapunov exponents are derived, respectively, based on the second-order perturbation method. The effects of different physical quantities on the stochastic stability of the coupled system are discussed under the given parameters. Correspondingly, the results are verified by means of the direct Monte Carlo simulation in Section 4. Finally, conclusions will be drawn in Section 5.


## 2. Problem Formulation

Considering a coupled viscoelastic system with nonviscous damped structure driven by white noise excitation, the governing equations can be expressed as(1)q¨1+ω12q1+ϵβ1∫0tγt-τq˙1τdτ+ϵH1q1-q2-ϵK1ξtq1=0,q¨2+ω22q2+ϵβ2∫0tγt-τq˙2τdτ+ϵH2q2-q1-ϵK2ξtq2=0,where q1 and q2 are generalized displacements. ω1 and ω2 are natural frequencies. ϵ is a small parameter. β1 and β2 are nonviscous damping coefficients. Hi and Ki are constants. ξt is a Gaussian white noise with zero mean and noise intensity σ2. γt is a relaxation function of nonviscous damping, which can be expressed by (2)γt=1αexp-tαin which α is the relaxation parameter. Note that the limit case α⟶0, the nonviscous damping force reduces to classical viscous damping.Letq1=x1, q˙1=ω1x2, q2=x3, q˙2=ω2x4; system (1) can be represented as(3)x˙1=ω1x2,x˙2=-ω1x1-ϵh1x1-x3-ϵζ1∫0tγt-τx˙1τdτ+ϵk1x1ξt,x˙3=ω2x4,x˙4=-ω2x3-ϵh2x3-x1-ϵζ2∫0tγt-τx˙3τdτ+ϵk2x3ξt,where(4)h1=H1ω1,h2=H2ω2,k1=K1ω1,k2=K2ω2,ζ1=β1ω1,ζ2=β2ω2Using a coordinate transformation,(5)x1=acos⁡φcos⁡θ1,x2=-acos⁡φsin⁡θ1,x3=asin⁡φcos⁡θ2,x4=-asin⁡φsin⁡θ2.where a represents the norm of the response. 0≤θ1≤2π and 0≤θ2≤2π are the angles process of the coupled stochastic system. 0≤φ≤π/2 denotes the coupling or exchange of energy between the coupled system. Equation (3) can be rewritten as(6)θ˙1=ω1+ϵζ1a1cos⁡φcos⁡θ1∫0tγt-τx˙1τdτ+ϵh121+cos⁡2θ1-2cos⁡θ1cos⁡θ2tan⁡φ-ϵk1cos2θ1ξt,θ˙2=ω2+ϵζ2a1sin⁡φcos⁡θ2∫0tγt-τx˙3τdτ+ϵh221+cos⁡2θ2-2cos⁡θ1cos⁡θ2cot⁡φ-ϵk2cos2θ2ξt,φ˙=ϵh1421-cos⁡2φsin⁡θ1cos⁡θ2-sin⁡2φsin⁡2θ1-ζ1asin⁡φsin⁡θ1∫0tγt-τx˙1τdτ+h24-21+cos⁡2φcos⁡θ1sin⁡θ2+sin⁡2φsin⁡2θ2+ζ2asin⁡φsin⁡θ2∫0tγt-τx˙3τdτ+ϵ4k1sin⁡2θ1-k2sin⁡2θ2sin⁡2φξt,a˙=ϵa4h11+cos⁡2φsin⁡2θ1-2sin⁡2φsin⁡θ1cos⁡θ2+ζ1acos⁡φsin⁡θ1∫0tγt-τx˙1τdτ+h21-cos⁡2φsin⁡2θ2-2sin⁡2φcos⁡θ1sin⁡θ2+ζ2acos⁡φsin⁡θ2∫0tγt-τx˙3τdτ-ϵ4ak11+cos⁡2φsin⁡2θ1+k21-cos⁡2φsin⁡2θ2ξtTherefore, the solution forms of θi can be expressed as θi=ωit+ϕit from (6), where θit is quickly varying process relative to the process ϕi due to ϵ≪1. Therefore, the following approximations will be obtained:(7)θit-s=ωit-s+ϕit-s≈θit-ωis,x˙1t-s=ω1x2t-s=-aω1cos⁡φsin⁡θ1t-s≈-aω1cos⁡φsin⁡θ1-ω1s=ω1x2cos⁡ω1s+ω1x1sin⁡ω1s,x˙3t-s=ω2x4t-s=-aω2cos⁡φsin⁡θ2t-s≈-aω2cos⁡φsin⁡θ2-ω2s=ω2x4cos⁡ω2s+ω2x3sin⁡ω2s.Substituting (7) into the nonviscous damping, one obtains(8)∫0tγt-τx˙1τdτ=1α∫0texp-1αsx˙1t-sds=M1x2+N1x1.Similarly, the following results also can be derived:(9)∫0tγt-τx˙3τdτ=M2x4+N2x3,where(10)M1=ω11+α2ω12,N1=αω121+α2ω12,M2=ω21+α2ω22,N2=αω221+α2ω22.Substituting (8) and (9) into (6) and letting P=ap=x12+x22+x32+x42p/2-∞<p<+∞, (6) can be written as the following Itô stochastic differential equations:(11)dθ1=m1θ1,θ2,φdt+σ11θ1,θ2,φdWt,dθ2=m2θ1,θ2,φdt+σ21θ1,θ2,φdWt,dφ=m3θ1,θ2,φdt+σ31θ1,θ2,φdWt,dP=Pm4θ1,θ2,φdt+Pσ41θ1,θ2,φdWt,where Wt is the standard Wiener process; the drift and diffusion coefficients can be rewritten as(12)m1θ1,θ2,φ=ω1-ϵζ12M1sin⁡2θ1-N11+cos⁡2θ1+ϵh121+cos⁡2θ1-2cos⁡θ1cos⁡θ2tan⁡φ-ϵk12σ282sin⁡2θ1+sin⁡4θ1,m2θ1,θ2,φ=ω2-ϵζ22M2sin⁡2θ2-N21+cos⁡2θ2+ϵh221+cos⁡2θ2-2cos⁡θ1cos⁡θ2cot⁡φ-ϵk22σ282sin⁡2θ2+sin⁡4θ2,m3θ1,θ2,φ=ϵh1421-cos⁡2φsin⁡θ1cos⁡θ2-sin⁡2φsin⁡2θ1+ζ14sin⁡2φM11-cos⁡2θ1-N1sin⁡2θ1-h2421+cos⁡2φcos⁡θ1sin⁡θ2-sin⁡2φsin⁡2θ2-ζ24sin⁡2φ×M21-cos⁡2θ2-N2sin⁡2θ2+σ264-4k121+2cos⁡2θ1+cos⁡4θ1sin⁡2φ+k121-cos⁡4θ1sin⁡4φ-4k1k2sin⁡2θ1sin⁡2θ2sin⁡4φ+4k221+2cos⁡2θ2+cos⁡4θ2sin⁡2φ+k221-cos⁡4θ2sin⁡4φ,m4θ1,θ2,φ=ϵp4h11+cos⁡2φsin⁡2θ1-2sin⁡2φcos⁡θ2sin⁡θ1-ζ11+cos⁡2φ×M11-cos⁡2θ1-N1sin⁡2θ1+h21-cos⁡2φsin⁡2θ2-2sin⁡2φcos⁡θ1sin⁡θ2-ζ21-cos⁡2φM21-cos⁡2θ2-N2sin⁡2θ2+σ2k123210+3p+16cos⁡2θ1-3p-2cos⁡4θ1+8+4p+16cos⁡2θ1-4p-2cos⁡4θ1cos⁡2φ+p-21-cos⁡4θ1cos⁡4φ+σ2k1k28p-2×1-cos⁡4φsin⁡2θ1sin⁡2θ2+σ2k223210+3p+16cos⁡2θ2-3p-2cos⁡4θ2-8+4p+16cos⁡2θ2-4p-2cos⁡4θ2cos⁡2φ+p-21-cos⁡4θ2cos⁡4φ,σ11θ1,θ2,φ=-ϵk1σcos2θ1,σ21θ1,θ2,φ=-ϵk2σcos2θ2,σ31θ1,θ2,φ=ϵσ4k1sin⁡2θ1-k2sin⁡2θ2sin⁡2φ,σ41θ1,θ2,φ=-ϵpσ4k11+cos⁡2φsin⁡2θ1+k21-cos⁡2φsin⁡2θ2.To obtain the eigenvalue problem for thepth moment Lyapunov exponents of a four-dimensional linear Itô stochastic system, a linear stochastic transformation [34] is introduced(13)S=Tθ1,θ2,φP,P=T-1θ1,θ2,φS,where S is a new norm process depending on the transformation function Tθ1,θ2,φ. The transformation function Tθ1,θ2,φ is defined on the independent stationary phase processes θ1, θ2, and φ in the ranges 0≤θ1≤2π, 0≤θ2≤2π, 0≤φ≤π/2. Applying Itô formula, one can obtain(14)dS=P12σ112Tθ1θ1′′+σ11σ21Tθ1θ2′′+σ11σ31Tθ1φ′′+12σ212Tθ2θ2′′+σ21σ31Tθ2φ′′+12σ312Tφφ′′+m1+σ11σ41Tθ1′+m2+σ21σ41Tθ2′+m3+σ31σ41Tφ′+m4Tdt+Pσ11Tθ1′+σ21Tθ2′+σ31Tφ′+σ41TdWt.If the transformation functionTθ1,θ2,φ is bounded and nonsingular, both processes P and S have the same stability behavior. Therefore, transformation function Tθ1,θ2,φ is chosen so that the drift term of the Itô differential equation (14) does not depend on the phase processes θ1, θ2, and φ, so one can obtain(15)dS=ΛpSdt+ST-1σ11Tθ1′+σ21Tθ2′+σ31Tφ′+σ41TdWt.Comparing (14) and (15), such transformation function Tθ1,θ2,φ can be written by the following equation:(16)L0+ϵL1Tθ1,θ2,φ=ΛpTθ1,θ2,φ.where L0 and L1 are the following first-order and second-order differential operators:(17)L0=ω1∂∂θ1+ω2∂∂θ2,L1=a1θ1,θ2,φ∂2∂θ12+a2θ1,θ2,φ∂2∂θ1∂θ2+a3θ1,θ2,φ∂2∂θ1∂φ+a4θ1,θ2,φ∂2∂θ22+a5θ1,θ2,φ∂2∂θ2∂φ+a6θ1,θ2,φ∂2∂φ2+b1θ1,θ2,φ∂∂θ1+b2θ1,θ2,φ∂∂θ2+b3θ1,θ2,φ∂∂φ+cθ1,θ2,φ,in which(18)
a
1
θ
1
,
θ
2
,
φ
=
1
2
k
1
2
σ
2
c
o
s
4
θ
1
,
a
2
θ
1
,
θ
2
,
φ
=
σ
2
k
1
k
2
c
o
s
2
θ
1
c
o
s
2
θ
2
,
a
3
θ
1
,
θ
2
,
φ
=
-
k
1
σ
2
4
k
1
sin
⁡
2
θ
1
-
k
2
sin
⁡
2
θ
2
c
o
s
2
θ
1
sin
⁡
2
φ
,
a
4
θ
1
,
θ
2
,
φ
=
1
2
k
2
2
σ
2
c
o
s
4
θ
2
,
a
5
θ
1
,
θ
2
,
φ
=
-
k
2
σ
2
4
k
1
sin
⁡
2
θ
1
-
k
2
sin
⁡
2
θ
2
c
o
s
2
θ
2
sin
⁡
2
φ
,
a
6
θ
1
,
θ
2
,
φ
=
σ
2
32
k
1
sin
⁡
2
θ
1
-
k
2
sin
⁡
2
θ
2
2
s
i
n
2
2
φ
,
b
1
θ
1
,
θ
2
,
φ
=
-
ζ
1
2
M
1
sin
⁡
2
θ
1
-
N
1
1
+
cos
⁡
2
θ
1
+
1
2
h
1
1
+
cos
⁡
2
θ
1
-
2
cos
⁡
θ
1
cos
⁡
θ
2
tan
⁡
φ
-
k
1
2
σ
2
8
×
2
sin
⁡
2
θ
1
+
sin
⁡
4
θ
1
+
k
1
σ
2
p
4
k
1
1
+
cos
⁡
2
φ
sin
⁡
2
θ
1
+
k
2
1
-
cos
⁡
2
φ
sin
⁡
2
θ
2
c
o
s
2
θ
1
,
b
2
θ
1
,
θ
2
,
φ
=
-
ζ
2
2
M
2
sin
⁡
2
θ
2
-
N
2
1
+
cos
⁡
2
θ
2
+
1
2
h
2
1
+
cos
⁡
2
θ
2
-
2
cos
⁡
θ
1
cos
⁡
θ
2
cot
⁡
φ
-
k
2
2
σ
2
8
×
2
sin
⁡
2
θ
2
+
sin
⁡
4
θ
2
+
k
2
σ
2
p
4
k
1
1
+
cos
⁡
2
φ
sin
⁡
2
θ
1
+
k
2
1
-
cos
⁡
2
φ
sin
⁡
2
θ
2
c
o
s
2
θ
2
,
b
3
θ
1
,
θ
2
,
φ
=
h
1
4
2
1
-
cos
⁡
2
φ
sin
⁡
θ
1
cos
⁡
θ
2
-
sin
⁡
2
φ
sin
⁡
2
θ
1
+
ζ
1
4
sin
⁡
2
φ
M
1
1
-
cos
⁡
2
θ
1
-
N
1
sin
⁡
2
θ
1
+
h
2
4
-
2
1
+
cos
⁡
2
φ
cos
⁡
θ
1
sin
⁡
θ
2
+
sin
⁡
2
φ
sin
⁡
2
θ
2
-
ζ
2
4
sin
⁡
2
φ
M
2
1
-
cos
⁡
2
θ
2
-
N
2
sin
⁡
2
θ
2
+
σ
2
64
-
4
k
1
2
1
+
2
cos
⁡
2
θ
1
+
cos
⁡
4
θ
1
×
sin
⁡
2
φ
+
k
1
2
1
-
cos
⁡
4
θ
1
sin
⁡
4
φ
-
4
k
1
k
2
sin
⁡
2
θ
1
sin
⁡
2
θ
2
sin
⁡
4
φ
+
4
k
2
2
1
+
2
cos
⁡
2
θ
2
+
cos
⁡
4
θ
2
sin
⁡
2
φ
+
k
2
2
1
-
cos
⁡
4
θ
2
sin
⁡
4
φ
-
σ
2
p
16
k
1
sin
⁡
2
θ
1
-
k
2
sin
⁡
2
θ
2
×
k
1
1
+
cos
⁡
2
φ
sin
⁡
2
θ
1
+
k
2
1
-
cos
⁡
2
φ
sin
⁡
2
θ
2
sin
⁡
2
φ
,
c
θ
1
,
θ
2
,
φ
=
p
4
h
1
1
+
cos
⁡
2
φ
sin
⁡
2
θ
1
-
2
sin
⁡
2
φ
cos
⁡
θ
2
sin
⁡
θ
1
-
ζ
1
1
+
cos
⁡
2
φ
M
1
1
-
cos
⁡
2
θ
1
-
N
1
sin
⁡
2
θ
1
+
h
2
1
-
cos
⁡
2
φ
sin
⁡
2
θ
2
-
2
sin
⁡
2
φ
cos
⁡
θ
1
sin
⁡
θ
2
-
ζ
2
1
-
cos
⁡
2
φ
M
2
1
-
cos
⁡
2
θ
2
-
N
2
sin
⁡
2
θ
2
+
σ
2
k
1
2
32
×
10
+
3
p
+
16
cos
⁡
2
θ
1
-
3
p
-
2
cos
⁡
4
θ
1
+
8
+
4
p
+
16
cos
⁡
2
θ
1
-
4
p
-
2
cos
⁡
4
θ
1
×
cos
⁡
2
φ
+
p
-
2
1
-
cos
⁡
4
θ
1
cos
⁡
4
φ
+
σ
2
k
1
k
2
8
p
-
2
1
-
cos
⁡
4
φ
sin
⁡
2
θ
1
sin
⁡
2
θ
2
+
σ
2
k
2
2
32
10
+
3
p
+
16
cos
⁡
2
θ
2
-
3
p
-
2
cos
⁡
4
θ
2
-
8
+
4
p
+
16
cos
⁡
2
θ
2
-
4
p
-
2
cos
⁡
4
θ
2
cos
⁡
2
φ
+
p
-
2
1
-
cos
⁡
4
θ
2
cos
⁡
4
φ
.In the following investigation, the perturbation theory will be used to obtain the solutionΛp of (16).


## 3. Moment Lyapunov Exponents

The method of deriving the eigenvalue problem for the moment Lyapunov exponents of a two-dimensional linear Itô stochastic system was first applied by Wedig [34]. The eigenvalue problem for a differential operator of three independent variables θ1, θ2, and φ will be identified from (16), in which Λp is the eigenvalue and Tθ1,θ2,φ is the associated eigenfunction. Meanwhile, the eigenvalue Λp is seen to be the Lyapunov exponent of the pth moment of system (1) from (15). Applying the perturbation theory, both the moment Lyapunov exponent Λp and the eigenfunction Tθ1,θ2,φ are expanded in power series of ϵ; that is, (19)Λp=Λ0p+ϵΛ1p+ϵ2Λ2p+⋯+ϵnΛnp+⋯,Tθ1,θ2,φ=T0θ1,θ2,φ+ϵT1θ1,θ2,φ+ϵ2T2θ1,θ2,φ+⋯+ϵnTnθ1,θ2,φ+⋯.Substituting (19) into (16) and equating terms of the equal powers of ϵ lead to the following equations:(20)ϵ0:L0T0θ1,θ2,φ=Λ0pT0θ1,θ2,φ,ϵ1:L0T1θ1,θ2,φ+L1T0θ1,θ2,φ=Λ0pT1θ1,θ2,φ+Λ1pT0θ1,θ2,φ,⋯ϵn:L0Tnθ1,θ2,φ+L1Tn-1θ1,θ2,φ=Λ0pTnθ1,θ2,φ+Λ1pTn-1θ1,θ2,φ+⋯+ΛnpT0θ1,θ2,φ,where Tiθ1,θ2,φ must be positive and periodic in the range 0≤θ1≤2π, 0≤θ2≤2π.


### 3.1. The Zeroth-Order Perturbation

Substituting (17) into (20), the zeroth-order perturbation equation is(21)ω1∂T0θ1,θ2,φ∂θ1+ω2∂T0θ1,θ2,φ∂θ2=Λ0pT0θ1,θ2,φ.Through Λ0=Λ00+ϵΛ10+ϵ2Λ20+⋯+ϵnΛn0=0, thus, Λn0=0, n=1,2… will be obtained. Then (21) does not contain p, and the eigenvalue Λ0p is independent of p. Therefore, Λ00=0 will lead to Λ0p=0 for arbitrary p. At this point, (21) will be reduced to(22)ω1∂T0θ1,θ2,φ∂θ1+ω2∂T0θ1,θ2,φ∂θ2=0.Using the method of separation of variables, so we can express T0θ1,θ2,φas(23)T0θ1,θ2,φ=Θ1θ1Θ2θ2ψφ,where Θ1θ1, Θ2θ2, and ψφ are functions about θ1, θ2, and φ, respectively.Substituting (23) into (22), (22) will become (24)ω1Θ1∂Θ1∂θ1+ω2Θ2∂Θ2∂θ2=0.Solving the above differential equation (24), one can easily obtain (25)Θ1θ1=C1ec1θ1,Θ2θ2=C2ec2θ2.In addition, since function T0θ1,θ2,φ is periodic about θ1 and θ2, the following boundary conditions can be obtained:(26)T0θ1+2π,θ2,φ=T0θ1,θ2+2π,φ=T0θ1,θ2,φ.Substituting (23) and (25) into (26), we can determine the constants c1=c2=0. Therefore, T0θ1,θ2,φ is a function about only φ to be determined; that is, T0θ1,θ2,φ=T0φ.The adjoint equation of (22) is (27)-ω1∂T0∗θ1,θ2,φ∂θ1-ω2∂T0∗θ1,θ2,φ∂θ2=0.Applying the method of separation of variables to solve (27), one obtains(28)T0∗=Fφ4π2,φ∈0,2π,where Fφ is an arbitrary function.


### 3.2. The First-Order Perturbation

From (20), the first-order perturbation equation is(29)L0T1θ1,θ2,φ+L1T0θ1,θ2,φ=Λ0pT1θ1,θ2,φ+Λ1pT0θ1,θ2,φ,Based on the analysis results in Section 3.1, one obtains(30)L0T1θ1,θ2,φ=Λ1pT0φ-L1T0θ1,θ2,φ.Then the solvability condition of (30) is(31)Λ1pT0φ-L1T0φ,T0∗=14π2∫0π/2Fφ∫02π∫02πΛ1pT0φ-L1T0φdθ1dθ2dφ=0,L1T0φ=cT0φ+b3∂T0φ∂φ+a6∂2T0φ∂φ2.Applying the sense of the arbitrary of function F(φ), (31) will become(32)∫0π/2Λ1p-pQφT0φ-μφ+pμ~φ∂T0φ∂φ-12R2φ∂2T0φ∂φ2dφ=0,where (33)Qφ=132126+pk12+k22σ2-8M1ζ1+M2ζ2+p+2k12-k22σ2-8M1ζ1-M2ζ2×cos⁡2φ+12-2+pk12+k22σ2cos22φ,μφ=-1322k12-k22σ2-8M1ζ1-M2ζ2sin⁡2φ-k12+k22σ2sin⁡2φcos⁡2φ,μ~φ=-132k12-k22σ2sin⁡2φ+k12+k22σ2sin⁡2φcos⁡2φ,R2φ=164k12+k22σ2sin22φ;then (32) is simplified as the following equation:(34)LpT0φ=Λ1pT0φ,Lp=12R2φ∂∂φ2+μφ+pμ~φ∂∂φ+pQφ.Then the boundary conditions of (34) are determined by considering the adjoint equation for the case of p=0,(35)L∗m~φ=0,L∗=12∂∂φ2R2φ-∂∂φμφ,where L∗ is the Fokker-Planck operator and φ=0 and φ=π/2 are the entrance boundaries [27]. The eigenfunction T0φ satisfies zero Neumann boundary condition, and Λ1p is the largest eigenvalue of (34) with zeros Neumann boundary. Therefore, the solution of (34) can be calculated from an orthogonal expansion [34]. Under zeros Neumann boundary conditions, T0φ can be expressed by a Fourier cosine series [27]; that is,(36)T0φ=∑n=0∞zncos⁡2nφ.Substituting (36) into (34), multiplying by cos⁡2nφ in both sides and integrating for φ, we can obtain(37)∑m=0∞am0zm=2Λ1pz0,∑m=0∞amnzm=Λ1pzn.n=1,2,…,where (38)amn=4π∫0π/2Lpcos⁡2mφcos⁡2nφdφ,m,n=0,1,2…,In order to guarantee the existence of the nontrivial solution for each zn, the coefficient matrix A=amn must equal zero. Therefore, the problem of calculating Λ1p is translated into calculating the leading eigenvalue of A. The sequence of approximations will be constructed by the eigenvalues of a sequence of the following submatrices: (39)12a00,12a00a0112a10a11,⋯,12a00a01a02⋯12a10a11a12⋯12a20a21a22⋯⋮⋮⋮⋱.Obviously, the set of approximate eigenvalues obtained by this method converges to the associated true eigenvalues asn⟶∞. In general, we can obtain the approximate eigenvalues by truncating n, which can be written as Λ1p=a00/2 and (40)Λp=ϵΛ1p+Oϵ2=1128ϵp10+3pk12+k22σ2-32M1ζ1+M2ζ2+Oϵ2.The comparison results of approximate analytical moment Lyapunov exponents obtained by truncatingn and the direct Monte Carlo simulation results are shown in Figure 1. It is easily found that the approximate results agree well with the simulation results. Therefore, the first-order perturbation will be reduced to(41)ω1∂T1θ1,θ2,φ∂θ1+ω2∂T1θ1,θ2,φ∂θ2+cθ1,θ2,φ=Λ1p.Without loss of generality, there is a relationship between two frequencies of the form m1ω1=m2ω2, in which m1 and m2 are integers. The second frequency can be rewritten as ω2=kω1. Then the function cθ1,θ2,φ in (41) can be rewritten as(42)cθ1,θ2,φ=Λ1p+f0θ1,θ2+f1θ1,θ2cos⁡2φ+f2θ1,θ2cos⁡4φ+f3θ1,θ2sin⁡2φ,where function frθ1,θ2 is the periodic function on θ1,θ2 and given as(43)f0θ1,θ2=p4h1+ζ1N1sin⁡2θ1+p82ζ1M1+k12σ2cos⁡2θ1+p4h2+ζ2N2sin⁡2θ2+p82ζ2M2+k22σ2cos⁡2θ2+3k12σ2128p2-pcos⁡4θ1+3k22σ2128p2-pcos⁡4θ2-k1k2σ232p2-psin⁡2θ1sin⁡2θ2,f1θ1,θ2=p4h1+ζ1N1sin⁡2θ1+p82ζ1M1+k12σ2cos⁡2θ1-p4h2+ζ2N2sin⁡2θ2-p82ζ2M2+k22σ2cos⁡2θ2+k12σ232p2-pcos⁡4θ1-k22σ232p2-pcos⁡4θ2+p328-ζ1M1+ζ2M2+k12-k22σ22+p,f2θ1,θ2=-p2-pσ2128k12+k22+p2-pσ232sin⁡2θ1sin⁡2θ2+p2-pσ2128k12cos⁡4θ1+p2-pσ2128k22cos⁡4θ2,f3θ1,θ2=-h1p2sin⁡2θ1cos⁡2θ2-h2p2sin⁡2θ2cos⁡2θ1.Note that the general solution T1θ1,θ2,φ of (41) can not be obtained except in some special cases. Considering the nature of the coefficients of (41), the series expansion of function T1θ1,θ2,φ can be presented in the following form:(44)T1θ1,θ2,φ=T10θ1,θ2+T11θ1,θ2cos⁡2φ+T12θ1,θ2cos⁡4φ+T13θ1,θ2sin⁡2φ.Figure 1
Variation of moment Lyapunov exponent withp. Solid-line in the analytical results of truncation of n, circle in the numerical results by Monte Carlo.Substituting (44) into (41) and equating terms of the equal trigonometry function to give a set of partial differential equations, one obtains (45)ω1∂T1rθ1,θ2∂θ1+ω2∂T1rθ1,θ2∂θ2+frθ1,θ2=0,r=0,1,2,3,where the function T1rθ1,θ2 can be written as(46)T10θ1,θ2=-p2ζ1M1+k12σ216ω1sin⁡2θ1+ph1+ζ1N18ω1cos⁡2θ1-p2ζ2M2+k22σ216kω1sin⁡2θ2+ph2+ζ2N28kω1cos⁡2θ2-3k12p2-pσ2512ω1sin⁡4θ1-3k22p2-pσ2512kω1sin⁡4θ2+kk1k2p2-pσ2641-k2ω1sin⁡2θ1cos⁡2θ2-k1k2p2-pσ2641-k2ω1sin⁡2θ2cos⁡2θ1+C0θ2-kθ1,T11θ1,θ2=-p8-ζ1M1+ζ2M2+k12-k222+pσ232ω1θ1-p2ζ1M1+k12σ216ω1sin⁡2θ1+ph1+ζ1N18ω1cos⁡2θ1+p2ζ2M2+k22σ216kω1sin⁡2θ2-ph2+ζ2N28kω1cos⁡2θ2-k12p2-pσ2128ω1sin⁡4θ1+k22p2-pσ2128kω1sin⁡4θ2+C1θ2-kθ1,T12θ1,θ2=p2-pk12+k22128ω1θ1-k12p2-pσ2512ω1sin⁡4θ1-k22p2-pσ2512kω1sin⁡4θ2-kk1k2p2-pσ2641-k2ω1sin⁡2θ1cos⁡2θ2+k1k2p2-pσ2641-k2ω1sin⁡2θ2cos⁡2θ1+C2θ2-kθ1,T13θ1,θ2=ph2-h1k21-k2ω1sin⁡θ1sin⁡θ2-ph1-h2k21-k2ω1cos⁡θ1cos⁡θ2+C3θ2-kθ1,in which Crθ2-kθ1 are arbitrary functions of two variables, and we make assumptions as follows:(47)Crθ2-kθ1=A1r+B1rsin2θ2-2kθ1+C1rsin4θ2-4kθ1,r=0,1,2,3.Here, the unknown constants A1r, B1r, C1r  r=0,1,2,3 will be determined by using the following conditions:(48)T1r0,0=T1r0,2π=T1r2π,0=T1r2π,2π=0,∂T1r0,0∂θ1=∂T1r2π,0∂θ1,∂T1r0,0∂θ2=∂T1r0,2π∂θ2,r=0,1,2,3.That is,(49)A10=-ph2+ζ2N2+h1k+kζ1N18kω1,A11=ph2+ζ2N2-h1k-kζ1N18kω1,A12=0,A13=ph1-h2k21-k2ω1,B10=0,B11=-pπ8-ζ1M1+ζ2M2+k12-k222+pσ2cos22kπ4sin⁡2kπ5cos⁡2kπ+cos⁡6kπω1,B12=pπk12+k222-pσ2cos⁡2kπ163sin⁡2kπ+sin⁡6kπω1,B13=0,C10=0,C13=0C11=pπ8-ζ1M1+ζ2M2+k12-k222+pσ232sin⁡2kπ5cos⁡2kπ+cos⁡6kπω1,C12=-pπk12+k222-pσ2644sin⁡4kπ+sin⁡8kπω1,


### 3.3. The Second-Order Perturbation

To further investigate the effects of parameterH2 on stochastic stability of system, the second-order perturbation equation of (20) will be used; that is,(50)L0T2θ1,θ2,φ+L1T1θ1,θ2,φ=Λ1pT1θ1,θ2,φ+Λ2pT0θ1,θ2,φ.Based on the solvability condition, one obtains(51)Λ2p=12π3∫0π/2∫02π∫02πa1θ1,θ2,φ∂2T1θ1,θ2,φ∂θ12+a2θ1,θ2,φ∂2T1θ1,θ2,φ∂θ1∂θ2+a3θ1,θ2,φ∂2T1θ1,θ2,φ∂θ1∂φ+a4θ1,θ2,φ∂2T1θ1,θ2,φ∂θ22+a5θ1,θ2,φ∂2T1θ1,θ2,φ∂θ2∂φ+a6θ1,θ2,φ∂2T1θ1,θ2,φ∂φ2+b1θ1,θ2,φ∂T1θ1,θ2,φ∂θ1+b2θ1,θ2,φ∂T1θ1,θ2,φ∂θ2+b3θ1,θ2,φ∂T1θ1,θ2,φ∂φ+cθ1,θ2,φ-Λ1pT1θ1,θ2,φdθ1dθ2dφ.The solution has the following form:(52)Λ2p=Λ21psin⁡4kπ2+cos⁡4kπ+Λ22pcos⁡4kπ2+cos⁡4kπ+Λ23p12+cos⁡4kπ,where the values Λ21p, Λ22p, Λ23p are given in Appendix A.Substituting (40) and (52) into (19), the second-order approximate solution of the moment Lyapunov exponent is obtained as(53)Λp=ϵΛ1p+ϵ2Λ2p+Oϵ3=ϵk12+k22σ2128p3p+10-p4ζ1M1+ζ2M2+ϵ2Λ2p+Oϵ3.The corresponding Lyapunov exponents can be expressed as(54)λ=dΛpdpp=0=ϵλ1+ϵ2λ2+Oϵ3=ϵ564k12+k22σ2-14ζ1M1+ζ2M2+ϵ2λ21sin⁡4kπ2+cos⁡4kπ+λ22cos⁡4kπ2+cos⁡4kπ+λ2312+cos⁡4kπ+Oϵ3,where the values λ21, λ22, and λ23 are given in Appendix B.


## 3.1. The Zeroth-Order Perturbation

Substituting (17) into (20), the zeroth-order perturbation equation is(21)ω1∂T0θ1,θ2,φ∂θ1+ω2∂T0θ1,θ2,φ∂θ2=Λ0pT0θ1,θ2,φ.Through Λ0=Λ00+ϵΛ10+ϵ2Λ20+⋯+ϵnΛn0=0, thus, Λn0=0, n=1,2… will be obtained. Then (21) does not contain p, and the eigenvalue Λ0p is independent of p. Therefore, Λ00=0 will lead to Λ0p=0 for arbitrary p. At this point, (21) will be reduced to(22)ω1∂T0θ1,θ2,φ∂θ1+ω2∂T0θ1,θ2,φ∂θ2=0.Using the method of separation of variables, so we can express T0θ1,θ2,φas(23)T0θ1,θ2,φ=Θ1θ1Θ2θ2ψφ,where Θ1θ1, Θ2θ2, and ψφ are functions about θ1, θ2, and φ, respectively.Substituting (23) into (22), (22) will become (24)ω1Θ1∂Θ1∂θ1+ω2Θ2∂Θ2∂θ2=0.Solving the above differential equation (24), one can easily obtain (25)Θ1θ1=C1ec1θ1,Θ2θ2=C2ec2θ2.In addition, since function T0θ1,θ2,φ is periodic about θ1 and θ2, the following boundary conditions can be obtained:(26)T0θ1+2π,θ2,φ=T0θ1,θ2+2π,φ=T0θ1,θ2,φ.Substituting (23) and (25) into (26), we can determine the constants c1=c2=0. Therefore, T0θ1,θ2,φ is a function about only φ to be determined; that is, T0θ1,θ2,φ=T0φ.The adjoint equation of (22) is (27)-ω1∂T0∗θ1,θ2,φ∂θ1-ω2∂T0∗θ1,θ2,φ∂θ2=0.Applying the method of separation of variables to solve (27), one obtains(28)T0∗=Fφ4π2,φ∈0,2π,where Fφ is an arbitrary function.


## 3.2. The First-Order Perturbation

From (20), the first-order perturbation equation is(29)L0T1θ1,θ2,φ+L1T0θ1,θ2,φ=Λ0pT1θ1,θ2,φ+Λ1pT0θ1,θ2,φ,Based on the analysis results in Section 3.1, one obtains(30)L0T1θ1,θ2,φ=Λ1pT0φ-L1T0θ1,θ2,φ.Then the solvability condition of (30) is(31)Λ1pT0φ-L1T0φ,T0∗=14π2∫0π/2Fφ∫02π∫02πΛ1pT0φ-L1T0φdθ1dθ2dφ=0,L1T0φ=cT0φ+b3∂T0φ∂φ+a6∂2T0φ∂φ2.Applying the sense of the arbitrary of function F(φ), (31) will become(32)∫0π/2Λ1p-pQφT0φ-μφ+pμ~φ∂T0φ∂φ-12R2φ∂2T0φ∂φ2dφ=0,where (33)Qφ=132126+pk12+k22σ2-8M1ζ1+M2ζ2+p+2k12-k22σ2-8M1ζ1-M2ζ2×cos⁡2φ+12-2+pk12+k22σ2cos22φ,μφ=-1322k12-k22σ2-8M1ζ1-M2ζ2sin⁡2φ-k12+k22σ2sin⁡2φcos⁡2φ,μ~φ=-132k12-k22σ2sin⁡2φ+k12+k22σ2sin⁡2φcos⁡2φ,R2φ=164k12+k22σ2sin22φ;then (32) is simplified as the following equation:(34)LpT0φ=Λ1pT0φ,Lp=12R2φ∂∂φ2+μφ+pμ~φ∂∂φ+pQφ.Then the boundary conditions of (34) are determined by considering the adjoint equation for the case of p=0,(35)L∗m~φ=0,L∗=12∂∂φ2R2φ-∂∂φμφ,where L∗ is the Fokker-Planck operator and φ=0 and φ=π/2 are the entrance boundaries [27]. The eigenfunction T0φ satisfies zero Neumann boundary condition, and Λ1p is the largest eigenvalue of (34) with zeros Neumann boundary. Therefore, the solution of (34) can be calculated from an orthogonal expansion [34]. Under zeros Neumann boundary conditions, T0φ can be expressed by a Fourier cosine series [27]; that is,(36)T0φ=∑n=0∞zncos⁡2nφ.Substituting (36) into (34), multiplying by cos⁡2nφ in both sides and integrating for φ, we can obtain(37)∑m=0∞am0zm=2Λ1pz0,∑m=0∞amnzm=Λ1pzn.n=1,2,…,where (38)amn=4π∫0π/2Lpcos⁡2mφcos⁡2nφdφ,m,n=0,1,2…,In order to guarantee the existence of the nontrivial solution for each zn, the coefficient matrix A=amn must equal zero. Therefore, the problem of calculating Λ1p is translated into calculating the leading eigenvalue of A. The sequence of approximations will be constructed by the eigenvalues of a sequence of the following submatrices: (39)12a00,12a00a0112a10a11,⋯,12a00a01a02⋯12a10a11a12⋯12a20a21a22⋯⋮⋮⋮⋱.Obviously, the set of approximate eigenvalues obtained by this method converges to the associated true eigenvalues asn⟶∞. In general, we can obtain the approximate eigenvalues by truncating n, which can be written as Λ1p=a00/2 and (40)Λp=ϵΛ1p+Oϵ2=1128ϵp10+3pk12+k22σ2-32M1ζ1+M2ζ2+Oϵ2.The comparison results of approximate analytical moment Lyapunov exponents obtained by truncatingn and the direct Monte Carlo simulation results are shown in Figure 1. It is easily found that the approximate results agree well with the simulation results. Therefore, the first-order perturbation will be reduced to(41)ω1∂T1θ1,θ2,φ∂θ1+ω2∂T1θ1,θ2,φ∂θ2+cθ1,θ2,φ=Λ1p.Without loss of generality, there is a relationship between two frequencies of the form m1ω1=m2ω2, in which m1 and m2 are integers. The second frequency can be rewritten as ω2=kω1. Then the function cθ1,θ2,φ in (41) can be rewritten as(42)cθ1,θ2,φ=Λ1p+f0θ1,θ2+f1θ1,θ2cos⁡2φ+f2θ1,θ2cos⁡4φ+f3θ1,θ2sin⁡2φ,where function frθ1,θ2 is the periodic function on θ1,θ2 and given as(43)f0θ1,θ2=p4h1+ζ1N1sin⁡2θ1+p82ζ1M1+k12σ2cos⁡2θ1+p4h2+ζ2N2sin⁡2θ2+p82ζ2M2+k22σ2cos⁡2θ2+3k12σ2128p2-pcos⁡4θ1+3k22σ2128p2-pcos⁡4θ2-k1k2σ232p2-psin⁡2θ1sin⁡2θ2,f1θ1,θ2=p4h1+ζ1N1sin⁡2θ1+p82ζ1M1+k12σ2cos⁡2θ1-p4h2+ζ2N2sin⁡2θ2-p82ζ2M2+k22σ2cos⁡2θ2+k12σ232p2-pcos⁡4θ1-k22σ232p2-pcos⁡4θ2+p328-ζ1M1+ζ2M2+k12-k22σ22+p,f2θ1,θ2=-p2-pσ2128k12+k22+p2-pσ232sin⁡2θ1sin⁡2θ2+p2-pσ2128k12cos⁡4θ1+p2-pσ2128k22cos⁡4θ2,f3θ1,θ2=-h1p2sin⁡2θ1cos⁡2θ2-h2p2sin⁡2θ2cos⁡2θ1.Note that the general solution T1θ1,θ2,φ of (41) can not be obtained except in some special cases. Considering the nature of the coefficients of (41), the series expansion of function T1θ1,θ2,φ can be presented in the following form:(44)T1θ1,θ2,φ=T10θ1,θ2+T11θ1,θ2cos⁡2φ+T12θ1,θ2cos⁡4φ+T13θ1,θ2sin⁡2φ.Figure 1
Variation of moment Lyapunov exponent withp. Solid-line in the analytical results of truncation of n, circle in the numerical results by Monte Carlo.Substituting (44) into (41) and equating terms of the equal trigonometry function to give a set of partial differential equations, one obtains (45)ω1∂T1rθ1,θ2∂θ1+ω2∂T1rθ1,θ2∂θ2+frθ1,θ2=0,r=0,1,2,3,where the function T1rθ1,θ2 can be written as(46)T10θ1,θ2=-p2ζ1M1+k12σ216ω1sin⁡2θ1+ph1+ζ1N18ω1cos⁡2θ1-p2ζ2M2+k22σ216kω1sin⁡2θ2+ph2+ζ2N28kω1cos⁡2θ2-3k12p2-pσ2512ω1sin⁡4θ1-3k22p2-pσ2512kω1sin⁡4θ2+kk1k2p2-pσ2641-k2ω1sin⁡2θ1cos⁡2θ2-k1k2p2-pσ2641-k2ω1sin⁡2θ2cos⁡2θ1+C0θ2-kθ1,T11θ1,θ2=-p8-ζ1M1+ζ2M2+k12-k222+pσ232ω1θ1-p2ζ1M1+k12σ216ω1sin⁡2θ1+ph1+ζ1N18ω1cos⁡2θ1+p2ζ2M2+k22σ216kω1sin⁡2θ2-ph2+ζ2N28kω1cos⁡2θ2-k12p2-pσ2128ω1sin⁡4θ1+k22p2-pσ2128kω1sin⁡4θ2+C1θ2-kθ1,T12θ1,θ2=p2-pk12+k22128ω1θ1-k12p2-pσ2512ω1sin⁡4θ1-k22p2-pσ2512kω1sin⁡4θ2-kk1k2p2-pσ2641-k2ω1sin⁡2θ1cos⁡2θ2+k1k2p2-pσ2641-k2ω1sin⁡2θ2cos⁡2θ1+C2θ2-kθ1,T13θ1,θ2=ph2-h1k21-k2ω1sin⁡θ1sin⁡θ2-ph1-h2k21-k2ω1cos⁡θ1cos⁡θ2+C3θ2-kθ1,in which Crθ2-kθ1 are arbitrary functions of two variables, and we make assumptions as follows:(47)Crθ2-kθ1=A1r+B1rsin2θ2-2kθ1+C1rsin4θ2-4kθ1,r=0,1,2,3.Here, the unknown constants A1r, B1r, C1r  r=0,1,2,3 will be determined by using the following conditions:(48)T1r0,0=T1r0,2π=T1r2π,0=T1r2π,2π=0,∂T1r0,0∂θ1=∂T1r2π,0∂θ1,∂T1r0,0∂θ2=∂T1r0,2π∂θ2,r=0,1,2,3.That is,(49)A10=-ph2+ζ2N2+h1k+kζ1N18kω1,A11=ph2+ζ2N2-h1k-kζ1N18kω1,A12=0,A13=ph1-h2k21-k2ω1,B10=0,B11=-pπ8-ζ1M1+ζ2M2+k12-k222+pσ2cos22kπ4sin⁡2kπ5cos⁡2kπ+cos⁡6kπω1,B12=pπk12+k222-pσ2cos⁡2kπ163sin⁡2kπ+sin⁡6kπω1,B13=0,C10=0,C13=0C11=pπ8-ζ1M1+ζ2M2+k12-k222+pσ232sin⁡2kπ5cos⁡2kπ+cos⁡6kπω1,C12=-pπk12+k222-pσ2644sin⁡4kπ+sin⁡8kπω1,


## 3.3. The Second-Order Perturbation

To further investigate the effects of parameterH2 on stochastic stability of system, the second-order perturbation equation of (20) will be used; that is,(50)L0T2θ1,θ2,φ+L1T1θ1,θ2,φ=Λ1pT1θ1,θ2,φ+Λ2pT0θ1,θ2,φ.Based on the solvability condition, one obtains(51)Λ2p=12π3∫0π/2∫02π∫02πa1θ1,θ2,φ∂2T1θ1,θ2,φ∂θ12+a2θ1,θ2,φ∂2T1θ1,θ2,φ∂θ1∂θ2+a3θ1,θ2,φ∂2T1θ1,θ2,φ∂θ1∂φ+a4θ1,θ2,φ∂2T1θ1,θ2,φ∂θ22+a5θ1,θ2,φ∂2T1θ1,θ2,φ∂θ2∂φ+a6θ1,θ2,φ∂2T1θ1,θ2,φ∂φ2+b1θ1,θ2,φ∂T1θ1,θ2,φ∂θ1+b2θ1,θ2,φ∂T1θ1,θ2,φ∂θ2+b3θ1,θ2,φ∂T1θ1,θ2,φ∂φ+cθ1,θ2,φ-Λ1pT1θ1,θ2,φdθ1dθ2dφ.The solution has the following form:(52)Λ2p=Λ21psin⁡4kπ2+cos⁡4kπ+Λ22pcos⁡4kπ2+cos⁡4kπ+Λ23p12+cos⁡4kπ,where the values Λ21p, Λ22p, Λ23p are given in Appendix A.Substituting (40) and (52) into (19), the second-order approximate solution of the moment Lyapunov exponent is obtained as(53)Λp=ϵΛ1p+ϵ2Λ2p+Oϵ3=ϵk12+k22σ2128p3p+10-p4ζ1M1+ζ2M2+ϵ2Λ2p+Oϵ3.The corresponding Lyapunov exponents can be expressed as(54)λ=dΛpdpp=0=ϵλ1+ϵ2λ2+Oϵ3=ϵ564k12+k22σ2-14ζ1M1+ζ2M2+ϵ2λ21sin⁡4kπ2+cos⁡4kπ+λ22cos⁡4kπ2+cos⁡4kπ+λ2312+cos⁡4kπ+Oϵ3,where the values λ21, λ22, and λ23 are given in Appendix B.


## 4. Stochastic Stability Analysis

Through the above-mentioned analysis, the analytic expressions of moment Lyapunov exponentsΛp and Lyapunov exponents λ have been derived. Next, we use two indices to go on stochastic stability analysis. The pth moment Lyapunov exponent is defined as(55)Λp=limt→∞1tlog⁡Εxp,where Ε· denotes the expected value and xp is a solution process of a random system. If Λp<0, Εxp→0 as t→∞, and it means that the pth moment is stable. Lyapunov exponent is the slope of Λp at the origin; that is,(56)λ=dΛpdpp=0=limt→∞1tlogx,If λ<0, it indicates that the system is almost-surely stable.Thepth moment and almost-sure stability boundary can be determined from Λp=0 and λ=0, respectively. By comparing the second-order perturbation with the first-order perturbation, we find that the moment Lyapunov exponents Λp and the Lyapunov exponents λ under the second-order perturbation contain all the parameter information of the coupled stochastic system. Therefore, we only consider the second-order perturbation in the subsequent research.The boundaries of thepth moment and the almost-sure stability with respect to nonviscous damping parameters β1 and β2 are shown in Figure 2. We find that the moment stability boundaries are more conservative than the almost-sure stability boundaries. This means that the almost-sure stability cannot assure the moment stability. And the boundaries of stability for pth moment would increase with the increase of the order p.Figure 2
Stability regions for almost-sure stability andpth moment stability with respect to nonviscous damping parameter β1, β2.The effects of relaxation parameterα on moment Lyapunov exponents are shown in Figure 3. It describes that the stability index and the root p of Λp=0 will decrease and the system becomes more unstable. Correspondingly, the analytical results and the direct Monte Carlo results are shown in Figure 3(b) for the given α. And Lyapunov exponent is negative; it means the system is almost-surely stable.Figure 3
Variation of the moment Lyapunov exponent for differentα. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)The moment Lyapunov exponentsΛp for different nonviscous damping parameter β2 are discussed in Figure 4. It is clear that the bigger values of β2 can enhance the system's stability for p>0. Correspondingly, the analytical results and the direct Monte Carlo results are shown in Figure 4(b) for the given different β2. It illustrates that the system also is almost-surely stable. However, when p becomes sufficiently large, Λp>0, the system becomes more unstable.Figure 4
Variation of the moment Lyapunov exponent for differentβ2. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)Next, the effects ofϵ on the moment Lyapunov exponents are studied in Figure 5. One easily finds the system is moment stability when 0<p<3, and the stability of the system increases with the increase of ϵ from 0.01 to 0.1. Afterward, the moment Lyapunov exponents are changed from negative to positive; the system is varied from moment stability to moment instability with the increase of p. Numerical results are verified by the direct Monte Carlo in Figure 5(b).Figure 5
Variation of the moment Lyapunov exponent for differentε. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)Then, in order to discuss the influences of the noise intensityσ on stochastic stability of coupled stochastic system, the moment Lyapunov exponents Λp are shown in Figure 6 with change of values of σ. They describe that the coupled stochastic system could be almost-surely stable. Meanwhile, the moment stability of the coupled system would reduce with the increase of the noise intensity σ.Figure 6
Variation of the moment Lyapunov exponent for differentσ. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)The variations of the moment Lyapunov exponents for differentH2 are also analyzed. From Figure 7, we find that the coupled stochastic system may be almost-surely stable due to the value of slope Λp less than 0 at the origin. However, the moment stability of system is not assured for p<0 or p is sufficiently large due to the moment Lyapunov exponents greater than 0. Besides, one can easily see that the moment stability of the system enhances with the increase of H2.Figure 7
Variation of the moment Lyapunov exponent for differentH2. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)Finally, the effects ofK2 on the moment Lyapunov exponents are investigated and shown in Figure 8. From those figures, the system may be almost-surely stable in terms of the Lyapunov exponents. It is also indicated that the system has the characteristic of the moment stability. And the moment stability of the coupled system will increase with the increase of K2 under given parameters. The approximate analytical results agree well with numerical simulation results by the direct Monte Carlo.Figure 8
Variation of the moment Lyapunov exponent for differentK2. Solid-line in the analytical results, circle in the numerical results by Monte Carlo.
(a) (b)


## 5. Conclusions

In this paper, the moment stability and almost-sure stability of coupled viscoelastic system with nonviscous damping subject to Gaussian white noise excitation are investigated. The nonviscously damped structure is assumed to follow the exponential integral relation, which is possible to model the rate of energy dissipation. Motion equations are transformed into the coupled Itô stochastic differential equations by means of the method of stochastic averaging. Then the linear transformation is introduced to derive the eigenvalue problem. In order to obtain the analytical expressions of eigenvalue, the second-order perturbation method is used. Finally, the effects of Gaussian noise, the relaxation parameter, the nonviscous damping coefficient, and two physical quantities, that isH2, K2, on the moment stability are discussed in detail. Then the corresponding results are verified through the direct Monte Carlo simulation technique. Based on the results above, the appropriate parameters will help to strengthen the stability of the coupled viscoelastic system with nonviscous damping.



---

*Source: 6532305-2018-09-16.xml*