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
to denote the ubiquitous time-continuous Markov process, which consists of shifts between three different levels, i.e.,
may jump between three values
,
and
, with
. 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
at different values are
thus
would be degenerated into the dichotomous form for
, the conditional transition probability between the three states is determined by the master equation [
21]
where
,
,
,
Combining the normalization condition
and the initial conditions
, one can use Equation (2) to give the a probability transition matrix of noise in regard to the time interval
, and the shift between the states
,
and
obey a pattern of Poisson distribution, which can be given as follows
All the characters of
are completely contained in (1) and (3), from which the mean value at arbitrary time
can be calculated as
and the autocorrelation function
is the switching rate, and correlation time
equals to the reciprocal of the noise switching rate
is the noise intensity. The mean square value and the counterpart mean quadruplicate value of the noise reads as
,
, thus one gets the flatness parameter [
22] of
:
,
. For different
,
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
and
, 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
wherein
,
and
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
is extremely short, say, modeled as GWN. The second assumption is that the damping coefficient
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
, the viscosity coefficient
should manifest as a time-dependent function
[
25] rather than a constant coefficient, the damping force at any time
is acquired by integrating the velocity at all times before
according to the weight function
. From the above consideration, it can be further turned into the generalized Langevin equation
where
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
The shape of damping kernel function with the friction coefficient
is shown in
Figure 1, it finds that
has a power-law attenuation trend with the speed of
, and it would degenerate into Delta function when
approach
. Inserting (6) into (5) and using the definition of Caputo derivative, Equation (5) can be further transformed to the following FDE
The well-known Brownian motion
is generally used in the traditional stochastic dynamics, to model the normal diffusion with mean-square displacement (MSD)
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
, 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
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:
, with
the Boltzmann constant and
the absolute temperature. Substituting (6) into the fluctuation dissipation equality gives the power-law form autocorrelation function of random force
Here we consider
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):
.
denotes the noise strength,
denotes the Hurst index with
.
is reduced to a standard Brownian motion for
, which has independent increments. For
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:
where the correlation function of FGN can be obtained by the covariance of FBM
When
,
corresponds to the subdiffusion movement with a negative correlation function,
corresponds the opposite situation, and it turns to the impulse function
for
, say, the normal diffusion case. Utilizing the fluctuation dissipation theorem, again one gets
The Hurst parameter should be constrained in to guarantee the mean-square value to be positive, thus the memory kernel is positive as well. Combining (10) and (8) one has , that is, the attenuation parameters of memory kernel function correspond exactly to the order of the fractional derivative. An intensity setting of FGN as 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
, disturbed by nonlinear colored noise and internal FGN, which is given by the following integro-differential equation:
and
denote the acceleration and velocity, respectively.
is external harmonic excitation,
the inherent frequency.
(
) represents the nonlinear trichotomous noise disturbance on damping intensity, which manifests itself as the polynomial form
.
is symmetric trichotomous colored noise with three states,
is zero-mean FGN with autocorrelation function given by (9) and memory damping kernel given by (6), respectively. Considering the potential field force
,
, where
represents delay intensity.
and
correspond, respectively, to non-delay and global item delay systems [
17],
corresponds to the system containing delay feedback, wherein
is also referred as the feedback gain [
29]. For computational convenience sake assume
, 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
3. Theoretical Analysis of Response
With the constructed equation
, and the arbitrary power of
can be expressed by
is a positive integer, with (13) the polynomial noise in Equation (12) can be simplified to
where “
” denotes just the value
rounded down to the nearest integer. For the case of
, 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 (
) 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 “
” in
and
will be dropped, and the simplified versions
and
will be used to describe the linear coefficient and nonlinear coefficient of the polynomial function
of the trichotomous colored noise, respectively.
As for the function
of the small time delay, when
no delay appeared, and this mainstream situation has attracted widespread attention. Meanwhile the global delay case happens for
and the stiffness term disappears [
15,
17]. The scenario of
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,
and can be expressed as a Taylor expansion series with respect to the time delay around
, with the convergence precision the same order of magnitude as
. In fact, such a Taylor expansion approximation of a function with respect to
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:
where
.
Combining Equation (15) and, as discussed below, Equation (14), we obtain a non-delayed stochastic differential equation, which is equivalent to Equation (12):
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
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
. For this purpose, the Shapiro–Loginov formula [
8] should be taken account for the exponential correlative trichotomous process:
and the fractional Shapiro-Loginov formula:
Here,
is considered as a function of
. Taking average of Equation (16), one gets
The prerequisite condition
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
are given by
and
It is apparent from the above autocorrelation function of the mean value
that
does not satisfy the conditions of the Shapiro–Loginov formula, however, the situation becomes gratifying for the adjustive variable
, that
Thus, the random variable
is applicable for various kinds of Shapiro–Loginov formulas (of course also including the fractional version (3.5))
(19) leads to the new correlation factor
appear in (18):
By substituting Equations (17) and (20) into (18), one gets the equation
To solve the new factor
, multiplying Equation (16) by
and averaging over the gained equation, one arrives at
in Equation (22) the uncorrelated postulate between
and the FGN
has already been adopted, i.e.,
,
Reusing the Shapiro–Loginov formula, we have
Substituting (17), (20) and (23) into Equation (22) yields
In order to solve the term
, multiplying Equation (16) by
and the averaging the gained results, one obtains:
the first averaged term in Equation (25) can also be derived by (19) and (23)
From (26) the averaging value
can be obtained as the expression
Inserting (17), (20) and (27) into Equation (25), meanwhile, replacing
,
for
and
, respectively, then (27) can be transformed to
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
and subtracting the gained equation from Equation (28), one can obtain:
The collection of (21), (24) and (29) truly consists of a close linear system of second-order differential equations for three variables
,
,
. By virtue of the Laplace transform technique, this equation set can be equivalently changed to a corresponding linear algebraic equation set, given as follows
Displacement’s first-order moment can be formally expressed formally as bellow:
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:
“” the convolution operation, , , , with the inverse Laplace operator the primitive function can be calculated by the integral of its corresponding Laplace transformation function : .
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
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
, by making use of the Tauberian theorem [
36], one gets the estimation:
,
,
, 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).
the complex susceptibility of the stochastic system given by Equation (12) can be given based on the form of (34) [
37]:
, where
denotes the imaginary unit with
,
and
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
. With the output amplitude
and the phase shift given as,
, it is worth nothing that
, thus
and
can also be calculated by
,
. 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
Through the expressions of the transfer function given in (31) and the coefficients below (29), after a series of simplification calculations, one can obtain
where
,
are listed as follows:
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.,
or
, then the corresponding system equation reads
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
after the Laplace transformation of Equation (37) with respect to the variable
, the mean value of steady-state for
can be calculated as
and
denote the corresponding amplitude and phase shift, respectively.
with
for any fractional order satisfying
, the steady-state output amplitude of the system (37) can be obtained as
In fact, according to the linear response theory, if one substitutes 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 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
, with system parameters
,
,
,
,
. It finds that the resonance of the system output signal is striking when
is small, while for the opposite situation
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
, the single-resonance point
gradually increases and
tends to be flat. Indeed, the critical resonant frequency
must satisfy
On the one hand, expression (38) implies
when the Hurst index is close to one. According to (39), the critical SR frequency of the system can be obtained as
, According to the form of
, single-peak SR will certainly emerge in the system for small delay, and
when
approaching to
. Therefore, in
Figure 2b
exhibits a saltation within a certain frequency range for
. On the other hand, consider the case
, then substituting
into Equation (39) give rise to the resonance peak
, with the occurrence condition for bona fide SR
and the optimal OAG
.
Figure 2a shows the critical frequency values of the peak points
and
are 1.38 and 0.67, respectively, which are consistent with the above analysis results.
When one considers the stochastic polynomial trichotomous noise
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
,
,
,
,
,
,
,
,
,
Figure 3 shows the bona fide SR in the disturbed system (12) with a small flatness parameter
under different noise transfer rates, three-peak resonance, double-peak resonance and single-peak resonance appear at
,
and
, 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.
Replacing (40) for
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
When
is close to
, the external excitation frequency corresponding to the three resonance extreme points can also be calculated in the same way, given as
As
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
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
. It should be pointed out that all the system parameters considered in this section meet the stability condition
.
The output amplitude of the system can also be approximately estimated for the fast-switching noise scenarios, as
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)
with the equivalent damping coefficient
. For the case of
(
), by solving
in Equation (41), one has
, the single-peak-resonance always exists in the system, which is consistent with the single peak point position
in
Figure 3c. For
, only single-resonance may occur at
, and the existence condition is
In the case of the slow-switching case
, and the bona fide SR phenomenon under different flatness parameters is shown in
Figure 4. For
, the bona fide SR occurs at three critical frequencies;
,
and
, as
increases, triple resonant peaks gradually converge to double and single peaks, which is similar to the situation in
Figure 3. The difference is that
decreases with the increase in
. When
is large enough, the resonance phenomenon disappears finally. Consider the other case
(we might as well set
),
,
,
Figure 5 shows SR phenomena under different polynomial noise coefficients. Since the condition (43) is not met for
, 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
corresponding to the resonant peak is almost monotonically dependent on the critical value of
. When
, 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
, under the polynomial trichotomous noise coefficient setting
,
. Considering the system parameters
,
,
,
,
, 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
. When
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
of the system output amplitude corresponds to the root of the equation
solving (44) and giving rise to the critical delay value
Inserting (45) into Equation (38), the optimal output amplitude at the critical delay
is given by
, 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
with the maximum output amplitude
, which is consistent with the critical delay analysis result in (45). Examining the disturbed situation
,
,
,
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
,
,
,
,
,
. 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
, 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
,
and
.
with
,
. 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
. When
is large enough, the critical delay of the single-peak GSR under the fast-switching noise case can be calculated as follows:
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 (
) in the system, so that the inherent resonance peak could only exist at the critical delay given by (45). When
,
is very close to
, while there is a large gap between
and
. Therefore, the peak position of GSR should be
and
in (46), and it is not difficult to get that to
and
, which is consistent with the resonance points of numerical results labeled in the figure. When
, the intensity of the slow switching noise is large enough, and the three equivalent damping coefficients
,
and
, 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
and
more carefully, different fractional orders are selected to give the joint dependence between GSR induced by delay and parameter
. It is found that SMR tends to appear for small
(cyan and green areas in
Figure 9a), the single-peak GSR region becomes larger and larger as
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
varied with the other system parameters that were fixed as
,
,
,
,
,
,
,
,
. 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
, 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
exists if the nonlinear part of colored noise is not considered (
), and the single-peak SR only occurs when the trichotomous noise switches to
. Indeed, by inserting
and (40) into Equation (38), one can obtain a real solution only if the equivalent damping coefficient is equal to
Two amplitude spikes appear in the system when
increases to
. 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
With
, two resonance points corresponding to
were calculated as 0.279 and 2.136, respectively, by using Equation (49). In addition, there was a minimum point at
, 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
as
, while the left peak also disappears for
, so no SR analogous to the traditional form appears in the system, which is due the fact that
Figure 10b describes the resonance under different delay intensity: the only peak of output amplitude appears for the none-delay case (
), which is due to
; 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
and
given in (49), and the minimum points of output amplitude all appear at
. 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 (
), delay feedback system with both
and
would exhibit more evident double- resonance.
Under low noise switching rate condition, examining the effect of damping order
and steady-state probability
of noise on the relationship curve of
, 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
,
,
,
,
,
,
,
,
. As the noise intensity increases from zero, the output amplitude reaches the resonant peak at
, followed by a sharp decrease in
with a slight increase in noise intensity, the first hypersensitive response occurs. The second hypersensitive response appears around
after reaching a minimum at
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
and before
, 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
and
exhibit a more gentle property than those in
Figure 11a, and the peak point
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
and
concretely, different resonance regions that may appear in
relation curve under all parametric conditions of
and
are given in
Figure 11d. The results suggest that two-time hypersensitivity response may occur only when parameters satisfy
and
simultaneously.
Considering the influence of fractional order
on the system output amplitude, it is difficult to give the explicit expression of the critical value
corresponding to the zero of
since the output amplitude
is a complex function of
. 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
,
,
,
,
,
,
,
. The system response manifests as a monotone increasing function of
under the linear trichotomous noise disturbance. When the colored noise is in nonlinear form (
), the output amplitude of the system peaks at
, 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
has a great influence on the GSR phenomenon induced by damping order. Moreover, the output amplitude decreases with the increase in damping order for
and
, while stochastic resonance occurs at
for
, and a minimal resonance point
emerges when
.
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 at time , i.e., , . At 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 that is falling in , which is then compared against the three values in the line of the conditional probability matrix. If , then we accept the value of the noise , i.e., , and if we accept the value or , otherwise, it will be located at , i.e., ;
2. Taking the site at time as the starting location again, marked as , another new uniformly distributed random number between is generated, and we conduct a comparison operation on with the elements in the line of (3). If or then the value of noise at time is confirmed as or , otherwise we accept . By repeating this two-step procedure one can generate a sequence of random numbers switching between three values , and .
It should be noted that the time interval between the two steps (
) is always fixed and be equal to
, which should be much smaller than the correlation time of the noise (
). Let
,
, and the total time
,
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
is often accompanied by a longer residence time of a particular state on average, with a larger
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
sample results; the first moment of generating time series is determined to satisfy the zero-mean characteristic of trichotomous colored noise (
). 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
under different characteristic relaxation time, as shown in
Figure 14, wherein red markers represent the simulation results for
samples.
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
, the original second-order FDE (12) can be reduced to first-order equation set:
According to the Caputo fractional derivative approximation method used in the previous work [
43], (50) can be discretized as follows:
,
,
denotes the time step,
the total time,
,
,
,
, the initial condition in
is set as
.
,
is the trichotomous noise state at
.
is the state value of FBM at
, 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 (
), the larger
H is, the smoother the sample path is. The situation of
is opposite to that of
, moreover, the range of FGN decreases with the increase in Hurst index.
Select
and
, 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
s, with parameters
,
,
,
,
,
,
,
,
,
. 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
show stronger periodic characteristics than those of
. (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
in both cases. In the case of
and
, the numerical results of the amplitude at the external excitation frequency are
and
, respectively, and the theoretical results obtained by the analytical expression (35) are
and
, 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
samples,
,
,
,
represents the total steps. The first-order moment of the system response at time
can be given by the following formula:
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
where FFT denotes the fast Fourier transform.
Firstly, examining the influence of time step on the time consumption of the program, fix
,
and select step in the range
. 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
increases exponentially with the decrease in the time step
, and the approximate fitting relationship reads
, the calculation time of a single sample increases greatly especially for
, 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
, the simulation results based on (53) and theoretical results based on (35) are compared in
Figure 17b. When
decreases, the numerical simulation amplitude gradually approaches the theoretical result
, and tends to 1.603 without significant fluctuation when
, 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
, so time step
is selected in the following part.
Under the same parameter setting as
Figure 16,
,
, the relationship between
and the total simulation time
for
and
was examined by (53), depicted by
Figure 18a,b, respectively. When
, the simulation results differ greatly from the theoretical results due to the transient influence of the initial condition. When
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
and the theoretical results given by (35)
where
represents the single sample amplitude estimation of the
experiment, given by
In virtue of (54), the mean square error of simulation results corresponding to
and
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,
for
and
for
, which reflects the reliability of the algorithm (51). Under the same parameter conditions, the influence of the number
of the Monte Carlo experiment on the steady-state output amplitude
and mean square error
was investigated, and the simulation time was fixed
. Comparison between
(
Figure 19a,b) and
(
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
, which indicates that the algorithm in this paper (51) is basically consistent with the theoretical analysis result given by (35).