Dynamics of simple earthquake model with time delay and variation of friction strength

Abstract. We examine the dynamical behaviour of the phenomenological Burridge–Knopoff-like model with one and two blocks, where the friction term is supplemented by the time delay τ and the variable friction strength c. Time delay is assumed to reflect the initial quiescent period of the fault healing, considered to be a function of history of sliding. Friction strength parameter is proposed to mimic the impact of fault gouge thickness on the rock friction. For the single-block model, interplay of the introduced parameters c and τ is found to give rise to oscillation death, which corresponds to aseismic creeping along the fault. In the case of two blocks, the action of c1, c2, τ1 and τ1 may result in several effects. If both blocks exhibit oscillatory motion without the included time delay and frictional strength parameter, the model undergoes transition to quasiperiodic motion if only c1 and c2 are introduced. The same type of behaviour is observed when τ1 and τ2 are varied under the condition c1 = c2. However, if c1, and τ1 are fixed such that the given block would lie at the equilibrium while c2 and τ2 are varied, the (c2, τ2) domains supporting quasiperiodic motion are interspersed with multiple domains admitting the stationary solution. On the other hand, if (c1, τ1) warrant oscillatory behaviour of one block, under variation of c2 and τ2 the system's dynamics is predominantly quasiperiodic, with only small pockets of (c2, τ2) parameter space admitting the periodic motion or equilibrium state. For this setup, one may also find a transient chaos-like behaviour, a point corroborated by the positive value of the maximal Lyapunov exponent for the corresponding time series.


