## Abstract

Complex fractional moment (CFM), which is defined as the Mellin transform of a probability density function (PDF), has been successfully employed to find the response PDF of a wide variety of integer-order nonlinear oscillators. In this paper, a CFM-based analysis is performed to determine the transient response PDF of nonlinear oscillators with fractional derivative elements under Gaussian white noise. First, an equivalent linear system is introduced for the purpose of deriving the Fokker–Planck (FP) equation for response amplitude. The equivalent natural frequency and equivalent damping coefficient of the system need to be determined, taking into account both the nonlinear and fractional derivative elements of the original oscillator. Moreover, to convert the FP equation into the governing equation of CFMs, these equivalent coefficients must be given in polynomial form of amplitude. This paper proposes formulas for appropriately determining the equivalent coefficients, based on an equivalent linearization technique. Then, applying stochastic averaging, the FP equation is derived from the equivalent linear system. Next, the Mellin transform converts the FP equation into coupled linear ordinary differential equations for amplitude CFMs, which are solved with a constraint corresponding to the normalization condition for a PDF. Finally, the inverse Mellin transform of the CFMs yields the amplitude PDF. The joint PDF of displacement and velocity is also obtained from the amplitude PDF. Three linear and nonlinear fractional oscillators are considered in numerical examples. For all cases, the analytical results are in good agreement with the pertinent Monte Carlo simulation results.

## 1 Introduction

Mechanical and structural systems in various engineering applications are exposed to random excitations, such as seismic ground motion, fluctuating pressure in turbulent flow, and hydraulic force of ocean waves. In order to grasp the response characteristics of dynamic systems subjected to random excitation and to assess their reliability, it is important to accurately determine the probability density function (PDF) of the system response. For a nonlinear system under Gaussian white noise, the system response can usually be regarded as a Markov process, and the time evolution of the response PDF is governed by a partial differential equation called Fokker–Planck (FP) equation. Unfortunately, exact solutions of the FP equation are not available except for the stationary response of a very limited class of nonlinear systems [1]. Solutions for transient/non-stationary responses are even more difficult to obtain, owing to their time dependence. Therefore, the development of approximate analytical or numerical techniques for finding system response PDFs is a significant issue in the field of stochastic dynamics.

In the last few decades, fractional calculus has attracted a great deal of attention due to its usefulness as a mathematical modeling tool for non-local or long-memory phenomena. Various topics on fractional calculus in science and engineering can be found in the books [2–7] and the papers [8–10]. One successful application of fractional calculus to engineering is the modeling of viscoelastic materials. Indeed, many experimental results have proven that fractional derivative models with a small number of parameters can describe viscoelastic behavior with a high degree of accuracy [11–18]. Therefore, in mechanical and structural dynamics, the fractional models have become widely used as models of viscoelastic dampers for seismic isolation and vibration control [16,19–23].

Under such circumstances, there has been a growing interest in the response determination of dynamic systems with fractional derivative elements under stochastic excitation. In this regard, Monte Carlo simulation (MCS) comes to mind as a general-purpose method. However, computing the fractional derivative operator requires considerable computational time and resources, since this operator is defined by a convolution integral, i.e., depends on the function’s entire past time history. Hence, direct numerical integration of system equations with fractional-order derivatives is much more computationally demanding than that of systems with only integer-order derivatives. For this reason, in recent years, a lot of efforts have been put into developing analytical techniques for determining the stochastic response of fractional-order systems. Meanwhile, machines and structures including viscoelastic materials have often been modeled as single-degree-of-freedom (SDOF) systems with fractional derivative elements [2,4,8,9,12–14,16]. Therefore, the development of analytical methods for SDOF systems with fractional derivatives is interesting and important from a practical point of view.

For moment analysis, an approach based on convolutional representation of the response was presented to obtain the response variance of SDOF linear systems with fractional derivatives [24,25]. More recently, approximate closed-form solutions of response moments were obtained for fractional linear oscillators to evolutionary stochastic excitation by expressing their impulse response with Prony series [26]. A moment equation approach has also been employed [27–29], where fractional-order systems are approximately converted into coupled integer-order differential equations. In addition, frequency-domain approaches have been developed for stationary moment analysis [30–32].

