## Abstract

Recently, methods have been developed to model low-dimensional chaotic systems in terms of stochastic differential equations. We tested such methods in an electronic circuit experiment. We aimed to obtain reliable drift and diffusion coefficients even without a pronounced time-scale separation of the chaotic dynamics. By comparing the analytical solutions of the corresponding Fokker–Planck equation with experimental data, we show here that crisis-induced intermittency can be described in terms of a stochastic model which is dominated by state-space-dependent diffusion. Further on, we demonstrate and discuss some limits of these modelling approaches using numerical simulations. This enables us to state a criterion that can be used to decide whether a stochastic model will capture the essential features of a given time series.

## 1. Introduction

Developing mathematical models and effective equations of motion for complex dynamical systems is one of the central issues in science and mathematics. The concept is fundamental to disciplines such as engineering, physics, chemistry and nowadays even biology and geology. While some communities highlighted such approaches about 50 years ago, subsumed under the heading adiabatic elimination, the issue has a much longer tradition, in particular when chaotic or random motion is at stake. To some extent, one can date the origin back to the foundation of non-equilibrium thermodynamics and statistical physics, i.e. to the pioneering works on Brownian motion by Einstein (1905) and Langevin (1908). Brownian motion is the classical textbook example to model a system with a distinct time-scale separation between a few slow and many fast degrees of freedom by a stochastic differential equation. The resulting effective equation of motion captures the evolution of the slow degrees of freedom, while the fast ones are replaced by a heat bath that leads to suitable stochastic forces.

One should keep in mind that a stochastic model is never able to reproduce a single trajectory of an underlying deterministic (or stochastic) system. As such, a comparison would depend on the particular realization of the noise process (e.g. Arnold (2003) on how to exploit such ideas in a purely mathematical framework). Of course, this does not imply by any means that stochastic models cannot give vital information about the underlying dynamics (cf. the seminal example of Brownian motion and the Wiener process as the appropriate stochastic description in, for example, Gardiner (2004)). The relevant quantities for estimating the quality of a stochastic model are thus given by the mean values and correlation functions. In our context, we will consider a model to be satisfactory if it is able to reproduce the main features of the stationary distribution and of some two-point temporal correlation function, i.e. it is able to model the main slow time-scale dynamics.

Nowadays, with growing computational power being available, stochastic modelling is widely used. For instance, in weather forecast and climate prediction, subgrid and fast processes are being replaced by stochastic processes and this may lead to an efficient probabilistic forecast (Palmer 2001; Bechtold *et al.* 2008). Replacing a high-dimensional deterministic system or parts of it by a low-dimensional stochastic model is tempting and sometimes the only way to make predictions on the future dynamics of a complex physical system. Naturally, the question arises about the limitations of these approaches: does a time-scale separation actually have to be several orders of magnitude and does one require a huge number of fast degrees of freedom to base the stochastic modelling on the variations of the central limit theorem?

We address this question by reporting on our investigation on the stochastic modelling of chaotic time-series data. Our investigation is motivated by the quite recent finding that stochastic modelling is applicable not only in systems with many fast degrees of freedom but also in low-dimensional chaotic systems (Beck 1995; Friedrich & Peinke 1997; Friedrich *et al.* 2000; Givon *et al.* 2004; Majda *et al.* 2006; Stemler *et al.* 2007). We test such methods using data from an experimental system that has only three degrees of freedom and shows crisis-induced intermittency.

Intermittency is a phenomenon which is normally associated with two time scales—here, a fast oscillation on two pre-critical attractors and a slow switching between both. This time-scale separation makes the phenomenon a good candidate for stochastic modelling. As in our experimental system the time-scale separation is not pronounced, it is far from obvious whether stochastic modelling works. Therefore, our investigation may be considered as a benchmark for the limits of stochastic modelling.

We report on the stochastic modelling of deterministic systems showing intermittency or a similar kind of switching dynamics. Time-series analysis was used to obtain the stochastic models’ parameters and we evaluated the models’ performance by comparison with the direct measurements of the systems. The overall aim behind these investigations was not to find a general stochastic model for intermittency but to demonstrate that time-series analysis enables us to build the stochastic models of intermittent systems that are able to capture certain aspects of the slow dynamics.

