Front. Phys. Frontiers in Physics Front. Phys. 2296-424X Frontiers Media S.A. 1238901 10.3389/fphy.2023.1238901 Physics Original Research Fractional stochastic vibration system under recycling noise Zhang et al. 10.3389/fphy.2023.1238901 Zhang Jian-Gang * Wang Fang Wang Hui-Nan School of Mathematics and Physics, Lanzhou Jiaotong University, Lanzhou, China

Edited by: Ji-Huan He, Soochow University, China

Reviewed by: Muhammad Nadeem, Qujing Normal University, China

Yajie Li, Henan University of Urban Construction, China

Ain Qura Tul, Guizhou University, China

 *Correspondence: Jian-Gang Zhang, Zhangjg7715776@126.com 02 08 2023 2023 11 1238901 12 06 2023 18 07 2023 Copyright © 2023 Zhang, Wang and Wang. 2023 Zhang, Wang and Wang

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The fractional stochastic vibration system is quite different from the traditional one, and its application potential is enormous if the noise can be deployed correctly and the connection between the fractional order and the noise property is unlocked. This article uses a fractional modification of the well-known van der Pol oscillator with multiplicative and additive recycling noises as an example to study its stationary response and its stochastic bifurcation. First, based on the principle of the minimum mean square error, the fractional derivative is equivalent to a linear combination of damping and restoring forces, and the original system is simplified into an equivalent integer order system. Second, the Itô differential equations and One-dimensional Markov process are obtained according to the stochastic averaging method, using Oseledec multiplicative ergodic theorem and maximal Lyapunov exponent to judge local stability, and judging global stability is done by using the singularity theory. Lastly, the stochastic D-bifurcation behavior of the model is analyzed by using the Lyapunov exponent of the dynamical system invariant measure, and the stationary probability density function of the system is solved according to the FPK equation. The results show that the fractional order and noise property can greatly affect the system’s dynamical properties. This paper offers a profound, original, and challenging window for investigating fractional stochastic vibration systems.

 van der Pol system fractional derivative recycling noise stochastic averaging method stochastic bifurcation section-at-acceptance Interdisciplinary Physics 1 Introduction

Fractional derivative [1, 2] is an extension of the theory of integer derivatives, and the study of fractional derivatives has a history of over 300 years. Some new materials have appeared, e.g., viscoelastic materials, nanomaterials, cement mortar, 3D-printed materials, and porous materials [3–8], which are different from either a solid or a fluid, and their constitutive relation is extremely difficult to be expressed correctly by the traditional calculus though much effort has been made to solve the problem, for example, using the fractal viscoelastic model [9] and the fractal rheological model [10]; the intractable constitutive relation has not yet been solved.

Considering its memory property, we consider that fractional calculus might be the best candidate for stochastic dynamical systems [11, 12]. Stochastic disturbances are widespread in nature, and fractional stochastic systems have become a hot spot in both mathematics and physics to deal with noise excitation. For example, energy-harvesting devices [13–15] are always subject to random excitation, and a fractional model can effectively reveal the bifurcation properties and multiple attractors of the energy-harvesting system, for example, Ref. [16]. Fractional models for Gaussian white noise also caught much attention [17–21], and the fractional convolution kernel neural network is a suitable mathematical tool for fault diagnosis [22–24]. Duffing oscillator [25, 26] is extended to its fractional partner under noise [27, 28]. Van der Pol oscillator [29] is another widely used model for the analysis of fractional stochastic P-bifurcation [30, 31].

In reality, noise exists in all aspects of practical applications, especially in nonlinear systems. The properties of stationary response, energy-harvesting efficiency, stability, and bifurcation will be greatly affected by noise excitation. At present, the research on the dynamic behavior of systems driven by recycling noise has attracted widespread attention from domestic and foreign scholars and achieved fruitful results, especially in the birhythmic biological system [32], stochastic resonance in asymmetric bistable systems [33], and the double entropic stochastic resonance phenomenon [34]. In this article, the fractional van der Pol model with recycling noise is adopted to investigate its dynamical properties.

 2 Model description

Balthazar van der Pol is a famous electronic engineer in the Netherlands. In 1927, he first deduced the famous van der Pol equation in order to describe the oscillation effect of triodes in electronic circuits, as shown below: x ¨ − μ 1 − x 2 x ˙ + x = 0

Afterward, as a classic nonlinear dynamic system, it is often used in mathematics and some nonlinear dynamic systems to demonstrate its dynamic behavior characteristics. In continuous research, the highest number of nonlinear terms considered is also constantly increasing, and there are also various methods for solving approximate solutions of such equations [35, 36]. From the classical van der Pol equation, changing the order of the equation can obtain systems with different dynamic behaviors, thereby better obtaining the dynamic behavior characteristics of the system. Therefore, we use the following equation to introduce the fractional generalized van der Pol model with multiplicative and additive recycling noise: x ¨ − − ε + α 1 x 2 − α 2 x 4 + α 3 x 6 − α 4 x 8 x ˙ + ω 2 x + D p x 0 c = η 1 t + x t η 2 t , where ε is the damping coefficient, α 1 , α 2 , α 3 , α 4 are nonlinear damping coefficients, ω is the frequency, η 1 t and η 2 t are independent recycling noises, i.e., D 1 ≠ D 2 , η i t = ξ i t + k ξ i t − τ , i = 1,2 . The power spectral density of recycling noise is obtained as: S i ω = 2 D i 1 + k 2 + 2 k c o s ω τ , i = 1 , 2 .

D p x t 0 c is the Caputo fractional derivative [1, 2] of p   0 ≤ p ≤ 1 order about x t defined as: D p 0 c x t = 1 Γ m − p ∫ 0 t x m u t − u 1 + p − m d u , m − 1 < p ≤ m , m ∈ N .

There are other definitions of fractional derivatives, for example, two-scale fractal derivative [37–41] and He’s fractional derivative [42]. The Caputo fractional derivative has memory property [43, 44], so it is used for the present study.

The D p 0 c x term in Eq. 1 can be expressed in a combination of spring stiffness and damping terms [45–48], hence, Eq. 1 becomes: x ¨ − − ε + α 1 x 2 − α 2 x 4 + α 3 x 6 − α 4 x 8 + C p x ˙ + ω 2 + K p x = η 1 t + x t η 2 t , where C and K are the equivalent damping and stiffness coefficients of fractional damping, respectively.

In order to identify C and K , we introduce an error function, which reads e = − C p x ˙ + K p x − D p x t 0 c ,

