Next Article in Journal
Synchronization of Fractional Stochastic Chaotic Systems via Mittag-Leffler Function
Next Article in Special Issue
The Tracking Control of the Variable-Order Fractional Differential Systems by Time-Varying Sliding-Mode Control Approach
Previous Article in Journal
Dynamics, Periodic Orbit Analysis, and Circuit Implementation of a New Chaotic System with Hidden Attractor
Previous Article in Special Issue
Modified Predictor–Corrector Method for the Numerical Solution of a Fractional-Order SIR Model with 2019-nCoV
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise

1
Division of Dynamics and Control, School of Mathematics and Statistics, Shandong University of Technology, Zibo 255000, China
2
Department of Applied Mathematics and Statistics, Technical University of Cartagena, Hospital de Marina, 30203 Cartagena, Spain
3
Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589, Saudi Arabia
4
State Key Laboratory of Mechanics and Control of Mechanical Structures, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(4), 191; https://doi.org/10.3390/fractalfract6040191
Submission received: 15 February 2022 / Revised: 24 March 2022 / Accepted: 28 March 2022 / Published: 30 March 2022

Abstract

:
A general investigation on the mechanism of stochastic resonance is reported in a time-delay fractional Langevin system, which endues a nonlinear form multiplicative colored noise and fractional Gaussian noise. In terms of theoretical analysis, both the expressions of output steady-state amplitude and that of the first moment of system response are obtained by utilizing stochastic averaging method, fractional Shapiro and Laplace methods. Due to the presence of trichotomous colored noise, the excitation frequency can induce fresh multimodal Bona fide stochastic resonance, exhibiting much more novel dynamical behaviors than the non-disturbance case. It is verified that multimodal pattern only appears with small noise switching rate and memory damping order. The explicit expressions of critical noise intensity corresponding to the generalized stochastic resonance are given for the first time, by which it is determined that nonlinear form colored noise induces much more of a comprehensive resonant phenomena than the linear form. In the case of slow transfer rate noise, a newfangled phenomenon of double hypersensitive response induced by a variation in noise intensity is discovered and verified for the first time, with the necessary range of parameters for this phenomenon given. In terms of numerical scheme, an efficient and feasible algorithm for generating trichotomous noise is proposed, by which an algorithm based on the Caputo fractional derivative are applied. The numerical results match well with the analytical ones.

1. Introduction

Based on the Markov process without after-effect, classical stochastic resonance theory grounded on ordinary Langevin equation and FPE tends to be perfect gradually. However, more and more systems under complex environments exhibit memory effects, such as long-range correlation and historical dependence in time. At this time, the traditional stochastic resonance theory with ODE is obviously no longer applicable as it often needs to be modeled by non-traditional mathematical models such as fractional order differential equation. Early scholars tried to use the research method corresponding to a normal diffusion system to approximate this type of issue, e.g., fractional relaxation equation [1], fractional diffusion equation [2] and fractional Fokker–Planck equations [3]. Due to the characteristics of the fractional calculus equation itself, the systematic description of stochastic resonance in the fractional system cannot be achieved by using analytical algorithm in most cases at present.
The concept of stochastic resonance (SR) was firstly proposed in paleometeorology to explain why the climatic alternating period of warm climate and glacial period can stay in sync with the offset period of the eccentric angle of earth revolution. The synergistic effect of nonlinearity, periodic signal imposed by sun and random forces of earth itself drive the response signal cross the potential barrier, and make the climate alternation possible. However, the nonlinearity is not a necessary condition for the occurrence of stochastic resonance, say, SR can exist in linear systems disturbed by multiplicative random noise [4,5]. As research on stochastic resonance in linear systems started later than that on nonlinear systems, the former is far inferior to the latter in terms of the diversity of research methods and application range. While studies of SR in fractional systems started much later than that in integer order systems, various bistable fractional order systems are gradually taken as the main research object, and became an important topic up to the beginning of the present century. Reports on SR in linear fractional systems has been gradually developed in recent years, Ren [6] considers fractional damping fluctuations manifests as the form of trichotomous noise for the second order linear FLE, confirming that the interaction between multiplicative noise and memory effect can lead to stochastic resonance. Kim [7] analyzes self-tuning stochastic resonance energy harvesting for rotating systems under modulated noise. For all this, the stochastic resonance of linear fractional-order systems is still at the preliminary exploratory stage.
In the studies on SR of fractional order systems induced by colored noise, dichotomous and trichotomous random telephone signal (RTS) are often presented in multiplicative and additive noise terms, to describe the disturbance on oscillator mass [8], inherent frequency [9], fractional damping [10], signal modulating noise [7], etc. Due to the one-sided understanding of the complex disturbance in practical application, noises in linear form are often used to simulate the random factors for formal convenience. However, nonlinear noise exists more widely in real world systems than linear ones [11]. Many nonlinear noises such as quadratic noise [12], nonlinear Markov noise [13] and nonlinear Ornstein–Uhlenbeck noise [14] really exist in the actual physical system. As for the nonlinear dichotomous and trichotomous noise, it has attracted more and more attention on account of its universal existence and potential application in practical nonlinear physical system. On the one hand, there are research works on a variety of nonlinear noises that mainly target integer order systems, and quite a few studies involving nonlinear color noise have been reported in fractional order systems. Zhong [15] studied underdamped and overdamped fractional order systems with nonlinear color noise natural frequency disturbance, respectively. As far as we are aware, no previous work on the research of SR in second-order fractional-order anomalous diffusion systems with fractional damping disturbance of polynomial form has been reported, ever. On the other hand, on account of different propagation mechanisms, both signal and noise have a time delay effect when transmitting through the system [16], and time delays exist in many natural or artificial systems. It is found that the coexistence of noise and time delay tend to change the dynamic behavior of stochastic systems that are described by the Langevin equation [17], and this kind of synergy can also have a significant effect on the resonance phenomena [18]. Nevertheless, most studies on the effect of time delay on stochastic resonance are limited to integer order systems, and, in a numerical simulation way, analytical studies on SR for fractional order systems with time delay are quite rare. Taking all these into account, we consider a class of anomalous diffusion systems in regard to fractional Brown particles, containing small delay, fractional Gaussian internal noise, and driven by nonlinear trichotomous colored noise and external periodic excitation simultaneously.
In this paper, stochastic resonance induced by fractional damping disturbance in the form of polynomial trichotomous noise is investigated, qualitatively and quantitatively. A new feasible numerical generating algorithm for trichotomous colored noise is proposed, by which a numerical scheme for the stochastic fractional Langevin system is established. The effect of the noise switching rate and fractional damping order on the occurrence of stochastic multi-resonance is expounded, which is an exploration in the category of fractional systems with stochastic damping disturbance. In addition, some novel phenomena of hypersensitivity response are discovered, which have not been reported before, revealing the complexity in a system driven by nonlinear telegraphic noise. By this, the resonance regimes in the complex system will be understood more deeply. The organization of this paper is as follows. In Section 2, starting with the traditional Brown motion, a particle diffusion model described by the fractional Langevin equation (FLE) is introduced by virtue of the generalized Langevin equation. Section 3 is devoted to the theoretical results of the exact solution of the response, from the first moment, and the steady-state amplitude by extending the stochastic averaging method. In Section 4, detailed analysis on the bona fide stochastic resonance under different system parameter conditions are investigated. Section 5 presents the dependence of generalized stochastic resonance phenomena on time delay, noise intensity and fractional order of the damping. In Section 6 the numerical verification of the theoretical predictions is provided. The conclusion is given in the Section 7.

2. Model

