Prediction of complex oscillations in the dynamics of coupled chaotic systems using transients

Journal Pre-proof Prediction of complex oscillations in the dynamics of coupled chaotic systems using transients O.N. Pavlova, A.N. Pavlov PII: DOI: Reference:

S0378-4371(19)32124-7 https://doi.org/10.1016/j.physa.2019.123818 PHYSA 123818

To appear in:

Physica A

Received date : 27 September 2019 Revised date : 3 November 2019 Please cite this article as: O.N. Pavlova and A.N. Pavlov, Prediction of complex oscillations in the dynamics of coupled chaotic systems using transients, Physica A (2019), doi: https://doi.org/10.1016/j.physa.2019.123818. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

*Highlights (for review)

Journal Pre-proof Highlights

Jo

urn

al

Pr e-

p ro

of

 The ability to recognize complex oscillations from transients is studied.  Prediction can be provided by analyzing numerical measures on the amount of samples.  Excluding the most nonstationary part of the data does not reduce the total number of samples used for forecasting.

*Manuscript Click here to view linked References

Journal Pre-proof

of

Prediction of complex oscillations in the dynamics of coupled chaotic systems using transients

b Yuri

p ro

O.N. Pavlovaa , A.N. Pavlova,b,∗

a Saratov State University, Astrakhanskaya Str. 83, Saratov 410012, Russia Gagarin State Technical University of Saratov, Politechnicheskaya Str. 77, Saratov 410054, Russia

Pr e-

Abstract

Forecasting complex oscillations based on transients is limited for stochastic and chaotic systems. Unlike restoring further samples, forecasting the emerging oscillatory mode is another important problem, especially near bifurcation points or in the case of multistability. Using the model of two diffusively coupled R¨ossler oscillators, we recognize various types of chaotic and hyperchaotic attractors, considering a relatively small amount of data, including transients. We show that different methods give similar estimates for a data set required, which is based on the behavior of several quantitative measures depending on the number of samples. The exclusion of the initial part of the data does not allow to reduce the total signal duration necessary for the forecast. Keywords: complex oscillations; transient process; detrended fluctuation analysis; return times; chaotic attractor PACS: 05.45.Tp

al

1. Introduction

Jo

urn

Predicting the behavior of a system based on available experimental data is a problem of major interest in various fields of human activity [1–10]. Many natural systems need forecasting methods to understand and possibly correct their further dynamics. Early prediction of pathological changes in the functioning of the body in medicine allows timely treatment that can save lives. In particular, predicting epileptic seizures remains an important issue and is still actively investigated [11, 12]. Prediction of failures in the operation of technical systems avoids accidents. Seismic waves, economic crises, or global temperature are just a few examples where forecasting is needed. Over the past decades, many approaches to forecasting have been proposed based, e.g., on machine learning or other techniques originating from computer science. Usually, methods for predicting time series include developing models that provide further values based on previous measurements. Linear models can be obtained using an autoregressive integrated moving average [13], which is a widely used tool. Nonlinear models are offered with generalized autoregressive conditional heteroskedasticity (GARCH) [14], artificial neural networks [15], etc. Despite the presence of many modeling methods, the advantages of nonlinear models are not always obvious [16], and the performance of linear models can be comparable. Nevertheless, recently developed approaches can outperform standard tools, and their application in solving many problems seems encouraging [17]. Time series forecasting is limited for stochastic and chaotic systems, and it becomes ineffective at time intervals exceeding the predictability horizon. Although the performance of different techniques varies, there is a limited duration of time for which a reasonable level of accuracy can be achieved. Taking into account the sensitive dependence ∗ Corresponding

author Email address: [email protected] (A.N. Pavlov)

1

Journal Pre-proof

2

Pr e-

p ro

of