According to the minimum mean square method [44], we have ∂ E e 2 / ∂ C p = 0 , ∂ E e 2 / ∂ K p = 0 .

Equation 6 leads to the following equations: E − C p x ˙ 2 + K p x x ˙ − x ˙ D p 0 c x = lim T → ∞ 1 T ∫ 0 T − C p x ˙ 2 + K p x x ˙ − x ˙ D p 0 c x d t = 0 , E − C p x x ˙ + K p x 2 − x D p 0 c x = lim T → ∞ 1 T ∫ 0 T − C p x x ˙ + K p x 2 − x D p 0 c x d t = 0 .

Assuming that x t = a t cos ⁡ φ t = a t cos ω t + θ and a ˙ t ≈ 0 , we have x ˙ t = − a t ω sin ⁡ φ t , x ¨ t = − a t ω 2 cos ⁡ φ t .

Considering Eq. 8, 9, we re-write Eq. 7 in the form lim T → ∞ 1 T ∫ 0 T − C p x ˙ 2 + K p x x ˙ − x ˙ D p x 0 c d t = lim T → ∞ 1 T ∫ 0 T − C p a 2 t ω 2 sin 2 ⁡ φ t − K p a 2 t ω φ t cos ⁡ φ t + a t ω sin ⁡ φ t D p x 0 c d φ ≈ − C p a 2 ω 2 + 1 Γ 1 − p lim T → ∞ 1 T ∫ o T a ω sin ⁡ φ ∫ 0 t x ˙ t − τ τ p d τ d φ = − C p a 2 ω 2 − 1 Γ 1 − p lim T → ∞ 1 T ∫ 0 T a 2 ω sin ⁡ φ × ∫ 0 t sin ⁡ φ cos ω τ − cos ⁡ φ sin ω τ τ p d τ d t = 0 ,

For the same reason, we have lim T → ∞ 1 T ∫ 0 T − C p x x ˙ + K p x 2 − x D p x 0 c d t = lim T → ∞ 1 T ∫ 0 T − C p a 2 t ω sin ⁡ φ t cos ⁡ φ t + K p a 2 t cos 2 ⁡ φ t − a t cos ⁡ φ t D p x 0 c d φ ≈ K p a 2 2 ω − 1 Γ 1 − p lim T → ∞ 1 T ∫ o T acos ⁡ φ ∫ 0 t x ˙ t − τ τ p d τ d φ = K p a 2 2 ω + 1 Γ 1 − p lim T → ∞ 1 T ∫ 0 T a 2 × cos ⁡ φ ∫ 0 t sin ⁡ φ cos ω τ − cos ⁡ φ sin ω τ τ p d τ d t = 0 .

Hence lim T → ∞ 1 T ∫ 0 T − C p x ˙ 2 + K p x x ˙ − x ˙ D p x 0 c d t = − C p a 2 ω 2 − 1 Γ 1 − p lim T → ∞ 1 T ∫ 0 T a 2 ω sin ⁡ φ ∫ 0 t sin ⁡ φ cos ω τ − cos ⁡ φ sin ω τ τ p d τ d t = 0 , lim T → ∞ 1 T ∫ 0 T − C p x x ˙ + K p x 2 − x D p x 0 c d t = K p a 2 2 ω + 1 Γ 1 − p lim T → ∞ 1 T ∫ 0 T a 2 cos ⁡ φ ∫ 0 t sin ⁡ φ cos ω τ − cos ⁡ φ sin ω τ τ p d τ d t = 0 .

To simplify Eq. 10 further, we use the following asymptotic integrals ∫ 0 t cos ω τ τ p d τ = ω p − 1 Γ 1 − p sin p π 2 + sin ω t ω t p + ο ω t − p − 1 , ∫ 0 t sin ω τ τ p d τ = ω p − 1 Γ 1 − p cos p π 2 − cos ω t ω t p + ο ω t − p − 1 .

In view of Eq. 11, the integral averaging of Eq. 10 with respect to φ results in C p = − ω p − 1 sin p π 2 , K p = ω p sin p π 2 .

Hence, the equivalent system (4) can be written in the form x ¨ − λ x + ω 0 2 x = η 1 t + x t η 2 t , where λ = − ε + α 1 x 2 − α 2 x 4 + α 3 x 6 − α 4 x 8 − ω p − 1 sin p π 2 , ω 0 2 = ω 2 + ω p cos p π 2 .

3 Model processing

Now the problem becomes relatively simple; we assume that the solution of Eq. 13 can be expressed as [49]. X = x t = a t cos ⁡ Φ t , Y = x ˙ t = − a t ω 0 sin ⁡ Φ t , Φ t = ω 0 t + θ t , where a t and θ t are the amplitude and initial phase of the system, respectively.

We re-write Eq. 13 in the form x ˙ = y , y ˙ = λ y − ω 0 2 x t + η 1 t + x t η 2 t .

By Eq. 15 and the stochastic averaging method [50], Eq. 16 becomes d a d t = F 11 a , θ + G 11 a , θ η 1 t + G 12 a , θ η 1 t , d θ d t = F 21 a , θ + G 21 a , θ η 1 t + G 22 a , θ η 1 t , where F 11 a , θ = a sin 2 ⁡ Φ − ε + α 1 a 2 cos 2 ⁡ Φ − α 2 a 4 cos 4 ⁡ Φ + α 3 a 6 cos 6 ⁡ Φ − α 4 a 8 cos 8 ⁡ Φ − ω p − 1 sin p π 2 , F 21 a , θ = sin ⁡ Φ cos ⁡ Φ − ε + α 1 a 2 cos 2 ⁡ Φ − α 2 a 4 cos 4 ⁡ Φ + α 3 a 6 cos 6 ⁡ Φ − α 4 a 8 cos 8 ⁡ Φ − ω p − 1 sin p π 2 , G 11 a , θ = − sin ⁡ Φ ω 0 , G 12 a , θ = − asin ⁡ Φ cos ⁡ Φ ω 0 , G 21 a , θ = − cos ⁡ Φ a ω 0 , G 22 a , θ = − cos 2 ⁡ Φ ω 0 .