In statistical thermodynamics, due to the anfractuous interaction of nonlinearity, noise and disturbances entrained in different sources, most natural systems manifest themselves as non-balanced. The modeling of stochastic systems usually focuses on part of relevant variables, and uses “random force” or “noise” with specific statistical characteristics to describe eliminated ones, with a characteristic time scale matching with the specific correlation time of the random noise [19]. When the time scale magnitude order of these eliminated variables is much smaller than that of the system’s characteristic response time, Gaussian white noise (GWN) with a pulse form correlation time can be employed for modeling. There can be practical problems, as the noises often manifest themselves into colored form with non-zero correlation time and nonconstant power spectrum. The dichotomous colored noise, which has a simple statistical property, is originally used to simulate real noise. As an extension, the trichotomous one is more flexible in modeling than the former one. In most instances, the system induced by trichotomous noise may exhibit more numerous and complicated dynamic behaviors than that induced by the dichotomous one [20].
Here we use ζ ( t ) to denote the ubiquitous time-continuous Markov process, which consists of shifts between three different levels, i.e., ζ ( t ) may jump between three values Δ , 0 and Δ , with Δ > 0 . The jumps of this stationary random telegraph process follow, in time, according to the pattern of a Poisson process, and the stationary probability for the symmetric trichotomous noise ζ ( t ) at different values are
P s ( ζ ( t ) = Δ ) = P s ( ζ ( t ) = Δ ) = p ,   P s ( ζ ( t ) = 0 ) = 1 2 p
thus ζ ( t ) would be degenerated into the dichotomous form for p = 1 / 2 , the conditional transition probability between the three states is determined by the master equation [21]
d d t P Δ i , t | x , s = v j = 1 3 A i j P Δ j , t | x , s
where Δ 1 = Δ , Δ 2 = 0 , Δ 3 = Δ , t > s
A i j = p 1 1 2 p p p 2 p p p 1 2 p p 1
Combining the normalization condition j = 1 3 P Δ j , t | x , s = 1 and the initial conditions P Δ j , t | x , t = δ Δ j x , one can use Equation (2) to give the a probability transition matrix of noise in regard to the time interval t s , and the shift between the states Δ , 0 and Δ obey a pattern of Poisson distribution, which can be given as follows
P = P i j t s i , j 1 , 2 , 3 = P ζ t = Δ j | ζ s = Δ i i , j 1 , 2 , 3 = 1 + p 1 1 e v t s 1 2 p 1 e v t s p 1 e v t s p 1 e v t s 1 2 p 1 e v t s p 1 e v t s p 1 e v t s 1 2 p 1 e v t s 1 + p 1 1 e v t s
All the characters of ζ ( t ) are completely contained in (1) and (3), from which the mean value at arbitrary time t can be calculated as
ζ ( t ) | x 0 , s = x Δ , 0 , Δ x P x , t | x 0 , s = 0 ,   t > s
and the autocorrelation function
ζ t ζ s = x 0 , x 1 Δ , 0 , Δ x 0 x 1 P t , x 1 | s , x 0 P s x 0 = 2 p Δ 2 e v t s = D v e v t s
v is the switching rate, and correlation time τ c o r equals to the reciprocal of the noise switching rate
D = 0 ζ ( t ) ζ ( t + u ) d u = 2 p Δ 2 v ,   τ c o r = 1 D 0 + ζ ( t ) ζ ( t + u ) u d u = 1 v
D is the noise intensity. The mean square value and the counterpart mean quadruplicate value of the noise reads as ζ 2 ( t ) = 2 Δ 2 p , ζ 4 ( t ) = 2 Δ 4 p , thus one gets the flatness parameter [22] of ζ ( t ) : κ = ζ 4 ( t ) / ζ 2 ( t ) = 1 / 2 p , 0 < p < 0.5 . For different p , κ could be any value from one to , and the trichotomous is reduced to a more concise case, say, the symmetric dichotomous noise and Gaussian colored noise for κ = 2 and κ = 3 , respectively. For the other degree of freedom, it proves itself to be more flexible when modelling practical fluctuations [23].
Since the phenomenon of pollen particles’ never-ending movement in liquid was first observed by the British botanist Brown in 1827, the theory of normal diffusion based on integer ODE has been developed gradually. The model of single Brownian particle’s motion subjected to potential field force and random force in a homogeneous fluid media can be described by
m d 2 d t 2 x t + m γ d d t x t = d d x U x + ξ t + F t
wherein m γ x ˙ t , d U x / d x and F t represent the resistance, the potential field force, and the external force, respectively. This classical model is established on two aspects of idealized assumptions: firstly, the volume of Brownian particles is assumed to be much larger than that of fluid molecules. Since the time scale of fluid molecules’ motion is much smaller than that of Brownian particles, the association time of random forces ξ t is extremely short, say, modeled as GWN. The second assumption is that the damping coefficient γ 0 should be constant, i.e., damping force at the current moment has nothing to do with t motion history. Under this Markov hypothesis, the fluid damping of Brownian particles can be determined by Stokes law, and thus can be described by the first-order derivative.
However, in practical problems many diffusion phenomena exhibit quite different pattern, e.g., solute migration in fractured lithosphere [23], anomalous diffusion in the cytoplasm of mammalian cells [24], etc. These processes can no longer be described by the standard method of statistical physics, are generally referred to as anomalous diffusion. Indeed, the abovementioned assumptions are often too idealistic, that does not correspond to the actual situation. On one hand, the size of the object particle may not be much larger than that of the impacting particles in a complex surrounding environment with nonhomogeneous media. In this case, the correlation time of collisions is no longer infinitesimal, occurrence of low-frequency collisions leads to the correlation between random forces, then the Gaussian white noise model is no longer applicable. On the other hand, in order to satisfy the characteristic equation v t v ˙ t = 0 , the viscosity coefficient γ 0 should manifest as a time-dependent function t 0 t γ t t v t d t [25] rather than a constant coefficient, the damping force at any time t is acquired by integrating the velocity at all times before t according to the weight function γ . From the above consideration, it can be further turned into the generalized Langevin equation
m d 2 d t 2 x t + m t 0 t γ t t x ˙ t d t = d d x U x + W H t
where γ t is the memory damping kernel, describing delayed effect of the friction.
In many physical and biological environments, a power-law form memory kernel is often used to describe the dependence between viscosity force and particles’ historical velocity, as follows
γ t = γ 0 t q Γ 1 q
The shape of damping kernel function with the friction coefficient δ = 1 is shown in Figure 1, it finds that γ t has a power-law attenuation trend with the speed of q , and it would degenerate into Delta function when q approach 1 . Inserting (6) into (5) and using the definition of Caputo derivative, Equation (5) can be further transformed to the following FDE
m d 2 d t 2 x t + m γ 0 D 0 , t q C x ( t ) = d d x U x + W H t
The well-known Brownian motion B ( t ) is generally used in the traditional stochastic dynamics, to model the normal diffusion with mean-square displacement (MSD) x 2 ( t ) = t induced by the random factor accompanied with all sorts of systems, it has further been confirmed to be quite serviceable in many situations. However, it is clunky when dealing with some other ubiquitous situations, i.e., the Brownian motion theory does not apply to the anomalous diffusion process [26], such as the long-range correlation behavior in history-dependent viscoelastic medium under the background of anomalous diffusion, with a magnitude of MSD the same as t α , the diffusion exponents α lies between zero and one for a subdiffusion situation and beyond one for a superdiffusion case, respectively [27]. By this reason, a noise with a long time memory effect is frequently adopted to describe the non-Markov dynamic behavior of anomalous diffused particles, and the correlation function of W H t tends to be in a power-law form, as in Equation (7) [26]. According to the fluctuation dissipation theorem [25], the friction produced by the microscopic forces of the system depends entirely on the correlation of the random forces: W H ( t ) W H ( s ) = k B T γ ( t s ) , with k B the Boltzmann constant and T the absolute temperature. Substituting (6) into the fluctuation dissipation equality gives the power-law form autocorrelation function of random force
W H ( t ) W H ( s ) = k B T γ 0 t s q Γ 1 q
Here we consider W H t to be the fractional Gaussian noise (FGN), and similar to the definition pattern of GWN, FGN is given by the generalized differential of the fractional Brownian motion (FBM): W H t = 2 D H d B H t / d t . D H denotes the noise strength, H denotes the Hurst index with 0 < H < 1 . B H t is reduced to a standard Brownian motion for H = 0.5 , which has independent increments. For H > 0.5 the increments are positively correlated, while they are negatively correlated for the opposite situation, and more details of FBM can be found in Ref. [28]. The mean value and autocorrelation function of FBM read as:
B H t = 0 ,   B H t B H s = t 2 H + s 2 H t s 2 H / 2 ,  
where the correlation function of FGN can be obtained by the covariance of FBM
W H ( t ) W H ( s ) = 2 D H H ( 2 H 1 ) t s 2 H 2 + 2 H t s 2 H 1 δ t s
When t s , 0 < H < 1 / 2 corresponds to the subdiffusion movement with a negative correlation function, 1 / 2 < H < 1 corresponds the opposite situation, and it turns to the impulse function W H ( t ) W H ( s ) = 2 D H δ t s for H = 1 / 2 , say, the normal diffusion case. Utilizing the fluctuation dissipation theorem, again one gets
γ ( t ) = 2 D H H ( 2 H 1 ) t 2 H 2 / k B T
The Hurst parameter should be constrained in ( 0.5 , 1 ) to guarantee the mean-square value W H 2 ( t ) to be positive, thus the memory kernel is positive as well. Combining (10) and (8) one has q = 2 2 H , that is, the attenuation parameters of memory kernel function correspond exactly to the order of the fractional derivative. An intensity setting of FGN as D H = k B T γ 0 / Γ 2 H + 1 can provide a result consistent with (6), which reveals the practical significance of considering FGN in an anomalous diffusion system.
In this paper, we consider a kind of anomalous diffusion system with small time delay x ( t τ ) , disturbed by nonlinear colored noise and internal FGN, which is given by the following integro-differential equation:
m x ¨ ( t ) + m γ 0 + Ψ ζ ( t ) , N t 0 t γ t t x ˙ t d t + ω 0 2 g x ( t ) , x ( t τ ) = A cos Ω t + W H ( t )
x ¨ ( t ) and x ˙ t denote the acceleration and velocity, respectively. A cos Ω t is external harmonic excitation, ω 0 2 the inherent frequency. Ψ ζ ( t ) , N ( N 2 ) represents the nonlinear trichotomous noise disturbance on damping intensity, which manifests itself as the polynomial form Ψ ζ ( t ) , N = k = 1 N γ k ζ k ( t ) . ζ ( t ) is symmetric trichotomous colored noise with three states, W H ( t ) is zero-mean FGN with autocorrelation function given by (9) and memory damping kernel given by (6), respectively. Considering the potential field force g ( x ( t ) , x ( t τ ) ) = ( 1 ρ ) x ( t ) + ρ x ( t τ ) , 0 ρ 1 , where ρ represents delay intensity. ρ = 0 and ρ = 1 correspond, respectively, to non-delay and global item delay systems [17], 0 < ρ < 1 corresponds to the system containing delay feedback, wherein ρ is also referred as the feedback gain [29]. For computational convenience sake assume m = 1 , then based on the generalized Langevin equation given by Equation (11) and the definition of Caputo fractional derivative [30], the system can be rewritten as the following fractional Langevin equation
x ¨ ( t ) + γ 0 + k = 1 N γ k ζ k ( t ) D 0 , t q C x ( t ) + ω 0 2 ( 1 ρ ) x ( t ) + ρ x ( t τ ) = A cos Ω t + W H ( t )

3. Theoretical Analysis of Response

