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

Abstract

This paper introduces higher-order solutions of the stochastic nonlinear differential equations with the Wiener-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. A model 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.

1. 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–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 𝐿(𝑥(𝑡))=−𝜀𝑥^𝑛(𝑡)(𝑡)𝑁(𝑡),(]+𝑓+𝑔𝑡∈0,𝑇 (1) 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.

Theorem 1. The solution of (1), if exists, is a power series in 𝜀; that is, 𝑥(𝑡)=∑^∞_𝑖=0𝜀^𝑖𝑥_𝑖(𝑡).

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 𝐿𝑥^(𝑘)(𝑡)=−𝜀𝑥^(𝑘−1)(𝑡)^𝑛+𝑓(𝑡)+𝑔(𝑡)𝑁(𝑡),𝑘≥1, (2) apply the inverse operator, 𝐿⁻¹, on both sides to get 𝑥^(𝑘)(𝑡)=−𝜀𝐿⁻¹𝑥^(𝑘−1)(𝑡)^𝑛+𝐿⁻¹𝑓(𝑡)+𝑔(𝑡)𝑁(𝑡). (3) Let 𝑥⁽⁰⁾(𝑡)=𝐿⁻¹(𝑓(𝑡)+𝑔(𝑡)𝑁(𝑡)) which is the solution at 𝜀=0.
Then the Picard’s 𝑘th approximation is computed as 𝑥^(𝑘)(𝑡)=𝑥⁽⁰⁾(𝑡)−𝜀𝐿⁻¹𝑥^(𝑘−1)(𝑡)^𝑛. (4) Now, we need to prove that the 𝑘th approximation is a power series in 𝜀; that is, it can be written as 𝑥^(𝑘)(𝑡)=^𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡), (5) 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 𝑥⁽¹⁾(𝑡)=𝑥⁽⁰⁾(𝑡)−𝜀𝐿⁻¹([𝑥⁽⁰⁾(𝑡)]^𝑛), 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 𝑥^(𝑘+1)(𝑡)=𝑥⁽⁰⁾(𝑡)−𝜀𝐿⁻¹𝑥^(𝑘)(𝑡)^𝑛, (6) or, after substituting the power series of 𝑥^(𝑘)(𝑡), we get 𝑥^(𝑘+1)(𝑡)=𝑥⁽⁰⁾(𝑡)−𝜀𝐿^−1𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡)^𝑛. (7) The second term in the right hand side can be expanded using the multinomial theorem to get ^𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡)^𝑛=_ℎ𝑐_ℎ^𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡)^𝑗𝑖ℎ, (8) where 𝑐_ℎ=𝑛!/∏^𝑀𝑘_𝑖=0𝑗^𝑖_ℎ! and the counter ℎ runs over all the 𝑛+𝑀_𝑘𝑛 combinations of the positive integers 𝑗⁰_ℎ,𝑗¹_ℎ,…,𝑗^𝑀𝑘_ℎ such that ∑^𝑀𝑘_𝑖=0𝑗^𝑖_ℎ=𝑛. In a more simplified form, we can write ^𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡)^𝑛=_ℎ𝑐_ℎ𝜀^{𝑤ℎ𝑀𝑘}_𝑖=0𝑥_𝑖^(𝑘)(𝑡)^𝑗𝑖ℎ, (9) where 𝑤_ℎ=∑^𝑀𝑘_𝑖=0𝑖𝑗^𝑖_ℎ. Let ∏^𝑀𝑘_𝑖=0[𝑥_𝑖^(𝑘)(𝑡)]^𝑗𝑖ℎ=𝑣_ℎ(𝑡); then the expansion takes the form ^𝑀𝑘_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)(𝑡)^𝑛=_ℎ𝑐_ℎ𝜀^𝑤ℎ𝑣_ℎ(𝑡). (10) Substitute in Picard’s (𝑘+1)th approximation to get 𝑥^(𝑘+1)(𝑡)=𝑥⁽⁰⁾(𝑡)−𝜀𝐿⁻¹_ℎ𝑐_ℎ𝜀^𝑤ℎ𝑣_ℎ(𝑡) (11) or 𝑥^(𝑘+1)(𝑡)=𝑥⁽⁰⁾(𝑡)−_ℎ𝑐_ℎ𝜀^1+𝑤ℎ𝐿⁻¹𝑣_ℎ(𝑡), (12) 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 Section 4. The simulations up to third order and fourth correction are shown in Section 5.