Introduction
Earthquakes represent complex deformation feature of the Earth's brittle crust.Their complexity reveals itself primarily in the power-law scaling describing the distribution of magnitudes (Turcotte et al., 2000), fractal spatial distribution of epicentres and fractal-like structure of faults (Okubo and Aki, 1987;Marsan et al., 2000).Nonetheless, earthquakes also involve complex patterns of temporal behaviour, which may be manifested by the chaotic dynamics in the recorded seismic time series (Beltrami and Mareshal, 1993).A plausible explanation for this complexity lies in that it is generated by nonlinear dynamics on a smooth fault, a point first raised by Horowitz and Ruina (1989).This line of research was followed by Carlson and Langer (1989) in the work on dynamics of Burridge-Knopoff arrays of spring-connected blocks with classical velocity-weakening friction laws.Also, Bak and Tang (1989) investigated a simple cellular-automata model of the failure at a critical stress.However, though all these approaches consider earthquakes from the phenomenological standpoint, there is still no general consensus on what constitutes the right model for describing the motions of earthquake faults.Therefore, it is useful to study a variety of models and, in doing so, learn how their various ingredients determine the behaviour that they exhibit (Langer et al., 1996).
Many different phenomenological models of fault dynamics are currently being investigated.The most prominent classes include the spring-block models (like the one in the present paper) and cellular automata.While the automata models facilitate rapid numerical computation, spring-block models have the advantage of allowing one to establish more direct associations between the parameters involved and the observations of the real seismic phenomena.The first continuum model of the latter type was proposed by Burridge and Knopoff (1967).This spring-block model, coupled with the appropriate friction law, is capable of generating powerlaw distribution of energy release, described by Gutenberg-Richter and Omori-Utsu laws (Burridge and Knopoff, 1967).The model consists of several blocks interconnected with springs of some elastic constant k, having the blocks also attached by a leaf spring to a driving plate, that causes the whole system to move along the rough surface in a stick-slip fashion (Fig. 1).
Apart from accounting for the statistical power-law distribution of released energy (Burridge and Knopoff, 1967), this model has been shown to exhibit rich dynamical behaviour (Carlson and Langer, 1989;De Sousa Vieira, 1995, 1999;Erickson et al., 2008).Carlson and Langer (1989) have demonstrated that chaotic dynamics may arise as a direct consequence of the friction law, which effectively causes small irregularities in the system to be amplified during the slipping motions.Their model produces three qualitatively distinct kinds of slipping events: microscopic events, large but localised events and delocalised great events.On the other hand, De Sousa Vieira (1999) analysed rich dynamical behaviour in a two block system, having identified the perioddoubling route to chaos.Erickson et al. (2008) reported that in a spring-block model with only one block coupled with the Dieterich-Ruina friction law chaotic behaviour also emerges via the Feigenbaum scenario.
All the research mentioned so far relies on the notion that the behaviour of the considered model is essentially controlled by the chosen friction law.Initially, Burridge and Knopoff (1967) in their original model used a simple friction law, where the friction force is dependent on the velocity of the block relative to the frictional surface.The main feature of this law is that the friction has to decrease more or less suddenly as the velocity deviates from zero, whereas for very small velocities, friction may attain any value less than the limiting friction.Afterwards, slip rate and state variable constitutive laws for rock friction were developed and introduced by Dieterich (1979), Ruina (1983) and Rice (Rice, 1983;Rice and Ruina, 1983).These early studies were designed to study frictional instability as a possible mechanism for the repetitive stick-slip failure and the seismic cycle.In this paper, we assume that the friction force at the contact of a block and a rough surface depends only on the velocity of the block.However, recalling that many experimental observations indicate that the rock friction also depends on the state of the contact surface (Marone, 1998), we introduce this memory effect by including the time delay parameter τ in the friction term.The role of this parameter is twofold.Firstly, it replaces the state variable θ incorporated in the Dieterich-Ruina friction law, which is usually interpreted as a function of history of sliding (Pomeau and Le Berre, 2011).Secondly, this time delay effect is directly observed both in laboratory experiments and along the real fault.According to the results of laboratory tests, upon the cessation of motion, static friction shows an initial period of retarded healing for a few hundred days, after which an increase in healing is observed (Marone, 1998).Moreover, it is determined that the length of this initial period of delayed healing varies with stiffness, which justifies our variation of the introduced time delay parameter τ .These laboratory results are further corroborated by seismic data, which indicate that the healing rate is reduced during the period immediately following earthquakes of similar size (less than 10-100 days after the last earthquake), with the small variations in the stress drop.
One should point out that this retarded initial period of fault healing is reminiscent of the refractory stage of the relaxation oscillators.Modelling the behaviour of spatially extended media comprised of relaxation oscillators often involves the delay-differential equations.Such an approach is widely accepted, particularly in the fields of mathematical biology (Golpasamy and Leung, 1996;Burić and Todorović, 2002) and life sciences (Smith, 2011).
The second analysed feature is the variation of the friction strength, which we model by introducing the parameter c.Under the real conditions in the Earth's crust, this strength variation of friction along the fault is usually induced by the formation and properties of the fault gouge material.The latter represents an unconsolidated and cohesionless type of fault rock, consisting almost entirely of the finely crushed material, originated as a final product of grinding when two sides of the fault zones move along each other (Sibson, 1977).There are some evidence that this material has significant effect on the stability of the natural faults (Scholz et al., 1969;Marone and Scholz, 1988) and rock friction (Byerlee, 1967;Engelder et al., 1975;Byerlee and Summers, 1976;Das et al., 1986), causing the change of the frictional resistance along the fault.This is significantly influenced by the clay minerals in a saturated state (Morrow et al., 2000;Behnsen and Faulkner, 2012), thickness of fault gouge (Marone et al., 1990) and fault depth (Mizoguchi et al., 2007).Having included the variability of friction strength in our analysis, we consider, at least qualitatively, the impact of thickness of fault gouge on the dynamical stability of motion along the fault.Also, by assuming different values of parameter c for the two blocks, the model acquires additional heterogeneity feature, which is certainly a common occurrence under the real conditions along the fault.
The paper is organised as follows.The second part provides a detailed description of the original model, together with the governing equations for the systems with one and two blocks.In Sect.3, we present the modified model.The dynamical behaviour of a single-block model is examined in Sect. 4 by considering its dependence on the variation of time delay τ and friction strength parameter c.In Sect.5, we analyse the model with two blocks, contingent on the values of four control parameters: c 1 , c 2 , τ 1 and τ 2 .For each setup, the different forms of dynamical behaviour are characterised by calculating the maximal Lyapunov exponent and the Fourier power spectrum.The final section contains the concluding remarks and some suggestions for future research.