With the constructed equation ζ ( t ) ( ζ ( t ) Δ ) ( ζ ( t ) + Δ ) = 0 , and the arbitrary power of ζ ( t ) can be expressed by
ζ 2 m + 1 ( t ) = Δ 2 m ζ ( t ) ,   ζ 2 m + 2 ( t ) = Δ 2 m ζ 2 ( t )
m is a positive integer, with (13) the polynomial noise in Equation (12) can be simplified to
Ψ ( ζ ( t ) , N ) = α N ζ ( t ) + β N ζ 2 ( t ) ,   α N = j = 1 ( N 1 ) / 2 + 1 γ 2 j 1 Δ 2 j 2 ,   β N = j = 1 N / 2 γ 2 j Δ 2 j 2
where “ x ” denotes just the value x rounded down to the nearest integer. For the case of N = 1 , the effect of the linear Markovian symmetric trichotomous noise on stochastic resonance has been investigated [6]. From the expansion expression (13) it is clear that a N-order ( N > 2 ) polynomial trichotomous noise can be simplified to an equivalent quadratic polynomial one, i.e., the first two orders play a post role in the nonlinear trichotomous noise with arbitrary (adequate) high-level power, thus for convenience, this paper adopts the quadratic form given in (14), and all the other cases with high-power exponents can be expressed by the first two orders, thus can be discussed similarly. Moreover, for notational convenience, hereinafter the subscript “ N ” in α N and β N will be dropped, and the simplified versions α and β will be used to describe the linear coefficient and nonlinear coefficient of the polynomial function Ψ ( ζ ( t ) , N ) of the trichotomous colored noise, respectively.
As for the function g ( x ( t ) , x ( t τ ) ) of the small time delay, when ρ = 0 no delay appeared, and this mainstream situation has attracted widespread attention. Meanwhile the global delay case happens for ρ = 1 and the stiffness term disappears [15,17]. The scenario of 0 < ρ < 1 refers to the time-delayed feedback in the system [29], the small time-delay pre-assumption enables us to expand the small delay approximation method (SDAM) [31] to the current fractional system case, say, g ( x ( t ) , x ( t τ ) ) and can be expressed as a Taylor expansion series with respect to the time delay around x ( t ) , with the convergence precision the same order of magnitude as O ( τ 2 ) . In fact, such a Taylor expansion approximation of a function with respect to x ( t ) has been proved workable when dealing with a small delay term, no matter whether in an integer-order differential system [32] or a fractional-order system [33]. Making use of this approximation method, one gets:
ω 0 2 g ( x ( t ) , x ( t τ ) ) ω ˜ 0 2 x ( t )
where ω ˜ 0 2 = ω 0 2 1 τ ρ ω 0 2 .
Combining Equation (15) and, as discussed below, Equation (14), we obtain a non-delayed stochastic differential equation, which is equivalent to Equation (12):
d 2 d t 2 x ( t ) + γ 0 + α ζ ( t ) + β ζ 2 ( t ) D q x ( t ) + ω ˜ 0 2 x ( t ) = A cos ( Ω t ) + η H ( t )
As mentioned above in Introduction section, there are many alternative assessment indexes to quantify the SR phenomenon in stochastic systems. Among them, the amplitude gain of the output signal is often considered as the one that could be utilized the most, from an analytical point of view. When describing the generalized SR and bona fide one, naturally many researchers have adopted this index in dealing with SR problems for its simplicity and effectiveness in practice. In this paper an analytical expression of the amplitude gain of the output signal in regard to the system given by Equation (12) will be discussed in detail, however, before doing this, we first calculate the exact solution of the first-order moment of the fractional system. This operation is feasible since the displacement process x ( t ) could be always stationary [33]. In fact, the fractional particle in Equation (12) will always be pulled back to the origin sooner or later due to the existence of the harmonic potential V ( x ) = ( ω ˜ 0 2 x 2 ) / 2 . For this purpose, the Shapiro–Loginov formula [8] should be taken account for the exponential correlative trichotomous process:
d d t ζ ( t ) x ( t ) = ζ ( t ) d d t x ( t ) v ζ ( t ) x ( t ) ,  
and the fractional Shapiro-Loginov formula:
ζ ( t ) D q x ( t ) = e v t D q ζ ( t ) x ( t ) e v t
Here, x ( t ) is considered as a function of ζ ( t ) . Taking average of Equation (16), one gets
d 2 d t 2 + ω ˜ 0 2 + γ 0 D q x ( t ) + α ζ ( t ) D q x ( t ) + β ζ 2 ( t ) D q x ( t ) = A cos ( Ω t )
The prerequisite condition η H ( t ) = 0 has been used in Equation (18), considering the properties of the trichotomous noise and using the transition probabilities and stationary probabilities, the mean-square value and the correlation function of ζ 2 ( t ) are given by ζ 2 ( t ) = 2 Δ 2 p and
ζ 2 ( t ) ζ 2 ( t + Δ t ) = z 1 { Δ , Δ } z 2 { Δ , Δ } Δ 4 P ( ζ ( t + τ ) = z 2 | ζ ( t ) = z 1 ) P s ( ζ ( t ) = z 1 ) = 4 p Δ 4 + 2 p 2 Δ 4 ( 1 2 p ) e v t
It is apparent from the above autocorrelation function of the mean value ζ 2 ( t ) that ζ 2 ( t ) does not satisfy the conditions of the Shapiro–Loginov formula, however, the situation becomes gratifying for the adjustive variable ζ 2 ( t ) 2 Δ 2 p , that
ζ 2 ( t ) 2 Δ 2 p = 0 ζ 2 ( t ) 2 Δ 2 p ζ 2 ( t + t ) 2 Δ 2 p = 2 p 2 Δ 4 ( 1 2 p ) e v t
Thus, the random variable ζ 2 ( t ) 2 Δ 2 p is applicable for various kinds of Shapiro–Loginov formulas (of course also including the fractional version (3.5))
( ζ 2 ( t ) 2 Δ 2 p ) D q x ( t ) = e v t D q ( ζ 2 ( t ) 2 Δ 2 p ) x ( t ) e v t = ζ 2 ( t ) D q x ( t ) 2 Δ 2 p D q x ( t )
(19) leads to the new correlation factor ζ 2 ( t ) D q x ( t ) appear in (18):
ζ 2 ( t ) D q x ( t ) = 2 p Δ 2 D q x ( t ) + e v t D q ζ 2 ( t ) x ( t ) 2 Δ 2 p x ( t ) e v t
By substituting Equations (17) and (20) into (18), one gets the equation
d 2 d t 2 + ω ˜ 0 2 x + γ 0 + 2 p Δ 2 β D q x + α e v t D q ζ x e v t + β e v t D q ζ 2 x e v t 2 p Δ 2 β e v t D q x e v t = A cos ( Ω t )
To solve the new factor ζ ( t ) x ( t ) , multiplying Equation (16) by ζ ( t ) and averaging over the gained equation, one arrives at
ζ d 2 d t 2 x + γ 0 ζ D q x + α ζ 2 D q x + β ζ 3 D q x + ω ˜ 0 2 ζ x = 0
in Equation (22) the uncorrelated postulate between ζ t and the FGN η H t has already been adopted, i.e., ζ k ( t ) η H s = 0 , k = 1 , 2 Reusing the Shapiro–Loginov formula, we have
ζ ( t ) d 2 d t 2 x ( t ) = d 2 d t 2 ζ ( t ) x ( t ) + 2 v d d t ζ ( t ) x ( t ) + v 2 ζ ( t ) x ( t )
Substituting (17), (20) and (23) into Equation (22) yields
d d t + v 2 + ω ˜ 0 2 ζ x + 2 p Δ 2 α D q x + γ 0 + β Δ 2 e v t D q ζ x e v t + α e v t D q ζ 2 x e v t 2 p Δ 2 α e v t D q x e v t = 0
In order to solve the term ζ 2 x , multiplying Equation (16) by ζ 2 ( t ) and the averaging the gained results, one obtains:
ζ 2 d 2 d t 2 x + γ 0 ζ 2 D q x + α ζ 3 D q x + β ζ 4 D q x + ω ˜ 0 2 ζ 2 x = 2 p Δ 2 A cos ( Ω t )
the first averaged term in Equation (25) can also be derived by (19) and (23)
ζ 2 2 p Δ 2 d 2 d t 2 x = d 2 d t 2 ζ 2 2 p Δ 2 x + 2 v d d t ζ 2 2 p Δ 2 x + v 2 ζ 2 2 p Δ 2 x = ζ 2 d 2 d t 2 x 2 p Δ 2 d 2 d t 2 x
From (26) the averaging value ζ 2 d 2 d t 2 x can be obtained as the expression
ζ 2 d 2 d t 2 x = d 2 d t 2 ζ 2 x + 2 v d d t ζ 2 x + v 2 ζ 2 x 4 p Δ 2 v d d t x 2 p Δ 2 v 2 x
Inserting (17), (20) and (27) into Equation (25), meanwhile, replacing ζ 3 ( t ) , ζ 4 ( t ) for Δ 2 ζ ( t ) and Δ 2 ζ 2 ( t ) , respectively, then (27) can be transformed to
d d t + v 2 + ω ˜ 0 2 ζ 2 x 4 p Δ 2 v d d t x 2 p Δ 2 v 2 x + 2 p Δ 2 γ 0 + 2 p Δ 4 β D q x + α Δ 2 e v t D q ζ x e v t + γ 0 + β Δ 2 e v t D q ζ 2 x e v t 2 p Δ 2 γ 0 + 2 p Δ 4 β e v t D q x e v t = 2 p Δ 2 A cos ( Ω t )
the harmonic terms on the right side of Equation (28) would not further be eliminated unless recurring to Equation (21), in fact, multiplying Equation (21) by 2 p Δ 2 and subtracting the gained equation from Equation (28), one can obtain:
d d t + v 2 + ω ˜ 0 2 ζ 2 x + γ 0 + β Δ 2 2 p Δ 2 β e v t D q ζ 2 x e v t + α Δ 2 1 2 p e v t D q ζ x e v t + 2 p Δ 2 β Δ 2 2 p 1 γ 0 e v t D q x e v t + 2 p Δ 4 β 1 2 p D q x 2 p Δ 2 v 2 + ω ˜ 0 2 x 4 p Δ 2 v d d t x 2 p Δ 2 d 2 d t 2 x = 0
The collection of (21), (24) and (29) truly consists of a close linear system of second-order differential equations for three variables m 1 ( t ) = x ( t ) , m 2 ( t ) = ζ ( t ) x ( t ) , m 3 ( t ) = ζ 2 ( t ) x ( t ) . By virtue of the Laplace transform technique, this equation set can be equivalently changed to a corresponding linear algebraic equation set, given as follows
a 11 M 1 + a 12 M 2 + a 13 M 3 = A s / s 2 + Ω 2 + m 1 ( 0 ) + b 11 m 1 ( 0 ) + b 12 m 2 ( 0 ) + b 13 m 3 ( 0 ) a 21 M 1 + a 22 M 2 + a 23 M 3 = m 2 ( 0 ) + b 21 m 1 ( 0 ) + b 22 m 2 ( 0 ) + b 23 m 3 ( 0 ) a 31 M 1 + a 32 M 2 + a 33 M 3 = m 3 ( 0 ) 2 p Δ 2 m 1 ( 0 ) + b 31 m 1 ( 0 ) + b 32 m 2 ( 0 ) + b 33 m 3 ( 0 )
a 11 = s 2 + ω ˜ 0 2 + ( γ 0 + 2 p Δ 2 β ) s q 2 p Δ 2 β ( s + v ) q ,   a 12 = α ( s + v ) q ,   a 13 = β ( s + v ) q ,   a 21 = 2 p Δ 2 α s q ( s + v ) q ,   a 22 = ( s + v ) 2 + ω ˜ 0 2 + ( γ 0 + β Δ 2 ) ( s + v ) q ,   a 23 = a 12 ,   a 31 = 2 p Δ 2 ω ˜ 0 2 + s + v 2 + γ 0 s + v q + β Δ 2 1 2 p s + v q s q a 32 = Δ 2 ( 1 2 p ) α ( s + v ) q ,   a 33 = ( s + v ) 2 + ω ˜ 0 2 + ( γ 0 + β Δ 2 2 p Δ 2 β ) ( s + v ) q ,   b 11 = s + ( γ 0 + 2 p Δ 2 β ) s q 1 2 p Δ 2 β ( s + v ) q 1 ,   b 12 = α ( s + v ) q 1 ,   b 13 = β ( s + v ) q 1 ,   b 21 = 2 p Δ 2 α s q 1 ( s + v ) q 1 ,   b 22 = s + 2 v + ( γ 0 + β Δ 2 ) ( s + v ) q 1 ,   b 23 = b 12 ,   b 31 = 2 p Δ 2 β Δ 2 2 p 1 γ 0 s + v q 1 + 2 p β Δ 4 1 2 p s q 1 4 p Δ 2 v 2 p Δ 2 s ,   b 32 = Δ 2 ( 1 2 p ) α ( s + v ) q 1 ,   b 33 = s + 2 v + ( γ 0 + β Δ 2 2 p Δ 2 β ) ( s + v ) q 1 . M k ( s ) = m k ( t ) = 0 m k ( t ) e s t d t ,   k 1 , 2 , 3 .  
Displacement’s first-order moment can be formally expressed formally as bellow:
M 1 ( s ) = H ˜ 10 ( s ) A s s 2 + Ω 2 + k = 1 3 H ˜ 1 k ( s ) m k ( 0 ) + k = 1 3 G ˜ 1 k ( s ) m k ( 0 )
H ˜ 10 ( s ) = a 22 a 33 a 23 a 32 / I s ,   H ˜ 11 ( s ) = b 11 a 22 a 33 b 11 a 23 a 32 + a 13 b 21 a 32 a 12 b 21 a 33 + a 12 a 23 b 31 a 13 a 22 b 31 / I s H ˜ 12 ( s ) = b 12 a 22 a 33 b 12 a 23 a 32 + a 13 b 22 a 32 a 12 b 22 a 33 + a 12 a 23 b 32 a 13 a 22 b 32 / I s H ˜ 13 ( s ) = b 13 a 22 a 33 b 13 a 23 a 32 + a 13 b 23 a 32 a 12 b 23 a 33 + a 12 a 23 b 33 a 13 a 22 b 33 / I s G ˜ 11 ( s ) = a 22 a 33 a 23 a 32 + 2 p Δ 2 ( a 13 a 22 a 12 a 23 ) / I ( s ) G ˜ 12 ( s ) = a 13 a 32 a 12 a 33 / I s ,   G ˜ 13 ( s ) = a 12 a 23 a 13 a 22 / I s I ( s ) = a 11 a 22 a 33 + a 12 a 23 a 31 + a 13 a 21 a 32 a 13 a 22 a 31 a 11 a 23 a 32 a 12 a 21 a 33
Applying the inverse Laplace transformation technique on Equation (31) and based on the corresponding uniqueness condition [34], the averaging mean value (the first-order moment) of the displacement of Equation (12) can be represented in the following form:
x ( t ) = A H 10 ( t ) cos ( Ω t ) + k = 1 3 H 1 k ( t ) m k ( 0 ) + k = 1 3 G 1 k ( t ) m k ( 0 )
” the convolution operation, H 10 ( t ) = 1 H ˜ 10 ( s ) , H 1 k ( t ) = 1 H ˜ 1 k ( s ) , G 1 k ( t ) = 1 G ˜ 1 k ( s ) , with the inverse Laplace operator the primitive function y ( t ) can be calculated by the integral of its corresponding Laplace transformation function Y ( s ) : y ( t ) = 1 Y ( s ) = 1 2 π i ε i ε + i Y ( s ) e s t d s .
Based on the stability theory [35], the stability of the solution expressed as (33) can be guaranteed under the circumstances that all the possible solutions of the equation I ( s ) = 0 have roots with a negative or zero real part. When considering the asymptotic behavior of the functions on the right hand of (33), i.e., at a long-time limit t , by making use of the Tauberian theorem [36], one gets the estimation: H 10 ( t ) = O ( t q ) , H 1 k ( t ) , G 1 k ( t ) = O ( t ( 1 + q ) ) , thus the memorizing effect of the initial conditions can be disregarded when the system has a long enough evolute, i.e., the stationary solution of the first-order moment of the displacement in a system described by Equation (12) can be obtained by (33).
x ( t ) a s = A 0 t H 10 ( t t ) cos ( Ω t ) d t
the complex susceptibility of the stochastic system given by Equation (12) can be given based on the form of (34) [37]: ( Ω ) = 0 H 10 ( t ) e i Ω t d t = H ( i Ω ) = 1 ( Ω ) + 2 ( Ω ) i , where i denotes the imaginary unit with i 2 = 1 , 1 and 2 representing the real and imaginary parts of the complex susceptibility, respectively. Moreover, with the help of the linear response theory [38], (34) can be rewritten in the following harmonic pattern as well x ( t ) a s = F cos ( Ω t + ϕ ) . With the output amplitude F = A / and the phase shift given as, ϕ = arctan ( 2 / 1 ) , it is worth nothing that H ˜ 10 ( i Ω ) = H ˜ 10 ( i Ω ) ¯ = ¯ = 1 2 i , thus F and ϕ can also be calculated by F = A H ˜ 10 ( i Ω ) , ϕ = arg ( H ˜ 10 ( i Ω ) ) . The response amplitude gain, say, the ratio of the amplitude of the output signal to the amplitude of input signal, then can be expressed as
R = F A = H ˜ 10 ( i Ω )
Through the expressions of the transfer function given in (31) and the coefficients below (29), after a series of simplification calculations, one can obtain
R = e 1 2 + e 2 2 e 3 2 + e 4 2 ,   ϕ = arctan ( e 2 e 3 e 1 e 4 e 1 e 3 + e 2 e 4 ) ,
where e k , k = 1 , , 4 are listed as follows:
e 1 = ω ˜ 0 4 + B 4 cos ( 4 θ ) + 2 ω ˜ 0 2 B 2 cos ( 2 θ ) + ω ˜ 0 2 J 1 B q cos ( q θ ) + J 1 B q + 2 cos q + 2 θ + J 2 B 2 q cos ( 2 q θ ) e 2 = B 4 sin ( 4 θ ) + 2 ω ˜ 0 2 B 2 sin ( 2 θ ) + ω ˜ 0 2 J 1 B q sin ( q θ ) + J 1 B q + 2 sin q + 2 θ + J 2 B 2 q sin ( 2 q θ ) + J 2 J 3 + ω ˜ 0 2 J 5 α 2 β J 1 J 6 β J 7 B 2 q cos ( 2 q θ ) + J 5 γ 0 α 2 + J 6 J 7 β J 2 B 3 q cos ( 3 q θ ) + J 1 J 3 B q + 2 cos q + 2 θ + J 1 J 4 α 2 Ω q J 5 β J 5 J 8 B q + 2 cos q + 2 θ + π 2 q ω ˜ 0 2 J 1 J 4 α 2 Ω q J 5 β J 5 J 8 B q cos θ + π 2 q + J 5 α 2 β J 1 J 6 β J 7 B 2 q + 2 cos 2 q + 2 θ + J 2 J 4 γ 0 α 2 Ω q J 5 + J 5 J 6 J 8 B 2 q cos 2 θ + π 2 q + 2 ω ˜ 0 2 J 4 B 2 cos ( 2 θ + π 2 q ) + J 4 B 4 cos ( 4 θ + π 2 q ) e 4 = 2 ω ˜ 0 2 B 2 J 3 sin ( 2 θ ) + B 4 J 3 sin ( 4 θ ) + ω ˜ 0 2 J 1 J 3 B q sin ( q θ ) + ω ˜ 0 4 J 4 sin ( π 2 q ) + J 2 J 3 + ω ˜ 0 2 J 5 α 2 β J 1 J 6 β J 7 B 2 q sin ( 2 q θ ) + J 5 γ 0 α 2 + J 6 J 7 β J 2 B 3 q sin ( 3 q θ ) + J 1 J 3 B q + 2 sin q + 2 θ + J 1 J 4 α 2 Ω q J 5 β J 5 J 8 B q + 2 sin q + 2 θ + π 2 q ω ˜ 0 2 J 1 J 4 α 2 Ω q J 5 β J 5 J 8 B q sin θ + π 2 q + J 5 α 2 β J 1 J 6 β J 7 B 2 q + 2 sin 2 q + 2 θ + J 2 J 4 γ 0 α 2 Ω q J 5 + J 5 J 6 J 8 B 2 q sin 2 θ + π 2 q + 2 ω ˜ 0 2 J 4 B 2 sin ( 2 θ + π 2 q ) + J 4 B 4 sin ( 4 θ + π 2 q ) B = v 2 + Ω 2 ,   θ = arctan ( Ω v ) , J 1 = 2 γ 0 β Δ 2 p β Δ 2 J 2 = γ 0 + β Δ 2 γ 0 + β Δ 2 2 p β Δ 2 Δ 2 1 2 p α 2 , J 3 = ω ˜ 0 2 Ω 2 ,   J 4 = γ 0 + 2 p Δ 2 β Ω q ,   J 5 = 2 p Δ 2 ,   J 6 = α 2 γ 0 β β 2 Δ 2 ,   J 7 = 2 p β Δ 2 β Δ 2 γ 0 ,   J 8 = β Δ 2 1 2 p Ω q .

