Fractional-Differential Models of the Time Series Evolution of Socio-Dynamic Processes with Possible Self-Organization and Memory

Abstract

This article describes the solution of two problems. First, based on the fractional diffusion equation, a boundary problem with arbitrary values of derivative indicators was formulated and solved, describing more general cases than existing solutions. Secondly, from the consideration of the probability schemes of transitions between states of the process, which can be observed in complex systems, a fractional-differential equation of the telegraph type with multiples is obtained (in time: β, 2β, 3β, … and state: α, 2α, 3α, …) using orders of fractional derivatives and its analytical solution for one particular boundary problem is considered. In solving edge problems, the Fourier method was used. This makes it possible to represent the solution in the form of a nested time series (one in time t, the second in state x), each of which is a function of the Mittag-Leffler type. The eigenvalues of the Mittag-Leffler function for describing states can be found using boundary conditions and the Fourier coefficient based on the initial condition and orthogonality conditions of the eigenfunctions. An analysis of the characteristics of time series of changes in the emotional color of users’ comments on published news in online mass media and the electoral campaigns of the US presidential elections showed that for the mathematical expectation of amplitudes of deviations of series levels from the size of the amplitude calculation interval (“sliding window”), a root dependence of fractional degree was observed; for dispersion, a power law with a fractional index greater than 1.5 was observed; and the behavior of the excess showed the presence of so-called “heavy tails”. The obtained results indicate that time series have unsteady non-locality, both in time and state. This provides the rationale for using differential equations with partial fractional derivatives to describe time series dynamics.

1. Introduction

When modeling the dynamics of various processes, three types of fractional-differential equations can be considered: equations with a fractional derivative by coordinate (in the case of one-dimensional space), equations with a fractional derivative by time, as well as mixed-type equations, including fractional operators both by coordinate and by time.

It should be noted that in different problems, depending on the nature of the processes under investigation, the concept of “coordinate” may have different meanings. For instance, in most physical studies, the coordinate, or in the case of three dimensions, coordinates, refer to spatial dimensions. When modeling the dynamics of time series, the coordinate is defined as the variable that describes the changes in the levels of the series and is measured in units corresponding to the nature of the observed process. For instance, if fractional-differential equations are used to analyze and model the dynamics of market indicators, the levels of the series represent the prices of stocks, commodities, currency ratios, etc.

Differential equation with fractional partial derivatives in the general case has the following form (this is a fractional diffusion equation):

$$\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}=D\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}$$

(1)

where α and β are the exponents of the fractional derivative (according to Caputo); $\rho (x,t)$ is the probability density function of observing state x at time t, depending on both time t and state x; and D is a constant coefficient. When describing diffusion, $\rho (x,t)$ represents the concentration of the diffusing substance, and D stands for the diffusion coefficient.

In general, the fractional derivative of order ν for the function ψ(x) is defined as follows [1]:

$$\frac{d^{\nu }\psi \left(x\right)}{d{x}^{\nu }}=\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }(-\nu )}·{\int }_a^x\frac{f\left(\xi \right)d\xi }{{\left\{x-\xi \right\}}^{\nu +1}},\nu <0\\ \frac{1}{\mathsf{\Gamma }(1-\nu )}·\frac{d}{dx}{\int }_a^x\frac{f\left(\xi \right)d\xi }{{\left\{x-\xi \right\}}^{\nu }},0\le \nu <1\\ \begin{array}{c}\frac{1}{\mathsf{\Gamma }(2-\nu )}·\frac{d^2}{d{x}^2}{\int }_a^x\frac{f\left(\xi \right)d\xi }{{\left\{x-\xi \right\}}^{\nu -1}},1\le \nu <2\\ \begin{array}{c}\frac{1}{\mathsf{\Gamma }(3-\nu )}·\frac{d^3}{d{x}^3}{\int }_a^x\frac{f\left(\xi \right)d\xi }{{\left\{x-\xi \right\}}^{\nu -2}}\\ \cdots \end{array},2\le \nu <3\end{array}\end{array}\right.$$

or

$$\frac{d^{\nu }\psi \left(x\right)}{d{x}^{\nu }}=\frac{1}{\mathsf{\Gamma }(1-\left[\nu \right])}·\frac{d^{\left(max\left\{\nu \right\}\right)}}{d{x}^{\left(max\left\{\nu \right\}\right)}}{\int }_a^x\frac{f\left(\xi \right)d\xi }{{\left\{x-\xi \right\}}^{\left[\nu \right]}}$$

where [ν] denotes the fractional part of the ν value in the exponent, and max{ν} represents rounding the fractional value of the ν to the nearest integer.

The definition of Caputo fractional derivatives differs from that of Riemann–Liouville. In the Caputo approach, the procedure initiates by differentiating the function with an integer order n, surpassing the non-integer order after rounding to the nearest integer (max{ν}). Subsequently, the obtained result experiences integration with an order of 1 − [ν].

At present, an analytical solution to the fractional diffusion Equation (1) has only been obtained for the case where $0<\beta \le 1$ and $1\le \alpha \le 2$ [2,3,4,5,6,7,8,9]. Solutions beyond the boundaries of these specified values of α and β are not presented in the literature.

For instance, in [4], a general solution to Equation (1) is described. Relying on the theory of Lie groups, the authors in [4] introduce a one-parameter family of scaling transformations: $\overline{x}={\lambda }^{a}x$; $\overline{t}={\lambda }^{b}t$; $\overline{\rho }(\overline{x},\overline{t})={\lambda }^c\rho (x,t)$, where a, b and c are constants, and $\lambda$ is a real parameter confined within an open interval $I$ containing $\lambda =1$.

By employing scaling transformations, it is possible to derive [4]

$$\frac{{\partial }^{\beta }\overline{\rho }(\overline{x},\overline{t})}{\partial t^{\beta }}-D\frac{{\partial }^{\alpha }\overline{\rho }\left(\overline{x},\overline{t}\right)}{\partial x^{\alpha }}={\lambda }^{c+b\beta }\frac{{\partial }^{\beta }\rho (\overline{x},\overline{t})}{\partial {\overline{t}}^{\beta }}-D{\lambda }^{c+a\alpha }\frac{{\partial }^{\alpha }\rho (\overline{x},\overline{t})}{\partial {\overline{x}}^{\alpha }}$$

If $a/b=\beta /\alpha$ and denoting $\gamma =c/b$, then by making the substitution $z=x{t}^{-\beta /\alpha }$ in Equation (1) and utilizing the Lie group method, it is possible to obtain a general scale-invariant solution ($\gamma \ge 0$) for this equation (given $0<\beta \le 1$ and $1\le \alpha \le 2$) [4]:

$$\rho \left(x,t\right)=t^{\gamma }\sum _{j=0}^n{C}_j·z^{\alpha -j}·W_{\left(-\beta ,1+\gamma -\beta +\frac{\beta }{\alpha }j\right),\left(\alpha ,1+\alpha -j\right)}(z^{\alpha }/D)$$

(2)

where $C_j$ are arbitrary real constants, $1\le j\le n$; $\gamma$ is the scaling coefficient; and $W_{\left(-\beta ,1+\gamma -\beta +\frac{\beta }{\alpha }j\right), \left(\alpha , 1+\alpha -j\right)}(z^{\alpha }/D)$ represents a Wright function of the form

$$W_{\left(-\beta ,1+\gamma -\beta +\frac{\beta }{\alpha }j\right),\left(\alpha ,1+\alpha -j\right)}\left(z^{\alpha }/D\right)=\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{z^{\alpha }/D\right\}}^k}{\mathrm{\Gamma }\left(1+\gamma -\beta +\frac{\beta }{\alpha }j-\beta k\right)\mathrm{\Gamma }\left(1+\alpha -j+\alpha k\right)}$$

(3)

$$\rho \left(x,t\right)=t^{\gamma }\sum _{j=0}^n{C}_j·z^{\alpha -j}·\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{z^{\alpha }/D\right\}}^k}{\mathrm{\Gamma }\left(1+\gamma -\beta +\frac{\beta }{\alpha }j-\beta k\right)\mathrm{\Gamma }\left(1+\alpha -j+\alpha k\right)}$$

(4)