For the analysis of PDFs, stochastic averaging has been widely utilized. The displacement and velocity responses of fractional systems are no longer a Markov vector due to the history dependence of the fractional terms, even if the excitation is a white noise process. Nevertheless, by applying stochastic averaging, the response amplitude can be approximately treated as a Markov process, resulting in the FP equation governing the amplitude PDF [33]. As an additional benefit, system dimension reduction is also achieved. In this manner, the amplitude PDF of the stationary response of various fractional nonlinear oscillators has been accurately obtained [33–41].

Compared to moment and stationary PDF analyses, the analysis of transient response PDF is a more difficult task. Although the number of studies is still limited at present, efforts are being made to seek an approximate solution of the transient amplitude PDF of fractional nonlinear systems using a path integral technique [42] or a Rayleigh distribution assumption [43]. However, these existing methods have some drawbacks in terms of computational cost and applicability. The path integration generally requires very short time-steps to obtain a highly accurate solution. The Rayleigh distribution assumption for amplitude PDF may not be appropriate for some types of nonlinear systems, such as oscillators possessing limit cycles. Therefore, it remains an important challenge to develop efficient and versatile analytical methods for determining the transient response PDF of nonlinear oscillators with fractional derivative elements.

Meanwhile, a novel statistical moment named complex fractional moment (CFM) was proposed for the characterization of random variables [44,45]. The CFM is a complex-order moment, which is calculated as the Mellin transform of a PDF. Similar to the Fourier transform relation between PDFs and characteristic functions, the inverse Mellin transform of the CFM reconstructs the original PDF. This striking feature of the CFM allows us to estimate PDFs via CFMs. Di Paola presented a CFM-based approach for solving the FP equation corresponding to first-order nonlinear systems excited by Gaussian white noise [46], in which by applying the Mellin transform to the FP equation, it is converted into a set of linear ordinary differential equations with respect to response CFMs. In other words, a partial differential equation is converted into linear ordinary differential equations much easier to solve. Then, solving the CFM equations and taking the inverse Mellin transform on the resulting CFMs, the response PDF is constructed. Through the above method, the PDFs of both transient and stationary responses were accurately determined in the whole domain, including the tail region. This method has also been applied successfully for white noise excitations of other types [47–49]. Jin et al. developed the CFM-based approach for SDOF oscillators with nonlinear damping under Gaussian white noise by combining stochastic averaging [50]. Stochastic averaging leads to the FP equation for response amplitude, which is then solved via CFMs. The CFM method combined with stochastic averaging yielded transient and stationary response PDFs with high accuracy, and the computational time was significantly reduced compared to the corresponding MCS. In the subsequent studies, this method has been extended to a wide class of integer-order nonlinear oscillators, including systems with nonlinearities in both stiffness and damping [51], nonlinear vibro-impact systems [52] and nonlinear systems with time delay [53]. Such great success of the CFM-based method suggests that it can also be a powerful tool in the transient response analysis of fractional-order nonlinear oscillators.

In this paper, a CFM-based analysis is performed to determine the transient response PDF of nonlinear oscillators with fractional derivative elements, whose nonlinearities are included in both stiffness and damping. This analysis is mainly based on an analytical method for integer-order nonlinear oscillators recently developed by some of the authors [51]. First, we aim to derive the FP equation for response amplitude using stochastic averaging. For this purpose, the equivalent natural frequency and equivalent damping coefficient of the oscillator need to be determined. As an important point in the analysis of fractional oscillators, these equivalent coefficients have to include not only the effects of nonlinearities, but also the contributions of the fractional derivative terms. Furthermore, the coefficients must be given in polynomial form of response amplitude to utilize the CFM-based approach. This paper proposes formulas for determining the equivalent coefficients that satisfy the above requirements, based on an equivalent linearization technique. Then, resorting to stochastic averaging, the FP equation is derived. The solution of the FP equation, i.e., the amplitude PDF, is sought through the CFM approach. The joint PDF of displacement and velocity is also obtained from the amplitude PDF. The accuracy and versatility of the present analytical technique are demonstrated through numerical examples of three linear and nonlinear oscillators with fractional derivatives.

## 2 Preliminaries

### 2.1 Mellin Transform.

*Fundamental Strip*(FS) of the Mellin transform [45].

### 2.2 Complex Fractional Moment.

## 3 Equation of Motion