The recycling noise is a stationary process and can be approximated by a 2-D diffusion process. After stochastic averaging, the drift and diffusion coefficients are as follows: m 1 = F 11 + ∫ − ∞ 0 − cos ⁡ Φ ω 0 − cos ⁡ Φ t + τ 1 a ω 0 + − cos ⁡ Φ sin ⁡ Φ ω 0 − asin ⁡ 2 Φ t + τ 1 2 ω 0 + − a cos 2 ⁡ Φ − sin 2 ⁡ Φ ω 0 − cos 2 Φ t + τ 1 ω 0 R τ 1 d τ 1 = F 11 + cos 2 ⁡ Φ a ω 0 2 S 1 1 + a cos 2 ⁡ Φ sin 2 ⁡ Φ + acos ⁡ 2 Φ cos 2 ⁡ Φ ω 0 2 S 2 1 , m 2 = F 21 + ∫ − ∞ 0 cos ⁡ Φ a 2 ω 0 − sin ⁡ Φ t + τ 1 ω 0 + sin ⁡ Φ a ω 0 − cos ⁡ Φ t + τ 1 a ω 0 + 2 cos ⁡ Φ sin ⁡ Φ ω 0 − cos ⁡ Φ 2 t + τ 1 ω 0 R τ 1 d τ 1 = F 21 − 2 cos ⁡ Φ sin ⁡ Φ a 2 ω 0 2 S 1 1 − 2 cos 3 ⁡ Φ sin ⁡ Φ ω 0 2 S 2 1 , B 11 = ∫ − ∞ + ∞ − sin ⁡ Φ ω 0 − asin ⁡ Φ t + τ 1 ω 0 R τ 1 d τ 1 = 2 sin 2 ⁡ Φ ω 0 2 S 1 1 , B 12 = ∫ − ∞ + ∞ − asin ⁡ 2 Φ 2 ω 0 − asin ⁡ 2 Φ t + τ 1 2 ω 0 R τ 1 d τ 1 = 2 a 2 cos 2 ⁡ Φ sin 2 ⁡ Φ ω 0 2 S 2 1 , B 21 = ∫ − ∞ + ∞ − cos ⁡ Φ a ω 0 − cos ⁡ Φ t + τ 1 a ω 0 R τ 1 d τ 1 = 2 cos 2 ⁡ Φ a 2 ω 0 2 S 1 1 , B 22 = ∫ − ∞ + ∞ − cos 2 ⁡ Φ ω 0 − cos 2 ⁡ Φ t + τ 1 ω 0 R τ 1 d τ 1 = 2 cos 4 ⁡ Φ ω 0 2 S 2 1 , where S i 1 is the value of power spectral density of η i t at ω = 1 . S i 1 = 2 D i 1 + k 2 + 2 k c o s τ , i = 1,2

For the deterministic averaging of φ t , we have m ¯ 11 = 1 2 π ∫ 0 2 π F 11 a , θ + cos 2 ⁡ Φ a ω 0 2 S 1 1 + a cos 2 ⁡ Φ sin 2 ⁡ Φ + acos ⁡ 2 Φ cos 2 ⁡ Φ ω 0 2 S 2 1 d Φ = − 1 2 a ε + ω p − 1 sin p π 2 + α 1 a 3 8 − α 2 a 5 16 + 5 α 3 a 7 128 − 7 α 4 a 9 256 + S 1 1 2 a ω 0 2 + 3 a S 2 1 8 ω 0 2 m ¯ 22 = 1 2 π ∫ 0 2 π F 21 a , θ − 2 cos ⁡ Φ sin ⁡ Φ a 2 ω 0 2 S 1 1 − 2 cos 3 ⁡ Φ sin ⁡ Φ ω 0 2 S 2 1 d Φ = 0 , B ¯ 11 = 1 2 π ∫ 0 2 π 2 sin 2 ⁡ Φ ω 0 2 S 1 1 d Φ = S 1 1 ω 0 2 , B ¯ 12 = 1 2 π ∫ 0 2 π 2 a 2 cos 2 ⁡ Φ sin 2 ⁡ Φ ω 0 2 S 2 1 d Φ = a 2 S 2 1 4 ω 0 2 , B ¯ 21 = 1 2 π ∫ 0 2 π 2 cos 2 ⁡ Φ a 2 ω 0 2 S 1 1 d Φ = S 1 1 a 2 ω 0 2 , B ¯ 22 = 1 2 π ∫ 0 2 π 2 cos 4 ⁡ Φ ω 0 2 S 2 1 d Φ = 3 S 2 1 4 ω 0 2 .

The corresponding Itô SDE is d a = m 1 a d t + σ 11 2 a d B 1 t + σ 12 2 a d B 2 t , d θ = m 2 a d t + σ 21 2 a d B 1 t + σ 22 2 a d B 2 t , where m 1 a = − 1 2 a ε + ω p − 1 sin p π 2 + α 1 a 3 8 − α 2 a 5 16 + 5 α 3 a 7 128 − 7 α 4 a 9 256 + S 1 1 2 a ω 0 2 + 3 a S 2 1 8 ω 0 2 , m 2 a = 0 , σ 11 2 a = S 1 1 ω 0 2 , σ 12 2 a = a 2 S 2 1 4 ω 0 2 , σ 21 2 a = S 1 1 a 2 ω 0 2 , σ 22 2 a = 3 S 2 1 4 ω 0 2 .

The one-dimensional Markov Process can be expressed as: d a = H 1 a 8 + α 1 a 3 8 − α 2 a 5 16 + 5 α 3 a 7 128 − 7 α 4 a 9 256 + H 2 2 a d t + H 2 1 2 d B 1 t + H 3 a 2 4 1 2 d B 2 t , where H 1 = − 4 ε + ω p − 1 sin p π 2 + 3 S 2 1 ω 0 2 , H 2 = S 1 1 ω 0 2 , H 3 = S 2 1 ω 0 2 .

4 Stochastic stability analysis 4.1 The local stochastic stability

Considering the case of α 1 = α 2 = α 3 = α 4 = H 2 = 0 and linear Itô stochastic stability, from Eq. 24, we obtain d a = H 1 8 a d t + H 3 4 a 2 1 2 d B 2 t , m a = H 1 8 a , σ a = H 3 4 1 2 a .

Therefore, it is obtained that m ˙ 0 = H 1 / 8 and σ ˙ 12 0 = H 3 / 4 1 / 2 , using Oseledec multiplicative ergodic theorem [51] and maximal Lyapunov exponents to judge local stability. According to Itô stochastic differential equation, the solution of Eq. 26 is a t = a 0 exp ∫ 0 t m ˙ 0 − σ ˙ 12 0 2 2 d s + ∫ 0 t σ ˙ 12 0 d B 2 s ,

Then the approximate solution of the Lyapunov exponent of Itô stochastic differential equation is obtained λ = lim t → + ∞ 1 t ln x t , t 0 = lim t → + ∞ 1 t ln a t 1 2 = m ˙ 0 − σ ˙ 12 0 2 2 / 2 = 1 2 H 1 8 − H 3 8 .