of chaotic systems on the choice of initial conditions, models describing natural phenomena can provide only a shortterm prognosis, while forecast errors increase dramatically for long-term prediction. Other reasons restricting the ability to predict the behavior of chaotic systems, which are highly unstable, are variations of system parameters and random disturbances. Even when dealing with scalar time series produced by benchmark models of nonlinear dynamics, the creation of a mathematical description based on reconstruction is rarely provided with appropriate accuracy [18] and does not allow for long-term forecasts. In the case of natural systems, such as, e.g., physiological systems, short-term prediction errors can significantly reduce the value of the forecast. Regardless, the exact reproduction of the time series is not always required to understand the behavior of the system. If an external force or sudden changes in the parameters take the system out of the steady-state mode, then it performs a transition process to a new (or return to the previous) dynamic mode. In the case of multistability or near bifurcation points, several attractors exist, and the system can reach one of them after the completion of the transient process. The duration of this process varies depending on several circumstances, such as, e.g., how close the parameters are to the bifurcation points, how far the phase space variables are from the attracting limit set, etc. Nevertheless, the choice of a particular attractor can be predicted prior the end of the transient process. In this study, we address the problem of forecasting an attracting limit set using the initial part of the time dependence of the measured variable. We examine whether the prediction can be done based on relatively short data sets associated with transients. Taking into account the nonstationary dynamics of the system, we perform signal processing with the detrended fluctuation analysis [19, 20], which is a commonly used approach for characterizing the correlation properties of natural systems. To confirm the results and conclusions, we consider other fairly universal numerical methods, namely, the wavelet-based multifractal formalism [21, 22] and multiresolution analysis based on the discrete wavelet-transform [23, 24]. As a test model, we use a system of two coupled R¨ossler oscillators [25] that produces a variety of complex dynamic modes, including synchronous and asynchronous chaotic oscillations, hyperchaotic and quasiperiodic dynamics [26]. The additional complexity of the model under consideration is associated with the phenomenon of phase multistability of chaotic limit sets. The manuscript is organized as follows: in Sec. 2, we briefly consider the model of two coupled R¨ossler oscillators and three data processing methods that will be applied to obtain information used for forecasting the attracting limit set from the transient process. The main results describing the recognition between several types of attractors based on the transient dynamics are provided in Sec. 3. Section 4 summarizes the conclusions of this study.

2.1. Coupled R¨ossler systems

al

2. Model and methods

urn

The model of two diffusively coupled non-identical R¨ossler systems [25] is selected due to a rich family of complex oscillations arising from changes in control parameters and initial conditions for phase variables [26]. This model is described by the following set of ordinary differential equations dx1,2 dt dy1,2 dt dz1,2 dt

= −ω1,2 y1,2 − z1,2 + γ(x2,1 − x1,2 ),

= ω1,2 x1,2 + ay1,2 ,

(1)

= b + z1,2 (x1,2 − c), ω1,2 = ω0 ± ∆

Jo

where the parameter values are selected as ω0 =1, a=0.15, b=0.2 and γ=0.02. The two remaining parameters, c and ∆, vary to cover areas of different types of complex oscillations. The details of the transitions between various periodic, quasiperiodic, chaotic and hyper-chaotic attractors, as well as the features of chaotic synchronization in the model (1), are described in Ref. [26]. Here we do not discuss all modes and paths from periodic to chaotic oscillations. In this regard, Fig. 1 shows a simplified bifurcation diagram with several types of complex dynamics, namely: synchronous chaotic attractors (S C0 and S C1 ), asynchronous chaotic attractor (AC), and hyper-chaotic attractor (HA). Besides, periodic oscillations of period 4 (4C0 and 4C1 ) and quasiperiodic oscillations (QA) are presented. Periodic oscillations 0 1 of a higher period are not shown; they occur in regions below the bifurcation lines lcr and lcr , where cascades of the period-doubling lead to the appearance of chaotic attractors. The subscripts 0 and 1 in the designation of periodic 2

Journal Pre-proof

3

(4C0 , 4C1 ) and chaotic attractors (S C0 , S C1 ) denote the families of in-phase and out-of-phase oscillations related to the phenomenon of phase multistability. Such oscillatory modes are characterized by different phase shifts between the variables of each system (0 and π, respectively), however, these attractors are quantified by almost identical correlation features, complexity measures, and other properties. Due to this, we will omit the subscripts when considering transitions between attractors with significantly stronger structural changes. 000000000000000 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 7.2111111111111111 000000000000000 111111111111111 0000000000000 1111111111111 111111111111111 000000000000000 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 000000000000000 111111111111111 0000000000000 1111111111111 7.0111111111111111 000000000000000 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 m 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 1 0000000000000 1111111111111 6.8 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0 1 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 6.6 0000000000000 1111111111111 0000000000000 1111111111111 0 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 cr 0000000000000 1111111111111 0000000000000 1111111111111 1 0000000000000 1111111111111 0000000000000 1111111111111 6.4 0000000000000 1111111111111 cr 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0

of

HA

L

L

c

SC

l

AC

p ro

SC

l

l+1

6.0 0.009

4C 0

4C 1

1

l +1

QA

Pr e-

6.2

0.0095

∆

0.010

Figure 1. A bifurcation diagram that contains the main types of attractors considered in this paper. A thorough description of the dynamics of model (1) is given in Ref. [26].

urn

al