The purpose of this paper is to obtain the transient response PDF of the oscillator expressed by Eq. (7). Some of the authors recently developed a CFM-based analytical approach for determining the transient response PDF of an integer-order nonlinear oscillator obtained by removing the fractional derivative term from Eq. (7) [51]. The present analysis for the fractional nonlinear oscillator (7) is also performed based on this analytical approach. An important point that must be addressed to accomplish this analysis is the appropriate determination of the equivalent natural frequency and equivalent damping coefficient for the fractional nonlinear oscillator.

## 4 Derivation of Fokker–Planck Equation for Response Amplitude

### 4.1 Pseudo-Harmonic Representation of Response.

### 4.2 Equivalent Linearization.

To solve the FP equation derived in the following section via the CFM-based approach, $\omega e2(A)$ and $\beta e(A)$ must be given as polynomials in amplitude $A$ [51]. Unfortunately, this requirement is not met, due to the presence of $\omega e\alpha (A)$ and $\omega e\alpha \u22121(A)$ in Eqs. (16) and (17) ($0<\alpha <1$). In order to overcome this problem, $\omega e\alpha (A)$ in Eq. (16) and $\omega e\alpha \u22121(A)$ in Eq. (17) are approximated by $\omega 0\alpha $ and $\omega 0\alpha \u22121$, respectively. This analysis assumes $\epsilon \u226a1$. In Eq. (16), the term to which the approximation is made includes $\epsilon $. On the other hand, Eq. (17) does not include $\epsilon $, but the damping term of the equivalent linear system in Eq. (11) is multiplied by $\epsilon $. Therefore, the effect of the approximation is considered sufficiently small.

### 4.3 Stochastic Averaging.

Finally, the magnitude of the parameter $\epsilon $ of the oscillator in performing the above analysis is mentioned. In Refs. [61] and [62], the stationary response PDFs and transient response statistics of various integer-order nonlinear systems were analyzed by combining the equivalent linearization and the stochastic averaging. It was shown that these PDFs and statistics can be obtained correctly when $\epsilon $ is less than $0.1$. Thus, the range of $\epsilon $ for which the combination of the equivalent linearization and the stochastic averaging works well is expected to be $\epsilon <0.1$.

## 5 Complex Fractional Moment Approach

### 5.1 Conversion of Fokker–Planck Equation to Complex Fractional Moment Equations.

### 5.2 Solving Complex Fractional Moment Equations.

In this way, the solution of the FP equation (22) is attained in terms of CFMs by solving Eqs. (34)–(36). The computation time is only a few seconds. When performing a stationary response analysis, linear algebraic equations derived by setting the time derivative term in Eq. (34) to zero are solved in conjunction with Eq. (35).

### 5.3 Calculation of Response Probability Density Function.

The displacement PDF $pX(x,t)$ and the velocity PDF $pX\u02d9(x\u02d9,t)$ are determined as the marginal PDFs of $pXX\u02d9(x,x\u02d9,t)$.

### 5.4 Summary of Analytical Procedure.

In Secs. 4 and 5, the derivation process for the equations used in the present analysis was described in detail. However, when applying the present approach, it is not necessary to perform all of the analyses in Secs. 4 and 5 in sequence. The actual analysis consists of fewer equations and procedures. The practical analytical procedure is summarized as follows:

Calculate the equivalent natural frequency $\omega e(A)$ and the equivalent damping coefficient $\beta e(A)$ using Eqs. (18) and (19). $\omega e2(A)$ is given in the form of Eqs. (25) and (29), while $\beta e(A)$ is given in the form of Eq. (30).

Substitute $\omega e(A)$ and $\beta e(A)$ into Eq. (34) to get the linear ordinary differential equations for amplitude CFMs $MpA(\gamma s\u22121,t)(s=\u2212m,\u2026,\u22121,1,\u2026,m)$.

Solve Eq. (34) in combination with Eq. (35) and the initial condition given by Eq. (36) to determine $MpA(\gamma s\u22121,t)$$(s=\u2212m,\u2026,m)$.

Apply the inverse Mellin transform to $MpA(\gamma s\u22121,t)$, i.e., substitute $MpA(\gamma s\u22121,t)$ into Eq. (6) to obtain the amplitude PDF $pA(a,t)$.

Calculate $a$ and $|J|$ in the joint PDF $pXX\u02d9(x,x\u02d9,t)$ of $X(t)$ and $X\u02d9(t)$ in Eq. (37) via Eqs. (38) and (39).

Substitute $pA(a,t)$, $a,$ and $|J|$ into Eq. (37) to obtain $pXX\u02d9(x,x\u02d9,t)$.

