Probabilistic analysis of random nonlinear oscillators subject to small perturbations via probability density functions: theory and computing

We study a class of single-degree-of-freedom oscillators whose restoring function is affected by small nonlinearities and excited by stationary Gaussian stochastic processes. We obtain, via the stochastic perturbation technique, approximations of the main statistics of the steady state, which is a random variable, including the first moments, and the correlation and power spectral functions. Additionally, we combine this key information with the principle of maximum entropy to construct approximations of the probability density function of the steady state. We include two numerical examples where the advantages and limitations of the stochastic perturbation method are discussed with regard to certain general properties that must be preserved.


Introduction
Many problems arising in physics and engineering lead to differential equations whose data (initial/boundary conditions, forcing term and/or coefficients) must be set from experimental information that involves, what is usually termed, epistemic (or reducible) randomness [1]. Although this type of uncertainty can possibly be reduced by improved measurements or improvements in the modelling process, there is another uncertainty source often met in mathematical modelling of real-world problems called aleatory (or irreducible) randomness, which comes from the intrinsic variability of the phenomenon to be modelled. This approach leads to formulate random/stochastic differential equations [2,3]. Apart from answering fundamental questions about existence, uniqueness, continuous dependence of the solution with respect to model parameter or stability, solving a random differential equation means not only to calculate, exact or approximately, its solution, which is a stochastic process, but also to determine its main statistical information like the expectation or the variance. However, a more ambitious goal is to calculate the finite distribution functions (usually termed the fidis) of the solution, being the first probability density function the main fidis since by integration a e-mail: jccortes@imm.upv.es (corresponding author) b e-mail: ellona1@upv.es c e-mail: jvromero@imm.upv.es d e-mail: drosello@imm.upv.es one can calculate any one-dimensional moment (so including the mean and the variance), and also the probability that the solution lies within an interval of specific interest [2]. In real-world applications of physics or engineering, this fact is a key point since it allows us to calculate relevant information such as, for example, the probability that the position of an oscillator lies within two specific values where it is under control, the probability of buckling of a tower subject to extreme loads, etc.
In the realm of stochastic vibratory systems, many problems can be directly formulated through differential equations of the form L[Y (t)] = Z (t), where L[·] is a linear operator and Z (t) is a stochastic external source, which acts upon the system producing random vibrations. For example, the model L[ , has been used, from the original contribution [4], to describe the effect on earthbound structures of earthquake-type disturbances being Y (t) the relative horizontal displacement of, for example, the roof of a building with respect to the ground [5]. Additionally, many vibratory systems are described by nonlinear equations, say N [Y (t)] = Z (t), where the nonlinear operator is defined in terms of a small perturbation , N [·; ]. For example, in the case of the foregoing model the following general equation where h, that is independent of , is a nonlinear function of the unknown, and Y (t) has also been extensively applied. In most contributions, the nonlinear term h has a polynomial form [6]. For example, for h(Y (t)) = Y 3 (t), model (1) corresponds to Duffing oscillator, which physically models an elastic pendulum whose spring's stiffness violates Hooke's law [7]. The goal of this paper is to tackle the stochastic study of oscillators of form (1) in the case that the nonlinear term h is a transcendental function using a polynomial approximation, based on Taylor's expansions, and then to apply the stochastic perturbation method to approximate the main statistical functions of the steady state. Afterwards, we take advantage of the principle of maximum entropy (PME), in order to determine approximations of the probability density function (p.d.f.) of the stationary solution. To conduct our study, we have chosen a general form of the pendulum equation where we will assume that β > 0, the external source, Z (t), is defined via zero-mean Gaussian stationary stochastic process, which corresponds to an important case in the analysis of vibratory systems [8,9]. In Sect. 2, we will obtain the main theoretical results throughout. Afterwards, in Sect. 3, we perform a numerical analysis through two examples, with criticism about the validity of the results obtained via the stochastic perturbation method. Conclusions are drawn in Sect. 4.

Probabilistic analysis
This section is devoted to study, from a probabilistic standpoint, the steady state of model (2). It is important to point out that the non-perturbated associated model has a steady-state solution provided β > 0 and regardless of the value of ω 2 0 . For the sake of completeness, in "Appendix A", we discuss about this issue in the general setting that Z (t) is a stationary Gaussian process and also in connection with the examples presented later. Firstly, in Sect. 2.1, we will apply the stochastic perturbation method to approximate its stationary solution, which is a random variable. Secondly, in Sect. 2.2, we will perform a stochastic analysis addressed to obtain the main statistical functions of the stationary solution. In particular, we will obtain the first one-dimensional moments and the correlation function.

Stochastic perturbation expansion
In the case of model (2), the method of perturbation consists in expanding the unknown, Y (t), in terms of a formal power series of the perturbative parameter , which is assumed to have a small value (| | 1), where the coefficients Y n (t) need to be determined. The idea of the perturbation method is to impose that expansion (4) satisfies the nonlinear differential equation (2), which does not have a closed-form solution, in order to obtain a set of exactly solvable equations whose unknowns are the coefficients Y n (t). As this reasoning leads to an infinite cascade system for Y n (t), in practice only a few terms of the expansion are considered, and so the method provides an approximation of the solution. This technique is usually applied by truncating expansion (4) to the first-order approximation since as is small, the terms associated with higher powers (Y n (t) n , n = 2, 3, . . .) in the series expansion can usually be neglected. This fact determines the legitimacy of the approximation only for a small range for the values of parameter , which is consistent with the initial assumption that is a small parameter [10]. It is important to point out that, as already reported in [2,Ch. 7], no proof is available to show that the stochastic process Y (t) given in (4) converges in the mean square sense (or other probabilistic convergences), while some restrictive convergence conditions have been established in the deterministic framework [11]. Therefore, the convergence of Y (t) can be formulated in terms of strong hypotheses about its sample behaviour, which, in practice, results in very restrictive assumptions. Based on these facts, we here will apply the stochastic perturbation method by checking its validity with regard that some important statistical properties are preserved. On the other hand, to derive a solvable family of differential equations after applying the perturbation method to model (2), we will also use a double approximation for the nonlinear term sin(Y (t)). Specifically, we first apply a truncation of its Taylor's series and secondly, we approximate Y (t) using (5), i.e.
Eur. Phys. J. Plus (2021) 136:723 Substituting expansions (7) and (5) into (2), and equating terms with the same power of , leads to the two following linear differential equations that can be solved in cascade. Although the method can be applied for any order of truncation associated with the Taylor's expansion (6), in practice the value of M is set when the approximations obtained with M and M + 1 are very closed with reference to a prefixed error or tolerance. In our subsequent analysis, we will consider M = 2 (that corresponds to a Taylor's approximation of order 5, sin( 5 ), since, as we shall show later, the approximations of the main statistics of the solution stochastic process do not significantly change with respect to M = 1 (that corresponds to a Taylor's approximation of order 3, sin( 3 ). Summarizing, we have performed a double truncation. The first one, based on Taylor expansion for the nonlinear term (sin(Y (t))), and the second one, when applying the stochastic perturbation expansion given in expression (5).
(This latter approximation is fixed in our analysis at order 1 in , as it is usually done in the literature.) Notice that the order M of the Taylor truncation is independent of the order of truncation applied for the perturbation method.
As previously indicated, we are interested in the stochastic analysis of the stationary solution or steady state. Based upon the linear theory of Laplace transformation [12], the solutions Y 0 (t) and Y 1 (t) are given by and where for the underdamped case ( β 2 Physically this situation corresponds to the case that the oscillator approaches zero oscillating about this value [13]. Finally, we point out that the integrals defining Y 0 (t) and Y 1 (t) in (9) and (10) must be interpreted in the mean square sense [2, Ch. 4].

Constructing approximations of the first moments of the stationary solution
It is important to point out that the solution of model (2) is not Gaussian, so in order to probabilistically describe the stationary solution, the approximations of the mean and the variance (or more generally, the correlation) functions are not enough. This fact motivates that this subsection is addressed to construct reliable approximations of higher moments of the stationary solution [represented by the first-order approximation (5)]. This key information will be used in the next section to obtain approximations of the p.d.f., which in turns permits to obtain any one-dimensional moment of the stationary solution. Specifically, in the subsequent development we will compute any statistical moments of odd order, E ( Y (t)) 2i+1 , i = 0, 1, 2, . . ., the second-order moment, E ( Y (t)) 2 , the correlation function, Γ Y Y (τ ), the variance, V Y (t) , and the spectral density function, S Y(t) ( f ).
As indicated in Sect. 1, we shall assume that the stochastic external source, Z (t), is a stationary Gaussian stochastic process centred at the origin, i.e. E [Z (t)] = 0, being Γ Z Z (τ ), its correlation function. Notice that the hypothesis E [Z (t)] = 0 is not restrictive since otherwise we can work, without loss of generality, with the processZ , whose mean is null.
The mean of the first-order approximation is calculated taking the expectation operator in (5) and using its linearity, To compute E [Y 0 (t)], we take the expectation operator in (9), then we first apply the commutation of the expectation and the mean square integral [2, Eq. (4.165), p. 104] and, secondly, To compute E [Y 1 (t)], we take the expectation operator in (10) (recall that we take M = 2) and we again apply the commutation between the expectation and the mean square integral as well as the integral representation of Y 0 (t) given in (9), ds 5 ds 4 ds 3 ds 2 ds 1 ds .
Now, observe that  (14) one gets Substituting (15) and (13) into (12), one obtains To obtain the second-order moment, E ( Y (t)) 2 , we square expression (5), but retaining up to the first-order approximation in the parameter , To calculate E (Y 0 (t)) 2 , we substitute expression (9) and apply Fubini's theorem, Observe that the correlation function is a deterministic function of a single variable, since (17), we substitute the expressions of Y 0 (t) and Y 1 (t) given in (9) and (10) (recall that we take M = 2) and apply Fubini's theorem, ] ds 6 ds 5 ds 4 ds 3 ds 2 ds 1 ds . (20) Now, we express the expectations that appear in the above integrals in terms of the correlation function, Γ Z Z (·), taking into account that Z (t) is stationary. For the first expectation, one gets To calculate the other two expectations, we will apply the symmetry in the subindexes and the Isserlis-Wick theorem [14], which provides a formula to express the higher-order moments of a zero-mean multivariate normal vector, say (Y 1 , . . . , Y n ), in terms of its correlation matrix The sum is over all distinct ways of partitioning the set of indexes {1, 2, . . . , n} into pairs {i, j} (the set of these pairs is denoted by P 2 n ), and the product is over these pairs. In our case, we apply this result taking into account that Z (t) is a zero-mean stationary Gaussian process. To determine the second expectation in (20), we will denote Then, the second integral in (20) can be computed as Notice that we have taken advantage of the symmetry of the correlation function, Γ Z Z (·), to express the last multidimensional integral. The third expectation in (20) can be analogously calculated Then, substituting (22) and (23) into (20), and taking again the advantage of the symmetry of the correlation function, Γ Z Z (·), to simplify the representation of the last multidimensional integral, one obtains So, substituting expressions (19) and (24) into (17), we obtain an explicit approximation of the second-order moment for the approximation Y (t), Now, we compute the third-order moment of Y (t), using the first-order approximation with respect to the perturbative parameter , By reasoning analogously as in the calculation of the foregoing statistical moments, we obtain and We here omit the details of these calculations since they are somewhat cumbersome and they can easily inferred from our previous developments. Then, substituting (27) and (28) into (26), we obtain In general, it can be straightforwardly shown that the statistical moments of odd order are null, The correlation function of Y (t), using the approximation of first-order with respect to the perturbative parameter , is given by The first term of (30) corresponds to the correlation function of Y 0 (t). It is determined, in terms of the correlation function Γ Z Z (·), by applying Fubini's theorem and the stationarity of Z (t), The two last expectations on the right-hand side of (30), correspond to the cross-correlation function of Y 0 (t) and Y 1 (t). They can be expressed explicitly in terms of the correlation function Γ Z Z (·),  and 6 ) ds 6 ds 5 ds 4 ds 3 ds 2 ds 1 ds .
From expressions (31)-(33), we observe that the correlation function, Γ Y Y (τ ), given by (30), only depends on the difference τ between two different instants, t and t + τ (as it has been anticipated in the notation in (30)). This fact together with E Y (t) = 0 (see (16)) allows us to say that the first-order approximation, Y (t), given in (5) and obtained via the perturbation technique, is a stationary stochastic process. Additionally, it is clear that the covariance function of the steady state coincides with the correlation function, i.e. Cov Y (t 1 ), Y (t 2 ) = Γ Y Y (τ ), where τ = |t 1 − t 2 | and, also the variance matches the second-order moment, which in turn can be calculated evaluating the correlation function at the origin, From the correlation function, Γ Y Y (τ ), it is straightforward to also determine the power spectral density function (or simply, the power spectrum) of Y (t) here, i = √ −1 is the imaginary unit and f is the angular frequency. Observe that, S Y(t) ( f ) is the Fourier transform of the correlation function of Y (t). This function plays a key role to describe the distribution of power into frequency components composing a signal represented via stationary stochastic process [15].

Numerical examples
In the previous section, we have obtained approximations of the main statistical moments of the steady state of problem (2). In this section, we combine this key information, together with the PME technique, to construct approximations of the p.d.f. of the stationary solution. We will show two examples where important stochastic processes play the role of the external source, Z (t), in problem (2). In the first example, we will take as Z (t) the white Gaussian noise, while in the second one, the Ornstein-Uhlenbeck stochastic process will be considered. Since the validity of perturbation method is restricted to small values of the perturbative parameter , the numerical experiments are carried out with criticism to this key point taking into account that the numerical approximations must retain certain universal properties such as the positivity of even-order statistical moments (for any stochastic process) and the symmetry of the correlation and power spectral functions; the correlation function reaches its maximum value at the origin and the positivity of the spectral function, in the case of stationary stochastic processes. Of course, additional properties about other statistical functions could also be added to further check the consistency of the numerical results obtained via the stochastic perturbation method.
For the sake of completeness, down below we revise the main definitions and results about the PME method and the power spectral function that will be used throughout the two numerical examples.
As detailed in [16], the PME is an efficient method to determine the p.d.f., f Y (y), of a random variable, say Y , constrained by the statistical information available (domain, expectation, variance, symmetry, kurtosis, etc.) on Y . The method consists in determining f Y (y) such that it maximizes the so-called Shannon's entropy (also referred to as differential entropy), which is given by y)) dy, subject to a number of constraints, which are usually defined by the N first moments, say a n , n = 1, 2, . . . , N , of the random variable Y together with the normalization condition, y 1 y n f Y (y) dy = a n , n = 0, 1, . . . , N }. Notice that n = 0 corresponds to the normalization condition (i.e. the integral of the p.d.f. is one) and that the rest of the restrictions are defined by the statistical moments, E Y n = a n , n = 1, 2, . . . , N . It can be seen, by applying the functional version of Lagrange multipliers associated with A, that In practice, the Lagrange multipliers, λ i , i = 0, 1, . . . , N , are determined by numerically solving the following system of N + 1 nonlinear equations defined by the constraints y 2 y 1 y n 1 [y 1 ,y 2 ] (y) e − N i=0 λ i y i dy = a n , n = 0, 1, . . . , N .
In words, the PME method selects, among all the p.d.f.'s that satisfy the constraints given by the available statistical information, the one that maximizes the Shannon or differential entropy as a measure of randomness. This can be naively interpreted as looking for the p.d.f. which maximizes the uncertainty from the minimal quantity of information [16].
In the examples, we will also calculate the power spectral density function S Y(t) ( f ) defined in (35). Using the Euler identity e ix = cos(x) + i sin(x), it is easy to check that, for any stationary process, the power spectral density is an even function, This property is also fulfilled by the correlation function, Moreover, the correlation function reaches its maximum at the origin, i.e. |Γ Y Y (τ )| ≤ Γ Y Y (0) [2, Ch. 3], while it can be proved that the power spectral function is non-negative, S Y(t) ( f ) ≥ 0 [15]. To reject the possible spurious approximations obtained via the stochastic perturbation method, we will check in our numerical experiments whether all these properties are preserved. We will also take advantage of the approximations of the power spectral density and of the correlation function to obtain the two following important parameters associated with a stationary stochastic process, the noise intensity (D) and the correlation time (τ ), respectively, defined by Notice that the value of D comes from evaluating S Y(t) ( f ) at f = 0 in (38). So, the noise intensity is defined as the area beneath the correlation function Γ Y Y (τ ), while the correlation time is a standardized noise intensity [17]. The parameters D andτ indicate how strongly the stochastic process is correlated over the time.
In both examples, the numerical results that we shall show correspond to third-order Taylor's approximations of the nonlinear term sin(Y (t)) in (7), since we have checked that no significant differences are obtained using fifth-order Taylor's approximations.
Example 1 Let us consider as external source the stochastic process Z (t) = ξ(t), where ξ(t) is a white Gaussian noise, i.e. is a stationary Gaussian stochastic process with zero mean, E [Z (t)] = 0, and flat power spectral density, S Z (t) ( f ) = N 0 2 , for all f . Then, its correlation function is given by Γ Z Z (τ ) = N 0 2 δ(τ ), where δ(τ ) is the Dirac delta function. We will take the following data for the parameters involved in model (2), ω 0 = 1, β = 5 100 and N 0 = 1 100 (that satisfy the conditions of "Appendix A", so ensuring the existence of the steady-state solution), so the nolinear random oscillator is formulated by We will now take advantage of the results derived in Sect. 2.2, to calculate the following statistical information of the first-order approximation, Y (t), obtained via the perturbation method and given in (5): (1) the moments up to order 3, i.e. E ( Y (t)) i , i = 1, 2, 3; (2) the variance, V Y (t) , and (3) the correlation function, Γ Y Y (τ ). We will use this information to compute approximations, first of the p.d.f. of Y (t), using the PME, and, secondly, of the spectral density function of Y (t).
From expression (29), in particular, we know that E Y (t) and E ( Y (t)) 3 are null. For the second-order moment, E ( Y (t)) 2 , using expression (25) we obtain This value is also the variance since E Y (t) = 0. On the other hand, as E ( Y (t)) 2 is always positive, we can obtain the following bound for the perturbative parameter, < 1.01258. Although this value bounds the validity of the perturbation method, down below we show that is a conservative bound.
To this end, we will compare the mean and standard deviation (sd) of Y (t) obtained via the perturbation method and the ones computed by Kloeden-Platen-Schurz algorithm [3]. The results are shown in Table 1. We can observe that both approximations are accurate for = 0 (which corresponds to the linearization of model (2)), = 0.01 and = 0.1. But significant differences in the standard deviation are revealed for = 0.5 and = 1.
According to expressions (30)-(33), the approximation of the correlation function is given by Notice that, in full agreement with expression (34), it is satisfied that In Fig. 1, we show the graphical representation of the correlation function, Γ Y Y (τ ), plotted from the expression (41) and using different values of the perturbative parameter . To emphasize this dependence on the parameter , hereinafter we will denote this function by Γ Y Y (τ, ). In Fig. 1a, we observe that the approximations of the Γ Y Y (τ, ) deteriorate as increases in full agreement with the results obtained in Table 1. The deterioration of the approximations as increases can also be confirmed by checking that the general property |Γ Y Y (τ )| ≤ Γ Y Y (0) for the correlation function [18] does not fulfil for = 0.5 and = 1 (see Fig. 1a). In contrast, for smaller values of = 0, 0.01, 0.1} this property holds. Notice that the correlation function for = 0 and = 0.01 is quite similar, as expected. For the sake of clarity, we have plotted these results in Fig. 1b. Finally, in Fig. 1c Now, we compute the approximation of the p.d.f., f Y (t) (y), of the steady state, using the PME. Therefore, according to the PME the p.d.f. is sought in the form where in our case, λ 0 , λ 1 , λ 2 and λ 3 are determined numerically solving system (37) with a 0 = 1, a 1 = 0, a 2 = 1 40 − 12641 512000 . In Table 2, we show the values of λ 0 , λ 1 , λ 2 and λ 3 and the corresponding domain [y 1 , y 2 ] for the following values of the perturbative parameter ∈ {0, 0.01, 0.1}. The domain has been determined using the Bienaymé-Chebyshev inequality [μ−kσ, μ+kσ ] (in our case μ = 0) with k = 10. This guarantees the 99% of the probability is contained in the above interval [μ − kσ, μ + kσ ] regardless of the distribution of the corresponding random variable [14]. In Fig. 2, we compare the graphical representations of the p.d.f, f Y (t) (y). From them, we can observe that the plots are quite similar.
To complete our numerical analysis, in Fig. 3 we show a graphical representation of the power spectral density for ∈ {0, 0.01, 0.1}. We observe that the approximation obtained via the stochastic perturbation method is able to retain the properties of symmetry and positivity of the power spectral density for ∈ {0, 0.01}; however, positivity begins to slightly fail for = 0.1, therefore restricting the validity of the results provided by the stochastic perturbation method. In Table 3, the noise intensity (D) and the correlation time (τ ) have been calculated for ∈ {0, 0.01}.

Example 2
In this second example, let us consider the Ornstein-Uhlenbeck stochastic process to play the role of the external source, Z (t). So, it is defined as the stationary solution of the Langevin equation where W (t) is the Wiener process [3]. Notice that α > 0 is a necessary and sufficient condition to have a stationary solution. Z (t) satisfies the hypotheses so that the stochastic perturbation method can be applied, i.e. is a zero-mean stationary Gaussian stochastic process, being Γ Z Z (τ ) = σ 2 e −α|τ | its correlation function. We take the following values for the parameters in Eq. (2), ω 0 = 1, β = 1/100, σ = 1/100 and α = 1/2 (thus, according to "Appendix A", it ensures the existence of the steady-state solution). So, in this case, the nonlinear random oscillator is given by We will now apply the same steps as in Example 1 to obtain approximations of the main statistical functions of the approximate stochastic solution Y (t). First, we will determine the three first statistical moments. From expression (29), E Y (t) and E Y 3 (t) are null. For the second-order moment, E Y 2 (t) , using expression (17) we obtain From this expression, we can obtain a rough bound for , since E Y 2 (t) > 0. In this case, < 0.558098. To check that our second-order moment approximation is consistent, we will compare it with the random linear oscillator obtained when → 0 [2, Example 7.2], Notice that, according to (34), expression (43) is also the variance of Y (t), V Y (t) . In Fig. 4, we can observe that for t large enough (corresponding in the limit as t → ∞ to the steady state), the second-order moment of the random linear equation approaches to our approximation for = 0. Observe in the plot that the E Y 2 (t) → 0.00206349 ≈ 13 6300 as t → ∞, in accordance with (43).
Once we have obtained the mean, E Y (t) and the standard deviation of Y (t) via the perturbation method, we compare them with the ones computed by the Kloeden-Platen-Schurz algorithm. The results are shown in Table 4. We can observe that both approximations are accurate for ∈ {0, 0.01, 0.1}; however, significant differences in the standard deviation are shown for = 0.5, thus showing the perturbation method does not provide acceptable approximations. The approximation of the correlation function, Γ Y Y (τ ), using (30) together with expressions (31)-(33), is given by Then, we can check that property (34) holds. In Now, we will obtain the approximation of the p.d.f. in the form f Y (t) (y) = e −1−λ 0 −λ 1 y−λ 2 y 2 −λ 3 y 3 using the PME. The values of λ 0 , λ 1 , λ 2 and λ 3 are shown in Table 5 together with the domain [y 1 , y 2 ] for the following values of the perturbative parameter ∈ {0, 0.01, 0.1}. As in Example 1, the intervals [y 1 , y 2 ] have been computed using the Bienaymé-Chebyshev inequality.
In Fig. 6, we compare the graphical representations of the p.d.f., f Y (t) (y). From them, we can observe that the approximations are very similar.
To complete our numerical example, we have calculated graphical representations of the power spectral density, S Y (t) ( f ), of the Y (t). In Fig. 7, we show two plots. Panel left corresponds to ∈ {0, 0.01}, and panel right to ∈ {0, 0.001}. We observe that the property of symmetry breaks down when increases, while positivity is retained. In Table 6, D andτ have been calculated for ∈ {0, 0.001}.

Conclusions
In this paper, we have studied a class of stochastic nonlinear oscillators, whose nonlinear term is defined via a transcendental function. We have assumed that the oscillator is excited by a zero-mean stationary Gaussian process. Since the nonlinear term is affected by a small parameter, to conduct our probabilistic analysis, we have approximated the nonlinear term using a Taylor's polynomial, and then, we have applied the stochastic perturbation method to obtain the main statistical moments of the stationary solution. After this theoretical analysis, we have carried out numerical examples where the stochastic excitation is driven by two important stochastic processes, the Gaussian white noise and the Ornstein-Uhlenbeck processes. Since a key point when applying the perturbation method is the accuracy of the approximations in terms of the size of the perturbative parameter, from the numerical results  rectly preserved. To better checking the accuracy of the approximations of the mean and the standard deviation via the perturbation method, we have compared them with the ones calculated by means of an accurate numerical scheme, showing good agreement for certain sizes of the perturbative parameter. This comparative analysis includes the linear case obtained when the perturbative parameter is null. In this limit case, the results are also fully consistent. Summarizing, our study shows, by means of a class of stochastic nonlinear oscillators, that the double approximation Taylor-perturbation method is able to approximate the statistics of the stationary solution. Additionally, we have taken advantage of the above computed statistics, in combination with the principle of maximum entropy, to construct reliable approximations of the density of the steady state. Our analysis has been performed with criticism with regard to the size of the perturbative parameter, as required when applying the stochastic perturbation method. Our approach can be useful to study other stochastic nonlinear oscillators whose small perturbations affect a transcendental functions with the additional advantage of computing the density of the stationary solution. In our future works, we plan to tackle this type of problems for other stochastic nonlinear oscillators that have not been analysed yet.
Acknowledgements This work has been supported by the Spanish Ministerio de Economía, Industria y Competitividad (MINECO), the Agencia Estatal de Investigación (AEI), and Fondo Europeo de Desarrollo Regional (FEDER UE) Grant PID2020-115270GB-I00. The authors express their deepest thanks and respect to the reviewers for their valuable comments.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Conflict of interest
The authors declare that they do not have any conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A Appendix: Existence of the steady-state solution
Notice that the linearized equation associated with (2) is and using the classical change of variable Y = (Y 1 (t), Y 2 (t)) = (Y (t), dY (t) dt ) can be written as i.e.
In the case that Z (t) = ξ(t) is the white noise (a zero mean, Gaussian and stationary process), i.e. ξ(t) =Ẇ (t) (W (t) is the Wiener process), equation (48) writes it is well known that it has a stationary or steady-state solution if the real parts of the eigenvalues of matrix F = 0 1 −ω 2 0 −2β are negative (in other words, F is a Hurwitz matrix). Assuming β > 0, it is easy to check that this condition fulfils since the spectrum os F is given by σ (F) = λ 1 = −β + β 2 − ω 2 0 , λ 2 = −β − β 2 − ω 2 0 , and Re(λ 1 ) = Re(λ 2 ) = −β < 0 (Re(·) denotes the real part). Observe that in the underdamped case (β 2 /ω 2 0 < 1), λ 1 and λ 2 are complex conjugate. This fact is used in Example 1 where β = 1 20 > 0 to guarantee the existence of the steady-state solution. In the more general case that Z (t) in (47) is such that satisfies a state-space stochastic differential equation (SDE) of the form model (47) together with (49) can be written as This system is of the form The well-known results about the existence of a steady-state solution of Itô SDEs then apply to study the larger dimensional system (50) [or equivalently (51) and (52)]. As a consequence, it is enough to check that the real parts of all the eigenvalues of matrixF are negative. In this case, σ (F) = λ 1 = r 1 , λ 2 = −β + β 2 − ω 2 0 , λ 3 = −β − β 2 − ω 2 0 .
Hence, since β > 0, it is sufficient that r 1 < 0. Notice that in Example 2, Z (t) is the Ornstein-Uhlenbeck process and r 1 = −α < 0 (since α > 0), so the existence of the steady state is also guaranteed. Finally, it is interesting to point out that in the general case that Z (t) is a (zero mean) stationary Gaussian process (like the white noise and Ornstein-Uhlenbeck processes in Examples 1 and 2, respectively), it is possible to use state-space representations of form (51) to reduce a (stationary) Gaussian process-driven ordinary differential equation of the form dY(t) dt = FY(t) + GZ (t), as (47) into a larger-dimensional ordinary differential equation driven by white noise, i.e. a linear Itô SDE. This can be done when the spectral density of the covariance function of Z (t) is a rational function.