When H 1 − H 3 < 0 , i.e., λ < 0 , Eq. 26 is stable in the sense of probability, and Eq. 16 is stable at the balance point. When H 1 − H 3 > 0 , i.e., λ < 0 , the effect is just the opposite.

 4.2 The global stochastic stability 4.2.1 Linear Itô stochastic stability

Judging global stability by the singularity theory, a = 0 is the first kind of singular boundary of Eq. 26. a = + ∞ is the second kind of singular boundary problem of Eq. 26. Calculating the diffusion index, drift indices, and characteristic value at boundary a = 0 and a = + ∞ , respectively, yields α a = 2 , β a = 1 , c a = lim a → 0 + 2 m a a − 0 α a − β a σ 12 2 a = lim a → 0 + 2 H 1 8 a 2 / H 3 4 a 2 = H 1 H 3 , α l = 2 , β l = 1 , c l = − lim a → + ∞ 2 m a a − 0 α l − β l σ 12 2 a = − lim a → + ∞ 2 H 1 8 a 2 / H 3 4 a 2 = − H 1 H 3 .

And the following conclusions are drawn, as shown in Table 1.

Global stability analysis.

Condition State Category Conclusion
H 1 / H 3 < 1 c a < 1 a = 0 Attract natural boundary The trivial solution of Eq. 26 is stable in the sense of probability, and the original system is probabilistically stable at the balance point
c l > − 1 a = + ∞ Exclude natural boundary
H 1 / H 3 > 1 c a > 1 a = 0 Exclude natural boundary The trivial solution of Eq. 26 is unstable in the sense of probability, and the original system is probabilistically unstable at the balance point
c l < − 1 a = + ∞ Attract natural boundary
H 1 / H 3 = 1 c a = 1 a = 0 Strict natural boundary The critical condition of system bifurcation
4.2.2 Stability of nonlinear Itô stochastic differential equation

When α 1 , α 2 , α 3 , α 4 , H 2 ≠ 0 , a = 0 is the first kind of singular boundary of Eq. 24. When a = + ∞ and m a = + ∞ , a = + ∞ is the second kind of singular boundary problem of Eq. 24. Calculating the diffusion index, drift indices, and characteristic value at boundary a = 0 and a = + ∞ , respectively, yields α a = 2 , β a = 1 , c a = lim a → 0 + 2 m a a − 0 α a − β a σ 11 2 a + σ 12 2 a = lim a → 0 + − a ε + ω p − 1 sin p π 2 + α 1 a 3 4 − α 2 a 5 8 + 5 α 3 a 7 64 − 7 α 4 a 9 128 + S 1 1 a ω 0 2 + 3 a S 2 1 4 ω 0 2 a / S 1 1 ω 0 2 + a 2 S 2 1 4 ω 0 2 = lim a → 0 + ω 0 2 − 128 a 2 ε + ω p − 1 sin p π 2 + 32 α 1 a 4 − 16 α 2 a 6 + 10 α 3 a 8 − 7 α 4 a 10 + 128 S 1 1 + 96 a 2 S 2 1 128 S 1 1 + 32 a 2 S 2 1 = 1 , lim α l = 2 , β l = 9 , c l = − lim a → + ∞ 2 m a a − 0 α l − β l σ 11 2 a + σ 12 2 a = − lim a → + ∞ − a ε + ω p − 1 sin p π 2 + α 1 a 3 4 − α 2 a 5 8 + 5 α 3 a 7 64 − 7 α 4 a 9 128 + S 1 1 a ω 0 2 + 3 a S 2 1 4 ω 0 2 a − 7 / S 1 1 ω 0 2 + a 2 S 2 1 4 ω 0 2 = − lim a → + ∞ ω 0 2 − 128 a 2 ε + ω p − 1 sin p π 2 + 32 α 1 a 4 − 16 α 2 a 6 + 10 α 3 a 8 − 7 α 4 a 10 + 128 S 1 1 + 96 a 2 S 2 1 128 S 1 1 a 8 + 32 S 2 1 a 10 = − 7 α 4 ω 0 2 32 S 2 1 .

Conclusion: when a = 0 and c a = 1 are a strict natural boundary; when a = + ∞ , c l > − 1 , and α 4 ω 0 2 / S 2 1 < 32 / 7 , the boundary is an exclude natural boundary; when c l < − 1 and α 4 ω 0 2 / S 2 1 > 32 / 7 , the boundary is an attract natural boundary; when c l = − 1 and α 4 ω 0 2 / S 2 1 = 32 / 7 , the boundary is a strict natural boundary. Therefore, c a = 1 is a critical condition of system bifurcation.

 5 Stochastic bifurcation analysis 5.1 D-bifurcation

If H 2 = H 3 = 0 , Eq. 13 becomes a deterministic system without a stochastic bifurcation phenomenon. Therefore, discussing the situation of H 3 ≠ 0 and α 1 = α 2 = α 3 = α 4 = H 2 = 0 , let σ 12 a = H 3 / 4 1 / 2 a and m a = H 1 / 8 − H 3 / 8 a , then the continuous dynamic system generated by Eq. 26 is ψ 1 t x = x + ∫ 0 t m ψ 1 s x d s + ∫ 0 t σ ψ 1 s x ∘ d B ,

Equation 27 is the only strong solution of Eq. 26 with x as the initial value. When m 0 = 0 and σ 12 0 = 0 , let m a be bounded, for all a ≠ 0 , the elliptic condition σ 12 0 ≠ 0 is satisfied, so there is only one stationary probability density. Therefore, the FPK equation corresponding to Eq. 26 is obtained. ∂ p ∂ t = − ∂ ∂ a H 1 8 a p − ∂ 2 ∂ a 2 H 3 4 a 2 p .

Let ∂ p / ∂ t = 0 get the stationary probability density corresponding to Eq. 28 p a = c σ 12 − 1 a exp ∫ 0 a 2 m u σ 12 2 u d u .

At this time, Eq. 27 has a non-trivial stationary state and a fixed-point equilibrium state. Assuming the invariant measures of these two kinds of stationary states are υ 1 and ϑ 1 , respectively, the density is Eq. 29 and ϑ 1 x , respectively. Hence, the solution of Eq. 28 is a t = a 0 exp ∫ 0 t m ˙ a + σ 12 a σ ¨ 12 a 2 d s + ∫ 0 t σ ˙ 12 a d B 2 .

The Lyapunov exponent of ψ 1 with respect to estimate u can be defined as follows λ ψ 1 u = lim t → + ∞ 1 t ln a t ,

