Probability density and stochastic stability for the coupled Van der Pol oscillator system

Abstract: In this paper, the probability density and almost sure asymptotic stability of the coupled Van der Pol oscillator system under the noise excitations are investigated. Through the stochastic averaging method and slow changing process theorem, averaged Fokker–Planck–Kolmogorov equation and exact solution of the dynamical system are obtained. Especially, the marginal density functions of the system excited by the additive white noises are derived. Then, the effects of coupled parameters and noise parameters on the marginal density functions are discussed through the numerical figures. In addition, the almost sure asymptotic stability under the parametric excitation by means of the maximal Lyapunov exponent is studied, and the stable demarcation points about noise intensity are presented.


Introduction
A great deal of attention has been paid to the investigation of stochastic dynamical system excited by the noises since Itô theory was constructed (Itô, 1951). Lots of stochastic methods were established, and one of which was the stochastic averaging method. The classic stochastic averaging *Corresponding author: Quanxin Zhu, School of Mathematical Sciences and Institute of Finance and Statistics, Nanjing Normal university, Nanjing 210023, China E-mail: zqx22@126.com Reviewing editor: Zudi Lu, University of Southampton, UK Additional information is available at the end of the article ABOUT THE AUTHORS Quanxin Zhu is a professor of Nanjing Normal University, and he has obtained the Alexander von Humboldt Foundation of Germany. Zhu is an associate editor of MPE, JAM, TJMAA, JAAM, RAM and he is a reviewer of Mathematical Reviews and Zentralblatt-Math. Also, Zhu is a senior member of the IEEE and he is the lead guest editor of "Recent Developments on the Stability and Control of Stochastic Systems" in MPE. Zhu has obtained 2011 Annual Chinese and the award of one Hundred "The Most Inuential International Academic Paper". Moreover, he has been one of most cited Chinese researchers in 2014-2016, Elsevier. Zhu is a reviewer of more than 50 other journals and he is the author or coauthor of more than 110 journal papers. His research interests include random processes, stochastic control, stochastic differential equations, stochastic partial differential equations, stochastic stability, non-linear systems, Markovian jump systems, stochastic complex networks and finance mathematics.