This paper is organized as follows: in the next section, we present details of the experimental system and the relevant range of control parameters where intermittency is observed. Section 3 describes the time-series analysis that allowed us to model the dynamics in terms of a Fokker–Planck equation. A detailed discussion of the Fokker–Planck equation and comparison between the theoretical and the experimental results can be found in §4. In §5 we discuss further applications of our method and provide some indicators that can be used to decide whether a stochastic model is able to describe the relevant dynamics of a deterministic system.

## 2. Experimental system

The time-series analysis is based on the data from an oscillator proposed by Shinriki *et al.* (1981). This circuit belongs to the Chua family of circuits (Wu 1987). A schematic of the circuit is given in figure 1*a*. The autonomous circuit consists of an RLC-element coupled through a nonlinearity with an additional RC element. The essential nonlinearity consists of a pair of Zener diodes *Z*_{1,2} (BZX85C3V3). The current–voltage characteristic (see inset in figure 1*a*) can be described by *I*_{D}(*V*_{1}−*V*_{2})=*I*_{D}(Δ*V*)=sgn(Δ*V*)*f*(|Δ*V* |−*V*_{Z}), where *V*_{Z}=(1.02±0.04)V and *f*(*x*)=(*Ax*^{2}+*Bx*^{3})Θ(*x*) denotes a third-order fit with parameter values *A*=(13.1±0.7) mA V^{−2} and *B*=(−1.59±0.15) mA V^{−3}. In parallel to the control parameter *R*, there is an operational amplifier circuit (AD711JN) which acts as a negative resistor with constant resistance, *R*_{N}, and provides a power supply for the circuit.

To measure the voltage *V*_{1}, we used a transient recorder card (Meilhaus ME2600). The digital output channels of this card were used for the online variation of the control parameter *R*, which was realized by several digital resistors (Xicor X9C102/4P). The fully computer-controlled experiment enabled us to measure long time series of the circuit for many different control parameter values.

The dynamics of the circuit is shown in the bifurcation diagram in figure 1*b*. The data are typical for a Chua-type circuit: for increasing *R*, a period doubling route to chaos is found. Several periodic windows occur for *R*≥53 kΩ before at *R*=*R*_{c}≈66 kΩ the mono-scroll attractor collides with its counterpart and a double-scroll attractor is observed. For *R*≥*R*_{c}, the circuit shows crisis-induced intermittency.

Our time-series analysis focuses on this intermittency regime where the dynamics shows two time scales: a fast oscillation on the chaotic saddles (the former stable mono-scroll attractors) and a slow jumping dynamics between them (see figure 2*a*). The aim of the stochastic model is to replace the fast oscillations by stochastic forces and to describe the slow jumping dynamics in terms of a one-dimensional stochastic differential equation. A detailed analysis on the crisis-induced intermittency can be found elsewhere (Werner *et al.* 2008). Here, we want to stress that (a) for a wide range of Δ*R*, the mean residence time on the saddles obeys a typical scaling law (Grebogi *et al.* 1983)*τ*∼Δ*R*^{−γ} with Δ*R*=*R*−*R*_{c} and *γ*≈0.7, and (b) for the parameter values of Δ*R* employed here, there are about 10–100 oscillations on a chaotic saddle before a jump occurs (see table 1). The famous Kramer’s problem, i.e. the motion of an over-damped particle in a double well potential driven by a Gaussian noise (Risken 1989), shows at a superficial level a qualitative similar dynamics. Furthermore, we find in intermittent systems similar residence time distribution as in Kramer’s model. Measured data from the circuit are shown in figure 2*b* for Δ*R*=1.1 kΩ. As in the stated bistable stochastic system, the residence time distribution decays exponentially. However, we will show that such an ad hoc stochastic model does not capture the essential features of crisis-induced intermittency.

## 3. Time-series analysis

Consider a nonlinear dynamical system giving rise to a time series via a scalar measurement *z*(*t*)=*g*[** x**(

*t*)]. The scalar time series might be, for example, a component of the state vector

**or some nonlinear function of the components. Our main goal is to explore whether such a time series can be reproduced by a Langevin equation 3.1where**

*x**h*(

*z*) is a deterministic force and

*k*(

*z*)

*ξ*(

*t*) is Gaussian white noise with state-dependent intensity