Further, we will mainly discuss the transitions between synchronous chaotic (S C) and asynchronous chaotic (AC) oscillations, as well as between hyperchaotic (HA) and asynchronous chaotic (AC) attractors. These modes are distinguished by a visual analysis of the projections of phase portraits (Fig. 2), however, the latter can be performed in stationary modes. In the case of transients, when the model (1) changes its behavior, such visual control becomes less suitable. Moreover, different types of oscillations are quite similar (e.g., S C and HA), and their visual recognition becomes doubtful if the phase trajectory includes a transient. Therefore, we need quantities that will allow us to separate the distinct modes. To perform recognition, we consider the sequences of return times in the Poincar´e section, which provide information on the multiscale structure of chaotic and hyperchaotic attractors. Let us note that these sequences are rather similar for several hundreds of return times (Fig. 2d). To discuss, how the choice of the secant plane can affect the results, we compare two Poincar´e sections: x1 = 0 and y1 + x2 = 0. 2.2. Detrended fluctuation analysis

Jo

The detrended fluctuation analysis (DFA) [19, 20] is chosen as the primary approach for signal processing in this study. It provides a characterization of power-law correlations in nonstationary experimental data and is easier to implement in comparison with the wavelet-transform modulus maxima (WTMM) method [21, 22], which performs a similar determination of the correlation properties. This approach was proposed for the analysis of nonstationary processes in biology and medicine. However, the advantages of DFA over standard correlation analysis offer a wide field of application [25–30]. In addition to the case of nonstationary dynamics, DFA provides a better characterization of long-range correlations, when the correlation function fluctuates near zero, and computational errors do not allow to obtain a reliable value of the scaling exponent. In DFA, the decreasing correlation function is replaced by a growing function computed as follows. The measured

3

Journal Pre-proof

4 12

12

SC

AC

4

4

y2

y2 −4

−4

(a)

−4

y1

4

12

−12 −12

(b)

of

−12 −12

−4

12

HA

12

4

p ro

SC

AC

x

y2 −4

−12 −12

−4

4

y1

12

0

(d)

200

400

HA

600

i

Pr e-

(c)

4

y1

Figure 2. Phase portraits of the attractors S C (a), AC (b) and HA (c) and fragments of the related sequences of return times (d).

time series x(i), i = 1, . . . , N is converted to a random walk y(k) =

k X i=1

[x(i) − hxi] , hxi =

N X

x(i).

(2)

i=1

al

The dependence y(k) is separated into non-overlapping fragments of length n, and the local trend yn (k) is estimated using the least-squares fitting. The simplest way is to use a piecewise linear function as yn (k). Further, the root-meansquare fluctuations around the function yn (k) are analyzed depending on n v u t N  1 X (3) F(n) = y(k) − yn (k) 2 ∼ nα . N k=1

urn

The function F(n) typically increases with n, showing the power-law dependence (3) for many natural processes, and the scaling exponent α is used as a measure of correlation properties. Earlier studies have established its relationship with an exponent describing the decay of the correlation function [19, 20]. Various applications of DFA are described in a number of publications, e.g., [31–35].

Jo

2.3. Wavelet-based multifractal analysis The WTMM [21, 22] approach is probably one of the most powerful tools for statistical analysis of nonstationary and inhomogeneous processes in various systems. It includes a continuous wavelet-transform of a function f (t) ! Z 1 ∞ t−b W(a, b) = f (t)ψ dt, (4) a −∞ a where a and b determine the scale and shift of the wavelet ψ. Usually, the distribution function is considered as f (t). Whereas numerical analysis is performed for a time series instead of a continuous function of time, the integral in Eq. (4) is replaced by the sum, and the estimate of the distribution function is similar to the construction of a random walk for processes with zero mean value. When computing the singularity spectrum, real-valued wavelets are applied. We used the MHAT wavelet being the second derivative of the Gaussian function. 4

Journal Pre-proof

5

A statistical analysis of the singularities of f (t) is performed based on the partition functions !q X ′ ′ Z(q, a) = sup | W(a , bl (a )) | ∼ aτ(q) , ′ l∈L(a) a ≤a

(5)

h(q) =

dτ(q) dq

p ro

and the singularity spectrum

of

where L is the full set of skeleton lines identified for each value of a, and bl denotes the position of the l-th line. The power-law behavior of the partition functions Z(q, a) enables computing the scaling exponents τ(q), which are used to estimate the local Hurst (or the H¨older) exponents

D(h) = qh − τ(q).

(6)

(7)