2. 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] 𝑥(𝑡;𝑤)=𝑥⁽⁰⁾(𝑡)+^∞_−∞𝑥⁽¹⁾𝑡;𝑡₁𝐻⁽¹⁾𝑡₁;𝑤𝑑𝑡₁+^∞_−∞𝑥⁽²⁾𝑡;𝑡₁,𝑡₂𝐻⁽²⁾𝑡₁,𝑡₂;𝑤𝑑𝑡₁𝑑𝑡₂+^∞_−∞𝑥⁽³⁾𝑡;𝑡₁,𝑡₂,𝑡₃×𝐻⁽³⁾𝑡₁,𝑡₂,𝑡₃;𝑤𝑑𝑡₁𝑑𝑡₂𝑑𝑡₃+⋯ (13) or after eliminating the parameters, for the sake of brevity, we get 𝑥(𝑡;𝑤)=𝑥⁽⁰⁾(𝑡)+^∞_{𝑘=1𝑅𝑘}𝑥^(𝑘)𝐻^(𝑘)𝑑𝜏_𝑘, (14) where 𝑑𝜏_𝑘=𝑑𝑡₁𝑑𝑡₂⋯𝑑𝑡_𝑘 and ∫_𝑅𝑘 is a 𝑘-dimensional integral over the variables 𝑡₁,𝑡₂,…,𝑡_𝑘. 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 𝑥^(𝑗)(𝑡;𝑡₁,𝑡₂,…,𝑡_𝑗) 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 𝐻^(𝑛)(𝑡₁,𝑡₂…,𝑡_𝑛) 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: 𝐻^(𝑛)𝑡₁,𝑡₂,…,𝑡_𝑛=𝐻^(𝑛−1)𝑡₁,𝑡₂,…,𝑡_𝑛−1𝐻⁽¹⁾𝑡_𝑛−^𝑛−1_𝑚=1𝐻^(𝑛−2)𝑡_𝑛1,𝑡_𝑛2,…,𝑡_𝑛𝑛−2×𝛿𝑡_𝑛−𝑚−𝑡_𝑛, (15) with 𝐻⁽⁰⁾=1 and 𝐻⁽¹⁾(𝑡₁)=𝑁(𝑡₁): 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, 𝐸𝐻^(𝑖)𝐻^(𝑗)=0∀𝑖≠𝑗. (16) The average of almost all Wiener-Hermite functionals vanishes, particularly, 𝐸𝐻^(𝑖)=0∀𝑖≥1. (17) The expectation and variance of the solution will be 𝐸[𝑥(𝑡)]=𝑥⁽⁰⁾(𝑡),Var[𝑥(𝑡)]=^𝑚_𝑘=1(𝑘!)_𝑅𝑘𝑥_𝑖^(𝑘)2𝑑𝜏_𝑘. (18) 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 𝑥^(𝑖)(𝑡;𝑡₁,𝑡₂,…,𝑡_𝑖), 𝑖=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].(1)Truncate the expansion (14) to contain only 𝑚+1 (𝑚≥1) kernels 𝑥^(𝑗), 0≤𝑗≤𝑚; that is, 𝑥(𝑡;𝑤)=𝑥⁽⁰⁾(𝑡)+∑^𝑚_𝑘=1∫_𝑅𝑘𝑥^(𝑘)𝐻^(𝑘)𝑑𝜏_𝑘.(2)Substitute into the stochastic oscillatory equation (1).(3)Use the multinomial theorem to expand the quadratic nonlinear term in (1), 𝑥^𝑛, 𝑛=2.(4)Multiply by 𝐻^(𝑗), 0≤𝑗≤𝑚, and then apply the ensemble average. This will lead to (𝑚+1) equations in the kernels 𝑥^(𝑗), 0≤𝑗≤𝑚.(5)For each kernel 𝑥^(𝑗), 0≤𝑗≤𝑚, apply the perturbation technique up to NC corrections; that is, 𝑥^(𝑗)=∑^NC_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑗). (6)Equate the coefficients of 𝜀^𝑘, 0≤𝑘≤NC, in both sides to get NC+1 equations for each kernel 𝑥^(𝑗), 0≤𝑗≤𝑚.