Background of the original model
Our numerical simulations of the spring-block model are based on the system of equations introduced by Carlson and Langer (1989): where dots indicate derivatives with respect to time t, j stands for the block index, X j are the displacements of blocks of mass m measured from their initial equilibrium positions, and υ represents the velocity of the upper plate.Parameter k c represents the strength of harmonic spring connecting the blocks, and k p is the strength of the leaf springs connecting the block and the upper driving plate (Fig. 1).Friction force F is given in the following form: where φ vanishes at large values of its argument and is normalised so that φ(0) = 1, while υ c represents the speed that characterises the velocity dependence of the friction force F (Carlson and Langer, 1989).For convenience, system (1) is transformed to a non-dimensional one by defining new variables: (3) The quantity 2π/ω p is the period of oscillation of a single block attached to a pulling spring in the absence of sliding friction (Carlson and Langer, 1989).Carrying out the nondimensionalisation, one arrives at the following system, suggested by De Sousa Vieira (1999): Dots denote differentiation with respect to T .System (4) is valid only when block j is moving.Parameters k 1 and k 2 represent the ratio of spring strength connecting the blocks, k c , and the spring strength connecting the blocks and the driving plate, k p , for the first and the second block, respectively.Variable U represents block displacement, U is the velocity of the block defined in the standing reference frame, ν is the dimensionless pulling speed, and t is the time variable.Parameters ν c 1 and ν c 2 stand for the dimensionless characteristic velocities.The corresponding friction law φ reads: (5) Note that the friction force is assumed to depend only on the velocity of the block.

Proposed modified model
Starting from model (4), in the case of a single block, we introduce time delay τ and friction strength parameter c in the following way: where the remaining parameters are the same as in Eqs. ( 4) and ( 5).In present analysis, in contrast to De Sousa Vieira (1999), we do not discard backwards motion.Time delay parameters τ 1 and τ 2 and the friction strength parameters c 1 and c 2 are also introduced into the model with two blocks: where meaning of the remaining parameters is the same as in Eq. ( 4).Within this paper, we assume for simplicity that k 1 = k 2 = 1, which corresponds to homogenous elastic properties of the medium surrounding the fault.In the analysis below, the time delay parameters τ 1 and τ 2 and the friction strength parameters c 1 and c 2 are varied in the range [0, 10], with the iteration step of 10 −1 .
The unstable equilibrium around which the orbits of block j move in the phase space, is given in the following way: which is determined by setting Üj = 0 and Uj = ν in Eq. ( 7).This stationary solution, describing the uniform motion of the blocks with constant velocity ν, corresponds to aseismic creep along the fault.The detailed derivation of Eq. ( 8) from Eq. ( 7) is provided in the Appendix.6) with c = 0.
In the present paper, we investigate the behaviour of the trajectories near the stationary solution.In order to study the asymptotic dynamics, we have made certain that all the transients are discarded.As for the fashion in which the delaydifferential equations are numerically solved, the initial function for U is selected such that its values within the interval [−τ, 0] are set by Eq. ( 6) with c = 0 for the single-block model or by Eq. ( 7) with c 1 = c 2 = 0 in case of the model with two blocks.
The results are obtained by varying the control parameters c and τ for the model with a single block and the parameter set (c 1 , c 2 , τ 1 , τ 2 ) for the setup involving two blocks.The observed forms of behaviour are characterised by calculating the Fourier power spectrum and the largest Lyapunov exponent.The latter's value has been obtained by the methods of Wolf et al. (1985) and Rosenstein et al. (1993), with the use of two distinct algorithms aimed at providing additional verification of the results.

Model with one block
We first consider the single-block model ( 6), for which we adopt ν = 0.1, consistent with (De Sousa Vieira, 1999).The analysis shows that, without the introduced c and τ, and under the different values of parameter ν c , the original model displays just oscillatory behaviour of various amplitudes.If only the parameter c is introduced, under the variation of ν c the behaviour of the model does not change, except for the very low values ν c < 10 −8 , when the motion of the block settles at the stationary solution.Next we consider the effects of including the time delay τ in the friction term.If the value of the friction strength parameter is kept at c =1, there is no change in the behaviour of the model.This is expected, since the friction remains the same (Burridge and Knopoff, 1967;Carlson and Langer, 1989).However, if we vary the value of c, time delay may cause an oscillation death, the effect which has previously been investigated in (Yamaguchi and Shimizu, 1984;Aronson et al., 1990;Reddy et al., 1998).In particular, as τ is increased, the system undergoes a transition from equilibrium state to periodic motion and back to equilibrium state.This feature could have significant practical implications.In an earthquake analogy, the occurrence of oscillation death indicates that the increase of time delay could suppress the seismic activity, and, consequently, the onset of earthquakes.
Figure 2 is intended as an illustration of the (τ, ν c ) parameter domains admitting equilibrium or periodic motion for the single-block model.In the particular instance, friction strength is set to c = 3. Figure 3 shows the time series and phase portraits corresponding to the equilibrium state (stationary solution) and oscillations, observed at the respective points 1 and 2 from Fig. 2. Note that qualitatively similar results are obtained if time delay is held constant, while parameter c is varied.
We have also considered the behaviour of the single-block model under variation of both c and τ for ν c = 1, whereby the latter value is consistent with De Sousa Vieira (1999).Similar to the results above, one may distinguish between the two well defined regions of distinct dynamical behaviour, though appearing for different parameter values, cf.Figs. 4  and 2. It turns out that the periodic motion may be observed for the values of c up to 10.