The value of H=h(0) is often close to the scaling exponent α of DFA [38]. Another important quantity is the width of the singularity spectrum, β = hmax − hmin , which can be treated as a measure of the inhomogeneity of the signal under study.

Pr e-

2.4. Multiresolution analysis The multiresolution analysis [23, 24] is another popular approach that allows studying signals and systems with time-varying properties. It provides the characterization of time series over a wide range of scales by decomposing the signal x(t) with low-pass ϕm,k and high-pass ψ j,k mirror filters X XX x(t) = sm,k ϕm,k (t) + d j,k ψ j,k (t), (8) k

ϕm,k

m/2

= 2

m

j ≥m k

ϕ(2 t − k), ψ j,k = 2 j/2 ψ(2 j t − k).

urn

al

In many applications, the Daubechies wavelets [24] are used for this decomposition, although the discrete wavelet transform offers different sets of filters that can illuminate various features of the signal x(t). Details of the decomposition and advantages of the wavelet bases applied in multiresolution analysis can be found in many related books and review papers. We used the D8 symmetric wavelet of the Daubechies family [24] although many other bases give similar results. It is usually preferable to achieve a compromise between regularity of function and support length. The D8 wavelet is an example of such a compromise. The coefficients d j,k are used to characterize the signal x(t) at several scales. A useful measure is the standard deviation of the decomposition coefficients v u t J−1 1X σ( j) = [d j,k − hd j,k i]2 (9) J k=0

Jo

which has been widely applied to diagnose the complex dynamics of natural systems [39]. Since it is evaluated as a function of the time scale, changes in σ( j) can be associated with mechanisms operating in distinct frequency domains. Here, we use the measure σ=σ(5) that was chosen after comparing different scales j and provided better separation between data sets under study. 3. Results and discussion

The scaling exponent of DFA clearly distinguishes between synchronous (S C) and asynchronous (AC) chaotic dynamics (Fig. 3a, c=6.8), as well as between hyperchaotic (HA) and chaotic (AC) modes (Fig. 3b, c=7.2). Recognition of synchronous chaos and hyperchaos is also performed using the appropriate range of lg n. In Figure 3, the ranges with the most pronounced distinctions are shown by dotted lines. These estimates are performed for relatively long 5

Journal Pre-proof

6

data sets containing 5,000 return times. Here, we consider not only stationary dynamics when the phase trajectory belongs to one of the attractors, but also a transient process from arbitrary initial conditions that are far enough apart from the attracting limit set. Usually, such amount of data is enough to complete the transient process and consider the features of steady-state dynamics. In the neighborhood of the bifurcation points, the duration of transient processes can essentially increase. 0.0

0.0

of

AC

−0.2

−0.3

AC

lg F

α=

0.

lg F

81

−0.4

−0.6

−1.0

1.0

1.5

SC

2.0

2.5

−1.2

3.0

lg n

(a)

HA

−0.9

α=0.01

−0.8

p ro

−0.6

(b)

1.0

1.5

0 0.8 α=

α=0.22 2.0

2.5

3.0

lg n

5000

5000

265

4000

τ

Pr e-

Figure 3. The dependencies F(n) in the log-log plot showing distinctions between the complex oscillations associated with the attractors S C and AC for c=6.8 (a), HA and AC for c=7.2 (b).

4000

τ

3000 205 0.00975

2000

∆

τ

0.00987

3000 2000

1

1

1000

1000

2

2 0.00979

(a)

∆

0.00983

0 0.00949

0.00987

al

0 0.00975

(b)

0.00954

∆

0.00959

0.00964

urn

Figure 4. A reducing of the time duration τ required to obtain the scaling exponent related to the asynchronous chaotic attractor AC, with the parameter ∆ changing from the bifurcation value. The cases of c=6.8 (a), and c=7.2 (b) are shown. The numbers 1 and 2 denote the results for the secant plane x1 =0 (1) and y1 + x2 =0 (2). The value of τ is measured by the number of return times.

Jo

Figure 4 illustrates the dependence of the time interval τ, which is necessary to achieve a state with the value of α associated with the asynchronous chaotic attractor (AC) for c=6.8 and c=7.2, respectively, if the phase trajectory starts from the initial conditions chosen rather far from AC (x j =y j =z j =20, j=1,2). Because the transient process and the limited amount of data affect the estimates, τ was computed as the time when the value of α falls into the range α∗ ±20%, where α∗ is the expected scaling exponent, computed for steady-state dynamics without restrictions of the length of the data set, and α does not leave such a range with an increase in the number of samples. According to Fig. 4, τ depends on the selection of the secant plane. Thus, it takes significantly larger values for the Poincar´e section x1 =0 (dependence 1) compared with y1 + x2 = 0 (dependence 2). This phenomenon was discussed in [40], where we showed that numerical measures of hyperchaotic dynamics, such as the second Lyapunov exponent, can be underestimated and, therefore, the oscillatory mode is wrongly identified if the secant plane does not take into account the dynamics of each of the interacting oscillators. The Poincar´e section x1 =0 means that the dynamics of the first R¨ossler system is mainly considered. The second system affects the behavior of the entire model (1), however, the sequence of return times may contain little information about its features in the case of weak coupling, and larger data sets are required to provide the appropriate characterization of the dynamics under study. When considering the 6