It should be noted that a series of solutions is obtained for the case of integer α and $0<\beta \le 1$ [4]. In the case of $\alpha =2$, scale-invariant solutions of Equation (1) for the entire real line are expressed [4] using the Wright function as follows:

$$\rho \left(\left|x\right|,t\right)=t^{\gamma }\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{\left|z\right|/\sqrt{D}\right\}}^k}{\mathrm{\Gamma }\left(1+\gamma -\frac{\beta }{2}k\right)k!}$$

(5)

$$\rho \left(\left|x\right|,t\right)=t^{\gamma }\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{\left|z\right|/D\right\}}^k}{\mathrm{\Gamma }\left(1+\gamma -\beta k\right)k!}$$

(6)

$$\rho \left(\left|x\right|,t\right)=t^{\gamma }\sum _{k=0}^{\mathrm{\infty }}\frac{{\left|z\right|}^{\alpha -1}·{\left\{{\left|z\right|}^{\alpha }/D\right\}}^k}{\mathrm{\Gamma }\left(1+\gamma +\frac{\beta }{\alpha }-\beta \left[k+1\right]\right)\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}$$

(7)

The obtained solutions coincide, for $\gamma =0$, with known solutions presented in [6,7,8,9] (the condition of scale invariance transformation of the probability density function $\overline{\rho }(\overline{x},\overline{t})=\rho (x,t)$, is related to the preservation of total probability under any transformations).

For practical application of Equation (4) in the analysis and modeling of time series behavior, it is necessary to determine the constants $C_j$. In practice, this depends on the specific nature of the problem being solved to describe the dynamics of the time series. The determination of constants $C_j$ is a separate issue for discussion, similar to the convergence of the series in Equation (4).

Due to the complexity of the problems and substantial challenges in obtaining analytical solutions, numerical methods are frequently employed in practice. For instance, in [10], numerical methods were utilized to investigate solutions over an interval for an equation of the form

$$\frac{\rho \left(x,t\right)}{\partial t}=D\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }} \mathrm{at} 1<\alpha <2$$

Meanwhile, in [11], the solution over an interval for an equation of this form was obtained:

$$\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}=D\frac{{\partial }^2\rho \left(x,t\right)}{\partial x^2} \mathrm{at} 1<\beta \le 1$$

The function $\rho (x,t)$ can, for instance, be interpreted as the probability density of observing a particular value level within a time series at time t, if these levels can randomly vary over time.

In the case of a normal distribution law, the dispersion ${\sigma }^2\left(t\right)$ exhibits a linear relationship with time (${\sigma }^2\left(t\right)~t$) If there is a slower growth of ${\sigma }^2\left(t\right)$ concerning t (${\sigma }^2\left(t\right)~\sqrt[\beta ]{t}$, where $0<\beta <1$), such a process is classified as subdiffusion [12,13].

If there is a faster growth observed for ${\sigma }^2\left(t\right)$ concerning t (${\sigma }^2\left(t\right)~t^{\beta }$, where $1<\beta <2$), such a process is classified as superdiffusion [14]. The distribution function describing subdiffusive processes is obtained from the solution of a fractional-differential equation of the form $\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}=D\frac{{\partial }^2\rho \left(x,t\right)}{\partial x^2}$, and superdiffusion from the equation $\frac{\partial \rho \left(x,t\right)}{\partial t}=D\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}$, where α and β are the exponents of fractional derivatives, and D is a certain coefficient (diffusion coefficient). Additionally, processes of anomalous diffusion [15,16] can be observed, which is described by a mixed-type fractional-differential equation (Equation (1)).

Currently, the most developed and studied fractional-differential models are the ones created to describe various kinds of physical problems and are mainly devoted to the study of the processes of physical kinetics, abnormal diffusion [17] and relaxation in various environments [18,19,20,21,22,23,24], as well as wave processes [25,26,27,28].

Since this article will further consider the conclusion of the stochastic fractional-differential equation with multiples (in time: β, 2β, 3β …… and status: α, 2α, 3α ……), orders of fractional derivatives and the analytical solution for one particular boundary problem (with an arbitrary value of β and α = 1) is presented; then, as a comparison and representation of the novelty of the results obtained, we give a description of the solution of similar wave-type fractional-differential equations. Work [25] deals with the solution of the fractional Zener wave equation:

$${\nabla }^2\mathcal{U}\left(x,y,z,t\right)-\frac{1}{c_0^2}\frac{{\partial }^2\mathcal{U}\left(x,y,z,t\right)}{\partial t^2}+{\tau }_{\sigma }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}{\nabla }^2\mathcal{U}\left(x,y,z,t\right)-\frac{{\tau }_{ϵ}^{\beta }}{c_0^2}\frac{{\partial }^{\beta +2}\mathcal{U}\left(x,y,z,t\right)}{\partial t^{\beta +2}}=0$$

where $\mathcal{U}\left(x,y,z,t\right)$, in the case of elastic waves, is the amount of displacement of medium particles from the equilibrium position at a point with coordinates x, y, z at time t; C₀ is the wave velocity; and ${\tau }_{ϵ}$ и ${\tau }_{\sigma }$ are positive time constants. In this case, α and β are different values of fractional derivatives over time and operator ${\nabla }^2$ is not fractional.

This equation can be derived from the fractional ratio of stress–strain of Zener in the propagation of elastic waves and accounting and in the linearized conservation of mass and momentum. The Ziner wave equation allows us to describe three different attenuation modes with characteristics of the power law [25]. Models based on this equation have very important applications in medical elastography and acoustics.

Work [27] provides a detailed description of the solution to the Zener equation for the case of α = β (α and β can be derivatives of fractional positive numbers specifying time derivatives t). This result is very important since solving any fractional-differential equations in analytic form has significant problems.

A further development of the wave equation with a fractional derivative in time was obtained in work [28] to describe the behavior of waves in non-Newtonian liquids.

An example of an application of fractional-differential equations of the diffusion type to describe the dynamics of social processes is presented in work [29].

Fractional derivation by the state variable x and time derivative t allow you to describe non-local processes in which the transition to a certain state of the system (or process) x at time t depends not only on the local characteristics of the process or the behavior of the system in the vicinity of the point x under consideration. but also on the values of x obtained throughout the entire interval under study, that is, it globally depends on the distribution over all states of x and on the history of the process (memory) throughout the time t. Using fractional equations according to both time and coordinates allows for memory effects. The effect of non-locality in time t and non-locality in state x on the probability density of detecting a system or process in state x is qualitatively different. Non-locality in time affects the probability density at the initial moment of time, which can lead to self-organization, and non-locality in x affects the asymptotic behavior of the probability density at very late moments of time [12,13,14,15,16].

2. Data Processing and Analysis of Observed Time Series

In several studies [30,31,32,33], it was demonstrated that observed time series of social processes exhibit fractality, while the systems they describe showcase memory and self-organization. For instance, analyzing the dependence of the mean and dispersion of amplitude changes in time series on the interval of calculating these amplitudes reveals complex relationships. For example, their dispersion is dependent on the size of the “sliding” window as a fractional root, significantly differing from, for instance, a normal distribution law.

An analysis of the observed data can be used to justify the possibility of using fractional-differential equations in modeling complex processes.

We analyzed the time series of the emotional attitude of users of the “RIA Novosti” portal toward the news published during the day, for a period of 1460 days from 1 January 2019 to 31 December 2022. The following emotions were chosen as the subjects of this study: “like” and “dislike”.