Substituting Eq. 30 into Eq. 31, here σ 12 0 = 0 and σ ˙ 12 0 = 0 , its Lyapunov exponent of the fixed-point reads λ ψ 1 ϑ 1 = lim t → + ∞ 1 t ln a 0 + m ˙ 0 ∫ 0 t d s + σ ˙ 12 0 ∫ 0 t d B 2 s = m ˙ 0 + σ ˙ 12 0 lim t → + ∞ B 2 t t = m ˙ 0 = H 1 8 − H 3 8 .

Invariant estimate υ 1 with Eq. 29 as density. Substituting Eq. 30 into Eq. 31. Assuming that σ ˙ and m ˙ + σ σ ˙ are bounded and integrable, respectively, the Lyapunov exponent can be obtained λ ψ 1 υ 1 = lim t → + ∞ 1 t ∫ 0 t m ˙ a + σ 12 a σ ¨ 12 a d s = ∫ R m ˙ a + σ 12 a σ ¨ 12 a 2 p a d a = − 2 ∫ R m a σ 12 a 2 p a d a m ˙ 0 = − 2 H 2 3 / 2 H 1 8 − H 3 8 exp 8 H 3 H 1 8 − H 3 8 .

Let α = H 1 − H 3 , when α < 0 and H 1 < H 3 , ϑ 1 is stable, υ 1 is unstable; when α > 0 and H 1 > H 3 , ϑ 1 is unstable, υ 1 is stable. So α is a D-bifurcation point of Eq. 13.

 5.2 P-bifurcation 5.2.1 Stochastic P-bifurcation under additive recycling noise

When additive noise just exists, D 1 ≠ 0 and D 2 = 0 . The following is an analysis of the stochastic P-bifurcation of the system in this case. Eq. 22, 23 show that the Itô stochastic differential equation corresponding to a t does not depend upon θ t , and it is a 1-D diffusion process; its corresponding FPK equation can be expressed as ∂ p a , t ∂ t = − ∂ ∂ a m 1 a p a , t + 1 2 ∂ 2 ∂ a 2 σ 11 2 a p a , t , the corresponding boundary conditions are p = c , c ∈ − ∞ , + ∞ , w h e n   a = 0 . p → 0 , ∂ p ∂ a → 0 , w h e n   a → ∞ .

In view of Eq. 35, the stationary probability density of the amplitude is p a = C σ 11 2 a exp ∫ 0 a 2 m 1 u σ 11 2 u d u , where C is the normalization constant, C = ∫ 0 ∞ 1 σ 11 2 a exp ∫ 0 a 2 m 1 u σ 11 2 u d u d a − 1 .

In view of Eq. 23, from Eq. 36, we obtain p a = C a ω 0 2 S 1 1 exp − a 2 ω 0 2 Δ 7680 S 1 1 , where Δ = 3840 ε + 3840 ω p − 1 sin p π 2 − 480 α 1 a 2 + 160 α 2 a 4 − 75 α 3 + 42 α 4 a 8 , ω 0 2 = ω 2 + ω p cos p π 2 , S 1 1 = 2 D 1 1 + k 2 + 2 k c o s τ .

The original system response meets a t = x 2 t + x ˙ 2 t , in view of Eq. 38, the joint probability density function of the system is p a = C x 2 t + x ˙ 2 t ω 0 2 S 1 t exp − x 2 t + x ˙ 2 t ω 0 2 3840 S 1 t Δ .

5.2.1.1 Influence of fractional order

As fractional damping is a combination of the equivalent stiff and equivalent damping, the fractional order is of paramount importance; its value can be calculated by He-Liu’s fractal formulation [52] for practical applications. According to Eq. 12, when p = 1 , fractional damping becomes a damping term, while when p = 0 , it is a stiff term.

Setting τ = 2 , k = 0.4 , ε = − 0.1 , α 1 = 1.51 , α 2 = 2.85 , α 3 = 1.693 , α 4 = 0.312 , and ω = 1 in Eq. 13 as that in Refs [30, 53], the stochastic P-bifurcation is studied hereby. Keeping D 1 = 0.005 constant, we draw the joint probability density function section and top view of Eq. 13 under the influence of different fractional orders.

When p = 0.06 , the joint probability density function diagram shows a crater shape; there is only one peak in the section, and there is only a large limit cycle. The response is shown as a vibration far beyond the origin (Figure 1).

Joint probability density function section and top view of Eq. 13 when p = 0.06 .

When p = 0.139 , from the section, it can be clearly seen that there are two peaks, but the second peak has a much larger amplitude. At this time, the system has a balance point and a large limit cycle; hence, the system response switches between two peaks, and the probability of a large amplitude vibration is high, as shown in Figure 2.

Joint probability density function section and top view of Eq. 13 when p = 0.137 .

When p = 0.141 , the section has three peaks, showing two peaks in addition to the origin. A balance point now coexists with a large and small limit cycle in the system, and the system response switches between the three peaks, which is a multimodal response. Due to the existence of the double limit point set, the relative heights of the joint probability density function peaks at the three peaks are different, implying that the system response peaks are different, as shown in Figure 3.

Joint probability density function section and top view of Eq. 13 when p = 0.14 .

When p = 0.145 , the section has two peaks, in contrast to Figure 2, the relative height of the peak changes, with the second peak being significantly smaller. At this time, the system has both a balance point and a small limit cycle; hence, the system response switches between two peaks, and the probability of a small amplitude vibration is high, as shown in Figure 4.

Joint probability density function section and top view of Eq. 13 when p = 0.143 .

Based on the above discussions, we conclude that the fractional order can cause stochastic P-bifurcation behavior in the system. From Figure 5, we find that an increasing fractional order will change the stationary response from a single mode to a dual mode and then to a tristable mode. The peak value changes from a single peak to two peaks and then to three peaks, so stochastic P-bifurcation occurs. Increasing the value of p to 0.145 again, the tristable disappears and the bistable appears; the peak value changes from three peaks to two peaks, so stochastic P-bifurcation occurs.

Stationary probability density function diagram of Eq. 13 when D 1 = 0.005 .

 5.2.1.2 Influence of noise intensity

Keeping the above parameters unchanged, and fixing p = 0.14 , we draw the joint probability density function section and top view of Eq. 13 under the influence of different noise intensity.

When D 1 = 0.03 , the joint probability density function diagram shows a crater shape, there is only one peak in the section, and there is only a large limit cycle. The response is shown as a vibration far beyond its origin, as shown in Figure 6.

Joint probability density function section and top view of Eq. 13 when D 1 = 0.03 .