Journal Pre-proof

7

secant plane y1 + x2 = 0, a comparable amount of information about the complex dynamics of individual oscillators is contained in a sequence of return times, which reduces the required number of samples for a reliable characterization of oscillatory modes. Regardless of the absolute value of τ, the duration of the transient process usually decreases with the distance from the bifurcation point. In Fig. 4, the value of ∆ begins with the bifurcation point separating the regions associated with the attractors S C and AC (Fig. 4a), and the attractors HA and AC. For other secant planes, the results are between the dependencies 1 and 2. 1.0

1.0

0.8

0.8

0.6

α

AC (1) 0.4

0.6

p ro

α

of

AC (2)

AC (2)

AC (1)

0.4

0.2

SC (2) 0.0 100

800

1500

2200

2900

3600

4300

0.0 100

5000

N

(a)

HA (1)

0.2

SC (1)

(b)

800

1500

HA (2)

2200

2900

3600

4300

5000

N

Pr e-

Figure 5. Typical dependencies of the scaling exponent on the number of samples for c=6.8 (a) and c=7.2 (b). The results are shown for the secant plane x1 =0 (marked with 1) and y1 + x2 = 0 (marked with 2). For short data sets, the maximum value of lg n used to estimate α and the slope of lg F(lg n) decrease, and the scaling exponent takes values approaching zero. The choice of initial conditions affects these dependencies, but the general features remain unchanged.

Jo

urn

al

Figure 5 shows typical dependencies of α on the number of samples of the data set under study (N), which includes a transient process. According to this Figure, the estimates of the scaling exponent tend to achieve the expected values with increasing amount of data for both secant planes x1 =0 (dependencies are marked with number 1) and y1 +x2 = 0 (marked with number 2), although the duration of the data sets required to obtain the values of α with the appropriate accuracy varies significantly between these two Poincar´e sections. For the second plane, a reduced number of samples can be used to get the scaling exponent. In the case of the AC-attractor, there is a “jump” in α(N), and the scaling exponent takes values near the expected exponent (α∗ ≃0.8). For the Poincar´e section x1 =0, we can predict the attracting limit set even if α is still far from the expected value α∗ . Thus, an almost monotonic increase in α for the attractor AC, and a decrease in α for the attractor S C are observed (see Fig. 5a). Similar distinctions in the behavior of α(N) occur for the attractors AC and HA in Fig. 5b. Such behavior means that an attracting limit set can be forecasted taking into account the local slopes r of the dependence α(N) estimated within some floating window. Figure 6 illustrates such local slopes for the dependencies shown in Fig. 5 for the cases c=6.8 and c=7.2, and the secant plane y1 + x2 = 0. They demonstrate that there is a number of samples N when r reaches maximums, that is, the scaling exponent shows the highest growth rate, approaching the expected value α∗ . A related amount of data is sufficient to recognize that the model (1) will generate asynchronous chaotic oscillations associated with the AC attractor, but not synchronous chaotic oscillations (the S C attractor) or hyperchaos oscillations (the HA attractor). The number of samples required takes about 230 to distinguish between the attractors AC and S C, and 380 to recognize the attractors AC and HA. The consideration of the dependencies α(N) for S C and HA in Fig. 5 also makes it possible to distinguish between these types of complex motions by comparing local slopes. This recognition requires a larger amount of the data set (about 700–800 samples), however, this amount is almost 3-fold less than the data set which allows us to clearly estimate the scaling exponent of the HA-attractor using the sequences of return times to the Poincar´e section y1 + x2 = 0 (Fig. 5b, the dependence HA(2)). Then we compared the capabilities of other approaches in recognizing the type of complex motions from a transient data set. Figure 7 illustrates the results of the wavelet-based multifractal formalism. In contrast to Fig. 5, we give here estimates only for the Poincar´e section y1 + x2 = 0, which provided stronger changes in the scaling exponent α with the number of samples. According to Fig. 7, estimates of the average H¨older exponent show a correlation with α, providing strong changes in H for the same amount of data sets: about 230 samples for c=6.8, and 380 samples for c=7.2. Thus, the possibilities of these two methods in predicting the stationary dynamics from transient processes 7