In addition to analyzing the activity of users commenting on news on social media, an analysis of the preferences of voters of presidential electoral campaigns in the United States in 2012 and 2016 was carried out (data taken from the resource: http://www.realclearpolitics.com/epolls/, accessed on 27 April 2018).

For a preliminary analysis of the time series dynamics and determining their characteristics (for example, the possible presence of memory), the normalized range method of Hurst [34] can be employed. This method allows for the determination of their fractal dimension and classification of behavior type.

It is possible to normalize the range of sequence levels, $R$, using the standard deviation of these values, S, and represent their ratio (normalized range) as an equation: $R⁄S=C\tau H \left(log\right[R⁄S]=H\ast log[\tau ]+log[C\left]\right)$, where C is a constant, $\tau$ is the number of observations (levels of the sequence) comprising the considered time series, and the H exponent is known as the Hurst coefficient or exponent.

The presence of breakpoints in the $R⁄S \left(\tau \right)$ dependency may indicate the existence of characteristic time scales and/or periodicities. The value of the Hurst coefficient $H$ allows for the classification of time series according to the nature of their behavior [35].

For the time series, the Hurst coefficient for the “like” emotion turned out to be 0.22 and for “not like”, it was 0.24. For the election campaign for the US presidential election in 2012, the Hurst indicator for the temporary series of the preferences of Obama voters turned out to be 0.29 and for Romney voters, it was 0.22. In 2016, the Hurst’s indicator for a temporary series of the preferences of Clinton (Hillary) voters turned out to be 0.36 and for Trump voters, it was 0.30.

In all cases, the H value is significantly less than 0.5 and, therefore, the observed time series are anti-persistent (ergodic). Since the values of the Hurst coefficient are significantly different from 0.5, it follows from this that the structure of these series has fractality, and the processes described by it can have short-term memory [35].

The possible presence of memory should be taken into account in the model describing the dynamics of processes observed in complex social systems, for example, such as user activity when commenting on the news of online mass media or a change in voter preferences during electoral campaigns.

In order to justify the possibility of using the dynamics of processes in complex systems of fractional-differential equations for modeling, it is necessary to check whether the conditions of non-locality of behavior in time (variable t) and variable (x) describing the level of the series (for example, the proportion of “likes” and “dislikes” when commenting on news resources, or the proportion of voters with a preference for one candidate over another) are met.

For this, there were studies on the dependence of the mathematical expectation, variance, asymmetry (third moment of distribution) and excess (fourth moment of distribution) of the amplitudes of the deviations of the levels of time series from the calculation time intervals (dimensions of the “sliding window”) of these amplitudes.

The study of the behavior of “excesses” allows you to show the presence of so-called “heavy tails” (when the graph of the distribution function lies above the graph of the Gaussian distribution function). In practice, “light tails” can also be observed. The presence of “tails” other than the Gaussian distribution indicates the non-locality of the process over the variable x. The study of the behavior of expectation can show the presence of non-stationary (or vice versa, stationary), and the study of the behavior of variance can show non-locality over time (variable t). The identified non-localities can be described using fractional derivatives.

The observed value of mathematical expectation, dispersion, asymmetry and excess can be calculated by the following equations:

$$\mu \left(t\right)=\frac{\sum _{j=1}^N{x}_j\left(t\right)}{{\sum }_{l=1}^M{n}_l}$$

$${\sigma }^2(t)=\frac{\sum _{j=1}^N{\left\{x_j\left(t\right)-\mathsf{\mu }(t)\right\}}^2}{{\sum }_{l=1}^M{n}_l}$$

$$As(t)=\frac{\sum _{j=1}^N{\left\{x_j\left(t\right)-\mathsf{\mu }(t)\right\}}^3}{{\sigma }^3{\sum }_{l=1}^M{n}_l}$$

$$Ex\left(t\right)=\frac{\sum _{j=1}^N{\left\{x_j\left(t\right)-\mathsf{\mu }\left(t\right)\right\}}^4}{{\sigma }^4{\sum }_{l=1}^M{n}_l}-3$$

where ${\sum }_{l=1}^M{n}_l$ is calculated from the number of $n_l$ amplitudes.

For example, when executing the normal distribution law for a stationary time series $\left(x,t\right)=\frac{1}{2\sqrt{\pi \mathrm{D}t}}{·e}^{-\frac{x^2}{4·D·t}}$, the amount of expectation $\mu (t)$ would have to be zero:

$$\mu (t)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x·\rho \left(x,t\right)dx=\frac{1}{2\sqrt{\pi \mathrm{D}t}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{x·e}^{-\frac{x^2}{4·D·t}}dx=0$$

where x is the value of the amplitude, and t is the interval of its calculation time (the value of the “sliding window”).

When fulfilling the normal law, a linear dependence on the “sliding window” value would have to be observed for amplitude dispersion:

$${\sigma }^2\left(t\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{x}^2\rho \left(x,t\right)dx=\frac{1}{2\sqrt{\pi \mathrm{D}t}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{x}^2{e}^{-\frac{x^2}{4·D·t}}dx=2Dt$$

Processing shows that

• Mathematical expectations of changes in the amplitudes of time series levels depend on the time interval for calculating these changes (the “sliding window”). This indicates the unsteadiness of the time series under consideration and the inability to describe their parameters using the normal distribution law.
• The values of the dispersion of the amplitudes of the change in user activity on news commentary depends on the time interval for calculating these amplitudes (the “sliding window”) in a complex way: it is proportional to the fractional degree from the time interval from which they are calculated. The fractional dependence on the time interval indicates that the studied processes have a non-locality in time t (i.e., have a consequence or memory).
• The studies of the excess distribution of amplitudes show the presence of the so-called “heavy tails”, which are significantly greater than the normal distribution (where the excess is 3). With significant positive deviation values, the distribution function decreases more slowly at a distance from the average than with small values. If the deviation is more than three, the distribution density plot will be higher than the normal distribution plot and lower than three. This indicates that the processes under consideration not only have non-locality in time t, but also non-locality in state—let us designate it as x.

3. Setting a Study Objective

Studies devoted to the analysis of the dynamics and forecasting of the development of processes in complex systems are of high importance both from the point of view of science and from a practical perspective. A variety of studies have been conducted in this area [36,37,38,39,40,41,42,43,44,45]. Of these, we highlight the use of the theory and methods of neural networks [36,37,38,39,40], the use of the fuzzy logic apparatus [41], and the creation of non-parametric models based on chaos theory and methods of supporting regression vectors [42]. Also, the development of sets of rules based on genetic algorithms [43,44], as well as the use of self-organizing adaptive models [45] should be distinguished.

However, these studies have not taken into account the problems associated with studying processes in complex systems. Complex systems such as socioeconomic systems can be defined as structures involving at least one human element. On the one hand, they are characterized by stochasticity due to the influence of various random factors. In addition, there is uncertainty associated with the sometimes irrational behavior of people. However, on the other hand, the presence of the human factor also creates conditions for self-organization in such systems and can enable the storage of a memory of previous states of the processes occurring in these systems.

The joint action of these factors leads to the emergence of organized complexity, or “emergence”. Its origin cannot be reduced to a simple addition of element characteristics; it is the result of the formation of system relationships and adaptive redistribution of functions between elements.

The presence of all the listed features of complex systems leads to the need to find new approaches and methods for analyzing and modeling the dynamics of the processes observed in them.

In conclusion, we note that non-local processes in which the transition to a given state x depends not only on the local characteristics or behavior of the system in the vicinity of the point in question at a given moment in time, but also from the values taken by it over the entire studied interval of the values of the series levels at any time, which can be described on the basis of differential equations with fractional partial derivatives along t and x (and not only of the diffusion type).

Taking into account all the above, we can conclude that the development of new methods for describing the dynamics of stochastic processes based on differential equations with fractional partial derivatives may allow us to take into account the presence of self-organization and memory (not Markov processes that take into account the consequences).

Sociodynamic processes are widespread and have a wide variety of manifestations. One of the most important objects where sociodynamic processes are observed is the Internet, particularly, news and blogs, under which, users of social networks and mass media leave their comments or express their emotional attitude towards them. These are one of the most important online phenomena and can act as indicators of public opinion and mood.

4. Theoretical Model

4.1. Setting and Solving the Boundary Problem for Arbitrary Values of Derivatives of Diffusion-Type Fractional-Differential Equations

Given that the processes observed in complex systems have features that vary over a wide range of values, for the construction of the model, it is impossible to consider only the non-locality of the state of the system x (variable describing the level of the time series) or only non-locality of time t. In addition, it is necessary to consider a more general case with arbitrary fractional values of α and β, and not just the case of $0<\beta \le 1$ and $1\le \alpha \le 2$. Or, we should consider differential equations containing multiple fractional derivatives, both in time and in state (this will be discussed later in one of the paragraphs in the section on the theoretical model).

To analyze and model the time series of processes observed in complex systems, you can consider the edge problem of the form

$$\rho \left(0,t\right)=\rho \left(L,t\right)=0$$

with the initial condition given by the delta function:

$$\rho \left(x,0\right)=\delta (x-0)$$

The choice of such boundary conditions is due to the fact that when considering, for example, the dynamics of the amplitudes of change in the levels of the time series (in economics, this is called volatility), the state of x can change within a certain permissible range of values from 0 to L, and at the time of the beginning of observations ($t=0$), any change must be equal to 0. The $\rho$ function $(x, t)$ can be considered as the probability density for the amplitude x over the time interval t.

Work [46] considers the solution of the boundary problem (on the segment [0, 1]) based on Equation (1) at $\beta =1$: $\rho \left(0,t\right)=\rho \left(1,t\right)=0$ with the initial condition given by the delta function $\rho \left(x,0\right)=\delta (x-0)$.

$$\rho \left(x,t\right)=\sum _{n=1}^{\mathrm{\infty }}{C}_n{e}^{{\lambda }_{n}D·t}{·x}^{\alpha -1}\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{x}^{\alpha }\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}$$

(8)

where $C_n$ are the coefficients of the Fourier series and ${\lambda }_n$ are the zeros of the Leffler-Mittag function: $E_{\alpha ,\alpha }({\lambda }_k)=\sum _{k=0}^{\infty }\frac{{\left\{{\lambda }_k\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}$.

To determine the $C_n$ coefficients [46],

$$\phi \left(x\right)=\rho \left(x,0\right)=\sum _{n=1}^{\mathrm{\infty }}{C}_n{x}^{\alpha -1}\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{x}^{\alpha }\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}$$

Note that the function system ${\theta }_n={\left\{x^{\alpha -1}\sum _{k=0}^{\infty }\frac{{\left\{{\lambda }_n{x}^{\alpha }\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}\right\}}_{k=1}^{\infty }$ forms a basis in L₂ (0, 1) [47]. Since it is not orthogonal, it is necessary to also create a system of functions

$${\varphi }_n(x)={\left\{(1-{x)}^{\alpha -1}\sum _{k=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{(1-x)}^{\alpha }\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}\right\}}_{n=1}^{\mathrm{\infty }}$$

that are biorthogonal to the ${\theta }_n$ system [48].

After that, unknown $C_n$ coefficients can be determined through the dot product ${\varphi }_n$ and $\phi \left(x\right)$ [46]: $\left(\phi \left(x\right)·{\varphi }_n\left(x\right)\right)$ using the initial condition.

It should be noted that the determination of the $C_n$ coefficients and zeros of the Leffler-Mittag function performed in work [46] was not carried out, but the existence of a solution to the boundary problem formulated by the authors was investigated in a general form.

Consider the solution of the formulated boundary problem for the case of arbitrary fractional values of α and β (not only the case of $\beta =1$ or $0<\beta \le 1$ and $1\le \alpha \le 2$).

In this case, the equations obtained in [2,3,4,5,6,7,8,9,10,11,46] cannot be used and it is necessary to look for other solutions. In this regard, we suggest the following. First, using the Fourier method, imagine the $\rho \left(x,t\right)$ as $\rho \left(x,t\right)=X\left(x\right)·T\left(t\right)$. After substituting into Equation (1) and separating the variables, we obtain

$$\frac{{\partial }^{\alpha }X\left(x\right)}{\partial x^{\alpha }}+\frac{\lambda }{D}X\left(x\right)=0$$

(9a)

$$\frac{{\partial }^{\beta }T(t)}{\partial t^{\beta }}+\lambda T\left(t\right)=0$$

(9b)

where λ is some constant that appears due to the fact that after substituting $\rho \left(x,t\right)=X\left(x\right)·T\left(t\right)$ into Equation (1), the left side will depend only on the variable t, and the right only on the variable x.

To define the X(x) function, we encountered an problem finding eigenvalues. The solution to this problem is presented in the works [21,49,50]. In these studies, a solution has only been proven for eigenvalues of ${\lambda }_n$ that are function zeros:

$$E_{\alpha ,\alpha }({\lambda }_n{x}^{\alpha })=\sum _{n=0}^{\mathrm{\infty }}\frac{{\left\{\frac{{\lambda }_n}{D}·x^{\alpha }\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}$$

There are eigenfunctions that are particular solutions of Equation (9a). Zeros can be found using the given boundary conditions $\rho \left(0,t\right)=\rho \left(L,t\right)=0$. The eigenfunctions are equal to

$$X_n={\left\{x^{\alpha -1}\sum _{n=0}^{\mathrm{\infty }}\frac{{\left\{\frac{{\lambda }_n}{D}·x^{\alpha }\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}\right\}}_{n=1}^{\mathrm{\infty }}$$

where ${\lambda }_n$ are the zeros of the Leffler-Mittag function.

Using the Laplace transform method for the function $T\left(t\right)$, one can write

$$p^{\beta }\overline{G\left(p\right)}-\frac{1}{p^{1-\beta }}-\lambda \overline{G\left(p\right)}=0$$

Next,

$$\overline{G\left(p\right)}=\frac{1}{p^{1-\beta }}·\frac{1}{p^{\beta }-\lambda }=\frac{1}{p}·\frac{1}{1-\frac{\lambda }{p^{\beta }}}=\frac{1}{p}\sum _{q=0}^{\mathrm{\infty }}{\left\{-\frac{\lambda }{p^{\beta }}\right\}}^q=\sum _{q=0}^{\mathrm{\infty }}\frac{{(-1)}^q{\lambda }^q}{p^{\beta +1}}$$

Let us make the inverse Laplace transform and go from p to t:

$$T_n(t)=\sum _{q=0}^{\mathrm{\infty }}\frac{{(-1)}^q{\lambda }_n^q{t}^{\beta q}}{\mathrm{\Gamma }\left(\mathsf{\beta }q+1\right)}$$

Thus, the general solution to the formulated edge problem can be written as

$$\rho \left(x,t\right)=x^{\alpha -1}\sum _{n=0}^{\mathrm{\infty }}{C}_n\sum _{q=0}^{\mathrm{\infty }}\frac{{(-1)}^q{\lambda }_n^q{t}^{\beta q}}{\mathrm{\Gamma }\left(\mathsf{\beta }q+1\right)}\frac{{\left\{\frac{{\lambda }_n}{D}·x^{\alpha }\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}$$

(10)

To find the zeros of the Mittag-Leffler function, you can use the boundary condition $\rho \left(L,t\right)=0$ and the results obtained in [51,52,53,54,55], which examined the behavior of the zeros of this function for various α values, which allows us, taking into account the boundary condition, to obtain the following result.

In the case of $\alpha <2$, $\alpha \in C$, zeros are ${\lambda }_n$ and determined using the equation

$${\lambda }_n^{\pm }=e^{\pm i\frac{\pi }{2}\alpha }{\left\{\frac{2\pi n}{L}\right\}}^{\alpha }\left\{1+O\left\{\frac{\log n}{n}\right\}\right\}$$

When α = 2,

$${\lambda }_n={\left\{\frac{\pi n}{L}\right\}}^2\left\{1+O\left\{\frac{1}{n}\right\}\right\}$$

In the case of $\alpha >2$, $\alpha \in \mathrm{C}$, all $z_n$ functions are large enough modulo zeros that $E_{\alpha ,\alpha }\left(z\right)=\sum _{n=0}^{\infty }\frac{{\left\{z\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}$ and are described by the equation

$$z_n={\lambda }_n{L}^{\alpha }=-{\left\{\frac{\pi }{\mathrm{s}\mathrm{i}\mathrm{n}\left[\frac{\pi }{\alpha }\right]}\left\{n+1/2+\frac{\alpha -1}{\alpha }\right\}+{\mathsf{\Omega }}_n\right\}}^{\alpha }$$

At 6 > α > 2, ${\mathsf{\Omega }}_n=O\left\{n^{-\alpha ·\tau }·e^{-\pi n·ctg\left[\frac{\pi }{\alpha }\right]}\right\}$

$$\tau =\left\{\begin{array}{c}\frac{1}{\alpha }, n\ne 0\\ \frac{1+\alpha }{\alpha }, n=0\end{array}\right.$$

At α = 6, ${\mathsf{\Omega }}_n=O\left\{e^{-\pi n·ctg\left[\frac{\pi }{\alpha }\right]}\right\}$

At α > 6, ${\mathsf{\Omega }}_n=O\left\{\frac{e^{-\pi n·\left[cos\left[\frac{\pi }{\alpha }\right]-cos\left[\frac{3\pi }{\alpha }\right]\right]}}{sin\left[\frac{\pi }{\alpha }\right]}\right\}$

Next, consider the definition of Fourier series coefficients $C_n$ through the dot product $\left({\phi }_n\left(x\right)·{\varphi }_m\left(x\right)\right)$, where

$${\phi }_n\left(x\right)=C_n{x}^{\alpha -1}·\sum _{n=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{x}^{\alpha }\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}$$

$${\varphi }_m(x)=(L-{x)}^{\alpha -1}\sum _{m=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{(L-x)}^{\alpha }\right\}}^m}{\mathrm{\Gamma }\left(\alpha \left[m+1\right]\right)}$$

Using the biorthogonality property of functions $\phi \left(x\right)$ and ${\varphi }_n$ in $L_2(0, L)$,

$$\underset{0}{\overset{L}{\int }}{\phi }_n\left(x\right)·{\varphi }_m\left(x\right)dx=\left\{\begin{array}{c}1, n=m\\ 0,n\ne m\end{array}\right.$$

Next, we assume that with a weight equal to 1 on the segment [0, L], the following condition should be met:

$$\underset{0}{\overset{L}{\int }}{\phi }_n\left(x\right)·{\varphi }_m\left(x\right)dx=\underset{0}{\overset{L}{\int }}\delta \left(x-0\right)·{\varphi }_m\left(x\right)dx$$

This allows you to define $C_n$:

$$C_n=\frac{L^{\alpha -1}\sum _{m=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{L}^{\alpha }\right\}}^k}{\mathrm{\Gamma }\left(\alpha \left[k+1\right]\right)}}{{\int }_0^L{x}^{\alpha -1}·(L-{x)}^{\alpha -1}\sum _{n=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{x}^{\alpha }\right\}}^n}{\mathrm{\Gamma }\left(\alpha \left[n+1\right]\right)}·\sum _{m=0}^{\mathrm{\infty }}\frac{{\left\{{\lambda }_n{(L-x)}^{\alpha }\right\}}^m}{\mathrm{\Gamma }\left(\alpha \left[m+1\right]\right)}dx}$$

4.2. Derivation of Telegraph-Type Fractional-Differential Equation with Multiple Orders of Fractional Derivatives from Consideration of Probability Schemes of Transitions between Process States

To derive the basic equation of our proposed model, we can use the approach to describe the stochastic dynamics of processes, taking into account memory, as well as possible self-organization, which we developed earlier and is described in works [30,31].

This method allows us, based on the schemes of probabilistic transitions between states, to formulate a boundary value problem regarding the probability of achieving any of the states x as a function of time t. One can then consider solving this problem (to obtain a theoretical approximating distribution function) based on a model that takes into account the memory about previous states and their potential self-organization.

The essence of this approach is as follows. Let us designate the current state as x_i (process state). Suppose the state change time interval is τ (extremely small). It is assumed that during this time period, the τ state of the system may increase by a value of ε (indicating an increasing trend) or decrease by a value of ξ (indicating a decreasing trend). Let us represent the entire set of states as X. The state observed at time t can be denoted as $x_i$ ($x_{i}X$), where x represents the level of values in the time series describing the observed process. We express the value of the current time as $t=h\tau$, where h is the number of transition steps between states (the transition process between states becomes quasi-continuous with an infinitesimal time interval τ), and h takes values of 0, 1, 2, 3, N. The current state $x_i$ at step h can increase by some amount of $\epsilon$ or decrease by an amount of $\xi$ after the transition in step ($h+1$), and accordingly, turn out to be equal to ($x_i+\epsilon$), or ($x_i-\xi$).

Consider the concept of the probability of finding a process in a certain state. Suppose, after a certain number of steps h in the described process, it can be argued that

• $P(x-\epsilon ,h)$—probability that it is in the state of $(x-\epsilon )$;
• $P(x,h)$—probability that it is in the state of $x$;
• $P(x+\xi ,h)$—probability that it is in the state of $(x+\xi )$.

Following each step, the state $x_i$ (hereinafter, the index i may be omitted for brevity) can change by an amount of $\epsilon$ or $\xi$.

The probability of $P(x, h+1)$—indicating the likelihood of the process being in the state x at the next $(h+1)$ step—will be influenced by multiple transitions (see Figure 1):

$$P(x,h+1)=P(x-\epsilon ,h)+P(x+\xi ,h)-P(x,h)$$

(11)

Let us explain the Expression (11) and the diagram shown in Figure 1. The probability of transition to state x in step h, denoted as $P(x, h+1)$, is determined by the sum of the probabilities of transition to this state from states $(x-\epsilon ): P(x-\epsilon ,h)$ and $(x+\xi ): P(x+\xi ,h)$ where the system was at step h. From this sum, we subtract the probability of the system transitioning (P(x, h)) from state x (where it was at step h) to any other state at step h + 1.

In this context, we consider a Markov continuous process devoid of state memory. However, in real conditions, it is possible to store information about the previous state. To take into account memory (not Markov processes), we will determine the probabilities $P(x–\epsilon , h)$, $P(x+\xi ,h)$ and $P(x,h)$ through the states in the previous h − 1 step. In this case, the following algebraic equation can be obtained for the probability of transition:

$$P(x,h+2)=P(x-2\epsilon ,h)+P(x+2\xi ,h)+P(x,h)+2\left\{P\right(x-[\epsilon -\xi ],h)-P(x-\epsilon ,h)-P(x+\xi ,h\left)\right\}$$

(12)

In this instance, the parameter h is augmented by the quantity m = 2.

Having carried out the necessary mathematical actions and considering that $t=h·\tau$ where t is the process time, h is the step number, $\tau$ is the duration of one step can be obtained for any arbitrary value m, the following recurrent expression is the probability $P (x,t+m\tau )$ that the process state, for some time t, at memory depth m, will be equal to x:

$$P\left(x,t+m\tau \right)=\left\{\begin{array}{c}{\displaystyle \sum _{k,l=0}^m}{\left(-1\right)}^{m-k-l}\frac{m!·\mathrm{P}\left(x-\left[k·\epsilon -l·\xi \right], t\right)}{k!·l!\left(m-k-l\right)!}, \mathsf{п}\mathsf{p}\mathsf{и} m-k-l\ge 0\\ 0, \mathsf{п}\mathsf{p}\mathsf{и} m-k-l<0\end{array}\right.$$

(13)

If we differentiate Equation (13) by the variable x, we obtain the equation for the probability density:

$$\rho \left(x,t+m\tau \right)=\left\{\begin{array}{c}{\displaystyle \sum _{k,l=0}^m}{\left(-1\right)}^{m-k-l}\frac{m!·\rho \left(x-\left[k·\epsilon -l·\xi \right], t\right)}{k!·l!\left(m-k-l\right)!}, \mathsf{п}\mathsf{p}\mathsf{и} m-k-l\ge 0\\ 0, \mathsf{п}\mathsf{p}\mathsf{и} m-k-l<0\end{array}\right.$$

(14)

In order to take into account the previously described properties of the observed processes, it is necessary to move from the obtained algebraic equation for probability density to a differential equation with derivatives of fractional orders (derivatives of Caputo). To accomplish this, we conduct the corresponding decompositions of the terms of this equation into a Taylor series over the derivatives of fractional order [16] in the vicinity of the transition point.

$$\mathsf{\rho }\left(x,t+m\tau \right)=\mathsf{\rho }\left(x,t\right)+\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\rho \left(x,t\right)}{\partial t^{2\alpha }}+\cdots$$

$$\mathsf{\rho }\left(x-\left[k·\epsilon -l·\xi \right],t\right)=\mathsf{\rho }\left(x,t\right)+{(-1)}^{\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}·\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}+{\left(-1\right)}^{2\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{2\alpha }}{\mathrm{\Gamma }\left(2\alpha +1\right)}·\frac{d^{2\alpha }\mathsf{\rho }\left(x,t\right)}{d{x}^{2\alpha }}+\cdots$$

Note that for fractional derivatives of β and α, the equation has multiples of orders β, 2β, 3β and α, 2α, 3α.

We substitute the corresponding expansions into Equation (14) and obtain, taking into account no more than the second derivatives (not to be confused with the order), the following differential equation for changing the density of the probability of detecting a process in some state x depending on the value of time t and memory depth m:

$$\mathsf{\rho }\left(x,t\right)+\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\mathsf{\rho }\left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\mathsf{\rho }\left(x,t\right)}{\partial t^{2\beta }}=\mathsf{\rho }\left(x,t\right)+\sum _{k,l=0}^m\frac{{\left(-1\right)}^{m-k-l}·m!}{k!·l!\left(m-k-l\right)!}\left\{\begin{array}{c}{(-1)}^{\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}·\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}+\\ +{(-1)}^{2\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{2\alpha }}{\mathrm{\Gamma }\left(2\alpha +1\right)}·\frac{d^{2\alpha }\mathsf{\rho }(x,t)}{d{x}^{2\alpha }}\end{array}\right\}$$

Or

$$\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\rho \left(x,t\right)}{\partial t^{2\beta }}=\sum _{k,l=0}^m\frac{{\left(-1\right)}^{m-k-l}·m!}{k!·l!\left(m-k-l\right)!}\left\{\begin{array}{c}{(-1)}^{\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}·\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}+\\ +{(-1)}^{2\alpha }\frac{{\left[k·\epsilon -l·\xi \right]}^{2\alpha }}{\mathrm{\Gamma }\left(2\alpha +1\right)}·\frac{d^{2\alpha }\rho (x,t)}{d{x}^{2\alpha }}\end{array}\right\}$$

(15)

The resulting equation contains two derivatives in time and two in state, which allows you to characterize it as a telegraph-type equation. It is interesting in that it allows you to take into account previous states due to different m values and the influence of the fractional derivatives of α and β of different orders.

For a member of a view equation $\frac{{\partial }^{\beta }\mathsf{\rho }(x,t)}{\partial t^{\beta }}$, if the rate of change of process states changes over time, then a term of the equation of the form $\frac{{\partial }^{2\beta }\mathsf{\rho }\left(x,t\right)}{\partial t^{2\beta }}$ can be, by analogy using physical kinetics, seen as acceleration. In this interpretation, the term $\frac{{\partial }^{2\alpha }\mathsf{\rho }(x,t)}{\partial x^{2\alpha }}$ of the equation takes into account random transitions (diffusion wandering of the state of the system), and the term $\frac{{\partial }^{\alpha }\rho (x,t)}{\partial x^{\alpha }}$ will describe ordered transitions (trend or demolition), for example, or when the value of the state increases ($\epsilon >\xi$), or decreases ($\epsilon <\xi$).

A significant difference between Equation (15) and fractional-diffusion and fractional wave equations is that the resulting equation, being stochastic (derived from the consideration of probability schemes of transitions between states), contains multiple orders of fractional derivatives, both in time (β, 2β) and state (α and 2α). This is a novelty that allows you to expand the class of fractional-differential equations and their application to describe the dynamics of complex systems. The Zener equation [25] contains only fractional time operators, and the fractional diffusion equations do not contain multiple derivatives.

Solving Equation (15) generally requires a separate study, and a simpler case with arbitrary values can be considered to begin with β, $\alpha =1$ and $\epsilon =\xi$.

Considering that

$$\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!}{k!·l!\left(m-k-l\right)!}=1$$

$$\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!}{k!·l!\left(m-k-l\right)!}·k=\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!}{k!·l!\left(m-k-l\right)!}·l=m$$

$$\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!·k^2}{k!·l!\left(m-k-l\right)!}=\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!·l^2}{k!·l!\left(m-k-l\right)!}=m^2$$

$$\sum _{k,l=0}^m{\left(-1\right)}^{m-k-l}\frac{m!·k·l}{k!·l!\left(m-k-l\right)!}=m\left(m-1\right)$$

Then,

$$\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\rho \left(x,t\right)}{\partial t^{2\beta }}=\frac{1}{2}\left\{m^2{\epsilon }^2-2m\left(m-1\right)\epsilon \xi +m^2{\xi }^2\right\}\frac{{\partial }^2\rho \left(x,t\right)}{\partial x^2}-m[\epsilon -\xi ]\frac{\partial \rho (x,t)}{\partial x}$$

