Stochastic stability and moment Lyapunov exponent for co-dimension two bifurcation system with a bounded noise

Shenghong Li1 , Jiancheng Wu2

1, 2School of Mathematics and Physics, Jiangsu university of Science and Technology, Zhenjiang 212003, Jiangsu Province, China

2State Key Lab of Mechanics and Control for Mechanical Structures, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, 29 Yudao Street, Nanjing 210014, China

1Corresponding author

Journal of Vibroengineering, Vol. 17, Issue 8, 2015, p. 4589-4602.
Received 22 April 2015; received in revised form 30 July 2015; accepted 11 August 2015; published 30 December 2015

Copyright © 2015 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

In this paper, the pth moment Lyapunov exponent of a co-dimension two bifurcation system that is parametrically excited by a real noise is investigated. By a linear stochastic transformation, the eigenvalue problem of moment Lyapunov exponent is obtained. Then through perturbation method, we deduce the joint probability density function of the phase processes and its eigenvalue problem, which is solved by a Fourier cosine series expansion. Thus, an infinite matrix yields and whose leading eigenvalue is the second order of the asymptotic expansion of the moment Lyapunov exponent. Because of the complexity of elements in matrix A, the eigenvalues of the low order sub-matrices of A are obtained by the truncation of n and the convergence of the eigenvalue sequence is numerically illustrated. Finally, the effects of the system and noise parameters on the moment Lyapunov exponent are discussed.

Keywords: moment stability, moment Lyapunov exponent, perturbation method, diffusion process.

1. Introduction

The moment Lyapunov exponent that describes the pth moment stability of stochastic dynamic system, is defined as:

(1)
Λ p , x 0 = l i m t 1 t l o g E x t , x 0 p ,

where x(t,x0) is a solution to a random dynamical system and E[] denotes expected value. If Λ(p,x0)<0, by definition, E[x(t,x0)p]0 as t, which represents that the pth moment of the stochastic system is stable. The relation between moment stability and almost-sure stability for an un-damped linear oscillator under real noise excitation was found firstly by Molcanov [1]. Later, the results were extended to an arbitrary d-dimensional system by Arnold in Ref. [2], meanwhile, the rigorous definition of the moment Lyapunov exponent was firstly defined and the following results were obtained. Under some conditions, the limit in Eq. (3) exists and is independent of the initial value x0. So the moment Lyapunov exponent can be written as Λ(p), which is a convex analytic function of pR. In addition, the maximal Lyapunov exponent is a derivative of moment Lyapunov exponent at p= 0, i.e.:

(2)
λ = Λ ( p ) p p = 0 = l i m t 1 t l o g x ( t , x 0 ) .

Then, Aronld investigated the linear stochastic systems excited by the real noise and white noise, and given a series of results of the moment Lyapunov exponent respectively in [3, 4]. Thus, the stochastic problem of linear system was resolved completely.

However, it is extremely difficult to investigate the stochastic stability of a nonlinear dynamical system except for several very special cases, especially for moment stability. So far, almost all the investigations on the moment Lyapunov exponents are confined to the approximate analytical methods. For a two dimensional stochastic system, L. Arnold [5] firstly applied a perturbation method to perform the asymptotic expansions of the pth moment Lyapunov exponents on a small noise intensity and a small value of p. Using the same method, Namachchivaya [6] obtained the small pth moment Lyapunov exponent for a system with two coupled oscillators that is excited by a real noise. For a linear conservative system with a white noise, Khasminskii and Moshchuk [7] proved that the finite pth moment Lyapunov exponent and the stability index can respectively be expressed only as the asymptotic expansions of small noise intensity. On the basis of the theorem given in [7], for the same system and random excitation as [6], Namachchivaya and Roessel [8] obtained the asymptotic expansion of the finite pth moment Lyapunov exponent. For the two dimensional system that was driven by the real noise and the bounded noise excitations respectively, Xie [9-10] implemented the weak noise expansions of the finite pth moment Lyapunov exponent, the maximal Lyapunov exponent and the stability index through the same procedure. In recent years, this topic is still a hotspot for the researchers in the field of random dynamical system. The stability properties of a Van der Pol-Duffing oscillator with a real noise was investigated by Liu [11]. Due to the complexity of approximate analytical methods, Higham [12] gave the numerical simulation of moment Lyapunov exponent in stochastic differential equations. Then, the moment Lyapunov exponent and the stochastic stability of a double-beam system under the compressive axial loading and moving narrow bands was studied by Kozic [13]. Shenghong Li [14, 15] presented the result of the moment Lyapunov exponent for a three dimensional stochastic system. For a binary airfoil system driven by non-Gaussian colored noise, DL Hu [16] obtained its moment Lyapunov exponent.

In this paper, the pth moment Lyapunov exponent of a co-dimension two bifurcation system that is parametrically excited by a real noise is investigated. For the same system, the maximal Lyapunov exponent was investigated by Li [17, 18]. In Section 2, the studied dynamical system is introduced and the eigenvalue problem of moment Lyapunov exponent is obtained. In Section 3, through a perturbation method, the probability density functions of differential operator about ε0, ε1 and ε2are solved respectively. In Section 4, via an expansion of orthogonal Fourier series, the eigenvalue problem of the differential operator leads to the eigenvalue problem for an infinite matrix, whose leading eigenvalue is the moment Lyapunov exponent. The numerical results are presented in Section 5 and a conclusion is given in Section 6.

2. Formulation

Consider a typical deterministic co-dimension two bifurcation system that is on a three-dimension central manifold and possesses one zero-eigenvalue and a pair of pure imaginary eigenvalues [19], i.e.:

