Approximating the Density of Random Differential Equations with Weak Nonlinearities via Perturbation Techniques

We combine the stochastic perturbation method with the maximum entropy principle to construct approximations of the first probability density function of the steady-state solution of a class of nonlinear oscillators subject to small perturbations in the nonlinear term and driven by a stochastic excitation. The nonlinearity depends both upon position and velocity, and the excitation is given by a stationary Gaussian stochastic process with certain additional properties. Furthermore, we approximate higher-order moments, the variance, and the correlation functions of the solution. The theoretical findings are illustrated via some numerical experiments that confirm that our approximations are reliable.


Introduction and Motivation
The analysis of stochastic perturbations in nonlinear dynamical systems is a hot topic in applied mathematics [1,2] with many applications in apparently different areas such as control [3], economy [4] and especially in dealing with nonlinear vibratory systems. The study of systems subject to vibrations is encountered, for example, in Physics (in the analysis of different types of oscillators) and in Engineering (in the analysis of road vehicles, response of structures to earthquakes' excitations or to sea waves). The nature of vibrations in this type of systems is usually random because they are spawned by complex factors that are not known in a deterministic manner but statistically characterized via measurements that often contain errors and uncertainties. Although, oscillators in Physics and Engineering systems have been extensively studied in the deterministic case [5,6], and particularly, in the nonlinear case [7][8][9], due to the above-mentioned facts the stochastic analysis is more suitable since provides better understanding of their dynamics.
Here, X(t) denotes the position (usually of the angle w.r.t. an origin) of the oscillatory system at the time instant t, the parameter β is given by β := 2ξω 0 , being ξ the damping constant and ω 0 > 0 the undamped angular frequency, and finally, is a small perturbation (| | 1) affecting a nonlinear function of the position, g(X(t)). The expression X(t) + g(X(t)) is referred to as the nonlinear restoring term. The right-hand side term, Y(t), stands for an external source/forcing term (vibration) acting upon the system. In the setting of random vibration systems, Y(t) is assumed to be a stochastic process, termed stochastic excitation, having certain characteristics that in the present study will be specified later.
Notice that the nonlinear restoring term in Equation (1) involves the parameter , which determines the magnitude of the nonlinear perturbation, whose shape is given by g(X(t)). When = 0, Equation (1) describes a random linear oscillator. In [10], authors analyze this class of oscillators considering two cases for the stochastic source term Y(t), first when is Gaussian and, secondly, when it can be represented via a Karhunen-Loève expansion. In the case that = 0, the inclusion of the nonlinear term makes more difficult (even simply impossible) to exactly solve Equation (1). An effective method to construct reliable approximations of Equation (1) in the case that represents a small parameter is the perturbation technique [11][12][13][14][15]. In the stochastic setting, this method has been successfully applied to study different type of oscillators subject to random vibrations. After pioneer contributions by Crandall [16,17], the analysis of random vibration systems has attracted many researchers (see, for instance, in [15,18,19] for a full overview of this topic). In [20], approximations of quadratic and cubic nonlinear oscillators subject to white noise excitations are constructed by combining the Wiener-Hermite expansion and the homotopy perturbation technique. The aforementioned approximations correspond to the first statistical moments (mean and variance) because, as authors indicate in the introduction section, the computation of the probability density function (PDF) is usually very difficult to obtain. In [21], the authors extend the previous analysis to compute higher-order statistical moments of the oscillator response in the case the nonlinearity is only quadratic. The previous methodology is extended and algorithmically automated in [22]. In [23], the author considers the interesting scenario of an harmonic oscillator with a random mass and analyses important dynamic characteristics such as the stochastic stability and the resonance phenomena. To conduct that study, a new type of Brownian motion is introduced. The perturbation technique has also been used to approximate the first moments, mainly the mean and the variance, of some oscillators subject to small nonlinearities. The computational procedures of this method often requires amendments to the existing solution codes, so it is classified as an intrusive method. A spectral technique that allows overcoming this drawback is non-intrusive polynomial chaos expansion (PCE) in which simulations are used as black boxes and the calculation of chaos expansion coefficients for response metrics of interest is based on a set of simulation response evaluations. In the recent paper [24], authors design an interesting hybrid non-intrusive procedure that combine PCE with Chebyshev Surrogate Method to analyze a number of uncertain physical parameters and the corresponding transient responses of a rotating system.
Besides computing the first statistical moments of the response or performing a stability analysis of systems under stochastic vibrations, we must emphasize that the computation of the finite distribution (usually termed "fidis") associated to the stationary solution, and particularly of the stationary PDF, is also a major goal in the realm of vibratory systems with uncertainties. Some interesting contributions in this regard include [25,26]. In [25], the authors first present a complete overview of methods and techniques available to determine the stationary PDF of nonlinear oscillators excited by random functions. Second, nonlinear stochastic oscilators excited by a combination of Gaussian and Poisson white noises are fully analyzed. The study is based on solving the forward generalized Kolmogorov partial differential equation (PDE) using the exponential-polynomial closure method. The theoretical analysis is accompanied with several illustrative examples. In the recent contribution [26], authors propose a new method to compute a closed-form solution of stationary PDF of single-degree-of-freedom vibro-impact systems under Gaussian white noise excitation. The density is obtained by solving the Fokker-Planck-Kolmogorov PDE using the iterative method of weighted residue combined with the concepts of the circulatory and potential probability flows. Apart from obtaining the density of the solutions, it is worth to pointing out that in some recent contributions one also determines the densities of key quantities, that belong to Reliability Theory, like the first-passage time for vibro-impact systems with randomly fluctuating restoring and damping terms (see [27] and references therein).
In this paper, we address the study of random cross-nonlinear oscillators subject to small perturbations affecting the nonlinear term, g, which depend on both the position, X(t), and the velocity,Ẋ(t), Here, the stochastic derivatives are understood in the mean square sense [28] (Chapter 4). In our subsequent analysis, we will consider the case that g(X(t),Ẋ(t)) = X 2 (t)Ẋ(t) and the excitation Y(t) is a mean square differentiable and stationary zero-mean Gaussian stochastic process whose correlation function, Γ YY (τ), is known. On the other hand, assuming that Y(t) is a stationary and Gaussian stochastic process is a rather intuitive concept, which has been extensively used in both theoretical and practical studies [29,30]. Stationarity means that the statistical properties of the process do not vary significantly over time/space. This feature is usually met in a number of modeling problems as the surface of the sea in both spatial and time coordinates, noise in time in electric circuits under steady-state operations, homogeneous impurities in engineering materials and media, for example [28] (Chapter 3). Now, we list the main novelties of our contribution.
• We combine mean square calculus and the stochastic perturbation method to study a class of nonlinear oscillators whose nonlinear term, g, involves both position, X(t), as velocity,Ẋ(t), specifically, we consider the case g = g(X(t),Ẋ(t)) = X 2 (t)Ẋ(t). This corresponds to the most complicated case, usually termed cross-nonlinearity.

•
The oscillator is subject to random excitations driven by a stochastic process, Y(t), having the following properties: Y(t) is mean square differentiable and stationary zero-mean Gaussian. • We compute reliable approximations, not only of the mean, the variance, and the covariance (as is usually done), but also of higher moments (including the asymmetry and the kurtosis) of the steady-state of the above-described nonlinear oscillator. • We combine the foregoing information related to higher moments and the entropy method to construct reliable approximations of the probability density function of the steady-state solution. The approximation is quite accurate as it is based on higher moments.
To the best of our knowledge, this is the first time that stochastic nonlinear oscillators with the above-described type of cross-nonlinearities is studied using our approach, i.e., combining mean square calculus and the stochastic perturbation method. In this sense, we think that our approach may be useful to extend our study to stochastic nonlinear oscillators having more general cross-nonlinearities, in particular of the form g(X(t),Ẋ(t)) = X n (t)Ẋ m (t), for n ≥ 3 and m ≥ 2.
The paper is organized as follows. In Section 2, we introduce the auxiliary stochastic results that will be used throughout the whole paper. This section is intended to help the reader to better understand the technical aspects of the paper. Section 3 is divided into two parts. In Section 3.1, we apply the perturbation technique to construct a first-order approximation of the stationary solution stochastic process of model (2) with g(X(t),Ẋ(t)) = X 2 (t)Ẋ(t). In Section 3.2, we determine expressions for the first higherorder moments, the variance, the covariance, and the correlation of the aforementioned first-order approximation. These expressions will be given in terms of certain integrals of the correlation function of the Gaussian noise, Y(t), and of the classical impulse response function to the linearized oscillator associated to Equation (2). In Section 4, we take advantage of the results given in Section 3 to construct reliable approximations of the PDF of the stationary solution using the principle of maximum entropy. In Section 5, we illustrate all theoretical findings by means of several illustrative examples. Our numerical results are compared with Monte Carlo simulations and with the application of Euler-Maruyama numerical scheme, showing full agreement. Conclusions are drawn in Section 6.

Stochastic Preliminaries
For the sake of completeness, in this section we will introduce some technical stochastic results that will be required throughout the paper.
Hereinafter, we will work on a complete probability space (Ω, F , P), i.e., Ω is a sample space; F is a σ-algebra of sets of Ω, usually called events; and P is a probability measure. To simplify, we will omit the sample notation, so the input and the solution stochastic processes involved in Equation (2) The following result will be applied to calculate some higher-order moments of the solution stochastic process, X(t), of the random differential Equation (2), since as it shall be seen later, X(t) depends on a product of the stochastic excitation, Y(t), evaluated at a finite number of instants, say t 1 , Proposition 1 (p. 28, [28]). Let the random variables Y 1 , Y 2 , . . . , Y n be jointly Gaussian with zero mean, E{Y i } = 0, 1 ≤ i ≤ n. Then, all odd order moments of these random variables vanish and for n even, The sum above is taken over all possible combinations of n/2 pairs of n random variables. The number of terms in the summation is 1 · 3 · 5 · · · (n − 3) · (n − 1).
The two following results permit interchange the expectation operator with the mean square derivative and the mean square integral. In [28] (Equation (4.130) in Section 4.4.2), the first result is established for n = 2 and then it follows straightforwardly by induction.

Proposition 2. Let {Y(t)
: t ≥ 0} be a mean square differentiable stochastic process. Then, provided the above expectations exists.
Proposition 3 (p. 104, [28]). Let {Y(t) : −∞ ≤ a ≤ t ≤ b ≤ +∞} be a second-order stochastic process integrable in the mean square sense and h(t) a Riemann integrable deterministic function on t ∈ (a, b). Then, The following is a distinctive property of Gaussian processes since they preserve Gaussianity under mean square integration. Proposition 4 (p. 112, [28]). Let {Y(t) : a ≤ t ≤ ∞} be a Gaussian process and let h(t) be a Riemann integrable deterministic function on (a, t) such that the following mean square integral, exists, then {X(t) : t ≥ a} is a Gaussian process.

Probabilistic Model Study
As it has been indicated in Section 1, in this paper we will study, from a probabilistic standpoint, the random cross-nonlinear oscillator The analysis will be divided into two steps. First, in Section 3.1 we will apply the perturbation technique to obtain an approximation, X(t), of the stationary solution stochastic process, X(t). Then, in Section 3.2 we will take advantage of X(t) to determine reliable approximations of the main statistical functions of X(t), namely, the first higher-order moments, E{X n (t)}, n = 1, . . . , 5; the variance, V{X(t)}; the covariance, Cov{X(t 1 ), X(t 2 )}; and the correlation, Γ XX (τ).

Perturbation Technique
Let us consider the Equation (3). The main idea of the stochastic perturbation technique is to consider that the solution X(t) can be expanded in the powers of the small parameter (| | 1), Replacing expression (4) into Equation (3), yields the following sequence of linear differential equations, with random inputs Notice that each equation can be solved in cascade. As usual, when applying the perturbation technique, we take the first-order approximation This entails that in our subsequent development we will only need the two first equations listed in (5).
As indicated in Section 1, now we will focus on the analysis of the steady-state solution. Using the linear theory, the two first equations in (5) can be solved using the convolution integral [31]: and where is the impulse response function for the underdamped case ζ 2 < 1. This situation corresponds to the condition in which damping of an oscillator causes it to return to equilibrium with the amplitude gradually decreasing to zero (in our random setting it means that the expectation of the amplitude is null); system returns to equilibrium faster but overshoots and crosses the equilibrium position one or more times. Although, they are no treated hereinafter, two more situations are also possible, namely, critical damping and overdamping. The former corresponds to ζ 2 = 1 and in that case the damping of an oscillator causes it to return as quickly as possible to its equilibrium position without oscillating back and forth about this position, while the latter corresponds to ζ 2 > 1, and in this situation damping of an oscillator causes it to return to equilibrium without oscillating; oscillator moves more slowly toward equilibrium than in the critically damped system [32].

Approximation of the Main Statistical Moments
This subsection is devoted to calculate the main probabilistic information of the stationary solution stochastic process, X(t), of model (3). As it has been previously pointed out, to this end, we assume that the input term Y(t) is a stationary zero-mean (E{Y(t)} = 0) Gaussian stochastic process whose correlation function, Γ YY (τ), is given. We will further assume that Y(t) is mean square differentiable. This additional hypothesis will be apparent later. At this point, it is convenient to recall that for any stationary stochastic process its correlation function is even, so Γ YY (τ) = Γ YY (−τ), (p. 47, [28]). This property will be extensively applied throughout our subsequent developments.
To compute the mean of the approximation, we first take the expectation operator in (6), Therefore, we now need to determine both E{X 0 (t)} and E{X 1 (t)}. To compute the E{X 0 (t)} we again use the expectation operator in (7), where we have applied Proposition 3 and that E{Y(t)} = 0. Now, we deal with the computation of E{X 1 (t)} in an analogous manner but using the representation of X 1 (t) given in (8), Notice that the assumption of mean square differentiability of the input process Y(t) appears naturally at this stage.
Let us justify the last step in expression (12). Let us denote u 1 = t − s − s 1 , u 2 = t − s − s 2 and u 3 = t − s − s 3 , then applying Propositions 2 and 1, both with n = 3, one gets Therefore, substituting (11) and (12) into (10), we obtain the expectation of the approximation is null, From the approximation (6) and neglecting the term 2 , the second-order moment for X(t) is given by The first addend can be calculated using expression (7) and Fubini's theorem, Notice that we have used that Y(t) is a stationary process, so Now, we calculate the second addend in (14). To this end, we substitute the expressions of X 0 (t) and X 1 (t) given in (7) and (8), respectively, Observe that in the step (I) of the above expression, we have first applied Proposition 2 and second Proposition 1. Indeed, let us denote by and now we apply Proposition 1, with n = 4, to the right-hand side. This yields In step (II) of expression (16) we have taken advantage of the symmetry of the indexes. Then, substituing (15) and (16) in (14) one gets Notice that E X 2 (t) does not depend on t. This is consistent with the fact that we are dealing with the stochastic analysis of the stationary solution. The same feature will hold when computing higher-order moments, E X n (t) , n > 2, later.
As E X(t) is null (see (13)), then the variance of the solution coincides with E X 2 (t) . Now, we calculate the third-order moment of X(t) keeping up to the first-order term of perturbation . Therefore, Reasoning analogously as we have shown before, we obtain where we have applied Proposition 1 in the last step. The second addend in (18) is calculated using Propositions 1 and 2, ·Ẏ(t − s − s 5 ) ds 5 ds 4 ds 3 ds 2 ds 1 ds = 0.
Indeed, we know that, On the one hand, let us observe that applying first Fubini's theorem and Proposition 3, and second Proposition 1 for n odd, one gets On the other hand, using the same reasoning as in (20), where first we have applied Proposition 2, in order to put the first derivative out of the expectation, and second, we have utilized that X n−1 0 (t), X 2 0 (t − s) andẊ 0 (t − s) depend upon n − 1, 2 and 1 terms of Y(·), respectively, together with Proposition 1 (notice that n + 2 is odd as n is odd).
To complete the information of statistical moments of the approximation, we also determine E X 4 (t) .
The fourth-order moment of X(t), based on the first-order approximation via the perturbation method, is given by Reasoning analogously as we have shown in previous sections, we obtain for the first addend and for the second addend ) ds ds 1 ds 2 ds 3 ds 4 ds 5 ds 6 . (25) Observe that in the last step of the above expression, first we have used Proposition 2, and second, Proposition 1. From this last proposition, we know that exist 15 combinations, but we can reduce the expression by the symmetry of involved indexes. Now we deal with the approximation of the correlation function of X(t) via (6), i.e., taking the first-order approximation of the perturbation expansion, The first addend in (26) corresponds to the correlation function of X 0 (t). It can be expressed as The two last addends in (26) represent the cross-correlation of X 0 (t) and X 1 (t). They are given, respectively, by Summarizing, As E{ X(t)} = 0, we observe that the covariance and correlation functions of

