Solution of Stochastic Van der Pol Equation Using Spectral Decomposition Techniques

In this paper, we introduce the study of the general form of stochastic Van der Pol equation (SVDP) under an external excitation described by Gaussian white noise. The study involves the use of Wiener-Chaos expansion technique (WCE) and Wiener-Hermite expansion (WHE) technique. The application of these techniques results in a system of deterministic differential equations (DDEs). The resulting DDEs are solved by the numerical techniques and compared with the results of Monte Carlo (MC) simulations. Also, we introduce a new formula that facilitates handling the cubic nonlinear term of van der Pol equations. The main results of this study are: 1) WCE technique is more accurate, programmable compared with WHE and for the same order, WCE consumes less time. 2) The number of Gaussian random variables (GRVs) is more effective than the order of expansion. 3) The agreement of the results with the MC simulations reflects the validity of the forms obtained through theorem 3.1.


Introduction
In recent years, vibration systems represented by oscillations have received considerable attention. One of these systems is the van der Pol (VDP) oscillator [1] [2] [3]. In 1920, the Dutch scientist named Van der Pol found stable oscillation in the electric circuit employing a vacuum tube [1]. This oscillation was modeled mathematically by a second-order differential equation. This differential equa- tion is used for modeling many biological and physical phenomena such as the oscillation of atoms, the contraction and the expansion of the human heart [4], etc. Furthermore, van der Pol discovered a new type of oscillation called the relaxation oscillation [5]. The relaxation oscillator is used in many applications such as an electronic beeper, blinking lights and others [6]. It is applied also to the nonlinear oscillation system such as gene activation systems [4] and earthquakes [3].
As it is known, accurate analytical solution of the VDP system is not obtained until now. So, the interest in studying the approximate solution of this system has increased lately, where the dynamic properties of such system have been investigated in [7]- [12], and the references therein. Using Liénard's theorem [13], it was proved that the deterministic VDP (DVDP) system has a unique limit cycle. The efforts made to study of this system are not just for the deterministic case, but also the stochastic case. Since the VDP system can be affected by an external force, this external force can be described by the Gaussian white noise. Therefore, the stochastic system can be modeled by a mathematical model and called the SVDP system. This system has been studied by many authors such as using the stochastic averaging method with the Fokker-Planck equation to find the probability density function of the stationary solution [14] [15] [16]. So, it is appropriate for studying this system using other numerical techniques such as the spectral decomposition techniques. The spectral decomposition techniques have been used for solving nonlinear stochastic differential equations (NSDEs). One of those techniques is the WHE technique which was suggested by Norbert Wiener [17]. The WHE was used by Meecham et al. [18] to study turbulence solution of Burger equation. Also, WHE was combined with the perturbation theory to solve the perturbed NSDE by El-Tawil, M. and his co-workers [19] [20] [21]. Also, El-beltagy et al. used this technique to study the higher order solution for many NSDEs [20] [21] [22] [23]. For, the second technique, it is called the WCE technique; it was developed by Cameron and Martin [24] in 1947. This technique is based on a discretization of the Gaussian white noise by using the Fourier expansion [25]. After applying the WCE to the stochastic differential equations (SDEs), a system of DDEs is obtained which is called the propagators [25] [27] [28]. Finally, using a simple deterministic numerical method such as the 5th order Rung Kutta [29] is an appropriate method for solving these propagators.
The aim of the current work is to study the dynamic behavior of SVDP equation and SVDP-Duffing equation using the stochastic spectral expansions; WCE and WHE. Some important theoretical results are deduced. The paper also studies the influence of the coefficients on the oscillatory motion and dynamic behavior limitation. The outline of this paper is: In Section 2 we investigate the mathematical equation of SVDP system and some basic theorems that are quite useful to deal with non-linear differential equations. In Section 3 we introduce the WCE method, its application on the SVDP equation and the error bound of the method, in addition to the analytic formula for the WCE of polynomial nonlinear term. In Section 4 we introduce the WHE method and its application on the SVDP equation. The numerical simulation for the results of the WCE method is shown in Section 5. The simulation of the comparison between the numerical results of WCE method, the WHE method and those of the MC simulations are shown in Section 6. In Section 7 the conclusions are summarized.