Journal Pre-proof

8 0.010

0.006

r

−0.002 150

250

350

450

550

650

p ro

N

of

0.002

Figure 6. The local slopes r of α(N) for the attractor AC and the secant plane y1 + x2 = 0 estimated for c=6.8 (solid line) and c=7.2 (dashed line). The values of r are obtained by averaging within the window of 100 samples.

1.5

1.5

H

1.0

H

H, β

β 0.5

0.0 150

Pr e-

H, β

1.0

β

0.5

300

450

600

750

N

900

1050

1200

0.0 150

300

450

600

750

900

1050

1200

N

al

Figure 7. The estimates of the mean H¨older exponent (H) and the width of the singularity spectrum (β) for the attractor AC and the values c=6.8 (a) and c=7.2 (b).

Jo

urn

are comparable, however, it should be noted that the choice of measure is an important problem. According to [38], the transition between synchronous and asynchronous chaotic modes in the dynamics of two coupled R¨ossler oscillators is accompanied by two important changes: a shift in the H¨older exponents, which reflect distinct correlation properties, and a change in the degree of multiscality of the sequence of return times quantified by the width β of the singularity spectrum. Figure 7a shows that the latter measure can exhibit a distinct behavior depending on the number of samples N, and the “jump” in β(N) occurs with a delay as opposed to H(N). Thus, WTMM does not improve the ability to predict the occurring oscillatory mode, but, if the numerical measure is selected improperly, more data may be required for reliable prediction of the attractor. The delay in the behavior of β as opposed to the Hurst exponent H is explained by a small amount of data. The β estimates require analysis of a wide range of scales that is difficult to perform for short time series. In general, restrictions on the number of samples are responsible for the behavior of H(N) and β(N). In Fig. 7b, both measures vary similarly, although the changes in the H¨older exponents are stronger, which is important for the forecast. Note that more complex dynamics require more data for prediction based on H(N). The application of multiresolution analysis also allows establishing the amount of data needed to predict the type of attractor from data sets, including the transient process. The difference in the application of this method is the use of a discrete wavelet-transform, where data sets with 2 j samples are considered. Figure 8 shows that the same types of attractors are well separated using N=256 samples for c=6.8, and N=512 samples for c=7.2, which corresponds to estimates with DFA and WTMM taking into account the above features of the multiresolution analysis. Thus, all 8

Journal Pre-proof

9 0.12

0.12

0.08

0.08

AC

σ

AC

σ 0.04

0.04

0.00

0

1024

2048

3072

4096

0.00

0

of

HA

SC

1024

N

2048

3072

4096

N

p ro

Figure 8. Estimates of σ(N) for the sequences of return times in the Poincar´e section y1 + x2 = 0 for the attractor AC and two values of the control parameter c: 6.8 (a) and 7.2 (b).

4. Conclusion

urn

al

Pr e-

the approaches considered confirm the ability to predict various types of complex motions in the model (1) using similar amounts of data. None of these approaches offers the advantage in reducing the available information for earlier recognition of the type of attractor, however, the data sets required increases if an inappropriate measure and Poincar´e section are chosen. In the case of WTMM analysis, the choice of the measure β does not mean that it cannot distinguish between different types of attractors, but this requires larger data sets. Another important issue is related to the choice of initial conditions. In this study, we used initial conditions that were far from the attracting limit set. The latter means the appearance of rather long transient processes and strongly nonstationary phase trajectories, which demonstrate oscillations with a clearly different magnitude in their initial parts. These parts provide higher variation of measures and, therefore, the question arises whether is is possible to improve the results if such parts are removed. We performed recognition procedures by exclusion of different number M of return times at the beginning of the data sets and processing next K return times. The latter allows us to use data sets of about K=90–100 samples for recognizing the type of attractor, but the total amount of the data set (M+K) does not reduce compared to the case when the entire sequence of return times is considered. These estimates are performed for different initial conditions and the duration of transients. Therefore, the exclusion of these parts is not necessary if the approaches used are able to ignore non-stationarity due to detrending (DFA) or the application of wavelet functions with several vanishing moments which ignore polynomial trends in the data sets. Further, we plan to study more complex forecasting problems, e.g. to consider partial synchronization patterns called chimera states which occur in networks of coupled oscillators. It has been shown that these states are chaotic transients [41] (transients towards synchronized regime) and they have been found for chaotic bistable systems [42]. Very often they coexist with the completely synchronized regime (i.e., there is multistability in the network).