Model with two blocks
The analysis similar to that for the single-block model has been carried out for the setup with two blocks represented by the system (7).In all the examined cases, we fix k 1 = k 2 = k =1 as in De Sousa Vieira (1999).Note that the original model, with friction strength parameter c 1 = c 2 = 1, as in De Sousa Vieira (1999), and without time delays, exhibits only the periodic motion.We first consider the effects induced by c 1 and c 2 , without introducing the delays.The quasiperiodic motion is obtained under the variation of both c 1 and c 2 while maintaining ν c 1 = ν c 2 = 1.It is found that the model exhibits the limit cycle oscillations only for c 1 = c 2 .
Next, if one introduces τ 1 and τ 2 and assumes homogeneous friction strength c 1 = c 2 , the model exhibits quasiperiodic motion for all the considered parameter values, except in the case τ 1 = τ 2 (when the motion is periodic).An instance of quasiperiodic motion for the parameter set and ν c 1 = ν c 2 = 1 is displayed in Fig. 5.The phase portrait is plotted in terms of U1 vs. U 1 − U e 1 , because the considered system ( 7) is nonautonomous, and our goal is to analyse the motion of two blocks in the vicinity of the stationary solution, which is explicitly time-dependent.The corresponding phase portrait in the ( U2 , U 2 −U e 2 ) plane is not shown, since the attractors of the second block are analogous to those of the first one, as implied in De Sousa Vieira (1999).Two peaks in the Fourier power spectrum (Fig. 6a) and approximately zero value of the maximal Lyapunov exponent (Fig. 6b) obtained for the time series in Fig. 5 confirm that the given motion is indeed quasiperiodic.
Next we consider the following setup: the first block is held at the equilibrium state (stationary solution) for the uncoupled case, c 1 = 0.1 and τ 1 = 3.0 (see Fig. 2), whereas c 2  and τ 2 for the second block are varied.From the results indicated in Fig. 7 one reads that the model mainly exhibits quasiperiodic motion, except for the two rather small regions admitting the equilibrium state.The latter scenario occurs only for the low values of the friction strength parameter (c 2 < 1).Note that the limit cycle oscillations (periodic solutions) appear only as a transient feature, resembling the weak quasiperiodic motion.
On the other hand, if (c 1 , τ 1 ) are selected such that the first block would show oscillatory behaviour in the uncoupled case (see Fig. 4), under the variation of c 2 and τ 2 the model primarily displays quasiperiodic motion.Periodic motion and equilibrium state occur as asymptotic solutions, but only within the isolated parameter domains (Fig. 8).Apart from the two latter asymptotic features, we have also captured a chaos-like behaviour, but only as a transient feature, see Fig. 9.The term "transient" refers to the fact that the chaos-like dynamics is observed within a relatively short time interval, cf.Fig. 9a, after which the system eventually converges to the quasiperiodic attractor.The occurrence of the transient chaos-like behaviour is confirmed by the continuous broadband noise in the power spectrum (Fig. 10).Moreover, we have calculated the value of the maximal Lyapunov exponent implementing two methods, that of Wolf et al. (1985) and the one of Rosenstein et al. (1993) (Fig. 11).Other parameter values are as in Fig. 5.The initial function for U is selected such that its values within the interval [−τ, 0] are set by Eq. ( 7) with c 1 = c 2 = 0.
In both cases, the obtained values of maximal Lyapunov exponent are positive and of the same order of magnitude, λ max = 0.0016 and 0.0019, respectively.
Note that the problem concerning the algorithm appropriate for calculating the maximal Lyapunov exponent in case of a transient chaos-like behaviour is still tentative and considered unresolved.In principle, the issue of qualifying certain transient motion as transient chaos should be treated by determining the finite-time Lyapunov exponent (Stefański et al., 2010).In the present paper, by using the methods of Wolf et al. (1985) and Rosenstein et al. (1993) we determined the Lyapunov exponent for the time series showing comparably large transients, whereby the standard procedures are complemented by performing additional averaging over a set of different initial conditions (Fig. 12).One notes that the approximately stationary values of the exponents are reached on the time scale significantly smaller than the length of the transient and that the values obtained for different initial conditions are quite similar.In particular, for all the examined cases, maximal Lyapunov exponents converge well to positive values (Fig. 12) of the same order of magnitude (10 −3 ).
The last stage of our analysis is focused on the issue of how selecting different initial conditions may affect the behaviour of the model with two blocks.Assuming the homogeneous initial conditions near the stationary solution ( U o 1 = U 0 2 = U 0 , U 0 1 = U 0 2 = U 0 ), we recover either the equilibrium state or periodic motion, cf.Fig. 13.In this case, quasiperiodic motion is not observed, since the initial conditions are chosen near the stationary solution.As apparent from Fig. 13, one is able to clearly distinguish between the domains supporting either of the two states, which may effectively provide an indication on the respective basins of attraction.This also implies that there possibly exists some distant limit cycle attractor which could not be captured by the analysis confined to the vicinity of the unstable stationary solution.Moreover, the obtained results suggest that the system is fairly sensitive to perturbations, meaning that even the small stress changes could induce motion along the fault, and, consequently, the onset of earthquakes.One should emphasise that the sensitivity of the block motion on initial conditions was already observed in the work of Szkutnik et al. (2003).In particular, the analysis there has shown that the character of the motion for the three-block model within a certain parameter range depends on the initial conditions.In other words, changing only the initial position of one of the blocks may induce transition from quasiperiodic and non-synchronized motion to the periodic solution where the two lateral blocks are synchronized.