*k*

^{2}. Equivalently, the stochastic dynamics can be described in terms of the corresponding Fokker–Planck equation for the evolution of the probability density

*ρ*3.2As we base our analysis entirely on the Fokker–Planck equation, we do not have to specify in detail the stochastic integrals, although within our context a Stratonovich approach would be the natural choice. There exists a standard recipe known as the Kramers–Moyal expansion (Risken 1989) to obtain the drift

*D*

_{1}and the diffusion coefficient

*D*

_{2}from a stochastic time series by evaluating the first and second moments 3.3a and 3.3b The (conditional) averages are taken for an ensemble with

*z*(

*t*)=

*Z*. Within the context of a Langevin equation subjected to Gaussian white noise, such expressions yield in the limit the correct values for

*D*

_{1}and

*D*

_{2}.

Obviously, the limit is not applicable in a deterministic set-up because on this time scale the dynamics is strongly correlated. A suitable Δ*t* value has to be large compared with the correlation time of the fast chaotic motion but apart from this constraint equations (3.3) have to be evaluated in the asymptotic limit of small Δ*t*. One cannot take it for granted that a stochastic model such as equation (3.2) will capture the essential features of a deterministic time series *z*(*t*). Equation (3.2) and the time-series analysis based on equations (3.3) rely on the possibility of describing *z*(*t*) by a Markovian process. Although it is not clear whether our circuit data fulfil this requirement, we will nevertheless first model the data and test for such properties later.

Our experimental data are the voltage *z*(*t*)=*V*_{1}(*t*). Each of the analysed time series consisted of 6×10^{5} data points measured with a sampling time of 80 *μ*s. Evaluation of equation (3.3) typically results in a dependence of the moment on the time lag Δ*t* that can be seen in figure 3: for intermediate values Δ*t*∈[2 and 8 ms] a plateau is found. Thus, the values in this range can be used to estimate the coefficients of equation (3.2). To estimate spatially resolved *D*_{1} and *D*_{2}, we choose Δ*t*=4.5 ms. The conditional averages in equation (3.3) have been computed on a spatial grid with a step size of Δ*Z*=50 mV. This grid leads to an ensemble size of approximately 5×10^{3} for each datum point. We made sure that the final results are stable against these particular choices.

Figure 4 shows *D*_{1} and *D*_{2} evaluated with Δ*t*=4.5 ms at Δ*R*=1.1 Ω. *D*_{1} is dominated by a linear decay with a superimposed oscillatory structure. *D*_{2} turns out to be unimodal: large values around *Z*=0 V and small values at the boundaries. These large-scale features are robust against the reasonable changes of the time-series analysis parameters Δ*t* and Δ*Z*. The fine structure depends slightly on the details of the analysis, in particular on the ensemble size. Thus, we conclude that the fluctuations visible in figure 4 are to some extent a finite size effect, caused, for example, by the spiky time series. For instance, large fluctuations of the diffusion *D*_{2} are visible for *z*∈±[0.7 and 1.4 V] while no such features are observed for the stationary distribution (shown as an inset in figure 4). But the stable large-scale results achieved lead to our first conclusion that, despite the lack of time-scale separation, the time-series analysis leads to reproducible drift and diffusion coefficients.

## 4. Properties of the stochastic model

A very simple analytical approach which captures the essential features of the coefficients shown in figure 4 is to approximate the drift by a linear function and the diffusion by an inverted parabola 4.1a and 4.1bSuch a model seems to be quite crude as it does not take the oscillating part into account. However, it is not within the scope of our approach to model all the fine structure, which in fact could be meaningless given that we approximate fast deterministic chaotic dynamics by stochastic white noise. First of all, the simple approximation (4.1) has been tailored such that the corresponding stochastic model can be solved by analytical means and thus allows a detailed theoretical analysis. In addition, if quite crude approximations of the drift and the diffusion are able to reproduce the main features of the intermittent dynamics, then we have some good indications that the modelling is a robust approach. Of course, it would be tempting to check by numerical means whether the inclusion of the small-scale oscillatory parts yields important changes. But our results demonstrate that it is apparently sufficient to take only the stable large-scale results into account that are stated above.