Approximating the PDF via the Maximum Entropy Principle
So far, we have calculated approximations of the moments E{ X n (t)}, n = 1, . . . , 5 to the first-order approximation, X(t), via the perturbation method, of the steady-state solution of the random nonlinear oscillator (3). Although this is an important information, a more ambitious goal is the approximation of the PDF, say f X(t) (x), as from it one can calculate key stochastic information as the probability that the output lies in a specific interval of interest, say [a 1 , a 2 ], for any arbitrary fixed time t. Furthermore, from the knowledge of the PDF one can easily compute confidence intervals at a specific confidence level α ∈ (0, 1), where µ X (t) = E{ X(t)} = 0 (see (13)) and σ X (t) = V{ X(t)}. Usually α is taken as α = 0.05 so that 95% confidence intervals are built, and k > 0 must be determined numerically.
As we have calculated the approximations E{ X n (t)}, n = 1, . . . , 5, a suitable method to approximate the PDF, f X(t) (x), is the Principle of Maximum Entropy (PME), [33]. For t fixed, the PME seeks for a PDF, f X(t) (x), that maximizes the so-called Shannon's Entropy, of random variable X(t) with support [a, b], defined via the following functional, Condition (29) guarantees f X(t) (x) is a PDF, and the M conditions given in (30) impose that the sampled moments, m n , match the moments, E X n (t) , obtained in our setting by the stochastic perturbation method. For each t fixed, the maximization of functional (28) subject to the constrains (29)- (30) can be solved via the auxiliary Lagrange function where m 0 = 1. It can be seen that the form of the PDF is given by [33] In Section 3, we have approximated, via the stochastic perturbation technique, the moments E X n (t) for n = 1, 2, . . . , 5. Therefore, to apply the PME we will take M = 5 in (30). Notice that, in practice, to calculate the parameters λ i , i = 0, 1, . . . , 5, we will need to numerically solve the system of nonlinear Equations (29) and (30).