Jo

Prediction of emerging oscillatory dynamics in the case of multistability or near bifurcation points can be complicated by rather long transient processes, when the system leaves the steady-state mode due to external forces or changing external conditions, and the phase trajectory tends to attract one of the existing stable oscillatory modes. If some of these modes are undesirable, control methods are used to avoid such types of motions. However, it is necessary to predict the further dynamics of the system in order to determine which mode will be established. In this paper we considered the possibility of forecasting using the model of two diffusively coupled R¨ossler oscillators and three approaches to processing complex signals, including nonstationary data. Main results of this study consist in the following: 1) We showed how early prediction of the emerging dynamical mode can be carried out by considering changes in numerical measures with an increase in the number of samples. If these measures are selected appropriately, all approaches require almost the same amount of data to separate between different types of attractors. Thus, there is no 9

Journal Pre-proof

10

p ro

of

preferred method, but there are characteristics, the assessment of which requires an increased amount of data, and it is better to compare several numerical measures to conclude about the predicted dynamics. 2) We compared two secant planes to demonstrate that the number of samples needed for the prediction depends on the Poincar´e section. If the sequence of return times contains enough information about the dynamics of interacting subsystems, a reduced data set is enough. Otherwise, the duration of the data set should be increased. 3) The number of samples needed to predict complex behavior can be reduced by eliminating the most nonstationary (initial) part of the data that appears when the initial conditions or parameter values change. However, the total amount of data sets used for forecasting will not be reduced using such procedure. We hope that the approach under consideration can be applied to the analysis of real data, e.g., for early prediction of pathological changes in physiology and medicine. Epileptic seizures and intracranial hemorrhages are just a few examples of potential applications. In fact, this approach has limitations. In particular, recognized types of complex dynamics must have substantially different properties. The latter limits the ability of forecasting dealing with the hidden stages of pathologies. Nevertheless, the adaptation of forecasting methods for such cases is an important area of further research. Acknowledgments

This work was supported by the Russian Science Foundation (Agreement 19-12-00037).

al

urn

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