## 6 Numerical Examples

The accuracy and versatility of the present analytical approach are assessed through numerical examples of three linear and nonlinear oscillators with fractional derivative elements. To the best of the author’s knowledge, no exact PDF solution has yet been found for oscillators with fractional derivatives, even for stationary responses. For this reason, the analytical results of response PDFs are compared with the results of the pertinent MCS.

### 6.1 Linear Oscillator With Fractional Derivative Element.

The reliability of the analytical method is validated with three different values of the fractional derivative order $\alpha $: $\alpha =0.3,0.5,$ and $0.7$. It is known that $0.3\u2264\alpha \u22640.7$ is frequently adopted as the best fitting value in modeling the dynamic properties of viscoelastic materials [11–15]. For this reason, the above three values are used for the fractional order. Other system and excitation parameters are set as $\omega 0=1$, $\epsilon =0.01$, $\lambda =7,$ and $S0=1$. Further, $\Delta \eta =0.3$, $\rho =3.7,$ and $m=220$ are adopted in solving Eq. (44).

For comparisons of the analytical results with simulation results, MCS is performed on Eq. (42). The difficulty inherent in simulations of fractional-order systems is that the fractional derivatives depend on the entire past time history of the system response, as can be seen from Eq. (8). As an efficient and accurate numerical technique, this study utilizes a method to approximately convert an equation of motion including fractional-order derivatives into a set of integer-order differential equations [63]. A brief summary of this methodology is provided in Appendix. Once the integer-order differential equations are obtained, standard numerical integration algorithms can be used to generate realizations of the system response. In this study, the fourth-order Runge–Kutta method is employed. The time increment size $\Delta t$ for the numerical integration is $\Delta t=10\u22124$. MCS results are computed from $5\xd7106$ realizations. $\Delta t$ and the number of realizations were chosen based on a preliminary study in order to obtain as accurate a solution as possible, since the MCS results are used as substitutes for the true response PDFs of the oscillator. The computer system used in the MCS has an Intel Core i9-10,900 K CPU with 32 GB RAM, and the operating system is Microsoft Windows 10. The source code for the MCS was prepared using the c language, which was compiled with the GNU compiler collection (version 12.3). The simulation method described herein is also adopted in the second and third examples shown later.

In Fig. 1, the PDFs $pA(a,t)$, $pX(x,t)$, and $pX\u02d9(x\u02d9,t)$ of amplitude $A(t)$, displacement $X(t)$, and velocity $X\u02d9(t)$ for the case of $\alpha =0.5$ are shown. The solid lines indicate the analytical results, whereas the discrete symbols denote the pertinent MCS results. The three times in the figure were chosen to show the early phase of the response, the stationary phase and a phase in between the two, while also taking into account the ease of viewing the figure. From the standpoint of evaluating system reliability, it is crucial to obtain the tails of response PDFs with high accuracy. Therefore, semi-log plots of $pX(x,t)$ and $pX\u02d9(x\u02d9,t)$ are also presented to confirm the agreement between the analytical and MCS results in the tail region. The results of $pA(a,t)$ demonstrate the validity of the assumptions about $pA(a,t)$ used in deriving Eq. (32) from Eq. (31). Slight deviations can be seen in the peaks of the displacement and velocity PDFs. However, the accuracy of the analytical results is overall quite satisfactory, including the tail region. Although not shown in the figure, good agreement was also achieved for the cases of $\alpha =0.3$ and $\alpha =0.7$ (some results are provided in Fig. 2). Additionally, the authors also carefully checked the validity for various other values of $\alpha $ in the range $0<\alpha <1$, and the analytical method accurately yielded transient response PDF solutions for any $\alpha $, including the cases of $\alpha =0.1$ and $0.9$.

Figure 2 shows a comparison of displacement and velocity PDFs between the three types of fractional order $\alpha $. The following tendencies can be observed in the response PDFs for larger $\alpha $:

The widths of the PDFs are narrower (i.e., their variance decreases).

They reach stationary PDFs faster.

Finally, the computation time for the present analytical method was $2.7$ s, whereas the MCS took about $25$ h (the computation times were similar in the second and third examples shown later as well). Thus, the computational time was significantly reduced compared to the corresponding MCS.

### 6.2 Oscillator With Fractional Derivative Element and Nonlinearities in Both Damping and Stiffness.