When D 1 = 0.015 , from the section, it can be clearly seen that there are two peaks, but the first peak is much smaller. At this time, the system has both a balance point and a large limit cycle; hence, the system response switches between two peaks, and the probability of a large amplitude vibration is high, as shown in Figure 7.

Joint probability density function section and top view of Eq. 13 when D 1 = 0.015 .

When D 1 = 0.005 , the section has three peaks, showing two peaks in addition to the origin. A balance point now coexists with a large and small limit cycle in the system, and the system response switches among the three peaks, which is a multimodal response. Due to the existence of the double limit point set, the relative heights of the joint probability density function peaks at the three peaks are different, implying that the vibration frequency of the system response peak is different, as shown in Figure 8.

Joint probability density function section and top view of Eq. 13 when D 1 = 0.005 .

When D 1 = 0.0013 , the section has two peaks, in contrast to Figure 7; the relative height of the peak changes, with the first peak being much larger. At this point, the system has a balance point and a small limit cycle, the system response switches between two peaks, and the probability of a small amplitude vibration is high, as shown in Figure 9.

Joint probability density function section and top view of Eq. 13 when D 1 = 0.0013 .

Based on the above discussions, it can be verified that changing the noise intensity affects greatly stochastic P-bifurcation property. From Figure 10, it can also be seen that with noise intensity being reduced, the stationary response of the system switches from a single mode to a dual mode and then to a tristable mode. The peak value of the stationary probability density function curve changes from a single peak to two peaks and then three peaks, so stochastic P-bifurcation occurs. Decreasing the value of D 1 to 0.0013 again, the tristable disappears and the bistable appears; the peak value changes from three peaks to two peaks, so stochastic P-bifurcation occurs.

Stationary probability density function diagram of Eq. 13 when p = 0.141 .

 5.2.2 Additive and multiplicative recycling noise

When D 1 ≠ 0 and D 2 ≠ 0 , the expression of the stationary probability density function of the amplitude of Eq. 13 is p a = C σ 11 2 a + σ 12 2 a exp ∫ 0 a 2 m u σ 11 2 u + σ 12 2 u d u , where C is the normalization constant, C = ∫ 0 ∞ 1 σ 11 2 a + σ 12 2 a exp ∫ 0 a 2 m 1 u σ 11 2 u + σ 12 2 u d u d a − 1 .

In view of Eq. 23, from Eq. 42, we have p a = 4 C a ω 0 2 4 S 1 1 + a 2 S 2 1 − Δ 1 S 2 5 1 exp Δ 2 768 S 2 4 1 , where Δ 1 = 2 ω 0 2 ε + ω p − 1 sin p π 2 S 2 4 + α 1 S 1 S 2 3 + 2 α 2 S 1 2 S 2 2 + 5 α 3 S 1 3 S 2 + 14 α 4 S 1 4 , Δ 1 = a 2 ω 0 2 384 α 1 S 2 3 + 768 α 2 S 1 S 2 2 + 1920 α 3 S 1 2 S 2 + 5376 α 4 S 1 3 + a 4 ω 0 2 − 96 α 2 S 2 3 − 240 α 3 S 1 S 2 2 − 672 α 4 S 1 2 S 2 + a 6 ω 0 2 40 α 3 S 2 3 + 112 α 4 S 1 S 2 2 − 21 a 8 ω 0 2 α 4 S 2 3 , ω 0 2 = ω 2 + ω p cos p π 2 , S 1 1 = 2 D 1 1 + k 2 + 2 k c o s τ , S 2 1 = 2 D 2 1 + k 2 + 2 k c o s τ .

Keeping the above parameters unchanged, we draw the joint probability density function section and top view of Eq. 13 under the influence of different fractional orders and noise intensity.

When p = 0.05 , let D 1 = 0.5 and D 2 = 1 . The joint probability density function diagram shows a crater shape; there is only one peak in the section, and there is only a large limit cycle. The response is shown as a vibration far away from the origin. At the same time, reducing the value of noise intensity reveals that the peak of the joint probability density does not change with only one peak. However, the system only has a limit cycle, which has been in a monostable. Therefore, there is no stochastic P-bifurcation phenomenon occurring, as shown in Figure 11.

Joint probability density function section and top view of Eq. 13 when p = 0.05 .

When p = 0.1 , let D 1 = 0.002 and D 2 = 0.01 . From the section, it can be clearly seen that there are two peaks, but the second peak has a much larger amplitude. At this time, the system has both a balance point and a limit cycle; hence, the system response switches between two peaks, and the large amplitude vibration has a higher probability. When the simultaneous improvement of the noise intensifies to D 1 = 0.008 and D 2 = 0.2 , the peak value of the stationary probability density function curve changes from two peaks to one peak. There is only a large limit cycle, and the system response becomes a vibration far from the origin. Therefore, increasing the noise intensity induces a stochastic P-bifurcation property, as shown in Figure 12.

Joint probability density function section and top view of Eq. 13 when p = 0.1 .

When p = 0.14 , let D 1 = 0.003 and D 2 = 0.01 . The section has a peak near the origin. There is only a balance point in the system at this time, and the response is shown as a vibration closer to the origin. When simultaneously improving its noise intensity to D 1 = 0.004 and D 2 = 0.1 , the peak value of the stationary probability density function curve changes from a single peak to two peaks. At this time, the system has both a balance point and a limit cycle; hence, the system response switches between two peaks, and the probability of a large amplitude vibration is small. Therefore, increasing the noise intensity induces the stochastic P-bifurcation phenomenon. Further increasing the noise intensity to D 1 = 0.006 and D 2 = 0.15 , the peak value of the section changes relatively, and the first peak is lower. The system response switches between two peaks, and the probability of a small one is small. Continuing to increase the value of noise intensity to D 1 = 0.04 and D 2 = 0.3 , the peak value of the curve changes from two peaks to a single peak. There is only a limit cycle in the system at this time, and the response is shown as a vibration far away from the origin. Therefore, increasing the noise intensity induces the second stochastic P-bifurcation phenomenon, as show in Figure 13.

Joint probability density function section and top view of Eq. 13 when p = 0.14 .

 6 Conclusion