(3)
r ˙ = μ 1 r + a 1 r z + a 2 r 3 + a 3 r 2 z + O r , z 4 ,
z ˙ = μ 2 z + b 1 r 2 + b 2 z 2 + b 3 r 2 z + b 4 z 3 + O r , z 4 ,
Θ ˙ = ω + O r , z 2 ,

where μ1 and μ2 are unfolding parameters, a1, a2, a3, b1, b2, b3 and ω is real constants. This normalized form arises in the classic fluid dynamic stability study of coquette flow.

According to Oseledec multiplication ergodic theorem in the theory of random dynamical system, both invariant manifold of the nonlinear system and the relevant invariant subspace of its linearized system are tangent at the equilibrium point, thus the almost asymptotic stabilities for the two systems at the equilibrium point are the same. Therefore, in order to study the stability for a nonlinear stochastic system at the equilibrium point, it is sufficient we only investigate the stochastic stability of its corresponding linearized system within the vicinity of the equilibrium point by the results.

Via the transformation of r=(x12+x22)1/2, z=x3, Θ=arctan(x1/x2) in the vicinity of equilibrium point x'=(x1,x2,x3)=(0, 0, 0), the linearization of the original system Eq. (3), which is subjected to a stochastic parametric perturbation, is the following as:

(4)
x ˙ = A 0 x - ε 2 A 1 x + ε c o s ξ t B x , d ξ t = μ d t + σ d W t ,

where:

A 0 = 0 ω 0 - ω 0 0 0 0 0 ,       A 1 = δ 1 0 0 0 δ 1 0 0 0 δ 2 ,       B = b 11 b 12 b 13 b 21 b 22 b 23 b 31 b 32 b 33 .

The symbol ‘’ indicates that the equations within Eq. (4) are Stratonovich stochastic differential equations, μ and σ are real constants, W(t) is a unit Wiener process and the bifurcation parameters μ1, μ2 are rescaled such that μ1=-εδ1, μ2=-εδ2.

The following spherical polar transformation from (x1,x2,x3) to (ρ,θ,ϕ):

(5)
x 1 = a c o s θ s i n ϕ ,       x 2 = a c o s θ c o s ϕ ,       x 3 = a s i n θ ,       ρ = l n a ,       P = a p ,
ϕ t = ω t + φ t ,         θ - π 2 , π 2 ,         ϕ ,   φ 0 ,   2 π ,

yields a set of equations of the arguments of ρ, θ, ϕ and the noise process ξ(t), i.e.:

(6)
d P = [ ε c o s ( ξ ) p P ρ 1 + ε 2 p P ρ 2 ] d t ,
d θ = [ ε c o s ( ξ ) θ 1 + ε 2 θ 2 ] d t ,
d ϕ = [ ω + ε c o s ( ξ ) ϕ 1 ] d t ,
d ξ = μ d t + σ d W t ,

where:

ρ 1 = 1 2 f 12 + f 31 s i n 2 θ + f 11 c o s 2 θ + f 32 s i n 2 θ ,      
ρ 2 = - δ 1 c o s 2 θ - δ 2 s i n 2 θ ,
θ 1 = 1 2 f 32 - f 11 s i n 2 θ + f 31 c o s 2 θ - f 12 s i n 2 θ ,      
θ 2 = 1 2 δ 1 - δ 2 s i n 2 θ ,
ϕ 1 = f 21 + f 22 t a n θ ,
f 11 = 1 2 k 1 + k 2 c o s 2 ϕ + k 3 s i n 2 ϕ ,      
f 12 = b 13 s i n ϕ + b 23 c o s ϕ ,
f 21 = 1 2 k 4 + k 3 c o s 2 ϕ - k 2 s i n 2 ϕ ,      
f 22 = b 13 c o s ϕ - b 23 s i n ϕ ,
f 31 = b 31 s i n ϕ + b 32 c o s ϕ ,       f 32 = b 33 ,
k 1 = b 22 + b 11 ,       k 2 = b 22 - b 11 ,      
k 3 = b 12 + b 21 ,       k 4 = b 12 - b 21 .

For the norm process P, a reversible linear stochastic transformation is introduced, i.e.:

(7)
S = T θ , ϕ , ξ P ,         P = T - 1 θ , ϕ , ξ S ,       - π 2 θ π 2 ,       0 ϕ 2 π ,       0 < ξ < 2 π ,

where the function T(θ,ϕ,ξ) is a scalar function of the phase processes (θ,ϕ,ξ). Thus the equation for the new norm process S is obtained through lemma:

(8)
d S = P μ T ξ + 1 2 σ 2 2 T ξ 2 + ω T ϕ + ε p ρ 1 T + θ 1 T θ + ϕ 1 T ϕ c o s ( ξ )
            + ε 2 p ρ 2 T + θ 2 T θ d t + σ T ξ P d W .

Since T(θ,ϕ,ξ) is bounded and non-singular, then both P and S have the same stability. Therefore, T(θ,ϕ,ξ) is selected such that the drift term of Eq. (11) is independent of the phase processes θ, ϕ and ξ, i.e.:

(9)
d S = Λ p S d t + σ T ξ T - 1 θ , ϕ , ξ S d W .

Comparing between Eq. (8) and Eq. (9) yields a fact that T(θ,ϕ,ξ) is presented by the following equation:

(10)
Λ p T = μ T ξ + 1 2 σ 2 2 T ξ 2 + ω T ϕ + ε p ρ 1 T + θ 1 T θ + ϕ 1 T ϕ c o s ξ
            + ε 2 p ρ 2 T + θ 2 T θ .

The above equation is written as:

(11)
L ε p T θ , ϕ , ξ = Λ p T θ , ϕ , ξ ,

where:

(12)
L ε p = L 0 p + ε L 1 p + ε 2 L 2 p ,
L 0 ( p ) = μ ξ + σ 2 2 2 ξ 2 + ω ϕ ,
L 1 p = θ 1 θ + ϕ 1 ϕ + p ρ 1 c o s ξ ,
L 2 p = θ 2 θ + p ρ 2 .

And its corresponding adjoint operator is:

(13)
L 0 * = - μ ξ + σ 2 2 2 ξ 2 - ω ϕ ,  
L 1 * = - θ 1 θ + ϕ 1 ϕ + p ρ 1 c o s ξ ,
L 2 * = - θ 2 θ + p ρ 2 .

We can see that Eq. (11) defines an eigenvalue problem with the second-order differential operator, in which, Λ(p) is the eigenvalue and T(θ,ϕ,ξ) is the eigenfunction. Based on Eq. (11), one can easily find that this eigenvalue is the pth moment Lyapunov exponent of system Eq. (4).

Consider the operator Lε(p) and its adjoint operator Lε*(p), according to the facts shown by Arnold [3, 4], Λ(p) is an isolated simple eigenvalue of Lε(p) with non-negative eigenfunction T(θ,ϕ,ξ) and T(θ,ϕ,ξ)=1. For Lε*p,T*(θ,ϕ,ξ) is the unique eigenfunction corresponding to Λ(p) with the property of T(θ,ϕ,ξ),T*(θ,ϕ,ξ)=1, i.e.:

(14)
L ε p T θ , ϕ , ξ = Λ p T θ , ϕ , ξ ,
L ε * p T * θ , ϕ , ξ = Λ p T * θ , ϕ , ξ ,
T ( θ , ϕ , ξ ) , T * ( θ , ϕ , ξ ) = 1 ,         p R .

3. Asymptotic analysis

Because it is practically impossible that the expression of Λ(p)is solved by Eq. (14), perturbation method is applied here. The following asymptotic expansions Λ(p) and T(θ,ϕ,ξ) are assumed in advance, i.e.:

(15)
Λ p = Λ 0 p + ε Λ 1 p + ε 2 Λ 2 p + + ε n Λ n p + ,
T θ , ϕ , ξ = T 0 θ , ϕ , ξ + ε T 1 θ , ϕ , ξ + ε 2 T 2 θ , ϕ , ξ + + ε n T n θ , ϕ , ξ + .

Substituting Eq. (15) into Eq. (14) and equating the terms of the equal powers of ε, then the following recursion equations are obtained:

(16)
ε 0 :       L 0 p T 0 θ , ϕ , ξ = Λ 0 p T 0 θ , ϕ , ξ ,
ε 1 :       L 0 p T 1 θ , ϕ , ξ + L 1 p T 0 θ , ϕ , ξ = Λ 0 p T 1 θ , ϕ , ξ + Λ 1 p T 0 θ , ϕ , ξ ,
ε 2 :       L 0 ( p ) T 2 ( θ , ϕ , ξ ) + L 1 ( p ) T 1 ( θ , ϕ , ξ ) + L 2 ( p ) T 0 ( θ , ϕ , ξ ) = Λ 0 ( p ) T 2 ( θ , ϕ , ξ )
                  + Λ 1 p T 1 θ , ϕ , ξ + Λ 2 p T 0 θ , ϕ , ξ .

3.1. ε0 order perturbation

According to Eq. (16), the zero order perturbation equation is equavelant to:

(17)
L 0 p T 0 θ , ϕ , ξ = Λ 0 p T 0 θ , ϕ , ξ ,

i.e.:

(18)
σ 2 2 2 T 0 ξ 2 + μ T 0 ξ + ω T 0 ϕ = Λ 0 p T 0 .

In order to make the problem tractable, we assume θ, ϕ and ξ are mutually independent. Applying the method of variable separation, i.e. T0(θ,ϕ,ξ)=F0(θ)Φ0(ϕ)H0(ξ), and dividing both sides of Eq. (18) by Φ0(ϕ)H0(ξ), we obtain:

(19)
σ 2 2 H ¨ 0 ( ξ ) H 0 ( ξ ) + μ H ˙ 0 ( ξ ) H 0 ( ξ ) - Λ 0 = - ω Φ ˙ 0 ϕ Φ 0 ϕ = k .

Solving the equation for Φ0 yields Φ0(ϕ)=ke-cωϕ, where k and c are constants. Since Φ0(ϕ) is a periodic function of ϕ, we obtain c=0. Hence Φ0(ϕ) can be chosen as 1. On the basis of the property of moment Lyapunov exponent, we know Λ(0)=0, which is substituted into Eq. (15), then leads to that Λ0(0)=0. Because the left side of Eq. (18) does not include p, Λ0(0)=0 yields Λ0(p)=0. Thus the equation for H0(ξ) becomes:

(20)
σ 2 2 H ¨ 0 ( ξ ) H 0 ( ξ ) + μ H ˙ 0 ( ξ ) H 0 ( ξ ) = 0 .

We easily obtain:

(21)
H 0 ξ = C 0 + C 1 e x p - 2 μ σ 2 .

Since H(ξ) is bounded, C1=0 is required. So H(ξ) is a constant, we choose H(ξ)=1.

Thus, we obtain:

(22)
T 0 θ , ϕ , ξ = F 0 θ ,         θ - π 2 ,   π 2 ,       ϕ 0 ,   2 π ,       ξ 0 ,   2 π .

It is the joint probability density function of the phase processes (θ,ϕ,ξ).

The corresponding adjoint differential equation of Eq. (17) is the following:

(23)
μ ξ T 0 * - σ 2 2 2 ξ 2 T 0 * - ω ϕ T 0 * = 0 .

Also T0*(θ,ϕ,ξ)=F0*θG0*ϕH0*(ξ), Eq. (23) becomes:

(24)
σ 2 2 H ¨ 0 * ( ξ ) H 0 * ( ξ ) - μ H ˙ 0 * ξ H 0 * ξ = ω Φ ˙ 0 * θ Φ 0 * θ = κ .

Solving the above equation, we obtain Φ0*(ϕ)=c1eκ/ω. Where κ and c1 are constants. Because Φ0(ϕ) is a periodic function of ϕ, then κ=0, and Φ0*(ϕ) is selected:

Φ 0 * ϕ = 1 2 π ,       ϕ 0 ,   2 π ,

which represents the probability density function of a uniformly distributed random variable ϕ[0, 2π].

Therefore, Eq. (24) is simplified as:

(25)
σ 2 2 H ¨ 0 * ( ξ ) H 0 * ( ξ ) - μ H ˙ 0 * ξ H 0 * ξ = 0 .

The solution of Eq. (25) is H0*(ξ)=B0+B1exp(2μ/σ2)

Since H*(ξ) is bounded, it is required that B1=0. ξ(t) is the angle of triangular function, and cosine is the periodic function with 2π, so we can choose H*(ξ)=1/2π, ξ[0, 2π]. Thus:

(26)
T 0 * θ , ϕ , ξ = 1 4 π 2 F 0 * θ ,         θ - π 2 ,   π 2 ,       ϕ 0 ,   2 π ,       ξ 0 ,   2 π ,

which is the joint probability density of the independent random variables (θ,ϕ,ξ).

3.2. ε1 order perturbation

From Eq. (16), the first order perturbation is the following:

(27)
L 0 p T 1 θ , ϕ , ξ + L 1 p T 0 θ , ϕ , ξ = Λ 0 p T 1 θ , ϕ , ξ + Λ 1 p T 0 θ , ϕ , ξ .

Substituting Λ0(p)=0 and T0(θ,ϕ,ξ)=F0(θ) into Eq. (27) results as follows:

(28)
L 0 p T 1 θ , ϕ , ξ = Λ 1 p T 0 θ , ϕ , ξ - L 1 p T 0 θ , ϕ , ξ .

The solvability condition of Eq. (28) is:

(29)
Λ 1 ( p ) T 0 ( θ , ϕ , ξ ) - L 1 ( p ) T 0 ( θ , ϕ , ξ ) , T 0 * = 0 ,

where T0* is shown in Eq. (26), ,  denotes the inner product that is defined as follows:

S 1 , S 2 = 0 2 π d ϕ - π / 2 π / 2 d θ 0 2 π S 1 ( θ , ϕ , ξ ) S 2 ( θ , ϕ , ξ ) d ξ .

Solving Eq. (29) obtains that the first order term of the moment Lyapunov exponent, i.e.:

(30)
Λ 1 p = L 1 p T 0 θ , ϕ , ξ , T 0 * .

Again T0(θ,ϕ,ξ)=F0(θ), so by calculating it is obtained that:

(31)
L 1 p T 0 θ , ϕ , ξ = c o s ξ θ 1 F 0 ' θ + p ρ 1 F 0 θ .

And T0*(θ,ϕ,ξ)=F0*(θ)/4π2, Eq. (31) is rewritten as:

(32)
Λ 1 p = 1 4 π 2 cos ξ θ 1 F 0 ' θ + p ρ 1 F 0 θ , F 0 * θ .

Integrating Eq. (32) for ξ from 0 to 2π, it is obtained that Λ1(p)=0.

Hence, Eq. (27) reduces to:

(33)
L 0 p T 1 θ , ϕ , ξ = - L 1 p T 0 θ , ϕ , ξ ,

i.e.:

(34)
μ ξ + σ 2 2 2 ξ 2 + ω ϕ T 1 θ , ϕ , ξ = - c o s ξ θ 1 F 0 ' θ + p ρ 1 F 0 θ .

For the convenience to write, let F(θ,ϕ)=θ1F0'(θ)+pρ1F0(θ).

In order to solve the measure T1(θ,ϕ,ξ), we introduce an auxiliary time t' in Eq. (34) and make it become as:

(35)
t ' + σ 2 2 2 ξ 2 + μ ξ + ω ϕ T 1 θ , ϕ , ξ , t ' = c o s ξ t F θ , ϕ .

Applying the linear transformation t'=ψ+s, ϕ=ω(ψ+s), then Eq. (35) is:

(36)
s + σ 2 2 2 ξ 2 + μ ξ T 1 θ , ψ , ξ , s = c o s ξ t F θ , ω ψ + s .

According to Duhamel’s principle [20], the solution for Eq. (36) is given by:

(37)
T 1 θ , ψ , ξ , s = 0 s f θ , ψ , ξ , s ; r d r ,

where f(θ,ψ,ξ,s;r) is the solution of the following homogeneous equations:

(38)
s + σ 2 2 2 ξ 2 + μ ξ f θ , ψ , ξ , s ; r = 0 ,         s > r ,
f θ , ψ , ξ , r ; r = c o s ξ t F θ , ω ψ + s ,             s = r .

In order to solve Eq. (38), the following equations are considered:

(39)
s + σ 2 2 2 ξ 2 + μ ξ P ξ , s ; z , t = 0 ,         s < t ,
P ξ , s ; z , t = l i m s t P ξ , s ; z , t = δ z - ξ .