(16)

If any state transitions have the same value, then $\epsilon =\xi$.

$$\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\rho \left(x,t\right)}{\partial t^{2\beta }}=m{\epsilon }^2\frac{{\partial }^2\rho \left(x,t\right)}{\partial x^2}$$

(17)

4.3. Setting and Solving the Boundary Problem for a Particular Case using Telegraph-Type Fractional-Differential Equation (Arbitrary $\beta$ and $\alpha =2$)

Consider for Equation (17) the solution of the previously formulated boundary problem $\rho \left(0,t\right)=\rho \left(L,t\right)=0$ with the initial condition given by the delta function: $\rho \left(x,0\right)=\delta (x-0)$.

Using the Fourier method, we present the $\rho \left(x,t\right)$ as $\rho \left(x,t\right)=X\left(x\right)·T\left(t\right)$. After substituting Equation (17) and separating the variables, we obtain

$$\frac{{\partial }^{2}X\left(x\right)}{\partial x^2}+\frac{\lambda }{m{\epsilon }^2}X\left(x\right)=0$$

(18a)

$$\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}·\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}·\frac{{\partial }^{2\beta }\rho \left(x,t\right)}{\partial t^{2\beta }}+\lambda T\left(t\right)=0$$

(18b)

where λ is some constant that appears as a result of substituting ρ (x,t) = X(x) ∙ T(t) into Equation (17). As a result of this substitution, the left side becomes dependent only on the variable t, and the right side only on the variable x. Then, you can find partial solutions of Equation (18a,b), and the general solution is represented in the form of a Fourier series.