Discussion and conclusion
In present paper, we have numerically investigated the dynamics of Burridge-Knopoff-like model with one and two blocks, having introduced the strength parameter c and time delay τ within the friction term.The analysis of the models' dynamics has been carried out in the vicinity of the stationary solution.Without the novel introduced parameters, the original models with one and two blocks exhibit only limit cycle oscillations.However, the coaction of c and τ may significantly alter the models' dynamics, giving rise to different forms of complex behaviour.In case of a single block, the The time series over a long simulation period is shown in (a), where the chaos-like region is marked with the rectangle.It is apparent that the system eventually converges to quasiperiodic behavior.co-effect of c and τ may induce the transition from periodic motion to equilibrium state, cf.Fig. 2, which is consistent with the delay-induced oscillation death.This phenomenon could have significant implications for the real earthquake dynamics, since it indicates the possibility that, under the certain conditions in the Earth's crust, motion along the fault could be suppressed, or reduced to aseismic creeping.In the model with two blocks, the coaction of c 1 , c 2 , τ 1 and τ 2 may give rise to the transition from periodic to quasiperiodic motion.From Figs. 7 and 8 one reads that such a transition  (Wolf et al., 1985).The method of of (Rosenstein et al., 1993) occurs for relatively small values of the friction parameter c(c < 1).If we consider friction parameter c and width of the fault zone as analogous parameters, low values of friction parameter c indicate relatively small gouge thickness.This is quite similar to Marone et al. (1990), where it has been shown that the friction coefficient increases with gouge  Wolf et al. (1985).Maximal Lyapunov exponents converge well to positive values of the order 10 −3 , the same as in Fig. 11.Note that time t is expressed in the units of iteration steps.layer thickness, indicating the change in friction coefficient µ for more than two times.It should also be stressed that for the model with two blocks, transient chaos-like behaviour is observed.It should be noted that the results of the conducted research set a solid base for the further investigation of complex dynamics of the presented models, including the global dynamical behaviour (far from the stationary solution) with heterogeneous initial conditions and different values of spring constants.

Fig. 1 .
Fig. 1.Burridge-Knopoff model with two blocks.Parameter m represents the mass of the block, k c and k p are the spring constants, and V denotes the velocity of the driving plate.