Eq. (39) are the Kolmogorov’s backward equations for the transition probability function P(ξ,s;z,t), which is the probability density function of random variable z(t) conditioned on ξ(s), t>s. The transition probability function with Eq. (4) is given:

(40)
P ξ , s ; z , t = 1 2 π t - s σ e x p - z - ξ + μ t - s σ 2 t - s .  

By Eq. (38) and Eq. (39), the solution for Eq. (38) is given by:

(41)
f ( θ , ψ , ξ , s ; r ) = F ( θ , ω ( ψ + r ) ) - + E c o s z r P ξ , s ; z , t d z ,

where:

(42)
E c o s z r = - + cos z r P ξ , s ; z , r d z = c o s ξ - μ r - s e x p - 1 2 σ 2 r - s .

After substituting Eq. (41) and Eq. (42) into Eq. (37) and via some calculation, T1(θ,ψ,ξ,s) is solved:

(43)
T 1 ( θ , ψ , ξ , s ) = e x p - 1 2 σ 2 t - s { c o s ( ξ ) 0 s F ( θ , ω ( ψ + r ) ) c o s [ μ ( r - s ) ] d r
            + s i n ( ξ ) 0 s F θ ,   ω ψ + r sin μ r - s d r } .

Meanwhile, the solution T1(θ,ψ,ξ) in Eq. (34) is obtained by inserting ϕ=ω(ψ+s) into Eq. (43) and calculating the limit s-.

3.3. ε2 order perturbation

According to Eq. (16), the second order perturbation is rewritten as:

(44)
L 0   T 2 = Λ 2 T 0 - L 1   T 1 - L 2   T 0   = Λ 2 T 0 - c o s ξ θ 1 θ + ϕ 1 ϕ 1 + p ρ 1 T 1
            - θ 2 F 0 ' θ - p ρ 2 F 0 θ .

The solvability condition of Eq. (44) is:

(45)
1 4 π 2 0 2 π d ϕ - π 2 π 2 d θ 0 2 π Λ 2 T 0 - c o s ξ θ 1 θ + ϕ 1 ϕ 1 + p ρ 1 T 1
            - θ 2 F 0 ' ( θ ) - p ρ 2 F 0 ( θ ) ] F 0 * ( θ ) d ξ = 0 .

Via an integral for ϕ on [0, 2π] and the massive calculations, Eq. (45) can be sorted into:

(46)
- π / 2 π / 2 L ( p ) - Λ 2 ( p ) F 0 ( θ ) F 0 * ( θ ) d θ = 0 ,
L p = 1 2 σ 2 θ d 2 d θ 2 + μ θ + p μ ^ θ d d θ + p q θ + 1 2 p 2 q ^ θ ,
σ 2 ( θ ) = 1 32 ( 2 α 1 2 σ 3 - α 2 σ 2 ) s i n 2 ( 2 θ )
            + 1 2 α 4 - α 5 c o s 2 θ + 1 2 α 3 c o s 4 θ σ 1 c o s 2 θ - c o s 4 θ ,
μ ( θ ) = 1 64 [ 9 α 6 μ 1 - ( 3 α 5 + 16 α 4 ) σ 1 + 4 ( α 2 σ 2 - ( k 1 2 - k 1 ) σ 3 ) + 32 Δ - ] s i n ( 2 θ )
            + 1 64 [ 2 α 1 2 σ 3 + 6 ( α 5 σ 1 + α 4 μ 1 ) - ( α 5 σ 1 + α 2 σ 2 ) + 4 α 4 σ 1 ] s i n ( 4 θ )
            + 1 64 α 6 μ 1 - 3 α 5 σ 1 s i n 6 θ + 1 128 α 3 σ 1 s i n 8 θ + 1 4 α 4 σ 1 t a n θ ,
μ ^ ( θ ) = 1 128 [ 2 ( α 2 σ 2 + 5 α 5 σ 1 ) - 8 ( k 1 2 - 4 b 33 2 ) σ 3
            - ( 11 b 32 2 - 16 b 23 2 - 6 b 31 2 + 16 b 13 2 ) σ 1 ] s i n ( 2 θ )
            - 1 64 [ 2 ( α 1 2 σ 3 + b 32 2 α 2 σ 1 ) + 4 ( α 4 + α 5 + b 23 2 ) σ 1 - α 2 σ 2 ] s i n ( 4 θ )
            + 1 128 2 α 3 + 2 b 31 2 + b 32 2 σ 1 s i n 6 θ ,
q ( θ ) = 1 2 [ α 6 μ 1 - ( α 4 + α 5 ) σ 1 ]
            + 1 4 [ α 1 2 σ 3 - α 2 σ 2 - 3 α 6 μ 1 + ( α 3 + 4 α 4 + 5 α 5 ) σ 1 - Δ - ] c o s 2 ( θ )
            + 1 8 α 2 σ 2 - 1 4 α 6 μ 1 + α 1 2 σ 3 + 5 α 5 - 3 α 3 - 2 α 4 σ 1 c o s 4 θ ,
q ^ ( θ ) = 1 2 b 33 2 σ 3 + 1 4 [ 2 α 1 2 σ 3 - ( α 3 + α 4 + 2 α 5 ) σ 1 ] c   o s 2 ( θ )
            + 1 16 [ 2 α 1 2 σ 3 - α 2 σ 2 + ( 4 α 3 + 4 α 4 + 8 α 5 ) σ 1 ] c o s 4 ( θ ) α 1
            = k 1 - 2 b 33 , α 2 = k 2 2 + k 3 2 ,