To determine the function $X\left(x\right)$, we must solve the eigenvalue problem (Sturm–Liouville problem). This solution under given boundary conditions is well known and has the form $X_n\left(x\right)=C_n\mathrm{s}\mathrm{i}\mathrm{n}\left({\lambda }_{n}x\right)$, where λ_n are eigenvalues (${\lambda }_n=\pi n/L, n=\mathrm{0,1},\mathrm{2,3}\dots$).

To find a solution for $T_n\left(t\right)$, we perform the Laplace transform according to t:

$$\frac{{\left\{m\tau \right\}}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}\left\{p^{\beta }\overline{G\left(p,x\right)}-\frac{1}{p^{1-\beta }}\right\}+\frac{{\left\{m\tau \right\}}^{2\beta }}{\mathrm{\Gamma }(2\beta +1)}\left\{p^{2\beta }\overline{G\left(p,x\right)}-\frac{1}{p^{1-2\beta }}\right\}+\lambda \overline{G\left(p,x\right)}=0$$

Next,

$$\overline{G\left(p,x\right)}=\frac{\frac{{\left\{m\tau \right\}}^{2\beta }}{p^{1-2\beta }·\mathrm{\Gamma }(2\beta +1)}+\frac{{\left\{m\tau \right\}}^{\beta }}{p^{1-\beta }·\mathrm{\Gamma }(\beta +1)}}{\frac{{\left\{m\tau \right\}}^{2\beta }{p}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}+\frac{{\left\{m\tau \right\}}^{\beta }{p}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}+\lambda }=\frac{1}{p}·\frac{a_1\left(p\right)}{a_1\left(p\right)+\lambda }=\frac{1}{p}·\frac{1}{1+\frac{\lambda }{a_1\left(p\right)}}$$