Fig. 2 .
Fig. 2. Parameter domains (τ, ν c ) admitting equilibrium state (ES) or periodic motion (PM) are illustrated for the single-block model at c = 3.Time series and phase portraits corresponding to points 1 and 2 are shown in Fig. 3.The initial function for U is selected such that its values within the interval [−τ, 0] are set by Eq. (6) with c = 0.

Fig. 4 .
Fig. 4. Parameter domains (τ, c) admitting equilibrium state (ES) or periodic motion (PM) at ν c = 1.Diagram refers to the single-block model, whereby the 0.1 step size is adopted for both c and τ .

Fig. 6 .
Fig. 6.Fourier power spectrum (a) and maximal Lyapunov exponent (b) for time series in Fig. 5: (a) Two peaks in the power spectrum indicate the appearance of quasiperiodic motion.(b) The maximal Lyapunov exponent converges approximately to λ ≈ 0. Note that in (b) time t is expressed in the units of iteration steps.

Fig. 7 .
Fig. 7. Parameter domains (τ 2 , c 2 ) admitting periodic motion (PM) or quasiperiodic motion (QM) of the second block.c 1 = 0.1 and τ 1 = 3.0 are fixed so that the first block lies in the equilibrium state.Diagram is constructed for the 0.1 step size for both c 2 and τ 2 .Other parameter values are as in Fig.5.The initial function for U is selected such that its values within the interval [−τ, 0] are set by Eq. (7) with c 1 = c 2 = 0.

Fig. 8 .
Fig. 8. Parameter domains (τ 2 , c 2 ) admitting equilibrium state (grey dots), periodic motion (black dots) and quasiperiodic motion (white area), for the fixed parameters c 1 = 0.2, τ 1 = 0.5.The latter values would warrant oscillatory motion for the first block in the uncoupled case.The diagram is constructed for the step size 0.1 for both c 2 and τ 2 .The remaining parameters are the same as in Fig. 5.

Fig. 9 .
Fig. 9. Temporal evolution of variable U1 (b) and corresponding phase portrait (c) for c 1 = 0.2, τ 1 = 0.5 c 2 = 0.2 and τ 2 = 3.In the uncoupled case, (c 1 , τ 1 ) would warrant the oscillatory motion of the first block, whereas (c 2 , τ 2 ) would hold the other block at the equilibrium.The remaining parameter values are:υ c 1 = υ c 2 = k 1 = k 2 = k = 1.The time series over a long simulation period is shown in (a), where the chaos-like region is marked with the rectangle.It is apparent that the system eventually converges to quasiperiodic behavior.

Fig. 10 .
Fig. 10.Fourier power spectrum for the time series in Fig. 9.The continuous broadband noise indicates relatively weak chaotic behaviour of the system.

Fig. 11 .
Fig. 11.Calculation of maximal Lyapunov exponent for the time series in Fig. 9. (a) indicates the value λ max = 0.0016, obtained by the method of(Wolf et al., 1985).The method of of(Rosenstein et al., 1993), illustrated in (b), implies λ max ≈ 0.0019.Note that in (a) time t is expressed in the units of iteration steps.In (b), effective expansion rate S( n) represents the average of the logarithm of D i ( n), defined as the average distance of all nearby trajectories to the reference trajectory as a function of the relative time n.The slope of dashed lines indicating the predominant slope of S( n) in dependence on ndt presents a robust estimate for the maximal Lyapunov exponent.The results are determined for 1000 reference points and neighbouring distance ε = 0.1-0.15.The obtained values of maximal Lyapunov exponent are of the same order of magnitude.
Fig. 11.Calculation of maximal Lyapunov exponent for the time series in Fig. 9. (a) indicates the value λ max = 0.0016, obtained by the method of(Wolf et al., 1985).The method of of(Rosenstein et al., 1993), illustrated in (b), implies λ max ≈ 0.0019.Note that in (a) time t is expressed in the units of iteration steps.In (b), effective expansion rate S( n) represents the average of the logarithm of D i ( n), defined as the average distance of all nearby trajectories to the reference trajectory as a function of the relative time n.The slope of dashed lines indicating the predominant slope of S( n) in dependence on ndt presents a robust estimate for the maximal Lyapunov exponent.The results are determined for 1000 reference points and neighbouring distance ε = 0.1-0.15.The obtained values of maximal Lyapunov exponent are of the same order of magnitude.