This will lead to the following (𝑚+1)(NC+1) equations: 𝑗!𝐿𝑥₀^(𝑗)=𝛿_𝑗0𝑓(𝑡)+𝛿_𝑗1𝑔(𝑡)𝛿𝑡−𝑡₁,𝐿𝑥0≤𝑗≤𝑚,𝑗!_𝑏^(𝑗)=−_𝑓𝑐_𝑓𝐷^(𝑗)_{𝑓,𝑏−1}𝐸^𝑗_𝑓,0≤𝑗≤𝑚,1≤𝑏≤NC, (19)(20) where 𝐷^(𝑗)_{𝑓,𝑏−1}=_𝑅𝑧Var𝑐⁻_𝑔^𝑚_𝑖=0^NC_𝑝=0𝑥_𝑝^{(𝑖)ℎ𝑝𝑔}𝑑𝜏_𝑧. (21) And the expectations 𝐸^𝑗_𝑓 are computed as 𝐸^𝑗_𝑓=𝐻^𝑚(𝑗)_𝑖=0𝐻^{(𝑖)𝑘𝑗𝑓}. (22) 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 𝛿_𝑗𝑙=1,𝑗=𝑙0,𝑗≠𝑙. (23) The counter 𝑓, in the summation in the right hand side of (20), runs over all of the (_𝑛^𝑛+𝑚) combinations of the positive integers 𝑘⁰_𝑓,𝑘¹_𝑓,…,𝑘^𝑚_𝑓 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≤𝑗≤𝑚; then they are used to compute the other components in (20). For 𝑗=0, the component 𝑢₀⁽⁰⁾ is obtained by solving 𝐿(𝑥₀⁽⁰⁾)=𝑓(𝑡) with the original initial and boundary conditions. For 𝑗=1, the component 𝑥₀⁽¹⁾ is obtained by solving 𝐿(𝑥₀⁽¹⁾)=𝑔(𝑡)𝛿(𝑡−𝑡₁) with zero initial and boundary conditions. The other components 𝑢₀^(𝑗), 𝑗≥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 𝑥_𝑖^(𝑗) is evaluated in terms of the previously computed components 𝑥_𝑘^(𝑝), 𝑝≤𝑗, 𝑘<𝑖. This means that the 1st corrections for all kernels, 𝑥₁^(𝑗),  0≤𝑗≤𝑚, are solved firstly and then the 2nd corrections for all kernels, 𝑥₂^(𝑗), 0≤𝑗≤𝑚,… up to the (NC)th corrections for all kernels 𝑥^(𝑗)_NC, 0≤𝑗≤𝑚.

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 𝐸[𝑥(𝑡)]=^NC_𝑖=0𝜀^𝑖𝑥_𝑖⁽⁰⁾,Var[𝑥(𝑡)]=^𝑚_𝑘=1(𝑘!)_𝑅𝑘^NC_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑘)2𝑑𝜏_𝑘. (24) If 𝑥^(𝑗)=∑^∞_𝑖=0𝜀^𝑖𝑥_𝑖^(𝑗), then it will be convergent if [16] |𝜀|≤𝑥_𝑖^(𝑗)𝑥^(𝑗)_𝑖+1 (25) for 𝑡∈[𝑡₀,𝑇]. This means that |𝜀| should obey an upper bound condition after which divergence is obtained.

3. The Equivalent Deterministic Equations

Apply the previous WHEP algorithm to get the following systems of equations of the quadratic (𝑛=2) nonlinear oscillatory equation and first-order (𝑚=1) Gaussian approximation and different number of corrections (NC). The initial conditions are assumed deterministic and hence only the zero-order and zero-correction kernel equation (𝐿(𝑥₀⁽⁰⁾)=𝑓(𝑡)) will have the initial conditions (𝑥₀ and 𝑥̇₀). Other kernels equations will have zero initial conditions.

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.