$$a_1(p)=\frac{{\left\{m\tau \right\}}^{2\beta }{p}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}+\frac{{\left\{m\tau \right\}}^{\beta }{p}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}$$

Let us designate $z=\frac{\lambda }{a_1\left(p\right)}$ and if $0\le z<1$, we will carry out binomial decomposition:

$$\frac{1}{1-\left(-z\right)}={\sum }_{q=0}^{\mathrm{\infty }}(-{z)}^q={\sum }_{q=0}^{\mathrm{\infty }}{\left\{-\frac{\lambda }{a_1\left(p\right)}\right\}}^q={\sum }_{q=0}^{\mathrm{\infty }}{\left(-1\right)}^q{\left\{\frac{\lambda }{a_1\left(p\right)}\right\}}^q=$$

$$=\sum _{q=0}^{\mathrm{\infty }}{(-1)}^q{\left\{\frac{\lambda }{\frac{{\left\{m\tau \right\}}^{2\beta }{p}^{2\beta }}{\mathrm{\Gamma }\left(2\beta +1\right)}+\frac{{\left\{m\tau \right\}}^{\beta }{p}^{\beta }}{\mathrm{\Gamma }\left(\beta +1\right)}}\right\}}^q=\sum _{q=0}^{\mathrm{\infty }}{(-1)}^q\frac{{\lambda }^q}{{\frac{{\left\{m\tau \right\}}^{2\beta q}{p}^{\beta q}}{{\left\{\mathrm{\Gamma }\left(2\beta +1\right)\right\}}^q}\left\{p^{\beta }+\frac{\mathrm{\Gamma }\left(2\beta +1\right)}{{\left\{m\tau \right\}}^{\beta }\mathrm{\Gamma }\left(\beta +1\right)}\right\}}^q}$$