In this paper, the stationary response and the stochastic bifurcation of the fractional van der Pol equation under multiplicative and additive recycling noise excitations are investigated. By the least square method, we obtain an equivalent integral nonlinear stochastic system. The Itô differential equation and One-dimensional Markov process are obtained according to the stochastic averaging method. We discuss the local and global stochastic stability and analyze the conditions for inducing D-bifurcation and P-bifurcation in the system. The analysis shows that when α < 0 and H 1 < H 3 , the point equilibrium state becomes stable, and the non-trivial stationary state becomes unstable; when α > 0 and H 1 > H 3 , the result is the opposite. So α is a D-bifurcation point of the original system. When only additive noise exists, the fractional order and the noise intensity will greatly affect the system’s property. It was found that reducing the order p or increasing the noise intensity D 1 can cause nonlinear jumping or significant oscillation in the system, leading to system instability. Through increasing the order p or reducing the noise intensity D 1 , the system response is in a monostable state or a small disturbance near the balance point. Similarly, when additive and multiplicative noise coexist, selecting appropriate parameters can maintain the system response at a monostable or small disturbance near the balance point. Therefore, in practical engineering, to avoid the potential adverse effects of high noise intensity on the system, the occurrence of stochastic bifurcation behavior can be controlled by changing the noise intensity or fractional order. In the future, we will combine theory with practice to explore the impact of recycling noise on the stationary response and stochastic bifurcation of systems in wind turbines. We will study the impact of changing noise intensity and fractional order on the system, and how to handle these adverse effects to achieve optimal system performance.

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

 Author contributions

J-GZ: supervision, writing–review and editing, funding acquisition, investigation, and project administration. FW: writing–original draft, writing–review and editing, and software. H-NW: writing–review and editing and software. All authors contributed to the article and approved the submitted version.

 Funding