The Proposed Stochastic VDP Model
The van der Pol equation is a non-conservative system, where adding and subtracting energy from the system is done in the periodic motion called the limit cycle [12]. So the nonlinear oscillatory motion of a unit mass particle under external stochastic force is described by the following mathematical model: where the oscillation behavior of the SVDP Equation (2.1) is affected by three forces: 1) The damping force ( ) 3) The spring force x t is the displacement; ζ is the coefficient of the linear damping, ε is the coefficient of the nonlinear damping; λ is the forcing intensity, is the Gaussian white noise and , b β are deterministic coefficients. The three forces act in the opposite direction of the displacement. Equation (2.1) can be converted into a system of first order stochastic differential equations as follows: Existence and Uniqueness Theorem Suppose we have the following stochastic differential equation system, f t x y f t x y g t x y g t x y x x f t x y f t x y g t x y g t x y y y , , , , then the condition is satisfied and results in which means that ( ) x t should be a second-order (finite-variance) process.
Furthermore, the derivatives of the functions , f g are: 0 .
Then the functions , 1,2 By checking the growth condition (2.4), for 2 i = , consider the following: , which results in ( ) where x represents the solution of the stochastic Van der Pol Equation (2.1) and y x =  is the velocity. The inequality (2.7) means that the velocity should be bounded and the bound is inversely proportional to the displacement. In fact, the existence and boundedness of the solution depend on the damping coefficients ζ and ε . Also, the solution is expected to be restricted in a local area respected by condition (2.7). Furthermore, it is found that if x is larger than ζ ε , the nonlinear term will be dominant and the damping force will be positive. If x is smaller than ζ ε , the damping force will be negative because the nonlinear term will become negligible. On the other hand, to prove the uniqueness and stability of the limit cycle for the mean solution of system (2.2), we need to check the satisfaction of Lineard's theorem [13]. This theorem is used to prove that the mean solution of the SVDP equation has unique and stable limit cycle. Theorem 2.3 The second order differential equations are called Lineard's equations and take the form: this equation is a generalization of the van der Pol equation Lineard's theorem states that system (2.9) has a unique, stable limit cycle under the following conditions on Z and Q.
is an odd function and has exactly one positive zero at . So, the mean solution of the SVDP system (2.2) has unique and stable limit cycle in the phase plane surrounding the origin point. This result will be shown in Section 5.