After substituting the result into the equation for $\overline{G\left(p,x\right)}$, we obtain

$$\overline{G\left(p,x\right)}=\sum _{q=0}^{\mathrm{\infty }}\frac{{\left(-1\right)}^q{\lambda }^q{\left\{\mathrm{\Gamma }\left(2\beta +1\right)\right\}}^q}{{\left\{m\tau \right\}}^{2\beta q}}·\frac{p^{-\beta q-1}}{{\left\{p^{\beta }+\frac{\mathrm{\Gamma }\left(2\beta +1\right)}{{\left\{m\tau \right\}}^{\beta }\mathrm{\Gamma }\left(\beta +1\right)}\right\}}^q}=\sum _{q=0}^{\mathrm{\infty }}\frac{{(-1)}^q{\lambda }^q{\left\{\mathrm{\Gamma }\left(2\beta +1\right)\right\}}^q}{{\left\{m\tau \right\}}^{2\beta q}}·\frac{p^{-\beta q-1}}{{\left\{p^{\beta }+\Omega \right\}}^q}$$

$$\frac{\mathrm{\Gamma }\left(2\beta +1\right)}{{\left\{m\tau \right\}}^{\beta }\mathrm{\Gamma }\left(\beta +1\right)}=\Omega$$

We perform the inverse Laplace transform and move from the image $\frac{p^{-\beta q-1}}{{\left\{p^{\beta }+\Omega \right\}}^q}$ by p to its original by t. Consider the integral for this:

$$\begin{gathered} {\int }_0^{\mathrm{\infty }}{e}^{-pt}·t^{b-1}·\sum _{j=0}^{\mathrm{\infty }}\frac{(-\Omega {t^{\mu })}^j}{\mathrm{\Gamma }\left(\mu j+b\right)}dt=\sum _{j=0}^{\mathrm{\infty }}\frac{{\left(-\Omega \right)}^j}{\mathrm{\Gamma }\left(\mu j+b\right)}·{\int }_0^{\mathrm{\infty }}{t}^{\mu j+b-1}·e^{-pt}dt= \\ \sum _{j=0}^{\mathrm{\infty }}\frac{{\left(-\Omega \right)}^j}{\mathrm{\Gamma }\left(\mu j+b\right)}·{\int }_0^{\mathrm{\infty }}\frac{1}{p}·\frac{y^{\mu j+b-1}}{p^{\mu j+b-1}}·e^{-y}dy=\sum _{j=0}^{\mathrm{\infty }}\frac{{\left(-\Omega \right)}^j}{p^{\mu j+b}}=\sum _{j=0}^{\mathrm{\infty }}\frac{1}{p^b}·{\left\{-\frac{\Omega }{p^{\mu }}\right\}}^j= \\ \frac{1}{p^b}·\frac{1}{1+\frac{\Omega }{p^{\mu }}}=\frac{p^{\mu -b}}{p^{\mu }+\Omega } \end{gathered}$$

(19)

where $\sum _{j=0}^{\infty }\frac{(-\Omega {t^{\mu })}^j}{\Gamma \left(\mu j+b\right)}$ is the Mittag–Leffler function, μ and b are some real numbers greater than 0, and t is a variable. When calculating the integral ${\int }_0^{\infty }{t}^{\mu j+b-1}·e^{-pt}dt$, replacement was used: pt = y. It should also be borne in mind that

$${\int }_0^{\mathrm{\infty }}{y}^{\mu q+b-1}·e^{-y}dy=\mathrm{\Gamma }\left(\mu q+b\right)$$

We differentiate the obtained Expression (19) on the right and left by Ω ν = (q − 1) times (here, the derivative of the integer is calculated, not a fractional order):

$${\int }_0^{\mathrm{\infty }}{e}^{-pt}·t^{b-1}·\left\{\frac{d^{\left(q-1\right)}}{d{\Omega }^{\left(q-1\right)}}\sum _{j=q-1}^{\mathrm{\infty }}\frac{(-\Omega {t^{\mu })}^j}{\mathrm{\Gamma }\left(\mu j+b\right)}\right\}dt={(-1)}^{q-1}\frac{\left(q-1\right)!·p^{\mu -b}}{{\left\{p^{\mu }+\Omega \right\}}^q}$$

Next,

$${\int }_0^{\mathrm{\infty }}{e}^{-pt}·t^{b-1}·\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+\nu \right)!(-{t^{\mu })}^{j+\nu }{\Omega }^j}{j!·\mathrm{\Gamma }\left(\mu j+\mu \nu +b\right)}dt={\int }_0^{\mathrm{\infty }}{e}^{-pt}·t^{b-1}·\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-{t^{\mu })}^{j+q-1}{\Omega }^j}{j!·\mathrm{\Gamma }\left\{\mu j+\mu \left(q-1\right)+b\right\}}dt={(-1)}^{q-1}\frac{\left(q-1\right)!·p^{\mu -b}}{{\left\{p^{\mu }+\Omega \right\}}^q}$$

Equating expressions $\frac{p^{\mu -b}}{{\left\{p^{\mu }+\Omega \right\}}^q}$ and $\frac{p^{-\beta q-1}}{{\left\{p^{\beta }+\Omega \right\}}^q}$, we find that at $\mu =\beta$ and $b=\beta \left(q+1\right)+1$, the following equality is performed:

$${\int }_0^{\mathrm{\infty }}{e}^{-pt}·t^{\beta \left(q+1\right)}·\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-{t^{\beta })}^{j+q-1}{\Omega }^j}{j!·\mathrm{\Gamma }\left\{\beta \left(2q+j\right)+1\right\}}dt={(-1)}^{q-1}\frac{\left(q-1\right)!·p^{-\beta q-1}}{{\left\{p^{\alpha }+\Omega \right\}}^q}$$

Next, ${(-1)}^{q-1}\frac{\left(q-1\right)!·p^{-\beta q-1}}{{\left\{p^{\alpha }+\Omega \right\}}^q}$ is equated to the expression

$$t^{\beta \left(q+1\right)}·\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-{t^{\beta })}^{j+q-1}{\Omega }^j}{j!·\mathrm{\Gamma }\left\{\beta \left(2q+j\right)+1\right\}}={\left(-1\right)}^{q-1}·t^{2\beta q}\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-\Omega {t^{\beta })}^j}{j!·\mathrm{\Gamma }\left\{\beta (2q+j)+1\right\}}$$

Next, we find

$${\left(-1\right)}^{q-1}·t^{2\beta q}\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-\Omega {t^{\beta })}^j}{j!·\mathrm{\Gamma }\left\{\beta (2q+j)+1\right\}}\risingdotseq {(-1)}^{q-1}\frac{\left(q-1\right)!·p^{-\beta q-1}}{{\left\{p^{\beta }+\Omega \right\}}^q}$$

Or, rewrite in the form