Table 1 

The equivalent deterministic system for first-order approximation with different number of corrections (NC).

Table 2 

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

Table 3 

The equivalent deterministic system for third-order approximation with different number of corrections (NC).

In case of zero initial conditions, 𝑓(𝑡)=0 and 𝑔(𝑡)=1; that is, RHS=−𝜀𝜔²𝑥²+𝑁(𝑡), and we will have the following reduced system of equations: 𝐿𝑥₀⁽⁰⁾=0,𝐿𝑥₀⁽¹⁾=𝛿𝑡−𝑡₁,𝐿𝑥₁⁽⁰⁾=−𝜔²_𝑅𝑥₀⁽¹⁾𝑡₁²𝑑𝑡₁,𝐿𝑥₁⁽¹⁾𝐿𝑥=0,₂₍₀₎=0,𝐿𝑥₂⁽¹⁾=−2𝜔²𝑥₁⁽⁰⁾𝑥₀⁽¹⁾𝑡₁,𝐿𝑥₃⁽⁰⁾=−𝜔²𝑥₁⁽⁰⁾²−2𝜔²_𝑅𝑥₀⁽¹⁾𝑡₁𝑥₂⁽¹⁾𝑡₁𝑑𝑡₁,𝐿𝑥₃⁽¹⁾𝐿𝑥=0,₄₍₀₎𝐿𝑥=0,₄₍₁₎=−2𝜔²𝑥₁⁽⁰⁾𝑥₂⁽¹⁾𝑡₁−2𝜔²𝑥₃⁽⁰⁾𝑥₀⁽¹⁾𝑡₁,𝐿𝑥₅⁽⁰⁾=−2𝜔²𝑥₁⁽⁰⁾𝑥₃⁽⁰⁾−2𝜔²_𝑅𝑥₀⁽¹⁾𝑡₁𝑥₄⁽¹⁾𝑡₁𝑑𝑡₁−𝜔²_𝑅𝑥₂⁽¹⁾𝑡₁²𝑑𝑡₁,𝐿𝑥₅⁽¹⁾=0. (26)

4. Oscillatory Equation and the Numerical Solver