Using the drift and diffusion coefficients (4.1), we can determine the stationary distribution of equation (3.2) as . For *α*<2*D*, the bimodal distribution reproduces the measured stationary distribution (see inset in figure 4), and we conclude that the model at least qualitatively captures the asymptotic spatial distribution of *z*.

The dynamical characteristics of the model are given by the eigenvalues of the Fokker–Planck operator governing the time evolution. Using a separation ansatz, they follow as
For *α*≪*D*, there appears a spectral gap between the two largest eigenvalues, Λ_{0} and Λ_{1}, and the rest of the spectrum. For a wide enough spectral gap, the mean residence time is determined by the leading non-trivial eigenvalue Λ_{1} and yields *τ*=2/*α*. Note that, in this approximation, *τ* is independent of the diffusion *D*.

We evaluated this analytical estimate by computing the slope of the linear trend *α* using a least-squares fit in the interval [−1.6V,1.6V]. This investigation is based on the time series from the circuit measured with a control parameter resolution of 100 Ω. A comparison between the theoretical estimate *τ*=2/*α* and the directly measured residence times is displayed in figure 5. Taking into account the simplicity of the model the agreement between the experimental and the theoretical results is quite accurate over a wide control parameter range.

Note that the experimental mean residence time has been evaluated from the statistics of the deterministic time series while the theoretical prediction is solely based on the drift and diffusion coefficients which themselves have been evaluated from experimental data. To some extent, such a comparison is a consistency check as both the experimental result and the theoretical prediction are based on the same experimental data. However, this consistency check is precisely what we are aiming for as the coincidence confirms the validity of the main modelling step, namely approximating the deterministic intermittent dynamics by an effective stochastic Markov model. The fact that the stochastic model reproduces as well qualitatively the correct stationary density would not suffice as a vindication because the mean residence time relates to higher eigenvalues of the corresponding Fokker–Planck operator while the stationary density just corresponds to the lowest eigenstate. In that respect, the stationary distribution is a spatial information and does not tell anything about the dynamics leading to it. From temporal eigenvalues, we gain insights of the dynamics but nothing about the spatial properties. Therefore, our results lead to the conclusion that the stochastic model captures the essential features of the intermittent dynamics on the slow time scale.

Although the accuracy of the mean residence time seems to confirm the validity of the Fokker–Planck model, we checked for higher order moments in the Kramers–Moyal expansion. At least the third- and fourth-order coefficients do not vanish. A strict test for Markovian property of the time series was not possible with the given ensemble size. But using conditional correlation functions, we were able to estimate that non-Markovian effects are visible at time scales of about 20 ms. Apparently, these shortcomings and the simplistic nature of the model that neglects the fine structure of drift and diffusion are of minor importance for the proper modelling of the measured intermittency, as demonstrated by the results of our analysis.

## 5. Estimating the validity of stochastic modelling

We are now going to explore some limits of the approach illustrated in the previous section. *A priori* there is of course no guarantee that the drift and diffusion coefficients computed from equation (3.3) yield a Fokker–Planck equation which models the underlying deterministic dynamics. In fact, it is trivial to construct counterexamples. For instance, consider a noise-driven particle that is moving in a double well potential in the *x*-direction while it is showing Brownian motion not subjected to any potential in the *y*-direction. Analysing *y*(*t*) will not lead to a valid description of the jump process. Thus, one issue to be addressed is the influence of the measured quantity on the modelling procedure. Furthermore, for deterministic dynamics the topology of the attractor can lead to a high level of convolution of the different states, when measurements are based on a single scalar quantity. This complexity can make even state identification in intermittent time series quite a challenging task (Jüngling *et al.* 2008). Altogether, such constraints can make stochastic modelling quite difficult.

We want to demonstrate such features by using the simulations of two different systems. The first system is the Duffing equation (Moon & Holmes 1979)
5.1with *δ*=1/4, *ω*=1 and *γ*=0.3. Similar to the Shinriki oscillator, the system undergoes for the increasing values of *γ* the classical route to chaos via period doubling bifurcations followed by a merging crisis. For *γ*=0.3, the system shows intermittent jumping dynamics between two symmetrical saddles (see the time-series sample *z*(*t*)=*u*(*t*) in figure 6*a*).

The second system is the famous Lorenz equations (Lorenz 1963)
5.2with *a*=10, *b*=8/3 and *r*=50. Although the Lorenz system does not show crisis-induced intermittency, the dynamics shows as well fast oscillations and slow jumping dynamics between two symmetrical saddles (see the time-series sample *z*(*t*)=*x*_{2}(*t*) in figure 6*b*).

We applied the time-series analysis method described in §3 to datasets *z*(*t*)=*u*(*t*) and *z*(*t*)=*x*_{2}(*t*) of appropriate length and time resolution. Using the adapted values of the spatial resolution Δ*Z*, the *Z* dependence of the conditional averages (3.3) can be analysed. For both datasets, a plateau similar to the one in figure 3 appears for the intermediate values of Δ*t* which are larger than the oscillation period but smaller than the mean residence time. Spatially resolved drift and diffusion coefficients were estimated using a constant Δ*t* value lying within the plateau.

For the Duffing system (5.1), the spatial dependence of *D*_{1}(*z*) and *D*_{2}(*z*) is similar to the result displayed in figure 4. The drift can be described by *D*_{1}(*z*)=−*αz*, the diffusion has a unimodal shape , where *α*<*D*. Consequently, the stochastic model reproduces at least qualitatively the measured invariant density shown in figure 7*a* and the overall findings are quite similar to the features reported for the Shinriki circuit.

For the Lorenz system (5.2), the findings are quite different. While the drift decays linearly, the diffusion coefficient is small for |*z*|≤10 and increases strongly for |*z*|≥10. Such findings are in fact in good agreement with the measured invariant density in figure 7*b*. For instance, the simple choice *D*_{1}(*z*)=−*αz* and in equation (3.2) results in a unimodal stationary distribution. It is tempting to analyse by analytical means the dynamical features of such kinds of stochastic models, i.e. the spectrum of the corresponding Fokker–Planck operator. At least it seems far from obvious that the estimates for the drift and diffusion coefficients can be used to describe the jumping dynamics shown in figure 6*b*.

The seeming shortcomings of describing the jumping dynamics of the Lorenz system in terms of a one-dimensional stochastic model are reflected by the invariant density of the data we wanted to model. A unimodal density can result even if the time-series data indicate bistable dynamics. In the Lorenz model, the invariant density reflects small values of around *x*_{2}=0 (Sparrow 1982). We want to stress that similar situations can occur in intermittent systems, e.g. the system described by Curry (1978). For the Lorenz system using a nonlinear measurement function, *z*(*t*)=*g*[**x**(*t*)]=*a*×sgn(*x*_{2})*x*_{3} results in a bimodal invariant density. Our analysis shows that for appropriate choices of *a* this dataset results in a cubic drift coefficient and a constant diffusion describing the jumping dynamics quite well.

## 6. Conclusion

Our investigation shows that, even without excessive time-scale separation, one may successfully develop a stochastic model for low-dimensional chaotic systems. We have demonstrated that reliable drift and diffusion coefficients can be extracted using the Kramers–Moyal expansion with a time lag that has been adapted to the deterministic nature of the data. The success of the approach can be obstructed by the nature of the deterministic dataset or by the features of the measured quantity. The empirical invariant density can be used to judge the quality of the stochastic modelling.

Although our simple stochastic model of the Shinriki oscillator approximates only the large-scale features of the drift and diffusion, the essential properties of the intermittency are captured quite well. The analytically determined mean residence times are in good agreement with the experimental data over a wide range of control parameters and so are the analytical and experimental stationary distribution. Both examples of crisis-induced intermittency, the circuit experiment and the Duffing system, yield a stochastic model where tunnelling is dominated by state-space-dependent diffusion, a mechanism which is in sharp contrast to the famous Kramer’s model.

The modelling of chaotic motion by suitable stochastic forces may have a wide range of applications even for low-dimensional systems with no pronounced time-scale separation. Stochastic models may reveal features about the underlying nonlinear dynamics. Thus, such an approach is by any means far more than just a practical tool to reduce a nonlinear system to a one-dimensional stochastic model.

## Acknowledgements

T.S. would like to thank Kevin Judd for the helpful discussions and acknowledges the support from ARC, grant no. DP 0662841. J.W. and H.B. acknowledge the support from DFG, grant no. BE 864/4-4.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society