Numerical Approximation of Higher-Order Solutions of the Quadratic Nonlinear Stochastic Oscillatory Equation Using WHEP Technique

This paper introduces higher-order solutions of the stochastic nonlinear differential equations with theWiener-Hermite expansion and perturbation (WHEP) technique. The technique is used to study the quadratic nonlinear stochastic oscillatory equation with different orders, different number of corrections, and different strengths of the nonlinear term. The equivalent deterministic equations are derived up to third order and fourth correction. Amodel numerical integral solver is developed to solve the resulting set of equations. The numerical solver is tested and validated and then used in simulating the stochastic quadratic nonlinear oscillatory motion with different parameters.The solution ensemble average and variance are computed and compared in all cases. The current work extends the use of WHEP technique in solving stochastic nonlinear differential equations.


Introduction
Analysis of the response of linear and nonlinear systems subjected to random excitations is of considerable interest to the fields of mechanical and structural engineering [1].Stochastic differential equations based on the white noise process provide a powerful tool for dynamically modeling complex and uncertain aspects.In many practical situations, it may be appropriate to assume that the nonlinear term affecting the phenomena under study is small enough; then its intensity is controlled by means of a small parameter, say  [2].
According to [3], the solution of stochastic partial differential equations (SPDEs) using Wiener-Hermite expansion (WHE) has the advantage of converting the problem to a system of deterministic equations that can be solved efficiently using the standard deterministic numerical methods.The main statistics, such as the mean, covariance, and higher-order statistical moments, can be calculated by simple formulae involving only the deterministic Wiener-Hermite coefficients.In WHE approach, there is no randomness directly involved in the computations.One does not have to rely on pseudorandom number generators, and there is no need to solve the stochastic PDEs repeatedly for many realizations.Instead, the deterministic system is solved only once.
The application of the WHE [4][5][6][7][8][9][10] aims at finding a truncated series solution to the solution process of a stochastic differential equation.The truncated series is composed of two major parts: the first is the Gaussian part which consists of the first two terms, while the rest of the series constitutes the non-Gaussian part.In nonlinear cases, there exist always difficulties in solving the resultant set of deterministic integrodifferential equations got from the applications of a set of comprehensive averages on the stochastic integrodifferential equation obtained after the direct application of WHE.Many authors introduced different methods to face these obstacles.Among them, the WHE with perturbation (WHEP) technique [4] was introduced using the perturbation technique to solve perturbed nonlinear problems.
The WHE was originally started and developed by Wiener in 1938 and 1958 [11].Wiener constructed orthonormal random bases for expanding homogeneous chaos depending on white noise and used it to study problems in statistical mechanics.Cameron and Martin [12] developed a more explicit and intuitive formulation for WHE (now known as Wiener chaos expansion, WCE).Their development is based on an explicit discretization of the white noise process through its Fourier expansion, which was missed in Wiener's original formalism.This approach is much easier to understand and more convenient to use and hence replaced Wiener's original formulation.Since Cameron and Martin's work, WHE has become a useful tool in stochastic analysis involving white noise (Brownian motion) [3].We will denote it by Imamura formulation [13].Another formulation was suggested and applied by Imamura and his coworkers [13,14].They have developed a theory of turbulence involving a truncated WHE of the velocity field.The randomness is taken up by a white-noise function associated, in the original version of the theory, with the initial state of the flow.The mechanical problem then reduces to a set of coupled integrodifferential equations for deterministic kernels.In [1], the WHE (Imamura formulation [13]) was used to compute the nonstationary random vibration of a Duffing oscillator which has cubic nonlinearity under white-noise excitation.Solutions up to second order are obtained by solving the equivalent deterministic system by an iterative scheme.El-Tawil and his coworkers [4-10, 15, 16] used the WHEP technique to solve a perturbed nonlinear stochastic diffusion equation and the quadratic and cubic nonlinear stochastic oscillatory equation with first-order approximations.
In [17], the analysis of nonlinear random vibration has been studied using several methods, such as equivalent linearization method [18], stochastic averaging method [19], the WHE approach with nonstationary excitations [1], the WHEP technique [15], eigenfunction expansions [20], and the method of detailed balance [21].All of the above methods are applied and used for nonlinear random oscillations of real systems subjected to random nonstationary (or stationary) excitations.
According to [5,6], quadrate oscillation arises through many applied models in applied sciences and engineering when studying oscillatory systems [22].These systems can be exposed to a lot of uncertainties through the external forces, the damping coefficient, the frequency, and/or the initial or boundary conditions.These input uncertainties cause the output solution process to be also uncertain.For most of the cases, getting the probability density function (p.d.f.) of the solution process may be impossible.So, developing approximate techniques (through which approximate statistical moments can be obtained) is an important and necessary work.There are many techniques which can be used to obtain statistical moments of such problems.The main goal of this paper is to introduce higher-order WHEP solutions and to suggest a numerical solver suitable to handle the equivalent deterministic system.
In [16], the WHEP technique is generalized to th nonlinearity, general order of WHE, and general number of corrections.Also, the extension to handle white noise in more than one variable and general nonlinearities are outlined.The generalized algorithm is implemented and linked to MathML [23] script language to print out the resulting equivalent deterministic system.
In the current work, the WHE formulation suggested by Meecham and his coworkers (Imamura formulation) is used to solve the stochastic nonlinear differential models of the form with the proper set of initial conditions which will be assumed to be deterministic.The differential operator  is a general linear operator.The nonlinearity is introduced as losses of degree  > 1 strengthened by a deterministic small parameter ().For the quadratic nonlinearity,  will be equal 2. The uncertainty is introduced through white noise () scaled by a deterministic envelope function ().The function () is a deterministic forcing function.Theorem 1 will be used in the derivation of the WHEP technique.
Proof.Using Picard's method which generates a sequence of approximations that converge to the solution; that is, the solution is obtained as () = lim  → ∞  () () where the th approximation is computed as apply the inverse operator,  −1 , on both sides to get Let  (0) () =  −1 (() + ()()) which is the solution at  = 0.
Then the Picard's th approximation is computed as Now, we need to prove that the th approximation is a power series in ; that is, it can be written as where   + 1 is the number of terms in the series.Using the mathematical induction, we need to prove that a power series solution will be obtained at  = 1 and at  + 1 provided that the th approximation is a power series solution.
At  = 1, the Picard's 1st approximation will be  (1) ), which is a power series in .Now we need to prove that  (+1) () is a power series in  given that  () () is a power series in .The Picard's ( + 1)th approximation is computed as or, after substituting the power series of  () (), we get The second term in the right hand side can be expanded using the multinomial theorem to get where  ℎ = !/∏   =0   ℎ ! and the counter ℎ runs over all the ( +   ) combinations of the positive integers  0 ℎ ,  1 ℎ , . . ., In a more simplified form, we can write where ; then the expansion takes the form Substitute in Picard's ( + 1)th approximation to get or which is a power series in .This completes the proof.We note that the previous theorem is applied also in case of deterministic force term; that is, () = 0. Also, the theorem is applied if the unknown function  is a function of more than one variable [16].
This paper is organized as follows.In Section 2, the WHEP technique is reviewed and the generalized WHEP derivation steps are outlined.The deterministic set of equations equivalent to the stochastic nonlinear oscillatory equation are tabulated in Section 3. The analytical and numerical solutions of the oscillatory equation are described in