K. Lehnertz, C.E. Elger, Phys. Rev. Lett. 80 (1998) 5019. P.A. Stott, J.A. Kettleborough, Nature 416 (2002) 723. M.P. Clements, P.H. Franses, N.R. Swanson, Int. J. Forecast. 20 (2004) 169. G.P. Zhang, D. Kline, IEEE. Trans. Neural. Netw. 18 (2007) 1800. K.J. Arrow, R. Forsythe, M. Gorham, R. Hahn, R. Hanson, J.O. Ledyard, S. Levmore, R. Litan, P. Milgrom, F.D. Nelson, G.R. Neumann, M. Ottaviani, T.C. Schelling, R.J. Shiller, V.L. Smith, E. Snowberg, C.R. Sunstein, P.C. Tetlock, P.E. Tetlock, H.R. Varian, J. Wolfers, E. Zitzewitz, Science 320 (2008) 877. M. Scheffer, Nature 467 (2010) 411. J.M. Drake, B.D. Griffen, Nature 467 (2010) 456. G. Sugihara, R. May, H. Ye, C.H. Hsieh, E. Deyle, M. Fogarty, S. Munch, Science 338 (2012) 496. C. Boettiger, A. Hastings, Nature 493 (2013) 157. B. Podobnik, A. Majdandzic, C. Curme, Z. Qiao, W.-X. Zhou, H.E. Stanley, B. Li, Phys. Rev. E 89 (2014) 042807. R.G. Andrzejak, C. Rummel, F. Mormann, K. Schindler, Scientific Reports 6 (2016) 23000. T. Chouzouris, I. Omelchenko, A. Zakharova, J. Hlinka, P. Jiruska, E. Sch¨oll, Chaos 28 (2018) 045112. G.E.P. Box, G.M. Jenkins, G.C. Reinsel, Time Series Analysis: Forecasting and Control, Prentice-Hall, Englewood Cliffs, 1994. T. Bollerslev, J. Econometrics 31 (1986) 307. J.S. Zirilli, Financial prediction using Neural Networks, International Thompson Computer Press, London, 1997. M.P. Clements, P.H. Franses, N.R. Swanson, Int. J. Forecast. 20 (2004) 169. A. Golestani, R. Gras, Scientific Reports 4 (2014) 6834. A.N. Pavlov, N.B. Yanson, V.S. Anishchenko, Journal of Communications Technology and Electronics 44(9) (1999) 999. C.-K. Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley, A.L. Goldberger, Phys. Rev. E 49 (1994) 1685. C.-K. Peng, S. Havlin, H.E. Stanley, A.L. Goldberger, Chaos 5 (1995) 82. J.-F. Muzy, E. Bacry, A. Arneodo, Phys. Rev. Lett. 67 (1991) 3515. J.-F. Muzy, E. Bacry, A. Arneodo, Int. J. Bifurcation Chaos 4 (1994) 245. S.G. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 1998. I. Daubechies, Ten Lectures on Wavelets, Philadelphia, S.I.A.M., 1992. O.E. R¨ossler, Phys. Lett. A 57 (1976) 397. D.E. Postnov, T.E. Vadivasova, O.V. Sosnovtseva, A.G. Balanov, V.S. Anishchenko, E. Mosekilde, Chaos 9 (1999) 227. H.E. Stanley, L.A.N. Amaral, A.L. Goldberger, S. Havlin, P.Ch. Ivanov, C.-K. Peng, Physica A 270 (1999) 309. P. Talkner, R.O. Weber, Phys. Rev. E 62 (2000) 150. P.Ch. Ivanov, L.A.N. Amaral, A.L. Goldberger, S. Havlin, M.G. Rosenblum, H.E. Stanley, Z.R. Struzik, Chaos 11 (2001) 641. K. Hu, P.Ch. Ivanov, Z. Chen, P. Carpena, H.E. Stanley, Phys. Rev. E 64 (2001) 011114. Z. Chen, P.Ch. Ivanov, K. Hu, H.E. Stanley, Phys. Rev. E 65 (2002) 041107. R.M. Bryce, K.B. Sprague, Sci. Rep. 2 (2012) 315. A.L. Goldberger, L.A.N. Amaral, J.M. Hausdorff, P.Ch. Ivanov, C.-K. Peng, H.E. Stanley, PNAS 99 (2002) 2466. S.V. Buldyrev, A.L. Goldberger, S. Havlin, R.N. Mantegna, M.E. Matsa, C.-K. Peng, M. Simons, H.E. Stanley, Phys. Rev. E 51 (1995) 5084. H.E. Stanley, S.V. Buldyrev, A.L. Goldberger, S. Havlin, C.-K. Peng, M. Simons, Physica A 273 (1999) 1. K. Ivanova, M. Ausloos, Physica A 274 (1999) 349.

Jo

[1] [2] [3] [4] [5]

Pr e-

References

10

Journal Pre-proof

11

Jo

urn

al

Pr e-

p ro

of

[37] N.S. Frolov, V.V. Grubov, V.A. Maksimenko, A. L¨uttjohann, V.V. Makarov, E. Sitnikova, A.N. Pisarchik, J. Kurths, A.E. Hramov, Scientific Reports 9 (2019) 7243. [38] A.N. Pavlov, O.V. Sosnovtseva, A.R. Ziganshin, N.-H. Holstein-Rathlou, E. Mosekilde, Physica A 316 (2002) 233. [39] S. Thurner, M.C. Feurstein, M.C. Teich, Phys. Rev. Lett. 80 (1998) 1544. [40] A.N. Pavlov, O.N. Pavlova, Y.K. Mohammad, J. Kurths, Phys. Rev. E 91 (2015) 022921. [41] M. Wolfrum, O.E. Omel´chenko, Phys. Rev. E 84 (2011) 015201(R). [42] I.A. Shepelev, A.V. Bukh, T.E. Vadivasova, V.S. Anishchenko, A. Zakharova, Commun. Nonlin. Sci. Numer. Simulat. 54 (2018) 50.

11

*Declaration of Interest Statement

Journal Pre-proof

Conflicts of Interest Statement Prediction of complex oscillations in the dynamics of coupled chaotic systems using transients

of

Manuscript title:

p ro

The authors whose names are listed immediately below certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Author names:

al

Pr e-

Pavlova O.N. Pavlov A.N.

Jo

Author names:

urn

The authors whose names are listed immediately below report the following details of affiliation or involvement in an organization or entity with a financial or non-financial interest in the subject matter or materials discussed in this manuscript. Please specify the nature of the conflict on a separate sheet of paper if the space below is inadequate.

Journal Pre-proof

This statement is signed by all the authors to indicate agreement that the above information is true and correct (a photocopy of this form may be used if there are more than 10 authors):

Author's name (typed)

Author's signature

Date

11/03/2019

p ro

of

Pavlova O.N.

Jo

urn

al

Pr e-

Pavlov A.N.

11/03/2019