4. Bona Fide Stochastic Resonance

Stochastic resonance phenomenon has attracted much attention in the past 40 years, and many achievements have been made in physics, biology, chemistry and engineering fields. The measurement methods used to determine the degree of stochastic resonance vary in different fields, and are usually manifest in three forms: signal-to-noise ratio (SNR) [39], which is defined as the ratio of the power spectrum at the corresponding signal frequency to the average power spectrum of background noise. However, this characterization method is no longer applicable in the case of aperiodic signal excitation systems. The second is the residence time characterization method [40], which measures stochastic resonance in a symmetric bistable system by calculating the average residence time of particles in a potential well of the system. The third one is output amplitude gain (OAG), which was proposed by Gitterman [41]. Traditional stochastic resonance refers to the optimal response of the system at a critical noise intensity when the system is subjected to both external and random forces. In recent years the concept has been expanded: SR in broad sense [41] proposed by Gitterman refers to the non-monotonic dependence of some functions of the output signal (such as SNR, autocorrelation function, amplitude gain, etc.) on noise or system parameters, showing that traditional SR is actually included in the concept of generalized stochastic resonance (GSR). Another extended concept is bona fide SR, which indicates the peak phenomenon of output signal amplitude induced by the fluctuation of signal frequency [42]. In this section, based on the theoretical output amplitude gain results in the previous section, different types of stochastic resonance cases will be analyzed for the fractional perturbation system with small hysteresis.
Before discussing the possible stochastic resonance phenomenon in the system (12), we first examine the case of no polynomial trichotomous noise disturbance, i.e., α = β = 0 or Δ = 0 , then the corresponding system equation reads
x ¨ ( t ) + γ 0 D q x ( t ) + ω 0 2 ( 1 ρ ) x ( t ) + ρ x ( t τ ) = A cos Ω t + W H ( t )
using the small time delay approximation method (15), Equation (36) is further transformed into a delay-free equation, taking the average of both sides of the obtained equation, one has
d 2 d 2 t x ( t ) + γ 0 D q x ( t ) + ω 0 ˜ 2 x t = A cos Ω t
after the Laplace transformation of Equation (37) with respect to the variable x ( t ) , the mean value of steady-state for t can be calculated as
x t a s = A cos Ω t H t = F cos Ω t + ϕ
A and ϕ denote the corresponding amplitude and phase shift, respectively. H t = L 1 H ˜ s ; t with
H ˜ s = 1 / s 2 + γ 0 s q + ω 0 ˜ 2
for any fractional order satisfying 0 < q < 1 , the steady-state output amplitude of the system (37) can be obtained as
F = A H ˜ i Ω = A ω 0 ˜ 2 Ω 2 + γ 0 Ω q cos q π / 2 2 + γ 0 Ω q sin q π / 2 2
In fact, according to the linear response theory, if one substitutes x t = F cos Ω t + ϕ into (37) directly, a same result as (38) can also be obtained by utilizing the undetermined coefficient method. In addition, the same situation happens when inserting α = β = 0 into the expression (35), showing that system (36) is actually a special case of system (12).
Figure 2 shows the bona fide SR phenomenon in system (36) that was induced by external excitation frequency Ω under different flatness parameter q , with system parameters A = 1 , ρ = 0.5 , γ 0 = 1 , ω 0 = 1 , τ = 0.2 . It finds that the resonance of the system output signal is striking when q is small, while for the opposite situation F changes slowly despite the occurrence of peak Figure 2b describes the combined effect of external excitation frequency and flatness parameters on OAG of the system. With the increase in q , the single-resonance point Ω S R gradually increases and F tends to be flat. Indeed, the critical resonant frequency Ω S R must satisfy
d F Ω d Ω = 0
On the one hand, expression (38) implies F A / ω 0 ˜ 2 Ω 2 + γ 0 when the Hurst index is close to one. According to (39), the critical SR frequency of the system can be obtained as Ω S R ω 0 ˜ 2 + γ 0 , According to the form of ω 0 ˜ 2 , single-peak SR will certainly emerge in the system for small delay, and F max when q approaching to 0 . Therefore, in Figure 2b F exhibits a saltation within a certain frequency range for q 0 . On the other hand, consider the case H 1 / 2 , then substituting F A / [ ω 0 ˜ 2 Ω 2 ] 2 + γ 0 2 Ω 2 into Equation (39) give rise to the resonance peak Ω S R ω 0 ˜ 2 γ 0 2 / 2 , with the occurrence condition for bona fide SR γ 0 < 2 ω 0 ˜ and the optimal OAG F max 2 A / γ 0 4 ω 0 ˜ 2 γ 0 2 . Figure 2a shows the critical frequency values of the peak points q = 0.05 and q = 0.95 are 1.38 and 0.67, respectively, which are consistent with the above analysis results.
When one considers the stochastic polynomial trichotomous noise Φ ( ζ ( t ) , N ) in the system, the critical result of analytical expression (35) obtained in Section 2 can be utilized to give the theoretical analysis of SR. Let α = 0.4 , β = 2 , A = 1 , ρ = 0.5 , γ 0 = 1 , ω 0 = 1 , Δ = 1 , p = 0.3 , τ = 0.2 , Figure 3 shows the bona fide SR in the disturbed system (12) with a small flatness parameter q = 0.05 under different noise transfer rates, three-peak resonance, double-peak resonance and single-peak resonance appear at v = 0.01 , v = 0.2 and v = 3 , respectively. It is called stochastic multiresonance (SMR) if SR manifests as more than one peak. Thus in Figure 3 one can say that the external excitation frequency induces bona fide SMR. The noise transfer rate is sometimes called the correlation rate, which reflects the memory effect of noise. A small correlation rate means a long correlation time, and the switching between different states of noise will be extremely slow, leading to excessive correlation time results in the resonance peak at three different frequencies. In fact, in this case a state transition of the colored noise will be accompanied by a long-time stagnation, during which the damping coefficient of the system is actually fixed. Therefore, the random damping disturbance coefficient of the original system (12) can be divided into three cases in a long time.
γ 0 + Ψ ( ζ ( t ) , N ) = γ 0 γ 0 α Δ + β Δ 2 γ 0 + α Δ + β Δ 2
Replacing (40) for γ 0 in the results without damping disturbance (38), and using resonance generation condition (39), it can be estimated that the positions of the three peak points that may appear in the slow-switching-rate system for Hurst index approaching to one (Figure 3a) are
Ω S R 1 ω 0 ˜ 2 + γ 0 ,   Ω S R 2 ω 0 ˜ 2 + γ 0 α Δ + β Δ 2 ,   Ω S R 3 ω 0 ˜ 2 + γ 0 + α Δ + β Δ 2
When H is close to 1 / 2 , the external excitation frequency corresponding to the three resonance extreme points can also be calculated in the same way, given as
Ω S R 1 2 ω 0 ˜ 2 γ 0 2 / 2 ,   Ω S R 2 ω 0 ˜ 2 γ 0 α Δ + β Δ 2 2 / 2 Ω S R 3 ω 0 ˜ 2 γ 0 + α Δ + β Δ 2 2 / 2
As v increases, the noise switching between the three states becomes more and more frequent, resulting in the three peak points gradually converging and eventually merging into one overlapping peak point. The influence of the transfer rate on the bona fide SR peak is given in Figure 3d, it can be seen that when v keep growing, the unique resonance peak of the system will also increase gradually. The system presents at least a single-peak stochastic resonance phenomenon within the range of the combined parameter pairs Ω , v . It should be pointed out that all the system parameters considered in this section meet the stability condition I 0 > 0 .
The output amplitude of the system can also be approximately estimated for the fast-switching noise scenarios, as
F A ω 0 ˜ 2 Ω 2 + γ 0 + 2 p Δ 2 β Ω q cos π q / 2 2 + γ 0 + 2 p Δ 2 β Ω q sin π q / 2 2
Through calculation, it is not difficult to find that (41) exactly corresponds to the amplitude of the steady-state solution of the noiseless system (12)
d 2 d t 2 x ( t ) + γ 1 D q x ( t ) + ω 0 ˜ 2 x t = A cos Ω t
with the equivalent damping coefficient γ 1 = γ 0 + 2 p Δ 2 β . For the case of H 1 ( q 0 ), by solving d F Ω / d Ω = 0 in Equation (41), one has Ω S R = ω 0 ˜ 2 + γ 0 + 2 p Δ 2 β , the single-peak-resonance always exists in the system, which is consistent with the single peak point position Ω c = 1.76 in Figure 3c. For q 1 , only single-resonance may occur at Ω S R = ω 0 ˜ 2 γ 0 + 2 p Δ 2 β 2 / 2 , and the existence condition is
2 ω 0 ˜ 2 γ 0 + 2 p Δ 2 β 2 > 0
In the case of the slow-switching case v = 0.01 , and the bona fide SR phenomenon under different flatness parameters is shown in Figure 4. For q = 0.01 , the bona fide SR occurs at three critical frequencies; 1.38 , 1.87 and 2.08 , as q increases, triple resonant peaks gradually converge to double and single peaks, which is similar to the situation in Figure 3. The difference is that F decreases with the increase in q . When q is large enough, the resonance phenomenon disappears finally. Consider the other case q 1 (we might as well set q = 0.99 ), α = 2 , v = 3 , Figure 5 shows SR phenomena under different polynomial noise coefficients. Since the condition (43) is not met for β = 2 , the variation in Ω cannot induce bona fide SR. In fact, due to the effect of polynomial nonlinear multiplicative trichotomous noise, the resonance conditions are determined by many parameters. When nonlinear noise coefficient and damping coefficient were taken as control parameters and other parameter values fixed in (43), the parameter distribution of the occurrence of bona fide SR under the influence of combined parameters was given in Figure 5b. The critical value of γ 0 corresponding to the resonant peak is almost monotonically dependent on the critical value of β . When γ 0 > 1.82 , any β will no longer meet the condition (43), and the external excitation frequency of the system will no longer induce the bona fide SR phenomenon.
Figure 6 shows the influence of noise intensity on bona fide SR of the system with respect to the parametric region q v , under the polynomial trichotomous noise coefficient setting α = 0.4 , β = 2 . Considering the system parameters ρ = 0.5 , γ 0 = 1 , ω 0 = 1 , p = 0.3 , τ = 0.2 , it can be seen from Equation (38) that the system always exhibits monostable bona fide SR phenomenon when the trichotomous noise is absent. By comparing Figure 6a–c, it is apparent that appropriate noise intensity is conducive to the occurrence of SMR. When the noise amplitude is small (Figure 6b), only single-resonance appears in the highest range of parameter values (light blue area in the figure), with non-resonance appearing around 0.88 . When q , v are close to zero, a small range of two-peak resonance (cyan) and three-peak resonance (green) regions appear, indicating that thedditionn of trichotomous colored noise has enriched the stochastic resonance phenomenon of the system. When Δ increases to one, the proportion of the double-resonance region and triple-resonance region increase significantly, which reflects the fact that the SMR phenomenon of the system can be enhanced by the appropriate noise amplitude. However, it is not the case that the stronger the noise is, the stronger the resonance of the system will be: excessive noise intensity (Figure 6c) results in the disappearance of SMR in the system. Moreover, excessive noise intensity will even lead to the non-stochastic resonance (white) with a wider range of parameters than the noise-free case.