The following parameter values are used: $\omega 0=1$, $\epsilon =0.01$, $\lambda =7$, $\alpha =0.5$, $\delta 1=8$, $\delta 3=8$, $\kappa 3=1$, $\kappa 5=1$, and $S0=1$. The parameters when solving Eq. (48) are given by $\Delta \eta =0.5$, $\rho =3.7$, and $m=220$.

Figure 3 displays the analytical results of transient and stationary response PDFs and the corresponding MCS results. There are small errors in the peaks of the displacement and velocity PDFs for the initial short time period, but otherwise the analytical results are in good agreement with the simulation results, including the tail region. The results of the amplitude PDF $pA(a,t)$ show that the assumptions about $pA(a,t)$ are valid for the system containing nonlinearities as well. Further, a comparison of Figs. 1 and 3 reveals that the response PDFs in this example are narrower than those of the linear oscillator with same $\alpha $. This is because the nonlinear damping and nonlinear stiffness of the system behave to suppress large responses. The analytical result accurately captures such a narrowed PDF shape caused by the system nonlinearities. Thus, the present analytical method is able to provide accurate solutions even for a system with nonlinear damping and nonlinear stiffness as well as a fractional derivative element.

### 6.3 Duffing-van der Pol Oscillator With Fractional Derivative Element.

The linear damping coefficient $\delta 1$ is set as $\delta 1=\u22125$ so that the oscillator (51) takes diffusive limit-cycle motion. The other parameter values are chosen as $\omega 0=1$, $\epsilon =0.01$, $\lambda =1$, $\alpha =0.5$, $\delta 3=1,$ and $S0=1$. Further, $\Delta \eta =0.5$, $\rho =3.7,$ and $m=220$ are used when solving Eq. (53).

The analytical results of response PDFs and the corresponding MCS results for comparison are depicted in Fig. 4. The oscillator behavior approaches diffusive limit-cycle motion with time. Accordingly, the displacement and velocity PDFs change from unimodal shapes to bimodal ones. It can be seen that the present analytical technique estimates such highly non-Gaussian PDFs and the time evolution of their shapes with high accuracy. Furthermore, due to the behavior of diffusive limit cycle, the amplitude PDF also has a unique shape that is different from the PDFs in the previous two examples. Nevertheless, the assumptions about the amplitude PDF are still valid in this example.

## 7 Conclusions

A CFM-based analysis has been performed for obtaining the transient response PDF of nonlinear oscillators with fractional derivative elements under Gaussian white noise. The oscillators have polynomial-type nonlinearities in both stiffness and damping.

First, the equivalent natural frequency and equivalent damping coefficient of the oscillator were determined. The fractional derivative element of the oscillator contributes to both elasticity and viscosity. For this reason, the equivalent coefficients need to be determined by taking account of not only the effects of nonlinearities but also the contributions of the fractional element. Furthermore, the coefficients must be given in polynomial form of response amplitude to utilize the CFM-based approach. This paper proposed formulas to determine the equivalent coefficients that satisfy the above requirements, based on an equivalent linearization technique. Then, the FP equation was derived using stochastic averaging. Further, by the use of the Mellin transform and the relation between CFMs with different real parts of order, the FP equation was converted into a set of linear ordinary differential equations for amplitude CFMs. Solving it combined with a constraint that corresponds to the normalization condition for a PDF, the amplitude CFMs were obtained. Finally, the inverse Mellin transform of the resulting CFMs yielded the amplitude PDF. The joint PDF of displacement and velocity was also attained from the amplitude PDF.

In numerical examples, three linear and nonlinear oscillators with fractional derivatives have been considered to demonstrate the accuracy and versatility of the present analytical technique. It has been shown that the present method accurately determines the shapes of response PDFs that vary depending on both the fractional order and nonlinearities of the oscillator.

The practical importance of acquiring response PDFs is that system reliability can be assessed from the PDF. Therefore, it is highly significant to perform reliability analysis of fractional nonlinear systems using response PDF solutions obtained by the present analytical technique. Such a study would also be an effective and practical way to quantitatively evaluate the accuracy of the present analytical method. This topic will be addressed in our future work.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant No. JP21K03944.

## Conflict of Interest

The authors have no competing interests to declare that are relevant to the content of this paper.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Appendix

This appendix provides a brief summary of a method to approximately convert an equation of motion including fractional derivative elements into a set of integer-order differential equations, which is used in MCS. Readers interested in more details of the methodology are referred to Ref. [63].