Numerical Examples
This section is devoted to illustrate the theoretical findings obtained in previous sections. We take the following data for the parameters of the random nonlinear oscillator (3), ζ = 0.05 (ζ 2 < 1) and ω 0 = 1.
Now we shall obtain approximations to the first moments, E{ X i (t)}, i = 1, . . . , 5, the correlation function and the variance, V{ X(t)}, of the approximate solution X(t) of random nonlinear oscillator (31).
This provides a stronger bound, < 0.000166641. Now, applying (27), we obtain the following approximation of the correlation function, In Figure 1, we show the graphical representation of the correlation function, Γ X X (τ), given in the expression (34) for different values of . We can see the higher the perturbation , the lower the variability. This graphical behavior is in full agreement with the physical interpretation of the oscillator dynamics. Indeed, let us rewrite Equation (31) as follows, X(t) + (0.1 + X 2 (t))Ẋ(t) + X(t) = ξ 1 cos(t) + ξ 2 sin(t).
As > 0 increases, the damped coefficient 0.1 + X 2 (t) does, so the mechanical system reduces its oscillations. It should be noted that = 0.0004 only satisfies the first bound ( < 0.0005); however, we can observe that the corresponding approximation preserves symmetry of correlation function. This might be due to the sample regularity of the random excitation, Y(t), which is differentiable. For the approximation of the PDF, f X(t) (x), we apply the results exhibited in Section 4 based on PME by taking = 0.00005, which satisfies the stronger bound previously determined ( < 0.000166641). We first compute the approximation based on the three first moments In Figure 2, we compare both graphical representations. From them, we can observe that both plots are quite similar, so giving evidence that computations are consistent. Finally, to check that our approximations are reliable, we compare the mean and standard deviation of the approximate solution obtained via the perturbation method against the ones calculated by Monte Carlo. The results are collected in Table 1. We can observe that both approximations agree. Example 2. To previously perform our theoretical analysis, we have required the stationary Gaussian stochastic excitation Y(t) be differentiable in the mean square sense (or equivalently, its correlation function, Γ Y (τ), be twice differentiable in the ordinary sense at τ = 0 [28] (Chapter 4)), so having differentiable sample trajectories [34]. If we carefully revise our previous development, we can notice this is an hypothesis coming from the fact the nonlinearity cross-term depends upoṅ X(t). In this second example, we shall show that using the general concept of differentiability, in the sense of distributions, we can still obtain good results via the perturbation techniques when the excitation is not differentiable. To this end, we have chosen, Y(t) = ξ(t), a Gaussian white-noise process with zero-mean and correlation function Γ YY (τ) = 1 2 Wδ(τ), where δ(τ) is the Dirac delta function and W is the noise power. This type of random noise has been extensively used in the literature since the earliest contributions [17]. Observe that Y(t) = ξ(t) is a stationary zero-mean Gaussian process but is not mean square differentiable (as its correlation function, given by the Dirac delta function, is not differentiable) and, consequently, its sample trajectories are not differentiable either. In this case, Equation (3) becomes X(t) + 0.1Ẋ(t) + X 2 (t)Ẋ(t) + X(t) = ξ(t). (35) As in the previous example, we are going to obtain approximations to the five first moments, E{ X i (t)}, i = 1, . . . , 5, the correlation function and the variance, V{ X(t)}, of the approximate solution X(t) of Equation (35). To implement the corresponding formulas derived throughout Section 3.2 saving computational time in Mathematica, we have taken into account the following properties of Dirac delta function, As mentioned in Example 1, the moments of odd order are null and using the positiveness of even order moments we can obtain some bounds for the perturbation parameter . First, using expression (17), the second-order moment is determined by so we obtain the bound < 4. Since E{ X(t)} = 0, expression (36) is also the variance of the first-order approximation. Second, using expression (23)- (25), This provides a stronger bound, < 1.59289. Now, applying (27), we obtain the following approximation of the correlation function, (38) In Figure 3, we show the plot of the correlation function, Γ X X (τ), given in the expression (38) for different values of satisfying the weaker and the stronger bounds previously determined. We can observe that for smaller values of the obtained approximation of the correlation function better preserves the symmetry as expected.  Applying the results presented in Section 4, we obtain the approximation of the PDF, f X(t) (x), for = 0.5, which satisfies the stronger bound 1.59289. We first compute the approximation based on the three first moments and, second, the approximation based on the five first moments f X(t) (x) = e −1+1.940−5.226·10 −11 x−17.837x 2 +1.580·10 −9 x 3 −42.679x 4 −7.904·10 −9 x 5 .
In Figure 4, we compare both graphical representations. We can observe, again, the similarity between them, thus showing full agreement in our numerical computations.
Finally, to check that our approximations are accurate, we compare the mean and standard deviation of X(t) obtained via the perturbation method against the ones computed by Euler-Maruyama numerical scheme [35]. The results are shown in Table 2. We can observe that both approximations are similar.

Conclusions
We have studied, from a probabilistic standpoint, a family of oscillators subject to small perturbations on the nonlinear term that depends both upon the position and the velocity (cross-nonlinearity) and whose forcing source is driven by a mean square differentiable stationary zero-mean Gaussian process. Despite the hypothesis of differentiability for the stochastic excitation, we have checked, via a numerical example, the method also provides good results when this hypothesis is not fulfilled, but involved computations are performed using the concept of general differentiability in the sense of distributions. We must point out that the majority of contributions dealing with this type of stochastic oscillators focus on the computation of the mean, the variance and correlation function. Our main contribution is the computation of reliable approximations of the probability density function of the stationary solution, by combining the stochastic perturbation method and the principle of maximum entropy. In this manner, we provide a fuller probabilistic description of the solution since from the density one can determine any one-dimensional moment as well as further probabilistic information of the steady-state. The proposed approach can be very useful to open new avenues in the analysis to other kind of nonlinear oscillators subject to small fluctuations and whose forcing term is a stochastic process that satisfies certain hypotheses. In our future research, we will work to continue contributing in this direction.