5. Parameter Induced GSR

Firstly, the influence of time delay τ is investigated. In the case of a system with deterministic fractional damping (42), the peak point τ S R of the system output amplitude corresponds to the root of the equation
ω 0 2 1 τ ρ ω 0 2 Ω 2 + γ 0 Ω q cos q π / 2 = 0
solving (44) and giving rise to the critical delay value
τ S R = ω 0 2 Ω 2 + γ 0 Ω q cos π q / 2 ρ ω 0 4
Inserting (45) into Equation (38), the optimal output amplitude at the critical delay τ S R is given by F = 1 / γ 0 Ω q sin q π / 2 , which can also be referred to as the time-delay induced intrinsic resonance peak occasionally. The dependence of the output amplitude of the system without trichotomous noise disturbance on the time delay is depicted in Figure 7. The single-peak of GSR appears at τ = 0.64 with the maximum output amplitude F max = 12.8 , which is consistent with the critical delay analysis result in (45). Examining the disturbed situation α = 0.5 , β = 1.2 , p = 0.35 , Δ 2 = 0.6 to be a comparison, based on the prediction result (35) of system (12), Figure 8 shows the relationship between system output amplitude and delay, and other system parameters read A = 1 , ρ = 0.8 , γ 0 = 0.5 , ω 0 = 1 , q = 0.1 , Ω = 0.99 . For a different noise correlation rate and noise intensity, variation in delay always causes GSR in the system, and can also induce the double- and triple- SMR phenomenon. When v = 0.01 , the damping coefficient of the system can be regarded as undisturbed over long periods, as the transfer rate of the trichotomous noise is very slow. By substituting the three cases of (40) into (45), one can obtain the possible critical time-delay of the resonance peak corresponding to ζ t = Δ , ζ t = 0 and ζ t = Δ .
τ S R 1 = τ S R ,   τ S R 2 = ω 0 2 Ω 2 + γ 1 Ω q cos π q / 2 ρ ω 0 4 ,   τ S R 3 = ω 0 2 Ω 2 + γ 2 Ω q cos π q / 2 ρ ω 0 4
with γ 1 = γ 0 α Δ + β Δ 2 , γ 2 = γ 0 + α Δ + β Δ 2 . Substituting the parameters into (46), the result is consistent with the data annotation value in Figure 8. Triple-peak GSR is converted to a double-peak GSR with the increase in v . When v is large enough, the critical delay of the single-peak GSR under the fast-switching noise case can be calculated as follows:
τ S R = ω 0 2 Ω 2 + γ 0 + 2 p Δ 2 β Ω q cos π q / 2 ρ ω 0 4
The fact that (47) is consistent with the peak position in Figure 8a. Figure 8b shows that despite the long correlation time of noise, the three equivalent damping coefficients corresponding to (45) will be very close due to the low noise intensity ( Δ 2 = 0.01 ) in the system, so that the inherent resonance peak could only exist at the critical delay given by (45). When Δ 2 = 0.3 , γ 0 is very close to γ 1 0.58 , while there is a large gap between γ 0 and γ 2 1.14 . Therefore, the peak position of GSR should be τ S R 1 and τ S R 3 in (46), and it is not difficult to get that to τ S R 1 = 0.642 and τ S R 3 = 1.453 , which is consistent with the resonance points of numerical results labeled in the figure. When Δ 2 = 0.6 , the intensity of the slow switching noise is large enough, and the three equivalent damping coefficients γ 0 , γ 1 and γ 2 , differ greatly from each other. At this time, three optimal peaks appear, respectively, at the three positions in (46), thus the system delay induces the three-peak SMR of the output amplitude. Based on the above analysis, it can be determined that when the system time delay is the only control parameter, a small noise transfer rate and an appropriate noise intensity would make the occurrence of multi-peak GSR phenomena more possible. In order to analyze the influence of parameters v and Δ more carefully, different fractional orders are selected to give the joint dependence between GSR induced by delay and parameter v , Δ 2 . It is found that SMR tends to appear for small q (cyan and green areas in Figure 9a), the single-peak GSR region becomes larger and larger as q increases, while the phenomenon of multi-peak GSR gradually disappears.
The aforementioned results confirmed that the noise intensity has a significant influence on both the bona fide SR induced by external excitation frequency (Figure 6) and GSR induced by time delay (Figure 8b). As such, we examined the change in the system output amplitude when Δ 2 varied with the other system parameters that were fixed as A = 1 , γ 0 = 0.5 , α = 1 , ω 0 = 1 , p = 0.35 , q = 0.5 , v = 0.01 , Ω = 0.99 , τ = 0.2 . As shown in Figure 10, variation in noise intensity leads to peak phenomenon of system output amplitude under the slow noise switching rate condition of v = 0.01 , analogous to the traditional form of SR that occurs in the system under the measurement index SNR. Figure 10a shows that the noise-induced stochastic resonance phenomenon presents different forms under the different nonlinear trichotomous noise factor β : One unique peak F max = 3.938 exists if the nonlinear part of colored noise is not considered ( β = 0 ), and the single-peak SR only occurs when the trichotomous noise switches to ζ t = Δ . Indeed, by inserting β = 0 and (40) into Equation (38), one can obtain a real solution only if the equivalent damping coefficient is equal to γ 2
Δ S R = γ 0 Ω q + ω 0 ˜ 2 Ω 2 cos q π / 2 α Ω q
Two amplitude spikes appear in the system when β increases to 0.5 . Similarly, by virtue of (38) and (40), it can be calculated that the system output amplitude may have a peak only when the trichotomous noise switches to Δ , and all possible extremum points are
Δ S R 1 = α Ω q κ 2 β Ω q ,   Δ S R 2 = α Ω q + κ 2 β Ω q ,   Δ S R 3 = α 2 β
With κ = α 2 4 β γ 0 Ω 2 q + 4 β Ω q Ω 2 ω 0 ˜ 2 cos q π / 2 , two resonance points corresponding to β = 0.5 were calculated as 0.279 and 2.136, respectively, by using Equation (49). In addition, there was a minimum point at Δ = α / 2 β , which was consistent with the numerical results of (12) (data marked in Figure 10a). This shows that, compared to the system with linear form trichotomous noise disturbed damping, appropriate nonlinear noise will induce new resonance phenomena. The resonant peak on the right side disappears when β = 1 as α Ω q < κ , while the left peak also disappears for β = 2 , so no SR analogous to the traditional form appears in the system, which is due the fact that κ < 0  Figure 10b describes the resonance under different delay intensity: the only peak of output amplitude appears for the none-delay case ( ρ = 0 ), which is due to κ < 0 ; when time delay exists, the phenomenon of double-peak stochastic resonance happens, with the positions of the two optimal points moving apart and decreasing gradually with the increase in delay intensity. Similar to the double-resonance in Figure 10a, the optimal points of SR are, respectively, given by Δ S R 1 and Δ S R 2 given in (49), and the minimum points of output amplitude all appear at Δ = α / 2 β = 1 . It means that on one hand the existence of time delay enriches the stochastic resonance induced by noise intensity; on the other hand, it also proved that compared with the full-delay case ( ρ = 1 ), delay feedback system with both x t and x t τ would exhibit more evident double- resonance.
Under low noise switching rate condition, examining the effect of damping order q and steady-state probability p of noise on the relationship curve of F Δ 2 , we find the hypersensitive response induced by noise intensity, i.e., the system output amplitude shows a sensitive dependence on the noise intensity. Analytical results are shown in Figure 11a, with system parameters set at A = 1 , ρ = 0.7 , γ 0 = 0.5 , α = 1 , β = 0.5 , ω 0 = 1 , v = 0.01 , Ω = 1 , τ = 0.2 . As the noise intensity increases from zero, the output amplitude reaches the resonant peak at Δ 2 = 0.1923 , followed by a sharp decrease in F with a slight increase in noise intensity, the first hypersensitive response occurs. The second hypersensitive response appears around Δ 2 = 2.174 after reaching a minimum at Δ 2 = 0.2765 and then slowly passes another peak point, as the output amplitude increases rapidly to the last peak point. The results obtained by substituting the parameters into (49) are consistent with the positions of the three resonance peak points in the figure, and the twice hypersensitive response phenomena appear after Δ S R 1 and before Δ S R 2 , respectively. Further comparison of Figure 11b,c shows that a two-time hypersensitive response only appears when the damping order and steady-state probability of noise are both small. When one of the two, or both of them, is larger, the amplitude curves before and after resonance positions Δ S R 1 and Δ S R 2 exhibit a more gentle property than those in Figure 11a, and the peak point Δ S R 3 = α / 2 β gradually disappears and eventually becomes the minimum point. At this juncture, noise-induced triple-SR is transformed to a double-SR, thus hypersensitive response disappears. In order to show the joint influence of parameter q and p concretely, different resonance regions that may appear in F Δ 2 relation curve under all parametric conditions of 0 < q < 1 and 0 < p < 1 / 2 are given in Figure 11d. The results suggest that two-time hypersensitivity response may occur only when parameters satisfy p < 0.219 and q < 0.276 simultaneously.
Considering the influence of fractional order q on the system output amplitude, it is difficult to give the explicit expression of the critical value q S R corresponding to the zero of d F q / d q since the output amplitude F is a complex function of q . Even so, dependence of the amplitude of the system output signal on the flatness parameter can be obtained by calculating the expression numerically, given in Figure 12 with parameters A = 1 , ρ = 0.8 , α = 1 , ω 0 = 1 , p = 0.35 , v = 0.01 , Ω = 0.99 , τ = 0.2 . The system response manifests as a monotone increasing function of q under the linear trichotomous noise disturbance. When the colored noise is in nonlinear form ( β = 0.5 ), the output amplitude of the system peaks at q S R = 0.42 , that is, the fractional damping order does make the system produce the phenomenon of GSR. When different damping intensity values are considered (Figure 12b), we determine that γ 0 has a great influence on the GSR phenomenon induced by damping order. Moreover, the output amplitude decreases with the increase in damping order for γ 0 = 0 and γ 0 = 0.75 , while stochastic resonance occurs at q = 0.445 for γ 0 = 0.5 , and a minimal resonance point F min = 1.844 emerges when γ 0 = 0.25 .