$$\frac{p^{-\beta q-1}}{{\left\{p^{\beta }+\Omega \right\}}^q}\risingdotseq \frac{t^{2\beta q}}{\left(q-1\right)!}\sum _{j=0}^{\mathrm{\infty }}\frac{\left(j+q-1\right)!(-\Omega {t^{\beta })}^j}{j!·\mathrm{\Gamma }\left\{\beta \left(2q+j\right)+1\right\}}=t^{2\beta q}\left\{\frac{1}{\mathrm{\Gamma }\left\{2\beta q+1\right\}}+\sum _{j=1}^{\mathrm{\infty }}\frac{(-\Omega {t^{\beta })}^j\prod _{i=1}^j(q+i-1)}{j!·\mathrm{\Gamma }\left\{\beta (2q+j)+1\right\}}\right\}$$

After all the necessary substitutions, we obtain

$$T_n\left(t\right)={\sum }_{q=0}^{\mathrm{\infty }}{(-1)}^q{\lambda }_n^q{\left\{\mathrm{\Gamma }\left(2\beta +1\right)\right\}}^q·{\left\{\frac{t}{m\tau }\right\}}^{2\beta q}\left\{\frac{1}{\mathrm{\Gamma }\left\{2\beta q+1\right\}}+{\sum }_{j=1}^{\mathrm{\infty }}\frac{(-\Omega {t^{\beta })}^j\prod _{i=1}^j(q+i-1)}{j!·\mathrm{\Gamma }\left\{\beta (2q+j)+1\right\}}\right\}$$

$${\lambda }_n=\pi n/L (n=\mathrm{0,1},\mathrm{2,3}), \Omega =\frac{\mathrm{\Gamma }\left(2\beta +1\right)}{{\left\{m\tau \right\}}^{\beta }\mathrm{\Gamma }\left(\beta +1\right)}$$

Accordingly, the general solution for $\rho \left(x,t\right)$ is

$$\mathrm{\rho }\left(x,t\right)=\sum _{n=0}^{\mathrm{\infty }}{C}_n{T}_n(t)·X_n(x)$$

Using the initial condition and orthogonality property of the function $X_n\left(x\right)=C_n\mathrm{s}\mathrm{i}\mathrm{n}({\lambda }_{n}x)$, we find that $C_n=2/L$.

5. Conclusions

The dynamics of the behavior of complex systems is very complex, including the presence of memory and the possibility of self-organization. It should be noted that different processes of nature can be observed at the same time, and the driving causes of which may be of a hidden nature. It should be noted that to date, the dynamics of the behavior of various physical systems have been the most well studied. In many works devoted, for example, to the study of physical kinetics, modeling is carried out without the use of fractional-differential equations. In this regard, we can mention [56], in which, a dynamic renormal group analysis of the Burgers equation was carried out, which is a nonlinear generalization of the diffusion equation taking into account the influence of random noise on the behavior of the processes studied. The obtained results were applied to the description of a growing interface with the media during polymerization and the appearance of transverse fluctuations in the directed growth of a polymer in a random medium. The occurrence of noise may not be related to the diffusion process, but may have correlations with it in space and/or time. Weak and strong noise, according to the results obtained by the authors of the article [56], lead to different scaling indicators and the appearance of correlations. For spatial correlations with sizes less than critical values, any amount of noise matters, resulting in a strong link. In the absence of temporal correlations, two modes can be observed, with sizes smaller than the critical size, either hydrodynamic behavior is determined by white noise and correlations are not important, or correlations dominate.

To some extent, the processes studied can be processes with memory and self-organization at certain scales in terms of coordinates and time.

The solution of the generalized diffusion equation containing Cardara–Parisi–Zhang nonlinearity taking into account the influence of spatially correlated noise in combination with the long-range nature of interactions was considered in works [57,58], where an approach based on renorms of groups and phase diagrams was also used.

It is possible that in terms of considering correlations in space and/or time, the use of the Burgers equation or other generalized equations of nonlinear diffusion may be one of the alternatives to fractional-differential equations.

However, models based on fractional-differential equations do not require the consideration of characteristic scales and the introduction of critical dimensions, which, for example, can be an advantage when considering processes in social and economic systems.

Studies of complex social processes, for example, electoral campaigns and user activity on social networks, show that the time series observed in practice have fractality, and the systems whose dynamics they describe have memory and show self-organization. For example, if we analyze the dynamics of changes in mathematical expectations and dispersions of amplitudes of time series levels depending on the time interval for calculating these amplitudes (using a “sliding window”), complex dependencies are observed.

For example, for the mathematical expectation, there is a root dependence of a fractional degree and for dispersion, there dependence on a power law for a fractional exponent greater than 1.5. Fractional time dependencies indicate the presence of non-locality for this variable.

The examination of the excess shows the presence of the so-called “heavy tails”, with its size significantly greater than that of the normal distribution. This excess behavior indicates the presence of non-locality in the states of the time series levels.

The obtained results indicate that the process under consideration has memory and the possibility of self-organization, and its time series have unsteady, as well as non-locality, both in time and state.

Non-local processes are characterized by the fact that the transition to a certain state of the system or process depends not only on the local characteristics of the process or the behavior of the system near the point in question at the current moment in time, but also on the values obtained throughout the studied interval at previous points in time, i.e., they are globally dependent on the distribution over all states and on the history of the process (memory). Non-locality over time affects the probability density at the initial point in time, which can lead to the phenomenon of self-organization, and non-locality in the state affects the asymptotic behavior of the probability density to detect a certain state x at time t at large intervals of time.

Various types of fractional partial derivative differential equations can be used to describe such processes. Currently, a fractional diffusion equation is used to build models of the dynamics of time series with non-locality over states and time: $\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}=D\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}$. Analytical or numerical process models based on this equation are obtained only for the case of $0<\beta \le 1$ and $1\le \alpha \le 2$. Solutions outside the limits of these values of α and β are not present in the literature.

Given that the processes observed in complex systems have features that vary over a wide range of values, it is impossible to only consider non-locality of the state of the system x (variable describing the level of the time series) or of time t for the construction of the model. In addition, it is necessary to consider a more general case with arbitrary fractional values of α and β, and not just the case of $0<\beta \le 1$ and $1\le \alpha \le 2$, or consider differential equations containing multiple fractional derivatives, both in time and in state.

Based on the fractional-diffusion equation, an edge problem with arbitrary values of derivatives was formulated and solved. An analytical solution was obtained that can be used in practice to analyze the dynamics of time series of processes observed in complex systems.

In addition, a fractional-differential equation of the telegraph type with multiples (in time: $\beta , 2\beta , 3\beta \dots$ and state: $\alpha , 2\alpha , 3\alpha \dots$) by orders of fractional derivatives and its analytical solution for one particular boundary problem was considered.

In solving edge problems, the Fourier method was used. This makes it possible to represent the solution in the form of nested time series (one in time t, the second in state x), each of which is a function of the Mittag-Leffler type. The eigenvalues of the Mittag-Leffler function for describing states can be found using boundary conditions and the Fourier coefficient based on the initial condition and orthogonality conditions of the eigenfunctions.

The resulting solutions can be tested based on the observed time series. It is necessary to calculate the parameters of the developed models ($\alpha , \beta , \epsilon , D$ values) and check the accuracy of the forecasts obtained.

The novelty of the models described in this article is that the differential stochastic equation obtained when considering the probability schemes of transitions between process states (time series levels) contains multiple orders of fractional derivatives, both in time (β, 2β) and state (α and 2α). Also, for the obtained equation, the solution of one of the edge problems is considered. This allows you to expand the class of fractional-differential equations and their use to describe the dynamics of complex systems. The Zener equation [25] contains only fractional time operators, and the fractional diffusion equations do not contain multiple derivatives.

In addition, it is new that the article presents a solution for one of the edge problems fractionally—a diffusion-type equation $\frac{{\partial }^{\beta }\rho \left(x,t\right)}{\partial t^{\beta }}=D\frac{{\partial }^{\alpha }\rho \left(x,t\right)}{\partial x^{\alpha }}$ for arbitrary α and β values—while existing solutions are limited to considering fractional derivative values of 0 < β ≤ 1 and 1 ≤ α ≤ 2.