The Wiener-Chaos Expansion Technique
The WCE is applied to the SVDP Equation ( where the multi-index set The application of WCE on the linear problem is rather easy than the nonlinear problem. So, the following theorem enables us to handle the cubic term in the Equation (3.1) using WCE. Theorem 3.1. Assume X and Y are two processes, they have Wiener chaos ≤ + ≥ , the above summation can be rewritten as: For simplicity, we denote , ; α α θ θ α θ γ so, the formula (3.4) will be: which can be written as: Consequently, Equation (3.6) can be rewritten as the following system: By truncating the system (3.7) for N order of polynomial chaos and K GRVs; then the set ℑ can be truncated as the multi-index set , K N ℑ which can be defined as:

The Wiener-Hermite Expansion Technique
The solution of Equation (2.1) can be expanded in terms of second order Wiener-Hermite functionals as, [17] [18] [19]: is the i th order Wiener-Hermite time-independent functional. The zero and first order terms x will account for the non-Gaussian part of the solution.
Higher order terms will be less dominant than the second order one and hence, can be neglected as in the current work.
Apply the second order WHE to Equation (2.1) to get: we can get the following deterministic differential equations:

t t x t t x t t t t x t t x t t x t t t t x t t x t t x t t t t x x t t x t t t x x t t x t t t x t t t x t t x t
For solving the deterministic differential Equations (4.3)-(4.5), we can use a numerical method such as the modified Euler method or the perturbation theory [19]. The perturbation technique has the disadvantage of small convergence interval and hence is not will applicable to many problems. Furthermore, the perturbation method is sensitive for the problems of the long time integral and needs complicated computations for higher order equations. So it is a bad method for solving the results of WHE for SVDP equation although it is appropriate for solving many other NSDEs. Consequently, we will solve the WHE resulting system using the modified Euler with ( ) 2 t ο ∆ .

Numerical Simulation Using WCE Technique
In this section, we analyze the results using WCE technique. Hence, we compare the results using different number of propagators and using MC simulations. For comparison, we will solve system (2.  [25]. So, we will choose a moderate intervals 50 T = and 30 T = . Accordingly, the oscillatory properties of SVDP system will be shown in these intervals through its statistical properties such as mean and variance.
From the comparison, we find that increasing the number of GRVs is more effective than increasing the order of polynomial. Also, the parametric values affect the dynamic behavior of the approximate solution through its statistical properties. The influence of the external force on the dynamic behavior of the system is greater than the influence of the damping force.
This influence needs to reduce the error (  tends to zero. Consequently, the oscillations of system (2.2) converge to the harmonic oscillation. Accordingly, the mean solution makes a periodic motion with a self-sustained limit cycle with quasi-absence of the external force. While the variance increases linearity with time. Furthermore, magnifying the value of ε reduces the amplitude x of the system and that makes the motion converges from the nature of damping with time. This is one of advantage of generalized the damping term, using two different coefficients for the linear damping term and non-linear damping term.   GRVs increases. We can note that, it is more effective to use WCE with higher number K of random variables than to increase the polynomial order N. Furthermore, they show the dynamic behavior of the oscillations at large value of linear damping coefficient ζ . With quasi-absence of the external force, the mean solution reacts as a negative damping system. Where, the system starts outside of the limit cycle and acquires energy. But, for 1 ζ  the system acts as a relaxation oscillator system as result of the sudden increase and decrease of the velocity x  , (see Figure 5(a), Figure  6(b)).
Also, Depending on conditions (2.5) and (2.7), we found that reducing the value of ε makes the nonlinear damping term 2 x x ε  tends to zero. Thence, the   . So, Figure 1(a) to Figure 6(a) have the same variation tendencies.
In Figures 7-9, the system acts as a positive damping system due to increase the influence of the external force. Therefore, the system dissipates the energy and converges to the origin. Depending on condition (2.7) where 2 2 1 0 x λ Furthermore, increasing the value of 1 λ  produces a deviation between the results of WCE and those of MC simulations. This deviation can be reduced by increasing MC samples, but at the other side, the computing time will increase and also the memory required. Especially we handle a second order differential equation. So, the oscillation of the results of MC simulations is been   chaotic. Figure 9 shows the limit cycles of the mean solution of the system (2.2).
Where the limit cycles dissipate the energy and converge rapidly to the origin point when the value of λ increases.

Comparison between the Results of WCE Technique and WHE Technique for SVDP Equation (2.1)
The following results discussed the comparison between the results of WCE method, WHE method with those of MC simulations. Although the existence of a small deviation of the results of WHE from the other methods, Figure 10 & Figure 11 show a good agreement with both MC simulations and WCE This deviation is due to using the numerical Euler method to solve the deterministic system of the WHE. This numerical method becomes better by using a very small time step by increasing the number of discretization of time. And, this is in some times impossible with the ability of the used machinery. In addition to, increasing the value of ζ makes the results of WHE deviate from the results of the others both techniques as in Figure 12. So, we need complicated computations for higher accuracy results of the WHE method. Therefore, we can deduce that the WCE gives good results than the WHE method. Also, we found that the dynamic behavior of the oscillation for the results of the three methods has the same parametric effect. Figure 13 shows the phase plane of the mean solution of the three methods for different values of ζ . Figure 14 shows the positive damping motion of the results of using both WCE and WHE as a result of increasing the effect of the stochastic external force. As it known, WCE method, MC simulations and WHE method are handling differently with the external force. Therefore, the deviation between the results of the three methods appears when λ is large.

Conclusions
In this paper, we investigated the SVDP equation under an external excitation described by Gaussian white noise. The accurate analytical solution of the SVDP oscillator was not obtained until now. Therefore we needed to use numerical techniques such as the spectral decomposition techniques to obtain an effective approximate solution. The dynamic behavior of this solution was studied through its statistical moments. Furthermore, depending on the study, it was found that: • WCE technique was more accurate if it is compared with WHE technique for the same order. • WCE technique was programmable if it is compared with WHE technique and for the same order, WCE consumed less time on the same processor.
• The accuracy of the WCE technique increases by increasing the number of Gaussian random variables (GRVs), much more than the order of the expansion. Hence, in creasing the number of GRVs makes the propagators capture more information about the stochastic solution. Especially, for equations need a long interval time.
• Depending on the growth condition, the solution of the SVDP system was locally existence as long as the values of both the linear damping coefficient and the external force coefficient were small. Furthermore, we deduced that as long as the influence of the external force is negligible, the results of the three techniques are converged. But, increasing the influence of the external force refers to the efficiency of each method in dealing with the randomness of external force. In that study, we found that the WCE method was more efficient since it depends on the discretization of the white noise using the Fourier expansion. This discretization was missed in the WHE method. For the MC simulations, it was based on the number of realizations, the increasing of this number may be exceeded the capacity of the used machinery.