This study was supported by the key project of the Gansu Province natural science foundation of China (No. 23JRRA882).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References He JH . A tutorial review on fractal spacetime and fractional calculus. Int J Theor Phys (2014) 53:3698–718. 10.1007/s10773-014-2123-8 He JH . Fractal calculus and its geometrical explanation. Phys (2018) 10:272–6. 10.1016/j.rinp.2018.06.011 He CH Liu C . Fractal approach to the fluidity of a cement mortar. Nonlinear Engineering-modeling Appl (2022) 11(1):1–5. 10.1515/nleng-2022-0001 Zuo YT Liu HJ . Fractal approach to mechanical and electrical properties of graphene/sic composites. Facta Universitatis: Ser Mech Eng (2021) 19(2):271–84. 10.22190/fume201212003z He CH Liu C . Fractal dimensions of a porous concrete and its effect on the concrete’s strength. Facta Universitatis Ser Mech Eng (2023) 21(1):137–50. 10.22190/FUME221215005H Jankowski P . Detection of nonlocal calibration parameters and range interaction for dynamic of FGM porous nanobeams under electro-mechanical loads. Facta Universitatis Ser Mech Eng (2022) 20(3):457–78. 10.22190/fume210207007j He CH Shen Y Ji FY He JH . Taylor series solution for fractal Bratu-type equation arising in electrospinning process. Fractals (2020) 28(1):2050011. 10.1142/s0218348x20500115 Liu FJ Zhang T He CH Tian D . Thermal oscillation arising in a heat shock of a porous hierarchy and its application. Facta Universitatis Ser Mech Eng (2022) 20(3):633–45. 10.22190/fume210317054l Liang YH Wang KJ . A new fractal viscoelastic element: Promise and applications to Maxwell-Rheological model. Therm Sci (2021) 25(2):1221–7. 10.2298/tsci200301015l Zuo YT . Effect of Sic particles on viscosity of 3-D print paste a fractal rheological model and experimental verification. Therm Sci (2021) 25(3):2405–9. 10.2298/tsci200710131z Long Y Xu BB Chen DY Ye W . Dynamic characteristics for a hydro-turbine governing system with viscoelastic materials described by fractional calculus. Appl Math Model (2018) 58:128–39. 10.1016/j.apm.2017.09.052 Wang L Xue LL Xu W Yue XL . Stochastic P-bifurcation analysis of a fractional smooth and discontinuous oscillator via the generalized cell mapping method. Int J Non-Linear Mech (2017) 96:56–63. 10.1016/j.ijnonlinmec.2017.08.003 He CH Amer TS Tian D Abolila AF Galal AA . Controlling the kinematics of a spring-pendulum system using an energy harvesting device. J Low Frequency Noise: Vibration Active Control (2022) 41(3):1234–57. 10.1177/14613484221077474 He CH Ei-Dib YO . A heuristic review on the homotopy perturbation method for non-conservative oscillators. J Low FrequCency Noise, Vibration Active Control (2022) 41(2):572–603. 10.1177/14613484211059264 He JH Ei-Dib YO . The reducing rank method to solve third-order Duffing equation with the homotopy perturbation. Numer Methods Differential Equations (2021) 37(2):1800–8. 10.1002/num.22609 Zhang WT Xu W Niu LZ Tang YN . Bifurcations analysis of a multiple attractors energy harvesting system with fractional derivative damping under random excitation. Commun Nonlinear Sci Numer Simulation (2023) 118:107069. 10.1016/j.cnsns.2022.107069 Hu DL Mao XC Han L . Stochastic stability analysis of a fractional viscoelastic plate excited by Gaussian white noise. Mech Syst Signal Process (2022) 177:109181. 10.1016/j.ymssp.2022.109181 Li YJ Wu ZQ Lan QX Cai YJ Xu HF Sun YT . Stochastic transition behaviors in a Tri-Stable van der Pol oscillator with fractional delayed element subject to Gaussian White Noise. Therm Sci (2022) 26(3):2713–25. 10.2298/tsci2203713l Li YJ Wu ZQ Lan QX Cai YJ Xu HF Sun YT . Transition behaviors of system energy in a bi-stable van Ver Pol oscillator with fractional derivative element driven by multiplicative Gaussian white noise. Therm Sci (2022) 26(3):2727–36. 10.2298/tsci2203727l Li YJ Wu ZQ Wang F Zhang GQ Wang YC . Stochastic P-bifurcation in a generalized Van der Pol oscillator with fractional delayed feedback excited by combined Gaussian white noise excitations. J Low Frequency Noise, Vibration Active Control (2021) 40(1):91–103. 10.1177/1461348419878534 Din A Ain QT . Stochastic optimal control analysis of a mathematical model: Theory and application to non-singular kernels. fractal and fractional (2022) 6:279. 10.3390/fractalfract6050279 Zhu R Wang MX Xu SY Li K Han QP Tong X Fault diagnosis of rolling bearing based on singular spectrum analysis and wide convolution kernel neural network. J Low Frequency Noise, Vibration Active Control (2022) 41(4):1307–21. 10.1177/14613484221104639 Kuo PH Tseng YR Luan PC Yau HT . Novel fractional-order convolutional neural network based chatter diagnosis approach in turning process with chaos error mapping. Nonlinear Dyn (2023) 111:7547–64. 10.1007/s11071-023-08252-w Kuo PH Chen YW Hsieh TH Jywe WY Yau HT . A thermal displacement prediction system with an automatic Lrgtvac-PSO optimized branch Structured bidirectional GRU neural network. IEEE Sensors Journa (2023) 23:12574–86. 10.1109/JSEN.2023.3269064 He JH Jiao ML Gepreel KA Khan Y . Homotopy perturbation method for strongly nonlinear oscillators. Mathematics Comput Simulation (2023) 204:243–58. 10.1016/j.matcom.2022.08.005 He JH Jiao ML He CH . Homotopy perturbation method for fractal Duffing oscillator with arbitrary conditions. Fractals (2022) 30. 10.1142/S0218348X22501651 Chen LC Zhu WQ . Stochastic jump and bifurcation of Duffing oscillator with fractional derivative damping under combined harmonic and white noise excitations. Int J Non-Linear Mech (2011) 46:1324–9. 10.1016/j.ijnonlinmec.2011.07.002 Chen LC Zhu WQ . Stationary response of duffing oscillator with fractional derivative damping under combined harmonic and wide band noise excitations. Chin J Appl Mech (2010) 3:517–21. He CH Tian D Moatimid GM Salman HF Zekry MH . Hybrid Rayleigh–van der pol–duffing oscillator: Stability analysis and controller. J Low Frequency Noise, Vibration Active Control (2022) 41(1):244–68. 10.1177/14613484211026407 Li YJ Wu ZQ . Stochastic P-bifurcation in a tri-stable Van der Pol system with fractional derivative under Gaussian white noise. J Vibroengineering (2019) 21:803–15. 10.21595/jve.2019.20118 Li YJ Wu ZQ Lan QX Hao Y Zhang XY . Stochastic P bifurcation in a tri-stable van der Pol oscillator with fractional derivative excited by combined Gaussian white noises. J Vibration Shock (2021) 40(16):275–93. 10.21595/jve.2019.20118 Chamgoué AC Yamapi R Woafo P . Bifurcations in a biorhythmic biological system with time-delayed noise. Nonlinear Dyn (2013) 73:2157–73. 10.1007/s11071-013-0931-7 Wu YZ Sun ZK . Residence-times distribution function in asymmetric bistable system driven by noise recycling. Acta Phys Sin (2020) 69(12):120501. 10.7498/aps.69.20201752 Guo XY Cao TQ . Phenomenon of double entropic stochastic resonance with recycled noise. Chin J Phys (2022) 77:721–32. 10.1016/j.cjph.2021.10.020 He KY Nadeem M Habib S Sedighi HM Huang DH . Analytical approach for the temperature distribution in the casting-mould heterogeneous system. Int J Numer Methods Heat Fluid Flow (2021) 32:1168–82. 10.1108/HFF-03-2021-0180 Fang JH Nadeem M Habib M Karim S Wahash HA . A new iterative method for the approximate solution of klein-gordon and sine-gordon equations. J Funct Spaces (2022) 2022:1–9. 10.1155/2022/5365810 He JH Kou SJ He CH Zhang ZW Gepreel KA . Fractal oscillation and its frequency-amplitude property. Fractals (2021) 29(4):2150105. 10.1142/s0218348x2150105x He JH Moatimid G Zekry M . Forced nonlinear oscillator in a fractal space. Facta Universitatis, Ser Mech Eng (2022) 20(1):001–20. 10.22190/fume220118004h Tian D Ain QT Anjum N He CH Cheng B . Fractal N/MEMS: From pull-in instability to pull-in stability. Fractals (2021) 29:2150030. 10.1142/S0218348X21500304 He CH . A variational principle for a fractal nano/microelectromechanical (N/MEMS) system. Intermational J Numer Methods Heat Fluid Flow (2023) 33(1):351–9. 10.1108/hff-03-2022-0191 Ain QT Sathiyaraj Karim S Nadeem M Mwanakatwe PK Kandege Mwanakatwe P . ABC fractional derivative for the alcohol drinking model using two-scale fractal dimension. Complexity (2022) 2022:1–11. 10.1155/2022/8531858 Wang Y An JY . Amplitude-frequency relationship to a fractional Duffing oscillator arising in microphysics and tsunami motion. J Low Frequency Noise Vibration Active Control (2019) 38(3-4):1008–12. 10.1177/1461348418795813 Mendes EM Salgado GH Aguirre LA . Numerical solution of Caputo fractional differential equations with infinity memory effect at initial condition. Commun Nonlinear Sci Numer Simulation (2019) 69:237–47. 10.1016/j.cnsns.2018.09.022 Zeng HJ Wang YX Xiao M Wang Y . Fractional solitons: New phenomena and exact solutions. Front Phys (2023) 11:1177335. 10.3389/fphy.2023.1177335 Chen LC Wang WH Li ZS Zhu WQ . Stationary response of Duffing oscillator with hardening stiffness and fractional derivative. Int J Non-Linear Mech (2013) 48:44–50. 10.1016/j.ijnonlinmec.2012.08.001 Li W Zhang MT Zhao JF . Stochastic bifurcations of generalized Duffing-van der Pol system with fractional derivative under colored noise. Chin Phys B (2017) 26:090501. 10.1088/1674-1056/26/9/090501 Chen LC Li ZS Zhuang QQ Zhu WQ . First-passage failure of single-degree-of-freedom nonlinear oscillators with fractional derivative. J Vibration Control (2013) 19:2154–63. 10.1177/1077546312456057 Chen LC Zhu WQ . Stochastic response of fractional-order van der Pol oscillator. Theor Appl Mech Lett (2014) 4:013010. 10.1063/2.1401310 Spanos PD Zeldin BA . Random vibration of systems with frequency-dependent parameters or fractional derivatives. J Eng Mech (1997) 123(3):290–2. 10.1061/(asce)0733-9399(1997)123:3(290) Zhu WQ . Nonlinear stochastic dynamics and control: Hamilton theoretical framework. Beijing, China: Science Press (2003). Oseledec VI . A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems. Trans Mosc Math. Soc (1968) 19(2):197–231. He CH Liu C . Fractal dimensions of a porous concrete and its effect on the concrete’s strength. Facta Universitatis Ser Mech Eng (2023) 21(1):137–50. 10.22190/FUME221215005H Zhang XY Wu ZQ . Bifurcations in tri-stable Duffing–Van der Pol oscillator with recycling noise. Mod Phys Lett B (2018) 32(20):1850228. 10.1142/s0217984918502287