α 3 = b 31 2 + b 32 2 ,         α 4 = b 13 2 + b 23 2 ,         α 5 = b 13 b 31 + b 23 b 32 ,
α 6 = b 13 b 32 - b 23 b 31 ,         Δ ± = δ 1 ± δ 2 ,
σ 1 = σ 2 σ 4 + 4 ( ω + μ ) 2 + σ 2 σ 4 + 4 ( ω - μ ) 2 ,         σ 2 = σ 2 σ 4 + 4 ( 2 ω + μ ) 2 + σ 2 σ 4 + 4 ( 2 ω - μ ) 2 ,
σ 3 = - 2 σ 2 σ 4 + 4 μ 2 ,         μ 1 = 2 ω + μ σ 4 + 4 ( ω + μ ) 2 + 2 ω - μ σ 4 + 4 ( ω - μ ) 2 .

Due to the arbitrariness of the function F0*(θ), the bracketed expression in Eq. (46) must vanish identically, which yields the eigenvalue problem of the operator L(p), i.e.:

(47)
L p F 0 θ = Λ 2 p F 0 θ ,         θ - π 2 ,       π 2 .

4. Solution of the eigenvalue problem

According to Namachchivaya [6, 8], for the eigenvalue problem defined in Eq. (47), at the boundariesθ=±π/2, the eigenfunction F0(θ) meets the zero Neumann boundary condition. Then based on Wedig [21] and Bolotin [22], F0(θ) is thought as an orthogonal expansion of a Fourier cosine series, i.e.:

(48)
F 0 θ = m = 0 z m cos 2 m θ .

By substituting Eq. (48) into Eq. (47), and multiplying both sides of the equation by cos(2nθ), then integrating for θ on [-π/2,π/2], the following equations can be obtained:

(49)
m = 0 a n m z m = Λ 2 p z n ,       a n m = - π / 2 π / 2 [ L ( p ) c o s ( 2 m θ ) ] c o s ( 2 n θ ) d θ ,         n = 0 ,   1 ,   2 , .

Eq. (49) can be written as:

(50)
A Z = Λ 2 p Z ,
(51)
A = a 00 a 01 a 02 a 10 a 11 a 12 a 20 a 21 a 22 ,         Z = z 1 z 2 z 3 .

From Eq. (50), it can be seen that Λ2(p) is the leading eigenvalue of sub-matrices of the matrix A, whose order is varied with n, then Λ2(p) is also an infinite sequence of the eigenvalues. Therefore, solving Λ2(p) is converted into evaluating the eigenvalue problem of matrix A that is defined by Eq. (50). Meanwhile, for Z in Eq. (50), to guarantee the existence of the non-trivial solution, the determinant of the coefficients must vanishes. Moreover, the eigenvalue sequence obtained by this method converges to the determined eigenvalue which is Λ2(p) as n. However, the amount of calculation increases drastically with the increase of n, so we obtain the approximate eigenvalue by the truncation of n.

For example, for n=0, Λ2(p)=a00. As n=1, the second order approximation of Λ2(p) is the eigenvalue of the second order sub-matrix. Likewise, for n= 2, the third order approximation of Λ2(p) is the eigenvalue of the third order sub-matrix. Until the curves of Λ2(p) are almost coincident for different n, the curve can be as the approximation of Λ2(p). Due to the complexity of expressions in the matrix A, here we only give the elements of the second order sub-matrix:

(52)
a 00 = π 64 { 2 α 1 2 σ 3 - 12 [ ( α 4 + α 5 ) σ 1 - α 6 μ 1 ] - 5 α 2 - 32 Δ } p
            + π 128 6 k 1 2 + 4 b 33 2 + 4 3 k 1 b 33 σ 3 - 4 α 3 + α 4 + 2 α 5 σ 1 - 3 α 2 σ 2 p 2 ,
a 01 = π 128 [ ( 15 α 5 + 16 α 4 - α 3 ) σ 1 - 4 α 2 σ 2 - 17 α 6 μ 1 - 32 Δ - ] p
            + π 64 2 k 1 2 - 4 b 33 2 σ 3 - α 2 σ 2 p 2 ,
a 10 = π 64 [ ( 5 α 5 - 32 α 4 ) σ 1 - 9 α 6 μ 1 - 4 α 2 σ 2 + 4 ( k 1 2 - k 1 ) σ 3 - 32 Δ - ]
            + π 128 { 8 ( k 1 2 - 4 b 33 2 ) σ 3 - 8 α 2 σ 2 + 5 ( α 4 + α 5 ) σ 1 - 17 α 6 μ 1 - 32 Δ - } p
            + π 64 2 k 1 2 - 4 b 33 2 σ 3 - α 2 σ 2 p 2 ,
a 11 = π 256 [ ( 112 α 4 - 20 α 5 + 3 α 3 ) σ 1 + 12 α 6 μ 1 - 4 ( α 1 2 - k 1 ) σ 3 + 4 α 2 σ 2 ]
            + π 256 [ 2 α 1 2 σ 3 - 11 α 2 σ 2 - 20 ( α 4 + α 5 + 2 b 32 2 ) σ 1 + 28 α 6 μ 1 ] p - 32 Δ - } p
            + π 512 14 k 1 2 + 4 b 33 2 + 4 7 k 1 b 33 σ 3 - 4 α 3 + α 4 + 2 α 5 σ 1 - 7 α 2 σ 2 p 2 .

5. Numerical results

Because of very complex expressions of the elements in matrix A, It is almost impossible to obtain the analytical solution of moment Lyapunov exponent through eigenvalue problem defined in Eq. (50), especially, for high order matrix A. Therefore, we present the numerical solutions to eigenvalue problems of the finite order sub-matrices, and obtain the approximate numerical results of moment Lyapunov exponent by the relationship between moment Lyapunov exponent Λ2(p) and Λ(p) in Eq. (15).