6. Numerical Verification Scheme

With these aforesaid preliminaries mentioned in Section 2, we are now in a position to generate realizations of a stochastic process for the trichotomous noise, which is done in the following way:
1. Without loss of generality, assume the noise is located initially at Δ i at time t n , i.e., ζ n = Δ i , i = 1 , 2 , 3 . At n + 1 th time point, to determine whether the process moves to the other two sites or merely remains at the same site, we consider the conditional probability matrix as given by (3). Then, it generates a uniformly distributed random number r n that is falling in 0 , 1 , which is then compared against the three values in the i th line of the conditional probability matrix. If r n < P i 1 h , then we accept the value of the noise Δ 1 , i.e., ζ n + 1 = Δ , and if r n > 1 P i 3 h we accept the value Δ 3 or ζ n + 1 = Δ , otherwise, it will be located at Δ 2 , i.e., ζ n + 1 = 0 ;
2. Taking the site at time t n + 1 as the starting location again, marked as ζ n + 1 = Δ k , another new uniformly distributed random number r n + 1 between 0 , 1 is generated, and we conduct a comparison operation on r n + 1 with the elements in the k th line of (3). If r n + 1 < P k 1 h or r n + 1 > 1 P k 3 h then the value of noise at time t n + 2 is confirmed as ζ n + 2 = Δ or ζ n + 2 = Δ , otherwise we accept ζ n + 2 = 0 . By repeating this two-step procedure one can generate a sequence of random numbers ζ t switching between three values Δ , 0 and Δ .
It should be noted that the time interval between the two steps ( h ) is always fixed and be equal to Δ t , which should be much smaller than the correlation time of the noise ( h < < τ c o r ). Let h = 10 4 , Δ = 5 , and the total time T = 10 , Figure 13 shows four realizations of the trichotomous noise under different steady-state probability parameters and different correlation times. Closer observation would reveal that a larger τ c o r is often accompanied by a longer residence time of a particular state on average, with a larger p value that would lead to a shorter residence time at the origin during the noise switching between Δ and Δ . Comparing the left and right figures, one finds that smaller relaxation time corresponds to a faster average switching rate between the three states. In order to examine the accuracy of the above numerical algorithm, we first calculate the arithmetic average value of 2000 sample results; the first moment of generating time series is determined to satisfy the zero-mean characteristic of trichotomous colored noise ( ζ t 10 3 ). In addition, to check whether the results satisfy the theoretical autocorrelation function, we give the theoretical and numerical results of the normalized auto-correlation function ζ t ζ t + s / ζ 2 t under different characteristic relaxation time, as shown in Figure 14, wherein red markers represent the simulation results for 2000 samples. τ c o r f i t is the exponential order of the fitting curves of discrete results. It can be seen that the relaxation time is very consistent with the theoretical setting, which in turn proves the feasibility (of this algorithm).
Introducing the velocity variable y t , the original second-order FDE (12) can be reduced to first-order equation set:
x ˙ t = y t y ˙ t + γ 0 + α ζ t + β ζ 2 t D 0 , t q C x t + ω 0 2 g x t , x t τ = A cos Ω t + W H t
According to the Caputo fractional derivative approximation method used in the previous work [43], (50) can be discretized as follows:
x n + 1 = x n + y n h y n + 1 = y n h γ 0 + α ζ n + β ζ n 2 n + ω 0 2 1 ρ x n + ω 0 2 ρ x n M A n + D H B H n + 1 B H n n = 1 h q Γ 2 q k = 0 n 1 n k 1 q n k 1 1 q x k + 1 x k
n = 0 , 1 , 2 N , N = T / h , h denotes the time step, T the total time, M = τ / h , t n = n h , x n = x t n , y n = y t n , the initial condition in 0 , τ is set as x t = y t = 0 . A n = A cos Ω t n , ζ n = ζ t n is the trichotomous noise state at t n . B H n = B H t n is the state value of FBM at t n , which is generated by the Wood–Chan algorithm [44]. Based on the discrete Fourier method, this method provides a high computational speed, even in the case of large time step. Figure 15 shows the simulated time series of the process of FBM and their corresponding FGN under different Hurst indices. It can be seen that compared with the classical White Gaussian noise ( H = 0.5 ), the larger H is, the smoother the sample path is. The situation of H < 0.5 is opposite to that of H > 0.5 , moreover, the range of FGN decreases with the increase in Hurst index.
Select h = 10 3 and T = 1000 , the simulation of discrete algorithm (51) for the output signal of the system under different intensity of nonlinear trichotomous noise are shown in Figure 16. Time domain and frequency domain simulation results is displayed for the first 1000 s, with parameters A = 1 , ρ = 0.8 , γ 0 = 0.5 , α = 1 , β = 1 , ω 0 = 1 , p = 0.35 , v = 0.1 , Ω = 0.99 , τ = 0.2 . By comparing (a) and (c), it can be seen that larger noise intensity does not drive the system to output signals with larger amplitude, and the results of Δ 2 = 0.1 show stronger periodic characteristics than those of Δ 2 = 1 . (b) and (d) show the amplitude of the system output signal at each frequency, and it is found that it peaks at the external excitation frequency Ω = 0.99 in both cases. In the case of Δ 2 = 1 and Δ 2 = 0.1 , the numerical results of the amplitude at the external excitation frequency are F n u m = 1.587 and F n u m = 2.545 , respectively, and the theoretical results obtained by the analytical expression (35) are F = 1.592 and F n u m = 2.545 , respectively, fully showing that the numerical results of a single sample are quite close to the theoretical prediction results.
In order to verify the feasibility of the numerical simulation algorithm (51), independent Monte Carlo repeated experiments should be conducted to obtain I samples, x i t n , i = 1 , 2 , I , n = 1 , 2 , N , N represents the total steps. The first-order moment of the system response at time t n can be given by the following formula:
x t n = i = 1 I x i t n / I
The Fourier coefficient of the response is taken as the measurement of the output amplitude, and the time series (52) is used to give the simulation results of the output amplitude as follows
F n u m = M a x 2 FFT i = 1 I x i t n / N I = F s 2 + F c 2
F c = 2 n = 1 N i = 1 I x i t n cos Ω t n / N I ,   F s = 2 n = 1 N i = 1 I x i t n sin Ω t n / N I
where FFT denotes the fast Fourier transform.
Firstly, examining the influence of time step on the time consumption of the program, fix T = 10 0 , I = 1000 and select step in the range 10 3 h 10 0 . The average time of 1000 Monte Carlo experiments using algorithm (651) was recorded, and results under different time steps are shown in Figure 17a. The time consuming T c o s t increases exponentially with the decrease in the time step h , and the approximate fitting relationship reads T c o s t 164 e 1253 h , the calculation time of a single sample increases greatly especially for h 10 2 , which is bound to lead to a sharp increase in the experimental time of a large number of independent samples. In order to verify the convergence of numerical simulation algorithms (51) and (53), taking the output amplitude as a function of h , the simulation results based on (53) and theoretical results based on (35) are compared in Figure 17b. When h decreases, the numerical simulation amplitude gradually approaches the theoretical result F = 1.592 , and tends to 1.603 without significant fluctuation when h 10 2 , the inherent error is caused by the numerical discretization of the fractional system (12). Numerical results are no longer significantly affected by step size when h 10 2 , so time step h = 0.01 is selected in the following part.
Under the same parameter setting as Figure 16, I = 1000 , h = 10 3 , the relationship between F n u m and the total simulation time T for Δ 2 = 0.1 and Δ 2 = 0.1 was examined by (53), depicted by Figure 18a,b, respectively. When T < 100 , the simulation results differ greatly from the theoretical results due to the transient influence of the initial condition. When T is large enough, the accuracy of the numerical results is no longer significantly affected by the initial condition, yet fluctuates in a small range near the theoretical ones rather than gradually approaching them. In fact, due to the existence of nonlinear trichotomous noise, the output amplitude is actually a random variable. The statistical result (53) with a fixed sample number is bound to bring certain error. To measure this error, we still use the Monte Carlo method to calculate the mean squared error (MSE) between F n u m and the theoretical results given by (35)
σ = 1 I i = 1 I F i F 2
where F i represents the single sample amplitude estimation of the i th experiment, given by
F i = 2 N M a x FFT x i t n = F i s 2 + F i c 2 F c = 2 n = 1 N x i t n cos Ω t n / N ,   F s = 2 n = 1 N x i t n sin Ω t n / N
In virtue of (54), the mean square error of simulation results corresponding to Δ 2 = 0.1 and Δ 2 = 1 are given in Figure 18b,d, respectively. It is found that the mean square error tends to be stable with the increase in simulation time, σ 7.3 × 10 2 for Δ 2 = 0.1 and σ 4.2 × 10 2 for Δ 2 = 1 , which reflects the reliability of the algorithm (51). Under the same parameter conditions, the influence of the number I of the Monte Carlo experiment on the steady-state output amplitude F and mean square error σ was investigated, and the simulation time was fixed T = 500 . Comparison between Δ 2 = 0.1 (Figure 19a,b) and Δ 2 = 1 (Figure 19c,d) shows that, when the number of experimental samples is large enough, the simulated output amplitude oscillates within a small range of the theoretical one, and the mean square error is around 6 × 10 2 , which indicates that the algorithm in this paper (51) is basically consistent with the theoretical analysis result given by (35).

