Local observability of state variables and parameters in nonlinear modeling quantified by delay reconstruction
Features of the Jacobian matrix of the delay coordinates map are exploited for quantifying the robustness and reliability of state and parameter estimations for a given dynamical model using an observed time series. Relevant concepts of this approach are introduced and illustrated for discrete and continuous time systems employing a filtered Hénon map and a Rössler system.
For many physical processes dynamical models (differential equations or iterated maps) are available but often not all of their variables and parameters are known or can be (easily) measured. In meteorology, for example, sophisticated large scale models exist, which have to be continuously adapted to the true temporal changes of temperatures, wind speed, humidity, and other relevant physical quantities. To obtain a model that “follows” reality, measured data have to be repeatedly incorporated into the model. In geosciences this procedure is called data assimilation, but the task to track state variables and system parameters by means of estimation methods occurs also in other fields of physics and applications. However, not all observables provide the information required to estimate a particular unknown quantity. In this article, we consider this problem of observability in the context of chaotic dynamics where sensitive dependance on initial conditions complicates any estimation method. A quantitative characterization of local observability employing delay coordinates is used to answer the question where in state and parameter space estimation of a particular state variable or parameter is feasible and where not.
To describe and forecast dynamical processes in physics and many other fields of science mathematical models are used, like, for example, ordinary differential equations (ODEs), partial differential equations (PDEs), or iterated maps. Some of these models are derived from first principles while others are the result of a general black-box modeling approach (e.g., based on neural networks). These models typically contain two kinds of variables and parameters: those that can be directly measured or are known beforehand (e.g., fundamental physical constants) and others whose values are unknown and very difficult to accessVTK04 (); RBKT10 (). To estimate the latter estimation methods have been devised that aim at extracting the required information from the dynamics, here represented by the model equations and the experimentally observed dynamical evolution of the underlying process. Different approaches for solving this dynamical estimation problem have been devised in the past, including (nonlinear) observer or synchronization schemes PJK96 (); NM97 (); HLN01 (); GB08 (); ACJ08 (); SO09 (); SRL09 (); K05 (); FMG05 (); CK07 (); Amritkar09 (), particle filters L10 (), a path integral formalism A09 (); QA10 (), or optimization based algorithms CGA08 (); B10 (); SBP11 ().
Before applying such an estimation method one may ask whether the available time series (observable) actually contains the required information to estimate a particular unknown value. In control theory this is called observability problem and it can for linear systems of ODEs be analyzed and answered by means of the so-called observability matrix Sontag (); NS90 (). Using derivative coordinates this approach can be generalized for nonlinear continuous systems HK77 (); Sontag (); N82 (). For state estimation of chaotic systems Letellier, Aguirre and Maquet LAM05 (); AL05 (); LAM06 (); LA09 (); FBL12 () considered continuous dynamical systems
that generate some observed signal where is the state of the system, is a smooth submanifold of , and denotes a measurement or observation function. Consider now -dimensional derivative coordinates Aeyels (); Takens (); SYC91 (); KS97 (); book_HDIA () of the observed signal
where stands for the -th temporal derivative of , is the reconstruction dimension, and is called derivative coordinates map SYC91 (). If this map is (at least locally) invertible, then we can uniquely determine the full state vector from the signal and its higher derivatives LAF09 (). Furthermore, small perturbations in should correspond to small perturbations in and vice versa. Therefore, we want the map to be an immersion, i.e. a smooth map whose derivative map is one-to-one at every point of , or equivalently, whose Jacobian matrix has full rank on the tangent space (here ). This does not imply that the map itself is one-to-one (), since the derivative coordinates (2) may provide states with two (or more) pre-images separated by a finite distance. Therefore, the observability analysis presented in the following is local, only, because it is based on analyzing (the rank) of the Jacobian matrix . The (global) one-to-one property of the map is not checked (what would be necessary, and for compact U also sufficient, to show that is an embedding SYC91 ()).
The Jacobian matrix can be computed by means of the vector field given in Eq. (1). In fact, for linear ODEs the Jacobian matrix conforms with the observability matrix known from (linear) control theory AL05 (). To estimate the rank of Letellier, Aguirre, and Maquet LAM05 (); AL05 () suggest to compute the eigenvalues of the - matrix
Nonzero eigenvalues indicate full rank of and thus local invertibility of at . To quantify the (local) invertibility of and thus the (local) observability of the full state Aguirre, Letellier, and Maquet LAM05 (); AL05 () introduced the observability index
where and denote the smallest and the largest eigenvalue of the matrix , respectively. Time averaging (along the available trajectory for ) yields
Instead of derivative coordinates we consider in the following delay coordinates Aeyels (); Takens (); SYC91 (); KS97 (); book_HDIA (). Furthermore, we extend the observability analysis to parameter estimation and compute a specific measure of uncertainty PSBL () for each state variable or parameter to be estimated. Last not least, we are not only interested in quantifying the average observability (like in Eq. (5)) but also in local variations that can be exploited during the state and parameter estimation process.
Ii Delay coordinates and observability
To motivate the concepts to be presented in the following we shall first consider a discrete time system (iterated map) where all model parameters are known and only state variables have to be estimated from the observed time series.
ii.1 Estimating state variables of a filtered Hénon map
For an dimensional discrete system
Again we assume that all states of interest lie within a smooth manifold . Here we consider delay coordinates forward in time. The function is therefore called forward delay coordinates map . It is also possible to use delay coordinates backward in time, or mixed forward and backward, and we shall address this issue in Sec. II.2.
As already discussed with derivative coordinates in the previous section a state is locally observable from the time series if is an immersion, i.e. if the Jacobian matrix has maximal (full) rank at . The corresponding Jacobian matrix can be computed using the iterated map Eq. (6). If the Jacobian matrix has maximal rank (assuming ) then is locally invertible (on ). “Local” means that still a delay vector could possess different pre images (separated by a finite distance).
To motivate and illustrate this analysis we consider the Hénon map
with parameters and . In the following we shall assume that the dynamics of this system is observed via a filtered signal provided by an FIR-filter
with filter parameter .
For two dimensional delay coordinates the delay coordinates map reads
The Jacobian matrix of the map is given by:
and its determinant
vanishes for all states on the singular line
For the FIR-filter is (asymptotically) deactivated and the critical line disappears (). For , however, the critical line crosses the chaotic attractor as shown for in Fig. 1a.
ii.2 Forward, backward, and mixed delay coordinates
Instead of using state space reconstruction based on forward delay coordinates (II.1) one could also use backward delay coordinates
or more general, a combination of forward and backward components
called mixed delay coordinates in the following, with reconstruction dimension . To obtain the backward components the inverse map and its Jacobian matrix are required (here we assume that the dynamics is time invertible). For discrete time systems (like the Hénon example) the underlying map (6) has to be inverted and for continuous time systems the inverse of the flow can in principle be computed by integrating the system ODEs (1) backward in time. In both cases, however, problems may occur in practice, because an explicit form of the inverse map may not exist and backward integration of dissipative systems results in diverging solutions and numerical instabilities (for longer integration times). Despite these difficulties inclusion of backward components turns out to be beneficial for the estimation task as will be demonstrated in the following for the Hénon examples and the Rössler system.
ii.3 Noisy observations and uncertainty
At states with the delay coordinates map is in principle invertible, but the inverse can be very susceptible to perturbations in like measurement noise. To quantify the robustness and the sensitivity of the inverse with respect to noise we consider the singular value decomposition of the Jacobian matrix of the delay coordinates map
where is an diagonal matrix containing the singular values and and are orthogonal matrices, represented by the column vectors and , respectively. is the transposed of coinciding with the inverse . Analogously, and the inverse Jacobian matrix reads
where . Multiplying by from the right we obtain or
This transformation is illustrated in Fig. 2 and it describes how small perturbations of in delay reconstruction space result in deviations from in the original state space. Most relevant for the observability of the (original) state is the length of the longest principal axis of the ellipsoid given by the inverse of the smallest singular value (see Fig. 2). Small singular values correspond to directions in state space where it is difficult (or even impossible) to locate the true state given a finite precision of the reconstructed state . For the filtered Hénon map we find that the closer the state is to the critical line given by (14) the more severe is this uncertainty. This is illustrated in Fig. 1a,b where at some points the ellipses spanned by the column vectors of the matrix are plotted. Figure 3 shows (color coded) the logarithm of the ratio smallest singular value (here: ) divided by the largest singular value vs. state variables and in a range of coordinates containing the chaotic Hénon attractor. In Fig. 3a dimensional forward delay (II.1) is considered where at the smallest singular value vanishes indicating the singularity (14) illustrated in Fig. 1a,b. If the reconstruction dimension is increased from to the singularity disappears as can be seen in Fig. 3b showing the ratio (color coded) for dimensional forward delay coordinates. For comparison, Figs. 3c,d show results obtained with mixed delay coordinates (II.2) and backward delay coordinates (II.2), respectively. The white areas in Fig. 3d correspond to ratios indicating poor observability (due to fast divergence of backward iterates of the Hénon map). Further increase of the reconstruction dimension ( or ) results in even larger values of (not shown here).
To assess the observability on the Hénon attractor we computed histograms of ratios at points. Figure 4 shows these histograms for the same coordinates used to generate the corresponding diagrams in Fig. 3. The best results (large ratios) provide mixed delay coordinates (Fig. 4c). We speculate that this is due to the fact that forward and backward components cover different directions in state space (similar to Lyapunov vectors corresponding to positive and negative Lyapunov exponents).
If the perturbations of the reconstructed state are due to normally distributed measurement noise they can be described by a symmetric Gaussian distribution centered at
where denotes the covariance matrix ( stands for the -dimensional unit matrix) and the standard deviation quantifies the noise amplitude. For (infinitesimally) small perturbations this distribution is mapped by the (pseudo) inverse of the delay coordinates map to the (non-symmetrical) distribution
centered at with the inverse covariance matrix
The marginal distribution of a given state variable centered at is given by
where the standard deviation is given by the square root of the diagonal elements of the covariance matrix
Using Eq. (22) the standard deviation of the marginal distribution can be written
and we consider in the following the factor
as our measure of uncertainty when estimating , because it quantitatively describes how the initial standard deviation is amplified when estimating variable PSBL ().
Figure 5 shows the uncertainties and vs. for different dimensional delay coordinates. In Fig. 4a,b results obtained with dimensional forward delay coordinates (II.1) are given. Large uncertainties occur mainly in a vertical stripe located near the singularity at (Eq. (14)) occurring for . Figures 5c,d,e,f show uncertainties of and obtained with mixed delay coordinates ((c),(d)) and backward delay coordinates ((e), (f)). For mixed delay coordinates (Fig. 5c,d) areas with very high uncertainties occur near the origin, but along the attractor and take only relatively low values. This is also confirmed by the -histograms on the attractor given in Fig. 6 for the same delay coordinates as used in Fig. 4. Again, the mixed delay coordinates turns out to be superior to the purely forward or backward coordinates. Furthermore, the dependance of the range of uncertainty values on the type of coordinates is different for different variables. While the uncertainty of increases when changing from forward to backward delay coordinates (Figs. 6a,e), the uncertainty of exhibits the opposite trend (Figs. 6b,f).
ii.4 State and parameter estimation
Until now only the state variables and are considered as unknowns to be estimated and the parameters and of the Hénon map and of the FIR filter are assumed to be known. We shall now consider the general case including unknown parameters of the dynamical system and unknown parameters of the measurement function . Let the dynamical model be a -dimensional discrete
or a continuous
dynamical system generating a flow
with discrete or continuous time. Furthermore, let’s assume that a time series of length is given observed via the measurement function
from a trajectory starting at .
This provides the -dimensional delay coordinates
with . Here the delay coordinates map is considered as a function of: (i) the state and the parameters of the underlying system, (ii) the parameters of the measurement function, (iii) the dimension parameters and , and (iv) the delay times and in backward and forward direction, respectively. The option to use different delay times, and for the backward and forward iterations is motivated by the fact that for dissipative systems backward solutions quickly diverge and therefore a choice may be more appropriate. For the same reason has typically to be smaller than . Since the reconstruction dimensions and the delay times are chosen a priori and are not part of the estimation problem they shall not be listed as arguments of to avoid clumsy notation. The Jacobian matrix of has the structure
and it can also be written as
and and denote the Jacobian matrices of the flow whose elements are derivatives with respect to the state variables and the parameters , respectively. For discrete dynamical systems (28) the Jacobians and can be computed using the chain rule and the recursion schemes
with ( unit matrix and . If backward iterations are required ( and ) similar recursion schemes exist based on the inverse map (providing ) and its Jacobian matrice and . Instead of recursion schemes one may also use symbolic or automatic differentiation AD (). For continuous systems (29) the required Jacobians can be obtained by simultaneously solving linearized systems equations as will be discussed in Sec. II.6. Inverse maps ( and ) may be computed via backward integration of the ODEs (at least for short periods of time before solutions diverge significantly). An extension for multivariate time series is straightforward.
ii.5 Parameter estimation for the Hénon map
We shall now extend the discussion to include not only state estimation but also parameter estimation. For better readability only forward delay coordinates are considered in the following, but all steps can also be done with mixed or backward delay coordinates, of course. We first consider the case where and are assumed to be known and only has to be estimated. In this case, unknown variables and unknown system parameter exists (while ). Therefore, delay coordinates with dimension or higher will be used. Figure 7 shows the ratio of singular values vs. in a plane in given by (and fixed parameters and ). For reconstruction dimension the ration is very small for extended subsets of the plane (white stripes in Fig. 7a). If the delay reconstruction dimension is increased to (Fig. 7b). these regions shrink or disappear. If the dimension is increased furthermore, the delay coordinates map is locally invertible in the full range of and values shown in Fig. 7 (results not shown here).
Now we include in the list of quantities to be estimated. Figure 8 shows the ratio of singular values for reconstruction dimensions and . For curves with very low singular value ratios exist crossing the Hénon attractor which disappear for .
Similar results are obtained if we include the remaining parameters in the estimation problem. Scanning the two-dimensional - subspace (plane) of the dimensional estimation problem for with fixed , , and indicates (almost) vanishing smallest singular values as long as . With dimensional delay coordinates the Jacobian matrix has clearly full rank almost everywhere within the chosen range as can be seen in Fig.. 10a. Figures 10b,c,d,e,f illustrate the uncertainties (Eq. (27)) of , respectively.
ii.6 Continuous dynamical systems
To compute the Jacobian matrix (34) of the delay coordinates map we have to compute the gradients (35) and (36) of the observation function and the Jacobian matrices and containing derivatives of the flow generated by the dynamical system (29) with respect to variables and parameters , respectively. The -matrix can be computed by solving the linearized dynamical equations in terms of a matrix ODE
where is a solution of Eq. (29) with initial value and is an matrix that is initialized as , where denotes the identity matrix. Similarly, the -matrix is obtained as a solution of the matrix ODE Kawakami ()
with . Solving (39) and (40) simultaneously with the system ODEs (29) we can compute , , etc. and use these matrices to obtain the Jacobian matrix of the delay coordinates map (34). For mixed or backward delay coordinates the required components can be computed by integrating the system ODE and the linearized ODEs backward in time.
ii.6.1 The Rössler system
To demonstrate the observability analysis for continuous systems we follow Aguirre and Letellier AL05 () and consider the Rössler system
with , , and .
Time series of different observables , or , or are considered, all of them consisting of values sampled with . Figure 11 shows the Rössler attractor where color indicates the uncertainty of estimating the variable (first column), or (second column), or (third column) using forward delay coordinates. The results in the first row are obtained when observing , while the diagrams in the second and third row show results for or time series. The reconstruction dimension equals and the delay time is . The bright yellow bullet indicates the state with the lowest uncertainty. This state and the following states plotted as thick red bullets underly the time series values that are used for the delay reconstruction. They span a window in time of length which is about one half of the mean period of the chaotic oscillations . The lowest uncertainties are obtained for states whose reconstruction involves trajectory segments following the vertical -excursion on the attractor. In contrast, trajectory segments starting from states with poor observability (large uncertainty) are located in the flat part of the Rössler attractor. Figures 11a and b show that using time series low values of occur on parts for the attractor where is high (and vice versa). Interestingly, this is not the case for delay reconstructions based on time series as can be seen in Figs. 11d and e, where low uncertainties of and occur in similar regions on the attractor.
In Fig. 12 distributions of uncertainty values of the Rössler system are shown that were obtained along an orbit of states sampled with . The distributions are shown as color coded histograms, estimated from the relative frequency of occurrence of the corresponding (in %). All diagrams show the dependance of the histograms on the delay time chosen for forward delay coordinates (horizontal axis). The reconstruction dimension is for all cases and all three parameters are assumed to be known (and are not part of the estimation task, i.e. ). In the first row (Figs. 12a,d,g) estimations are based on a time series from the Rössler system, and in rows two and three, and time series are used, respectively. The uncertainties of the given observable (Figs. 12a,e,i) mostly equal one () or are smaller (due to the additional information provided by the delay coordinates). In general, lowest uncertainties for all variables are obtained when using time series (Figs. 12a,d,g) while data provide highest uncertainties (Figs. 12c,f,i).
Figure 13 shows the same diagrams but now computed using four dimensional mixed delay coordinates with and . Similar to the results obtained with the Hénon map the uncertainties computed for mixed delay coordinates are typically smaller than those obtained with forward coordinates. Furthermore, the histograms shown in Fig. 13 suggest that for mixed delay coordinates using as measured variable provides the best results, followed by the time series. This is in contrast to forward coordinates (Fig. 12) where data yield the smallest uncertainties for the other variables ( and ). Similar results have been obtained with three dimensional forward or mixed coordinates.
The fact that mixed delay coordinates provide the lowest uncertainties when using time series is consistent with results for derivative coordinates obtained by Letellier et al. LAM05 () who found a ranking (for a different set of model parameters). For better comparison with their results we computed the (attractor) average
of the ratio
that provides the delay reconstruction analog of the observability index (4). Figure 14 shows vs. the delay time for different delay coordinates (rows) and different measured time series (columns). While for forward delay coordinates the largest values of the observability index occur if is measured (Fig. 14c), time series provide best observability if mixed delay coordinates are used (Fig. 14d). Note that in most cases high observability occurs for which is very close to the first zero of the autocorrelation function (that is often used as preferred value for delay reconstruction).
Fig. 15 shows similar histograms but now for the full estimation problem ( variables and parameters). Forward delay coordinates are used and the reconstruction dimension is increased to and an time series of length is used (with sampling time ). For delay times that are an integer multiple of half of the mean period relatively high uncertainties of occur, in particular for , , , and . This is due to the well known fact that for these delay times the attractor reconstruction results in points scattered near a straight (diagonal) line (an effect that also occurs when considering delay coordinates of a sinusoidal signal).
Starting from the question “Does some particular (measured) time series provide sufficient information for estimating a state variable or a model parameter of interest” we revisited the observability problem for nonlinear (chaotic) dynamical systems. In particular we considered delay coordinates and the ability to recover not measured state variables and parameters from delay vectors. This requires to “invert” the delay coordinates construction process which is at least locally possible, if the Jacobian matrix of the delay coordinates map has maximum (full) rank. Furthermore, we investigated how states near the delay vector are mapped back to the state and parameter space of the systems. In this way it is possible to quantify the amplification of small perturbations in delay reconstruction space in different directions of the state and parameter space. This reasoning gave rise to the concept of uncertainties of estimated variables and parameters. Both, observability and uncertainties may vary considerably in state space and on a given (chaotic) attractor. This feature was demonstrated with a discrete time system (filtered Hénon map) and a continuous system (Rössler system). Local observability and uncertainties also depend on the available measured variable (time series) and the type of delay coordinates. Best results were obtained with mixed delay coordinates, containing at least a one step backward in time.
The obtained information about (local) uncertainties in state and parameter estimation can be used in several ways for subsequent analysis. First of all, it may help to decide whether the planned estimation task is feasible at all or whether another observable has to be measured instead. For continuous time systems relevant time scales (delay times) can be identified where uncertainties are minimal.
The strong variations of local uncertainty values in state space (along a trajectory) occurring with the examples shown here are typical and should be taken into account by any estimation method. If the system is in a state where, for example, the uncertainty of the first variable is high then it might be better not to try to estimate this variable in this state or close to it, because the estimate might be poor and may spoil the overall results. Instead it makes more sense to wait until the trajectory enters a region of state space where can be estimated more reliably from the given time series.
The concrete implementation of such an adaptive approach depends on details of the estimation algorithm. For Newton-like algorithms, for example, it may consist of a simple strategy decreasing correction step sizes. Another potential application of uncertainty analysis is the identification of redundant parameters, i.e., parameter combinations that provide the same dynamical output.
Acknowledgements.The research leading to these results has received funding from the European CommunityÕs Seventh Framework Program FP7/2007-2013 under grant agreement no HEALTH-F2-2009-241526, EUTrigTreat. We acknowledge financial support by the German Federal Ministry of Education and Research (BMBF) Grant No. 031A147, by the Deutsche Forschungsgemeinschaft (SFB 1002: Modulatory Units in Heart Failure), and by the German Center for Cardiovascular Research (DZHK e.V.).
- (1) H.U. Voss, J. Timmer, J. Kurths, Int. J. Bif. Chaos 14(6), 1905-1933 (2004).
- (2) A. Raue, V. Becker, U. Klingmüller, J. Timmer, Chaos 20, 045105 (2010).
- (3) H. Nijmeijer and I.M.Y. Mareels, IEEE Trans. Circuits Syst. I 44, 882-890 (1997).
- (4) H.J.C. Huijberts, T. Lilge, H. Nijmeijer, Int. J. Bif. Chaos 11(7), 1997-2006 (2001).
- (5) U. Parlitz, L. Junge and L. Kocarev, Phys. Rev. E 54, 6253-6529 (1996).
- (6) D. Ghosh and S. Banerjee, Phys. Rev. E 78, 056211 (2008).
- (7) H. D. I. Abarbanel, D.R. Creveling, J.M. Jeanne, Phys. Rev. E 77, 016208 (2008).
- (8) F. Sorrentino and E. Ott, Chaos 19, 033108 (2009).
- (9) I. G. Szendro, M. A. Rodríguez, and J. M. López, J. Geophys. Res. 114, D20109 (2009).
- (10) R. Konnur, Phys. Lett. A 346, 275-280 (2005).
- (11) U.S. Freitas, E.E.N. Macau, C. Grebogi, Phys. Rev.E 71, 047203 (2005).
- (12) M. Chen and J. Kurths, Phys. Rev. E 76, 027203 (2007).
- (13) R.E. Amritkar, Phys. Rev. E 80, 047202 (2009).
- (14) P.J. van Leeuwen, Q. J. R. Meteorol. Soc. 136, 1991-1999 (2010).
- (15) H.D.I. Abarbanel, Phys. Lett. A 373, 4044-4048 (2009).
- (16) J.C. Quinn, H.D.I. Abarbanel, Q. J. R. Meteorol. Soc. 136, 1855-1867 (2010).
- (17) D.R. Creveling, P.E. Gill, H. D. I. Abarbanel, Phys. Lett. A 372, 2640-2644 (2008).
- (18) J. Schumann-Bischoff and U. Parlitz, Phys. Rev. E 84, 056214 (2011).
- (19) J. Bröcker, Q. J. R. Meteorol. Soc. 136, 1906-1919 (2010).
- (20) E.D. Sontag, Mathematical Control Theory: Deterministic Finite Dimensional Systems, (Springer, New York,1998) 2nd ed.
- (21) H. Nijmeijer and A.J. van der Schaft, Nonlinear Dynamical Control Systems, (Springer,New York, 1990).
- (22) R. Hermann and A. J. Krener, IEEE Trans. Autom. Contr. AC-22(5), 728-740 (1977).
- (23) H. Nijmeijer, Int. J. Control 36(5), 867-874 (1982).
- (24) C. Letellier, L.A. Aguirre, and J. Maquet, Phys. Rev. E 71, 066213 (2005);
- (25) L. A. Aguirre and C.Letellier, J. Phys. A: Math. Gen. 38, 6311 - 6326 (2005).
- (26) C. Letellier, L.A. Aguirre, and J. Maquet, Comm. Nonl. Sci. Num. Sim. 11, 555-576 (2006).
- (27) C. Letellier and L.A. Aguirre, Phys. Rev. E 79, 066210 (2009).
- (28) M. Frunzete, J.-P. Barbot, and C. Letellier, Phys. Rev. E 86, 026205 (2012).
- (29) H. Kantz and T. Schreiber, Nonlinear Time Series Analysis, Cambridge Nonlinear Science Series 7, (Cambridge University Press, Cambridge, 1997).
- (30) H.D.I. Abarbanel, Analysis of Observed Chaotic Data, (Springer Verlag, 1997), 2nd ed.
- (31) T. Sauer, J.A. Yorke, and M. Casdagli J. of Stat. Phys. 65(3/4), 579-616 (1991).
- (32) D. Aeyels, SIAM J. Contr. Optimiz. 19(5), 595603 (1981).
- (33) F. Takens, Lect. Notes Math. 898, 366-381 (1981).
- (34) C. Letellier, L. A. Aguirre, and U. S. Freitas, Chaos 19, 023103, (2009).
- (35) U. Parlitz, J. Schumann-Bischoff, S. Luther, Phys. Rev. E 89, 050902(R), (2014).
- (36) J. Schumann-Bischoff, S. Luther, U. Parlitz, Commun. Nonlin. Sci. Num. Sim. 18(10) 2733Ð2742 (2013).
- (37) H. Kawakami, IEEE Trans. Circ. Syst. CAS-31(3), 248-260 (1984).