In Fig. 1 and Fig. 2, the curves of moment Lyapunov exponent Λ(p) for the different n in two cases are given respectively. It can be seen that the deviation of the curves of approximate moment Lyapunov exponent becomes less and less with the increase of n, which represents the order of sub-matrix. In particular, in cases that as n= 3 and n= 4, the curves of Λ(p) are nearly overlap. Thus, we conclude that the results are convergent with increase of n and it is sufficient for us to calculate the four order approximate of Λ2(p). In addition, we obtain the moment Lyapunov exponent of the system Eq. (4) using Monte Carlo simulation [23] in order to verify the effectiveness of the results. As can be seen in Figs. 1-2, asymptotic analytical result of the moment Lyapunonv exponent is nearly consistent with its numerical simulation.

Fig. 1. Variation of moment Lyapunov exponent with n and p for the case: k1= 2, k2=k3= 0, k4= 2, b13=b23= 1, b31=b32= –2, b33= 1, μ=σ=ω= 1, δ1=δ2= –1

 Variation of moment Lyapunov exponent  with n and p for the case: k1= 2, k2=k3= 0,  k4= 2, b13=b23= 1, b31=b32= –2,  b33= 1, μ=σ=ω= 1, δ1=δ2= –1

Fig. 2. Variation of moment Lyapunov exponent with n and p for the case: k1= 2, k2= 0, k3= 3, k4=1, b13=b23=b31=1, b32=0, b33=1, μ=σ=ω=1, δ1=δ2= 1

 Variation of moment Lyapunov exponent  with n and p for the case: k1= 2, k2= 0, k3= 3, k4=1, b13=b23=b31=1, b32=0, b33=1,  μ=σ=ω=1, δ1=δ2= 1

Fig. 3. Variation of moment Lyapunov exponent with parameters p and μ for the case: k1= 2, k2= 0, k3= 3, k4=1, b13=b23=b31=1, b32=0, b33=1, σ=ω= 1, δ1=δ2= 1

 Variation of moment Lyapunov exponent  with parameters p and μ for the case: k1= 2,  k2= 0, k3= 3, k4=1, b13=b23=b31=1,  b32=0, b33=1, σ=ω= 1, δ1=δ2= 1

Fig. 4. Variation of moment Lyapunov exponent with parameters p and σ for the case: k1= 2, k2= 0, k3= 3, k4=1, b13=b23=b31=1, b32=0, b33=1, σ=ω= 1, δ1=δ2= 1

 Variation of moment Lyapunov exponent  with parameters p and σ for the case: k1= 2,  k2= 0, k3= 3, k4=1, b13=b23=b31=1,  b32=0, b33=1, σ=ω= 1, δ1=δ2= 1

In addition, it can be shown from the expressions of the elements in matrix A that the moment Lyapunov exponents are relevant to the parameters of the system and the noise excitation. Fig. 3 and Fig. 4 describe the effect of noise parameters on the moment Lyapunov exponent. One can easily find that the apexes of the moment Lyapunov exponents move to the left and rise as both μ and σ increase, which illustrates the stability of the system become poor with the increased μ and σ. Moreover, by comparing two figures, we can see the effect of μ is bigger than that of σ. The trend of the curves in Fig. 5 is similar to Fig. 3, but the varied range is less than the one in Fig. 3. Fig.6 displays that the curves of moment Lyapunov exponent transform from the right side of vertical axis to the left with the increase of δ1, and the bigger δ1 is, the larger the influence on the moment Lyapunov exponent. Finally, in view of Fig. 7, we see that the peak of the moment Lyapunov exponents moves to the left and ascends rapidly with increase of δ2, and the vary of the curves is also alike to Fig. 3. Therefore, it can be seen from Fig. 3-7 that the parameter μ, δ1 and δ2 play the significant influence on the moment Lyapunov exponent.

Fig. 5. Variation of moment Lyapunov exponent with parameters p and ω for the case: k1= 2, k2= 0, k3= 3, k4= 1, b13=b23=b31= 1, b32= 0, b33= 1, μ=ω= 1, δ1=δ2= 1

 Variation of moment Lyapunov exponent  with parameters p and ω for the case: k1= 2,  k2= 0, k3= 3, k4= 1, b13=b23=b31= 1,  b32= 0, b33= 1, μ=ω= 1, δ1=δ2= 1

Fig. 6. Variation of moment Lyapunov exponent with parameters p and δ1 for the case: k1= 2, k2= 0, k3= 3, k4= 1, b13=b23=b31= 1, b32= 0, b33= 1, μ=σ=ω= 1, δ2= 1

 Variation of moment Lyapunov exponent  with parameters p and δ1 for the case: k1= 2,  k2= 0, k3= 3, k4= 1, b13=b23=b31= 1,  b32= 0, b33= 1, μ=σ=ω= 1, δ2= 1

Fig. 7. Variation of moment Lyapunov exponent with parameters p and δ2 for the case: k1= 2, k2= 0, k3= 3, k4= 1, b13=b23=b31= 1, b32= 0, b33= 1, μ=σ=ω= 1, δ1= 1

 Variation of moment Lyapunov exponent with parameters p and δ2 for the case: k1= 2, k2= 0, k3= 3, k4= 1, b13=b23=b31= 1, b32= 0, b33= 1, μ=σ=ω= 1, δ1= 1

6. Conclusions