WHEP Technique
As a consequence of the completeness of the Wiener-Hermite set [1], any arbitrary stochastic process can be expanded in terms of the Weiner-Hermite polynomial set, and this expansion converges to the original stochastic process with probability one.The solution function (; ) can be expanded in terms of Wiener-Hermite functionals as [4] (1) (;  1 )  (1) or after eliminating the parameters, for the sake of brevity, we get where integral over the variables  1 ,  2 , . . .,   .The first term in the expansion ( 14) is the nonrandom part or ensemble mean of the function.The first two terms represent the normally distributed (Gaussian) part of the solution.Higher-order terms in the expansion depart more and more from the Gaussian form [16].The Gaussian approximation is usually a bad approximation for nonlinear problems, especially when high-order statistics are concerned [18].The components  () (;  1 ,  2 , . . .,   ) are called the deterministic kernels of the WHE of ().They are functions of time and space variables and fully account for the time dependence of () as well as for its statistical properties [3]. is a random output of a triple probability space (Ω, , ), where Ω is a sample space,  is a -algebra associated with Ω, and  is a probability measure.For simplicity,  will be dropped later on.
The functional  () ( 1 ,  2 . . .,   ) is the th order Wiener-Hermite time-independent functional.The WH functionals form a complete set [1], and they satisfy the following recurrence relation for  ≥ 2: with  (0) = 1 and  (1) ( 1 ) = ( 1 ): the white noise.By construction, the Wiener-Hermite functionals are symmetric in their arguments and are statistically orthonormal with respect to the weighting function  −(1/2) ∑  =1  2 (  ) ; that is, The average of almost all Wiener-Hermite functionals vanishes, particularly, The expectation and variance of the solution will be The WHE method can be elementarily used in solving stochastic differential equations by expanding the solution as well as the stochastic input processes via the WHE.The resultant equation is more complex than the original one due to being a stochastic integrodifferential equation.Taking a set of ensemble averages together with using the statistical properties of the WHE functionals, a set of deterministic integrodifferential equations are obtained in the deterministic kernels  () (;  1 ,  2 , . . .,   ),  = 0, 1, 2, . ... To obtain approximate solutions of these deterministic kernels, one can use perturbation theory in the case of having a perturbed system depending on a small parameter .Expanding the kernels as a power series of , another set of simpler iterative equations in the kernel series components are obtained.This is the main algorithm of the WHEP algorithm [4].
The technique was successfully applied to several nonlinear stochastic equations; see for example [4,5,10].
The WHEP technique for general nonlinear exponent (), general order (), and general number of corrections (NC) follows the following steps [16].
This will lead to the following ( + 1)(NC + 1) equations: (!)  ( where And the expectations    are computed as It was explained in [16] how to get    in terms of the Dirac delta functions and then use them to reduce the integrals that appear in  () ,−1 .The summation ∑ Var means that all variations ℎ   , 0 ≤  ≤ , 0 ≤  ≤ NC, that satisfy the equality  = ∑  =0 ∑ NC =0 ℎ   are selected.This can be done be a searching technique.For these variations, the factors   =    !/∏ NC =0 ℎ   ! will be multiplied by each other to get  −  ; that is,  −  = ∏ Var   .The Kronecker delta functions are defined as The counter , in the summation in the right hand side of (20), runs over all of the ( +  ) combinations of the positive integers  0  ,  1  , . . .,    such that ∑  =0    = .Equations ( 19) and (20) can always be solved using the proper sequence.The first  + 1 equations of (19) are solved independently to get  () 0 , 0 ≤  ≤ ; then they are used to compute the other components in (20).For  = 0, the component  (0)  0 is obtained by solving ( (0) 0 ) = () with the original initial and boundary conditions.For  = 1, the component  (1)  0 is obtained by solving ( (1)  0 ) = ()( −  1 ) with zero initial and boundary conditions.The other components  () 0 ,  ≥ 2, will be zeros due to zero right hand side and zero initial and boundary conditions.Equation (20) specifies the solution sequence to be followed.The component  These results are consistent with the known results obtained using WHE.In WHE, higher-order kernels are driven by lower-order kernels, and at the bottom, the Gaussian kernels are driven by the random forcing directly.So, the lower-order kernels are usually dominant in magnitude [3].
The statistical properties of the solution will be calculated as , then it will be convergent if [16] || ≤             approximation and different number of corrections (NC).The initial conditions are assumed deterministic and hence only the zero-order and zero-correction kernel equation (( (0) 0 ) = ()) will have the initial conditions ( 0 and ẋ 0 ).Other kernels equations will have zero initial conditions.

The Equivalent Deterministic Equations
For the quadratic nonlinear oscillatory stochastic equation, the application of the WHEP technique will result in the following set of equations (Tables 1, 2, and 3).The equations will be written for the first, second, and third orders.For a certain correction level, the deterministic system will include also the equations from previous levels.

Oscillatory Equation and the Numerical Solver
For the linear oscillatory equation, the linear operator  will be This means that the linear oscillatory equation can be written as () = ().The parameter  is the undamped angular frequency of the oscillator and  is the damping ratio.The initial conditions will be considered as (0) =  0 and ẋ 0 (0) = ẋ 0 .In case of zero initial conditions, the exact solution that can be obtained using different methods such as the theory of linear differential equations or the Laplace transform will be the convolution where ℎ() = (1/  ) − sin(  ) with   = √1 −  2 , assuming underdamping ( < 1).

Results
The following output is simulated using the previous developed numerical solver.The solution (34) of the model equation ( 31) is used to get all kernels with the proper right hand side.The mean response and the response variance are then calculated from the kernels using (24).
Figures 3 and 4 show the first-order response mean and variance for different values of .The simulations are done for the case of zero initial conditions, zero deterministic excitation, and unit envelope function multiplied by the white noise.The angular frequency is  = 1 and the damping ratio is  = 0.5.For the response mean, the 3rd and 4th corrections are coincident.For the variance, the 2nd and 3rd corrections are coincident and also the 4th and 5th corrections.
As it is shown in the figures, the nonlinearity strength  greatly affects the amplitudes of the mean and variance.It should not be increased after a certain value to obtain a convergent solution.This value depends on the different parameters of the problem.As the nonlinearity strength  increases, we need higher number of corrections.For  = 0.1 and at  = 10 seconds, the ratio between the response mean of each correction and the proceeding one is around Figures 12 and 13 show a comparison, at  = 0.3, between the response mean and variance for different orders and different number of corrections.As it was described earlier, the solution is convergent at  = 0.3.As it is shown in the figures, the response variance is more sensitive to the approximation order than the response mean.The variance converges as the approximation order increases.This means that, for a convergent solution, higher-order approximations are required for accurate prediction of the stochastic response.
Figures 14 and 15 show a comparison, at  = 1.0, between the response mean and variance for different orders and different number of corrections.In this case, the solution is divergent.The response variance diverges as the approximation order increases and also as the correction level increases.
It is worth to note that the WHEP technique used in the current work can be extended to solve stochastic PDEs with white noise in multiple dimensions and of different colors as described in [16].

Summary and Conclusions
The mean response of the quadratic nonlinear oscillatory system subjected to nonstationary random excitation was investigated using WHEP technique.The equivalent deterministic equations are derived up to third order.The solution is approximated up to fifth correction for the first order and up to fourth correction for the second and third orders.Numerical integral solution of the equivalent deterministic set of equations was applied using the FVM.The numerical treatment is validated after comparing the results with the analytical solution.The numerical solver is utilized in simulating the mean and variance of the nonlinear stochastic oscillatory motion with higher order, higher number of corrections, and different strengths of the nonlinear term.The values of the nonlinearity strength required for convergent solution are estimated.It was found that the numerical solution is efficient, and higher-order approximations are required for accurate prediction of the stochastic response.

Figure 1 :Figure 2 :
Figure 1: Comparison between the exact and the numerical solutions at Δ = 0.1.

Figure 9 :
Figure 9: The second-order response mean and variance for  = 0.5 and time interval of 20 seconds.

Figure 12 :
Figure 12: Comparison between the response mean and variance for the first, second, and third orders with number of corrections 1 and 2,  = 0.3.

Figure 13 :
Figure 13: Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 3 and 4,  = 0.3.

Table 2 :
The equivalent deterministic systems for second-order approximation with different number of corrections (NC).

Table 4 :
The computing time for each level and for the response mean and variance compared with the total time.