PUBLIC INTEREST STATEMENT
The stochastic stability is an important topic of non-linear dynamical field. At present, we mainly focus on the study of the stochastic differential equation, stochastic dynamics and stability. Through the invariant measure, the stochastic P-bifurcation behaviours are analysed. By means of the Lyapunov exponent, Stochastic D-bifurcation, stability and stable boundaries are discussed. method has constituted a rigorous theoretical system and all the result has been concluded in (Lin & Cai, 1995). Then, Zhu promoted the theories for the Hamilton system excited by the white noise, bounded noise and wide band noise (Zhu, Huang, & Suzuki, 2001a, 2002Zhu, Huang, & Yang, 1997). Stochastic averaging method is a combination of stochastic averaging principle and Fokker-Planck-Kolmogorov (FPK) equation. Its merits are that after applying stochastic averaging method, the FPK equation is often greatly simplified. Moreover, the dimension of obtained non-linear system is less half than that of original one, which decreases the difficulty of solving FPK equation. The exact stationary solution of FPK equation were firstly generalized to some non-linear damping systems and multiple degree of freedom non-linear systems by Caughey (Caughey & Ma, 1982a, 1982bCaughey & Payne, 1967), then, was extended to some general non-linear dissipated and partially integrable Hamilton systems under external and parametric excitation of Gaussian white noises (Zhu, Cai, & Lin, 1990Zhu, Huang, & Suzuki, 2001b). Lately, the stochastic averaging method was improved to apply a non-linear oscillator system (Zhao, 2015), particularly, extended to the system excited by the fractional Gaussian noise excitation (Deng & Zhu, 2015). But so far, the exact solution of high-dimensional FPK equation is still a very difficult problem.
The almost sure asymptotic stability of a random system can be determined according to the sign of the maximal Lyapunov exponent. If the value of the maximal Lyapunov exponent is less than zero, the system is stable, otherwise, unstable (Arnold, 1998). A method for exact evaluation of the maximal Lyapunov exponent in a linear Itô stochastic differential equation was first created by Khasminskii (1967). Subsequently, many results for two-dimensional systems excited by the white noises or real ones were presented by some researchers (Arnold & Papanicolaou, 1986;Mitchell & Kozin, 1974;Nishoka, 1976). By applying the stochastic averaging method, Ariaratnam and Xie (1992) obtained the asymptotic expansions of maximal Lyapunov exponent for two coupled oscillators excited by a real noise. Through a perturbation method, a more general four-dimensional linear stochastic system was investigated by Doyle and Namchchivaya (1994). Then, the maximal Lyapunov exponents for a co-dimension two bifurcation system under a white noise and a bounded noise excitation were given by Li and Liu (2012) and Liu and Liew (2004). However, there have been no any results about the research of almost sure stability for above four-dimensional dynamical systems up to date.
In practice, due to the role of various modes, the dynamical systems mostly present as coupled multi-degree of freedom oscillators whose dynamical behaviours are also very complex, have aroused the wide attention of many scholars. Ariaratnam andXie (1992, 1994), Doyle and Namachchivaya (1994), Namachchivaya, Vedula, and Onu (2011) and Liu and Zhu (2013, 2014, 2015 researched this kind of the system. In this paper, the probability response and stochastic stability for coupled Van der Pol oscillator system excited by the white noises are investigated. In Section 2, the formulation of the problem is presented through the stochastic averaging method. The marginal probability densities of the system excited by the additive white noises are obtained, and the effects of system parameters on probability density are discussed in Section 3. In Section 4, the analytical expression of the maximal Lyapunov exponent is derived and the almost sure stability is given. The matter of the article is summarized in Section 5.

Model and formulation
Consider the following coupled Van der Pol oscillator system subjected to the noise excitations. The motion equation of the system is of the form: where ω i (i = 1, 2) is an irreducible natural frequency, and α j and β j (j = 0, 1, 2) are coupled parameters. ξ i (t)(i = 1, 2, 3, 4) is independently Gaussian white noises with zero mean and the intensity 2D i (D i ≥ 0). ɛ is a very small positive parameter. Thus, Equation (1) defines a quasi-Hamiltonian system with light damping. Let x 1 = x,ẋ 1 = x 2 , y 1 = y,ẏ 1 = y 2 , and introduce the Wong-Zakai correction terms, thus, Equation (1) is rewritten as the Itô stochastic differential equations, i.e.
where B i (t)(i = 1, 2, 3, 4) is a unit Wiener process. Obviously, the Hamiltonian system governed by Equation (2) is quasi-periodic and integrable as ɛ = 0 according to Zhu et al. (1997), and is supposed non-resonant. Thus, the corresponding Hamiltonian function of the Equation (2) is: It can be seen from the above equations that the model (2) is a quasi-integrable Hamilton system. Via the following transformation: and according to Itô's formula, the Itô stochastic differential equations about z i , i (i = 1, 2) are obtained where In Equation (5), z i (i = 1, 2) is slowly changing process and θ i (i = 1, 2) is fast changing process. According to Khasminskii's theorem (1968), z i (i = 1, 2) weakly converges to a two-dimensional diffusion process as → 0 in the time interval of the order ɛ −1 . If the two-dimensional limit process is still marked as z i (i = 1, 2), the averaged Itô stochastic differential equation for the process is in the following: and the corresponding FPK equation is:

The probability density of the coupled system
Because of the complexity of Equation (8), it is almost impossible that the exact solution can be deduced. Consequently, we only consider the coupled system is additively excited by the noises. When there is only external excitation, i.e. D 1 = D 3 = 0, one can calculate the exact stationary solution of Equation (8).
Under the compatibility condition of 2 1 ∕D 2 2 = 1 2 ∕D 4 1 = , the exact stationary solution of Equation (8) is derived as: where C is the normalized constant, ϕ(z 1 , z 2 )is a probability potential and its expression is as follows: Through integrating over z 2 from −∞ to +∞ about Equation (10), the marginal density function of z 1 is obtained. Then substituting z 1 = ( 2 1 x 2 1 + x 2 2 )∕2 1 , the stationary joint probability density of x 1 and x 2 is given as: Similarly, the stationary joint probability density of y 1 and y 2 is also obtained where z 2 = ( 2 2 y 2 1 + y 2 2 )∕2 2 .
In order to intuitively observe the variety of the density function, Figures 1-6 are displayed. Figures  1-2 give the change of the joint density function of x 1 and x 2 with the increase of parameters α 1 and β 2 , respectively. It can be seen from Figure 1 that the top of joint probability density is like a crater shape, which is stochastic analogue of limit cycle and the most probability value lies on the circle. Simultaneously, both the height and the diameter of crater decline with the increase of parameter α 1 . In Figure 2, the joint probability density appears the shape of mountain peak, and whose height falls with the increased β 2 .
The varieties of joint density function of y 1 and y 2 with the increase of parameters α 1 and β 2 are presented in Figures 3-4, respectively. Since the system (1) is coupled, the curved surfaces of joint probability densities are similar to Figures 1-2. In detail, Figure 3 is alike Figures 2 and 4 is close to Figure 1, which is because the status of α 1 and β 2 is identical in the coupled oscillator system.
The effects of noise intensities on the joint probability density function are described in Figures  5-6. Figures 5-6 show that the altitude of peaks and crater raises as the intensities of noises change from D 2 = 0.1, D 4 = 0.12 to D 2 = 0.12, D 4 = 0.144.

Almost sure asymptotic stability
When there is only parametric excitation in the system (1), thus, D 2 = D 4 = 0. In this section, we discuss the almost sure asymptotic stability of system (2). According to the theory of stochastic dynamics, the invariant manifold of non-linear dynamical system is tangent with the invariant subspace of its corresponding linear system in equilibrium point, and so almost sure asymptotic stability of the two systems is consistent, which is central conclusion of Oseledec multiplicative ergodic theorem in the theory of stochastic dynamics. Therefore, in order to investigate the stability of coupled Van der (10) (z 1 , z 2 ) = (4 0 1 z 1 + 1 z 2 1 )∕4D 2 + (4 0 2 z 2 + 2 z 2 2 )∕4D 4 + z 1 z 2 p(x 1 , x 2 ) = 2C √ D 2 ∕ 1 exp (D 2 z 2 + 2 0 1 2 ) 2 ∕D 2 1 − ( 2 z 2 2 + 4 0 2 z 2 )∕4D 4 . Pol oscillator system in the vicinity of equilibrium point, we only need to study the stability of its corresponding linear system. Consider Itô stochastic linearized differential equation of Equation (7) Introduce the transformation Obviously, i ∈ [0, 1], i = 0, 1. The stochastic differential equations for ρ and ϕ 1 can be obtained by

(d)
It can be seen from Equation (16) the process ϕ 1 is independent of variable ρ, and so ϕ 1 is a diffusion process. Thus, the corresponding FPK equation of (15) is as follows: There are two first kind singular boundaries ϕ 1 = 0, 1, and no other singular point exists at 0 < ϕ 1 < 1. Via a simple calculation, we can easily obtain the diffusion index of Equation (16)   index is 1. ϕ 1 = 0, 1 are both repulsively natural boundaries under the conditions 2c 1 > c 2 − , 2c 2 > c 2 − , thus, ϕ 1 (t) is ergodic in the interval (0,1) and the stationary probability density of Equation (17) exists. Through the direct integral, one obtains According to Khasminskii' theory (Deng & Zhu, 2015), the maximal Lyapunov exponent of system (14) is: where c 1 , c 2 , c ± are shown in Equation (17). The plan of maximal Lyapunov exponent and noise intensities D 1 and D 3 is given in Figures 7-8. The maximal Lyapunov exponents rise with the increased D 1 and D 3 , and are less than zero when D 1 < 0.5585 and D 3 < 0.3997, at this moments, the system (1) is almost sure asymptotic stable. Moreover, the maximal Lyapunov exponent increases more quickly in Figure 8 than in Figure 7 as noise intensities are increased.

Conclusion
In this paper, through applying the stochastic averaging method, the averaged FPK equation of coupled Van der Pol oscillator system under the noise excitation is obtained. Moreover, its exact marginal densities are deduced under the compatibility conditions. Then, the effects of the parameters in the system (1) on the marginal densities are discussed. From the Figures 1-6, it can be seen that the curved surfaces of density functions decline with the increase of parameters α 1 and β 2 , however, the influences of the noise intensities D 1 and D 3 are contrary to them. Finally, the almost sure asymptotic stability of the system (1) is researched through the maximal Lyapunov exponent, and whose analytical expression is derived. Furthermore, the stable boundary point about noise strengths D 1 and D 3 are presented in Figures 7-8.