7. Conclusions

In this paper, a general investigation into stochastic resonance in the fractional Langevin system with internal FGN is reported theoretically and numerically, wherein the fractional Brown particle is driven by an external periodic excitation and multiplicative polynomial colored trichotomous noise. The theoretical analysis of the steady-state amplitude of the system output signal and the first moment of the system response, relying on a stochastic vibrational mechanism, are obtained by extending the stochastic average method and utilizing a fractional Shapiro–Loginov formula and fractional Laplace transformation law. Bona fide stochastic resonance induced by the frequency of periodic signal and generalized stochastic resonance induced by delay τ , noise intensity Δ 2 and damping order q are analyzed in detail. An efficient and feasible algorithm for generating trichotomous noise is proposed, by which a numerical discretization algorithm for the fractional Langevin system is presented, and the numerical results agree well with the theoretical prediction.
The theoretical treatment provides a convenient and shortcut approach to reveal the mechanism of two kinds of stochastic resonance. For bona fide SR, new SMR phenomena are found in the system with trichotomous noise (especially, the one with a slow switching rate), which does not appear in the system without colored noise. Under appropriate noise intensity, bimodal and trimodal bona fide SR are more likely to occur with a slow switching colored noise and small fractional damping order (Figure 6). In other ways, dependence of GSR on joint parameters v , Δ 2 (Figure 9) indicates that multimodal GSR induced by delay τ happens only for the colored noise with a slow switching rate. Moreover, the two-peak GSR phenomenon induced by noise intensity implies that compared to the linear trichotomous noise, non-linear noise with appropriate coefficient β can induce a new resonance pattern. It is determined that the damping intensity coefficient has a significant effect on the GSR induced by the damping order. In addition, explicit expressions of critical noise intensity corresponding to the peak value of GSR are obtained for the first time. In the case of slow transfer rate noise, a newfangled phenomenon of double hypersensitive responses, induced by a variation in noise intensity, is discovered, by which a novel fact has been validated and confirmed that the nonlinear colored noise induces a new dynamic behavior, which has not been reported in previous work. It can be further predicted that double hypersensitivity response only occurs under the condition p < 0.219 and q < 0.276 . This allows us to control a system effectively by selecting appropriate parameters. On account of the significance of the SR analysis, the methods in this paper might be useful in dealing with some real stochastic fractional-order issues. It is believed that these results shed new light on the studies of SR induced by non-linear colored noises, and make further exploration of dynamical behaviors of anomalous diffusion systems with delayed feedback possible.

Author Contributions

Methodology, J.L.G.G., H.C. and T.S.; validation, J.L.G.G., H.C. and T.S.; formal analysis, J.L.G.G., H.C. and T.S.; investigation, J.L.G.G., H.C. and T.S.; resources, J.L.G.G., H.C. and T.S.; data curation, J.L.G.G., H.C. and T.S.; writing—original draft preparation, Z.Y. and X.L.; writing—review and editing, Z.Y. and X.L.; visualization, Z.Y. and X.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by National Natural Science Foundation of China (No.12002194; No.12072178); Project (No.ZR2020QA037; No.ZR2020MA054; No.2020KJI003) supported by Shandong Provincial Natural Science Foundation.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

This paper is not based on any data.

Conflicts of Interest

There is no conflict of interest.

References

  1. Ardjouni, A.; Djoudi, A. Positive solutions for first-order nonlinear Caputo-Hadamard fractional relaxation differential equations. Kragujev. J. Math. 2021, 45, 897–908. [Google Scholar] [CrossRef]
  2. Li, Y.; Liu, F.; Turner, I.W.; Li, T. Time-fractional diffusion equation for signal smoothing. Appl. Math. Comput. 2018, 326, 108–116. [Google Scholar] [CrossRef] [Green Version]
  3. Singh, J.; Jassim, H.K.; Kumar, D. An efficient computational technique for local fractional Fokker Planck equation. Phys. A Stat. Mech. Its Appl. 2020, 555, 124525. [Google Scholar] [CrossRef]
  4. Gitterman, M. Mean first passage time for anomalous diffusion. Phys. Rev. E 2000, 62, 6065. [Google Scholar] [CrossRef] [PubMed]
  5. Barzykin, A.V.; Seki, K. Stochastic resonance driven by Gaussian multiplicative noise. EPL (Europhys. Lett.) 1997, 40, 117. [Google Scholar] [CrossRef]
  6. Ren, R.; Luo, M.; Deng, K. Stochastic resonance in a fractional oscillator subjected to multiplicative trichotomous noise. Nonlinear Dyn. 2017, 90, 379–390. [Google Scholar] [CrossRef]
  7. Kim, H.; Tai, W.C.; Parker, J.; Zuo, L. Self-tuning stochastic resonance energy harvesting for rotating systems under modulated noise and its application to smart tires. Mech. Syst. Signal Process. 2019, 122, 769–785. [Google Scholar] [CrossRef]
  8. Ren, R.; Deng, K. Noise and periodic signal induced stochastic resonance in a Langevin equation with random mass and frequency. Phys. A Stat. Mech. Its Appl. 2019, 523, 145–155. [Google Scholar] [CrossRef]
  9. Zhong, S.; Wei, K.; Gao, S.; Ma, H. Stochastic resonance in a linear fractional Langevin equation. J. Stat. Phys. 2013, 150, 867–880. [Google Scholar] [CrossRef]
  10. Lai, L.; Zhang, L.; Yu, T. Collective behaviors in globally coupled harmonic oscillators with fluctuating damping coefficient. Nonlinear Dyn. 2019, 97, 2231–2248. [Google Scholar] [CrossRef]
  11. Tian, Y.; Zhong, L.F.; He, G.T.; Yu, T.; Luo, M.K.; Stanley, H.E. The resonant behavior in the oscillator with double fractional-order damping under the action of nonlinear multiplicative noise. Phys. A Stat. Mech. Its Appl. 2018, 490, 845–856. [Google Scholar] [CrossRef]
  12. Li, P.; Tang, Z.; Zhang, Q.; Zhang, W. Dynamic behaviors of a star-coupled underdamped system with multiplicative quadratic noise and periodic excitation. Phys. Scr. 2021, 96, 085218. [Google Scholar] [CrossRef]
  13. Wang, R.; Wang, B. Random dynamics of p-Laplacian lattice systems driven by infinite-dimensional nonlinear noise. Stoch. Processes Appl. 2020, 130, 7431–7462. [Google Scholar] [CrossRef]
  14. Gonçalves, P.; Jara, M. Quadratic fluctuations of the symmetric simple exclusion. ALEA 2019, 16, 605–632. [Google Scholar] [CrossRef]
  15. Zhong, S.; Wang, H.; Ma, H.; Luo, M.; Zhang, L. Nonlinear effect of time delay on the generalized stochastic resonance in a fractional oscillator with multiplicative polynomial noise. Nonlinear Dyn. 2017, 89, 1327–1340. [Google Scholar] [CrossRef]
  16. Töpfer, J.D.; Sigurdsson, H.; Pickup, L.; Lagoudakis, P.G. Time-delay polaritonics. Commun. Phys. 2020, 3, 1–8. [Google Scholar] [CrossRef]
  17. Du, L.C.; Mei, D.C. Stochastic resonance in a bistable system with global delay and two noises. Eur. Phys. J. B 2012, 85, 75. [Google Scholar] [CrossRef]
  18. Song, X.; Wang, H.; Chen, Y. Coherence resonance in an autaptic Hodgkin–Huxley neuron with time delay. Nonlinear Dyn. 2018, 94, 141–150. [Google Scholar] [CrossRef]
  19. Sarkar, M. Synchronization transition in the two-dimensional Kuramoto model with dichotomous noise. Chaos Interdiscip. J. Nonlinear Sci. 2021, 31, 083102. [Google Scholar] [CrossRef]
  20. Peng, H.; Ren, R.; Li, P.; Yu, T. Trichotomous noise induced resonance behavior of two coupled harmonic oscillators with fluctuating mass. Phys. Scr. 2020, 95, 075214. [Google Scholar] [CrossRef]
  21. Jin, Y.; Wang, H. Noise-induced dynamics in a Josephson junction driven by trichotomous noises. Chaos Solitons Fractals 2020, 133, 109633. [Google Scholar] [CrossRef]
  22. Bobryk, R.V. Stochastic multiresonance in oscillators induced by bounded noise. Commun. Nonlinear Sci. Numer. Simul. 2021, 93, 105460. [Google Scholar] [CrossRef]
  23. Evangelista, L.R.; Lenzi, E.K. Fractional Diffusion Equations and Anomalous Diffusion; Cambridge University Press: Cambridge, UK, 2018. [Google Scholar]
  24. Sabri, A.; Xu, X.; Krapf, D.; Weiss, M. Elucidating the origin of heterogeneous anomalous diffusion in the cytoplasm of mammalian cells. Phys. Rev. Lett. 2020, 125, 058101. [Google Scholar] [CrossRef] [PubMed]
  25. Cui, B.; Zaccone, A. Generalized Langevin equation and fluctuation-dissipation theorem for particle-bath systems in external oscillating fields. Phys. Rev. E 2018, 97, 060102. [Google Scholar] [CrossRef] [Green Version]
  26. Kou, S.C.; Xie, X.S. Generalized Langevin equation with fractional Gaussian noise: Subdiffusion within a single protein molecule. Phys. Rev. Lett. 2004, 93, 180603. [Google Scholar] [CrossRef]
  27. Fulinski, A. Fractional brown motions. Acta Phys. Pol. B 2020, 51, 1097–1129. [Google Scholar] [CrossRef]
  28. Mishura, I.U.S.; Mišura, J.S.; Mishura, Y.; Mišura, Û.S. Stochastic Calculus for Fractional Brownian Motion and Related Processes; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  29. Milton, J.G.; Insperger, T.; Cook, W.; Harris, D.M.; Stepan, G. Microchaos in human postural balance: Sensory dead zones and sampled time-delayed feedback. Phys. Rev. E 2018, 98, 022223. [Google Scholar] [CrossRef] [Green Version]
  30. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998. [Google Scholar]
  31. Guillouzic, S.; L’Heureux, I.; Longtin, A. Small delay approximation of stochastic delay differential equations. Phys. Rev. E 1999, 59, 3970. [Google Scholar] [CrossRef] [Green Version]
  32. Ricardo, H.J. A Modern Introduction to Differential Equations; Academic Press: Cambridge, MA, USA, 2020. [Google Scholar]
  33. Gao, S.L. Generalized stochastic resonance in a linear fractional system with a random delay. J. Stat. Mech. Theory Exp. 2012, 2012, P12011. [Google Scholar] [CrossRef]
  34. Verma, D. Applications of Laplace Transformation for Solving Various Differential Equations with Variable Coefficients. Int. J. Innov. Res. Sci. Technol. 2018, 4, 124–127. [Google Scholar]
  35. Mainardi, F. Fractional calculus: Theory and applications. Mathematics 2018, 6, 145. [Google Scholar] [CrossRef] [Green Version]
  36. Balcerzak, M.; Leonetti, P. A Tauberian theorem for ideal statistical convergence. Indag. Math. 2020, 31, 83–95. [Google Scholar] [CrossRef]
  37. Déjardin, P.M.; Titov, S.V.; Cornaton, Y. Linear complex susceptibility of long-range interacting dipoles with thermal agitation and weak external ac fields. Phys. Rev. B 2019, 99, 024304. [Google Scholar] [CrossRef] [Green Version]
  38. Heller, M.P.; Serantes, A.; Spaliński, M.; Svensson, V.; Withers, B. Hydrodynamic gradient expansion in linear response theory. Phys. Rev. D 2021, 104, 066002. [Google Scholar] [CrossRef]
  39. Zeng, X.; Lu, X.; Liu, Z.; Jin, Y. An adaptive fractional stochastic resonance method based on weighted correctional signal-to-noise ratio and its application in fault feature enhancement of wind turbine. ISA Trans. 2021, 120, 18–32. [Google Scholar] [CrossRef]
  40. Escotet-Espinoza, M.S.; Moghtadernejad, S.; Oka, S.; Wang, Z.; Wang, Y.; Roman-Ospino, A.; Schäfer, E.; Cappuyns, P.; Assche, I.V.; Futran, M.; et al. Effect of material properties on the residence time distribution (RTD) characterization of powder blending unit operations. Part II of II: Application of models. Powder Technol. 2019, 344, 525–544. [Google Scholar] [CrossRef]
  41. Gitterman, M. Classical harmonic oscillator with multiplicative noise. Phys. A Stat. Mech. Its Appl. 2005, 352, 309–334. [Google Scholar] [CrossRef]
  42. Yang, S.; Deng, M.; Ren, R. Stochastic resonance of fractional-order Langevin equation driven by periodic modulated noise with mass fluctuation. Adv. Differ. Equ. 2020, 2020, 81. [Google Scholar] [CrossRef] [Green Version]
  43. Yan, Z.; Liu, X. Fractional-order harmonic resonance in a multi-frequency excited fractional Duffing oscillator with distributed time delay. Commun. Nonlinear Sci. Numer. Simul. 2021, 97, 105754. [Google Scholar] [CrossRef]
  44. Wood, A.T.A.; Chan, G. Simulation of stationary Gaussian processes in [0, 1]. J. Comput. Graph. Stat. 1994, 3, 409–432. [Google Scholar] [CrossRef]