For the linear oscillatory equation, 𝑥̈+2𝜔𝜁𝑥̇+𝜔²𝑥=𝑓(𝑡), (27) the linear operator 𝐿 will be 𝑑𝐿=²𝑑𝑡²𝑑+2𝜔𝜁𝑑𝑡+𝜔². (28) 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)=𝑥₀ and 𝑥̇₀(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 𝑥(𝑡)=ℎ(𝑡)∘𝑓(𝑡), (29) where ℎ(𝑡)=(1/𝜔_𝑑)𝑒^{−𝜔𝜉𝑡}sin(𝜔_𝑑𝑡) with 𝜔_𝑑=𝜔1−𝜉², assuming underdamping (𝜉<1).

For 𝑓(𝑡)=𝑒^−𝑡, the solution will be 𝑥(𝑡)=∫^𝑡₀ℎ(𝑡−𝜏)𝑓(𝜏)𝑑𝜏, which results in 𝑥(𝑡)=11−2𝜔𝜉+𝜔²×𝑒^−𝑡−𝑒^{−𝜔𝜉𝑡}cos𝜔_𝑑𝑡+1−𝜔𝜉𝜔_𝑑𝑒^{−𝜔𝜉𝑡}sin𝜔_𝑑𝑡. (30) The numerical solution can be obtained for a model equation and then used for all kernels in the proper sequence. The model equation in this case will take the form 𝑥̈+𝑎₁𝑥̇+𝑎₀𝑥=𝑓(𝑡), (31) where 𝑎₁ and 𝑎₀ are assumed constant. We can use any difference scheme, but as we are working with a white noise, the Dirac delta function 𝛿(𝑡−𝑎) is expected to appear in the equations of some kernels. So, an integral numerical scheme such as FEM or FVM will be more suitable as the integration of the Dirac delta function is easier to be handled. In our case, the finite volume method (FVM) will be considered. The time axis where 𝑡∈[0,𝑇] will be discretized into 𝑛 equal intervals of size Δ𝑡. The interval extending from 𝑡_𝑖−1 to 𝑡_𝑖 is taken as the control volume. This technique in one dimension will be equivalent to the trapezoidal integration rule. The second-order linear equation (31) will be decomposed into two simultaneous first-order equations. This can be done by substituting 𝑥̇=𝑦, and then we will have 𝑥̇=𝑦with𝑥(0)=𝑥₀,𝑦̇+𝑎₁𝑦+𝑎₀𝑥=𝑓(𝑡)with𝑦₀=𝑥̇(0). (32) Integrate (32) along the control volume to get 𝑥̇^𝑡𝑖_𝑡𝑖−1𝑑𝑡=^𝑡𝑖_𝑡𝑖−1𝑦𝑑𝑡. (33) Approximate the integral with the trapezoidal rule which is of accuracy proportional to (Δ𝑡)² to get 𝑥_𝑖=𝑥_𝑖−1+Δ𝑡2𝑦_𝑖−1+𝑦_𝑖. (34) Also, integrate (32) along the control volume and substitute with 𝑥_𝑖 from (34) to get 𝑦_𝑖=4𝐹_𝑖−1+𝑦_𝑖−14−2𝑎₁Δ𝑡−𝑎₀(Δ𝑡)²−4𝑎₀Δ𝑡𝑥_𝑖−14+2𝑎₁Δ𝑡+𝑎₀(Δ𝑡)², (35) where 𝐹_𝑖−1=∫^𝑡𝑖_𝑡𝑖−1𝑓(𝑡)𝑑𝑡=0.5(𝑓_𝑖−1+𝑓_𝑖)Δ𝑡. If the Dirac delta function appears in the right hand side, a special treatment for 𝐹_𝑖−1 is considered. In this case 𝐹_𝑖−1=∫^𝑡𝑖_𝑡𝑖−1𝑞(𝑡)𝛿(𝑡−𝑡₁)𝑑𝑡=𝑞(𝑡_𝑖) when 𝑡_𝑖=𝑡₁ and 𝐹_𝑖−1=0 when 𝑡_𝑖≠𝑡₁. Equation (35) will be solved at each node 𝑡_𝑖 to get 𝑦_𝑖 and then (34) computed to get 𝑥_𝑖.

The numerical solution can be validated with the exact solution in case of 𝑓(𝑡)=𝑒^−𝑡 (see (30)). With 𝑥₀=𝑥̇₀=0, 𝜔=1, 𝜉=0.5, 𝑇=10, and Δ𝑡=0.1, the comparison in Figure 1 shows that the numerical solution has sufficient accuracy with a relative error of 0.03% at Δ𝑡=0.1. The overall convergence rate of the developed scheme was tested, and it has a convergence order of 1.5 as shown in Figure 2.

Figure 1 

Comparison between the exact and the numerical solutions at Δ𝑡=0.1.

Figure 2 

The convergence order of the developed numerical solver.

5. 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.

(a)

(b)

(c)

Figure 3 

First-order response mean and variance for (a) 𝜀=0.1, (b) 𝜀=0.3, and (c) 𝜀=0.5.

(a)

(b)

Figure 4 

First-order response mean and variance for (a) 𝜀=0.7 and (b) 𝜀=1.0.

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 0.98, which is greater than 𝜀. This means that the condition of convergence is satisfied and hence the solution converges in this case. Also, the convergence condition is satisfied for 𝜀=0.3 and 𝜀=0.5. In case 𝜀=0.7, the ratio between the 4th and 5th corrections is 0.6 and between the 2nd and 3rd corrections is 0.58. Both ratios are lesser than 0.7, which means that we will have a divergent solution for 𝜀=0.7. Also, 𝜀=1.0 will produce a divergent solution.

Figures 5 and 6 show the first-order response mean and variance using the same parameters in Figures 3 and 4, but the envelope function 𝑔(𝑡) is taken as 𝑒^−0.5𝑡. As it is clear from the figures, the effect of the nonlinearity strength 𝜀 on the response variance is negligible. The multiplication of 𝑒^−0.5𝑡 attenuates the effect of the white noise and it becomes negligible with the time increase. The variance vanishes with the time and the solution becomes nearly deterministic. Also, the response variance is nearly invariant with the number of corrections.

(a)

(b)

(c)

Figure 5 

The first-order response mean and variance for (a) 𝜀=0.1, (b) 𝜀=0.3, and (c) 𝜀=0.5. Case of exponential function 𝑒−0.5𝑡 multiplied by the white noise.

(a)

(b)

Figure 6 

The first-order response mean and variance for (a) 𝜀=0.7 and (b) 𝜀=1.0. Case of exponential function 𝑒−0.5𝑡 multiplied by the white noise.

Figures 7 and 8 show the second-order response mean and variance using the same parameters in the first-order simulations, Figures 3 and 4. Figure 9 shows the second-order response mean and variance for longer time interval, 𝑇=20 seconds. This was done to ensure that 𝑇=10 seconds is a sufficient interval for the response mean and variance to reach their steady state values. For the second-order approximation, the response variance is computed practically as [𝑥(𝑡=Var)]_{𝑅1NC𝑖=0}𝜀^𝑖𝑥_𝑖⁽¹⁾²+2_𝑅2^NC_𝑖=0𝜀^𝑖𝑥_𝑖⁽²⁾²𝑑𝑡₂𝑑𝑡₁. (36) Figures 10 and 11 show the third-order response mean and variance with the same parameters used in the first- and second-orders simulations. For the third-order approximation, the response variance is computed practically, for time and memory saving, as [𝑥(𝑡=Var)]_𝑅NC𝑖=0𝜀^𝑖𝑥_𝑖⁽¹⁾²+2_𝑅^NC_𝑖=0𝜀^𝑖𝑥_𝑖⁽²⁾²+3_𝑅^NC_𝑖=0𝜀^𝑖𝑥_𝑖⁽³⁾²𝑑𝑡₃𝑑𝑡₂𝑑𝑡₁. (37) Table 4 shows a comparison between the computing times for different levels of correction for the third-order approximation. The computing times for the response mean and variance are also shown. The table displays also the order of integrals in each level. The time step was 0.2 seconds and the computations are done on Intel Core i5, 2.4 GHz, machine with W7, 32 bits. We can note that around 90% of the computing time is consumed in level 4 (4th correction level). This is due to the higher orders of integrals computed in this level. As the level of correction increases, the order of integral increases and hence the computational time will also increase.

Table 4 

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

(a)

(b)

(c)

Figure 7 

The second-order response mean and variance for (a) 𝜀=0.1, (b) 𝜀=0.3, and (c) 𝜀=0.5.

(a)

(b)

Figure 8 

The second-order response mean and variance for (a) 𝜀=0.7 and (b) 𝜀=1.0.

(a)

(b)

Figure 9 

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

(a)

(b)

(c)

Figure 10 

The third-order response mean and variance for (a) 𝜀=0.1, (b) 𝜀=0.3, and (c) 𝜀=0.5.

(a)

(b)

Figure 11 

The third-order response mean and variance for (a) 𝜀=0.7 and (b) 𝜀=1.0.

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.

(a) 1st correction (mean of 1st, 2nd, and 3rd orders is coincident) (variance of 2nd and 3rd orders is coincident)

(b) 2nd correction (mean of 1st, 2nd, and 3rd orders is coincident)

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.

(a) 3rd correction (mean of 2nd and 3rd orders is coincident)

(b) 4th correction (mean of 2nd and 3rd orders is coincident)

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.

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.

(a) 1st correction (mean of 1st, 2nd, and 3rd orders is coincident) (variance of 2nd and 3rd orders is coincident)

(b) 2nd correction (mean of 1st, 2nd, and 3rd orders is coincident)

Figure 14 

Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 1 and 2, 𝜀=1.0.

(a) 3rd correction (mean of 2nd and 3rd orders is coincident)

(b) 4th correction (mean of 2nd and 3rd orders is coincident)

Figure 15 

Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 3 and 4, 𝜀=1.0.

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].

6. 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.