In this paper, the moment Lyapunov exponent of a co-dimension two bifurcation system, that is on a three dimensional center manifold and is parametrically excited by a bounded noise, is investigated. By a reversible linear transformation and perturbation method, a corresponding eigenvalue problem and each order perturbation of the moment Lyapunov exponent are obtained. Zero, one and two order perturbation are solved through differential equation and stochastic process theory, thus, the eigenvalue problem of probability density is yielded. Finally, the eigenvalue problem is solved by Fourier cosine series expansion, and an infinite matrix yields and whose leading eigenvalue is the second order perturbation of the moment Lyapunov exponent. Furthermore, the convergence of procedure is numerically verified in two cases. Finally, the effects of system parameters and noise parameters on the expression of the moment Lyapunov exponent are discussed.

References

  1. Molcanov S. A. The structure of eigenfunctions of one-dimensional unordered structures. Mathematics of the USSR – Izvestiya, Vol. 12, 1978, p. 69-101. [Search CrossRef]
  2. Arnold L. A formula connecting sample and moment stability of linear stochastic systems. SIAM Journal of Applied Mathematics, Vol. 44, 1984, p. 793-802. [Search CrossRef]
  3. Arnold L., Kliemann W., Oejeklaus E. Lyapunov Exponents of Linear Stochastic Systems. Lyapunov Exponents Lecture Notes in Mathematics, Vol. 1186, 1986, p. 85-125. [Search CrossRef]
  4. Arnold L., Oejeklaus E. Almost Sure and Moment Stability for Linear Equations. Lyapunov Exponents, Lecture Notes in Mathematics, Vol. 1186, 1986, p. 129-159. [Search CrossRef]
  5. Arnold L., Doyle M., Namachchivaya Sri N. Small noise expansion of moment Lyapunov exponents for two-dimensional systems. Dynamics and Stability of Systems, Vol. 12, 1986, p. 187-211. [Search CrossRef]
  6. Namachchivaya Sri N., Van Roessel H., Doyle M. Moment Lyapunov exponent for two coupled oscillators driven by real noise. SIAM Journal of Applied Mathematics, Vol. 56, 1996, p. 1400-1423. [Search CrossRef]
  7. Khasminskii R., Moshchuk N. Moment Lyapunov exponent and stability index for linear conservative system with small random perturbation. SIAM Journal of Applied Mathematics, Vol. 58, 1988, p. 245-256. [Search CrossRef]
  8. Namachchivaya Sri N. Moment Lyapunov exponent and stochastic stability of two coupled oscillators driven by real noise. Journal of Applied Mechanics, Vol. 68, 2001, p. 903-914. [Search CrossRef]
  9. Xie W.-C. Moment Lyapunov exponents of a two dimensional system under real noise excitation. Journal of Sound and Vibration, Vol. 239, 2001, p. 139-155. [Search CrossRef]
  10. Xie W.-C. Moment Lyapunov exponents of a two dimensional system under bounded noise parametric excitation. Journal of Sound and Vibration, Vol. 263, 2003, p. 593-616. [Search CrossRef]
  11. Liu X. B., Liew K. M. On the stability properties of a Van Der Pol-Duffing oscillator that is driven by a real noise. Journal of Sound and Vibration, Vol. 285, 2005, p. 27-49. [Search CrossRef]
  12. Higham D., Mao X., Yuan C. Almost sure and moment exponential stability in the numerical simulation of stochastic differential equations. SIAM Journal on Numerical Analysis, Vol. 45, 2007, p. 592-609. [Search CrossRef]
  13. Kozic P., Janevski G., Pavlovic R. Moment Lyapunov exponent and stochastic stability of a double-beam system under compressive axial loading. International Journal of Solids and Structures, Vol. 47, 2010, p. 1435-1442. [Search CrossRef]
  14. Li S. H., Liu X. B. Moment Lyapunov exponent for a three dimensional stochastic system. IUTAM Symposium on Nonliear Stochastic Dynamics and Control, Vol. 29, 2011, p. 191-200. [Search CrossRef]
  15. Li S. H., Liu X. B. Moment Lyapunov exponent for three-dimensional system under real noise excitation. Applied Mathematics and Mechanics, Vol. 34, 2013, p. 613-626. [Search CrossRef]
  16. Hu D. L., Huang Y., Liu X. B. Moment Lyapunov exponent and stochastic stability of binary airfoil driven by Non-Gaussian colored noise. Nonliear Dynamics, Vol. 70, 2012, p. 1847-1859. [Search CrossRef]
  17. Liew K. M., Liu X. B. The maximal Lyapunov exponent for a three-dimensional stochastic system. Journal of Applied Mechanics, Vol. 71, 2004, p. 677-690. [Search CrossRef]
  18. Li S. H., Liu X. B. The maximal Lyapunov exponent for a co-dimensional two-bifurcation system excited by bounded noise excitation. Acta Mechanica Sinica, Vol. 28, 2012, p. 511-519. [Search CrossRef]
  19. Guckenheimer G., Holmes P. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Field. Spring-Verlag, New York, 1983 [Search CrossRef]
  20. Zauderer E. Partial Differential Equations of Applied Mathematics. 2nd Edition, Wiley, New York, 1989. [Search CrossRef]
  21. Wedig W. Lyapunov Exponent of Stochastic Systems and Related Bifurcation Problems. Stochastic Structural Dynamics Progress in Theory and Applications. Elsevier Applied Science, New York, 1988, p. 315-327. [Search CrossRef]
  22. Bolotin V. V. The Dynamic Stability of Elastic Systems. Holden-Day, San Francisco, 1964. [Search CrossRef]
  23. Xie W. C., Huang Q. H. On the Monte Carlo simulation of moment Lyapunov exponents. Advances in Engineering Structures, Mechanics and Construction, Vol. 140, 2006, p. 627-636. [Search CrossRef]