Figure 1. (Color online) The shape of damping kernel function γ t under different q .
Figure 1. (Color online) The shape of damping kernel function γ t under different q .
Fractalfract 06 00191 g001
Figure 2. (a) Output amplitude of system (36) versus excitation frequency Ω under different fractional order. (b) Analytical result of F * versus Ω and q according to Equation (38).
Figure 2. (a) Output amplitude of system (36) versus excitation frequency Ω under different fractional order. (b) Analytical result of F * versus Ω and q according to Equation (38).
Fractalfract 06 00191 g002
Figure 3. (ac) Triple-, double- and single-peak Bona fide SR under different noise switching rate. (d) Analytical result of output steady-state amplitude of system (12) versus Ω and v .
Figure 3. (ac) Triple-, double- and single-peak Bona fide SR under different noise switching rate. (d) Analytical result of output steady-state amplitude of system (12) versus Ω and v .
Fractalfract 06 00191 g003
Figure 4. (a) Output steady-state amplitude versus excitation frequency with different fractional order under slow switching noise disturbance. (b) Dependence of output steady-state amplitude on excitation frequency and fractional damping order.
Figure 4. (a) Output steady-state amplitude versus excitation frequency with different fractional order under slow switching noise disturbance. (b) Dependence of output steady-state amplitude on excitation frequency and fractional damping order.
Fractalfract 06 00191 g004
Figure 5. (a) Output steady-state amplitude versus excitation frequency with different noise parameter β for the case of fast switching noise. (b) Different bona fide patterns occur for different parameters β and γ 0 .
Figure 5. (a) Output steady-state amplitude versus excitation frequency with different noise parameter β for the case of fast switching noise. (b) Different bona fide patterns occur for different parameters β and γ 0 .
Fractalfract 06 00191 g005
Figure 6. Dependence of bona fide SR on combined parameters v , q under different noise intensity. (a) Δ 2 = 1 , (b) Δ 2 = 0.1 , (c) Δ 2 = 6 .
Figure 6. Dependence of bona fide SR on combined parameters v , q under different noise intensity. (a) Δ 2 = 1 , (b) Δ 2 = 0.1 , (c) Δ 2 = 6 .
Fractalfract 06 00191 g006
Figure 7. Delay induced GSR under noiseless situation in system (42) according to Equation (38).
Figure 7. Delay induced GSR under noiseless situation in system (42) according to Equation (38).
Fractalfract 06 00191 g007
Figure 8. GSR induced by a delay in the system (12) under different correlation rate and intensity of noise (a) Δ 2 = 0.6 , (b) v = 0.01 .
Figure 8. GSR induced by a delay in the system (12) under different correlation rate and intensity of noise (a) Δ 2 = 0.6 , (b) v = 0.01 .
Fractalfract 06 00191 g008
Figure 9. Different GSR patterns relying on combined parameters v , Δ 2 under (a) q = 0.1 , (b) q = 0.2 , (c) q = 0.25 , (d) q = 0.9 . Single-peak SR happens in the grey area, double-peak SR happens in the cyan area, and triple-peak SR happens in the green area, respectively.
Figure 9. Different GSR patterns relying on combined parameters v , Δ 2 under (a) q = 0.1 , (b) q = 0.2 , (c) q = 0.25 , (d) q = 0.9 . Single-peak SR happens in the grey area, double-peak SR happens in the cyan area, and triple-peak SR happens in the green area, respectively.
Fractalfract 06 00191 g009
Figure 10. Influence of the variation in noise intensity on the output steady-state amplitude of system (16) under different nonlinear noise coefficient and delay intensity (a) ρ = 0.8 , (b) β = 0.5 .
Figure 10. Influence of the variation in noise intensity on the output steady-state amplitude of system (16) under different nonlinear noise coefficient and delay intensity (a) ρ = 0.8 , (b) β = 0.5 .
Fractalfract 06 00191 g010
Figure 11. (a) Phenomenon of double hypersensitive response for p = 0.01 , q = 0.01 . (b) effect of q on hypersensitive response for p = 0.01 , (c) effect of p on hypersensitive response for q = 0.01 , (d) parametric regions I, II and III correspond, respectively, to double, single and non- hypersensitive response.
Figure 11. (a) Phenomenon of double hypersensitive response for p = 0.01 , q = 0.01 . (b) effect of q on hypersensitive response for p = 0.01 , (c) effect of p on hypersensitive response for q = 0.01 , (d) parametric regions I, II and III correspond, respectively, to double, single and non- hypersensitive response.
Fractalfract 06 00191 g011
Figure 12. GSR induced by the order of fractional damping under different noise coefficients and damping intensity. (a) γ 0 = 0.5 , (b) β = 0.5 .
Figure 12. GSR induced by the order of fractional damping under different noise coefficients and damping intensity. (a) γ 0 = 0.5 , (b) β = 0.5 .
Fractalfract 06 00191 g012
Figure 13. Generation of trichotomous colored noise for (a) τ c o r = 0.3 , p = 1 / 4 , (b) τ c o r = 0.06 , p = 1 / 4 , (c) τ c o r = 0.3 , p = 2 / 5 , (d) τ c o r = 0.06 , p = 2 / 5 .
Figure 13. Generation of trichotomous colored noise for (a) τ c o r = 0.3 , p = 1 / 4 , (b) τ c o r = 0.06 , p = 1 / 4 , (c) τ c o r = 0.3 , p = 2 / 5 , (d) τ c o r = 0.06 , p = 2 / 5 .
Fractalfract 06 00191 g013
Figure 14. Plots of normalized correlation function for p = 0.4 , Δ = 5 , for three different noise profiles with given correlation times τ = 0.5 , 2 and 4 . Colored markers are numerical results, and bold curves correspond to fitted results.
Figure 14. Plots of normalized correlation function for p = 0.4 , Δ = 5 , for three different noise profiles with given correlation times τ = 0.5 , 2 and 4 . Colored markers are numerical results, and bold curves correspond to fitted results.
Fractalfract 06 00191 g014
Figure 15. Realizations of FBM and FGN under three values of the Hurst index. (a,d) H = 0.2 , (b,e) H = 0.5 , (c,f) H = 0.8 .
Figure 15. Realizations of FBM and FGN under three values of the Hurst index. (a,d) H = 0.2 , (b,e) H = 0.5 , (c,f) H = 0.8 .
Fractalfract 06 00191 g015
Figure 16. Simulation results of x t in terms of time and frequency domain for (a,b) Δ 2 = 1 , (c,d) Δ 2 = 0.1 .
Figure 16. Simulation results of x t in terms of time and frequency domain for (a,b) Δ 2 = 1 , (c,d) Δ 2 = 0.1 .
Fractalfract 06 00191 g016
Figure 17. Effect of step size on simulated results of output steady-state amplitude and time cost (a) T = 10 0 , (b) Δ 2 = 1 .
Figure 17. Effect of step size on simulated results of output steady-state amplitude and time cost (a) T = 10 0 , (b) Δ 2 = 1 .
Fractalfract 06 00191 g017
Figure 18. Dependence of output steady-state amplitude and MSE on total simulation time (a,b) Δ 2 = 0.1 , (c,d) Δ 2 = 1 .
Figure 18. Dependence of output steady-state amplitude and MSE on total simulation time (a,b) Δ 2 = 0.1 , (c,d) Δ 2 = 1 .
Fractalfract 06 00191 g018
Figure 19. Effect of samples on numerically simulated results of output steady-state amplitude and MSE. (a,b) Δ 2 = 0.1 , (c,d) Δ 2 = 1 .
Figure 19. Effect of samples on numerically simulated results of output steady-state amplitude and MSE. (a,b) Δ 2 = 0.1 , (c,d) Δ 2 = 1 .
Fractalfract 06 00191 g019
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yan, Z.; Guirao, J.L.G.; Saeed, T.; Chen, H.; Liu, X. Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise. Fractal Fract. 2022, 6, 191. https://doi.org/10.3390/fractalfract6040191

AMA Style

Yan Z, Guirao JLG, Saeed T, Chen H, Liu X. Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise. Fractal and Fractional. 2022; 6(4):191. https://doi.org/10.3390/fractalfract6040191

Chicago/Turabian Style

Yan, Zhi, Juan L. G. Guirao, Tareq Saeed, Huatao Chen, and Xianbin Liu. 2022. "Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise" Fractal and Fractional 6, no. 4: 191. https://doi.org/10.3390/fractalfract6040191

APA Style

Yan, Z., Guirao, J. L. G., Saeed, T., Chen, H., & Liu, X. (2022). Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise. Fractal and Fractional, 6(4), 191. https://doi.org/10.3390/fractalfract6040191

Article Metrics

Back to TopTop