

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


MODELING CHAOTIC SYSTEMS By Khaled Marzoug AlMughadhawi Bachelor of Science in Electrical Engineering King Fahad University of Petroleum and Minerals Dhahran, Saudi Arabia 1989 Master of Science in Electrical Engineering Southern Methodist University Dallas, Texas 1999 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY In Electrical Engineering July, 2005 ii COPYRIGHT By Khaled Marzoug AlMughadhawi July, 2005 iii MODELING CHAOTIC SYSTEMS Thesis Approved: Dr. Martin T. Hagan Thesis Adviser Dr. Gary E. Young Dr. Rao. K. Yarlagadda Dr. Carl D. Latino Dr. A. Gordon Emslie Dean of the Graduate College iv ACKNOWLEDGMENTS First and foremost, I would like to thank Allah (the God) for helping me and making it possible for me to reach this point of knowledge. Peace and blessing of Allah are on his messenger Muhammad, his last prophet and the seal of all messengers. My thanks extend to my parents: Marzoug, and Mohelah for their support, patience, and extreme love. My father grew up as an orphan in a small town known as Abar Ali (south of Madina, Saudi Arabia). When I was born, in 1966, he had to stop his education in the Madina teaching institute. He did that in order to support his new born child. He use to say to me “Khaled, I want you to compensate my loss of education by bringing me the highest degree in your major”. “I hope I did that Dad.” I would like to thank my advisor Dr. Martin Hagan for his support, patience, and help. He use to arrange weakly seminars for me, often lasting for three hours. He listened very carefully to my presentations and corrected me, suggesting new ways to solve the research problems. My thanks extend to all his graduate students who attended these seminars. I would like also to thank Dr. Yong Hu, and Dr. Tagi Menhaj for their valuable suggestions and comments. My thanks extend also to the dissertation committee members for their comments and suggestions. v Three professors in the Mathematics department gave me a much needed help and support. Dr. Saleh Mehdi gave me three lectures in differential geometry which helped me in understanding the embedding theorem. He and I corresponded often regarding the mathematical challenges I faced. Dr. Roger Zierau was very helpful to me. We spent long hours for many weeks discussing the linear Oseledec theorem. His support was very valuable to me. Dr. Leticia Barchini’s suggestions and comments made the proof of the theorem possible as well. My friends in the Islamic society of Stillwater made the difficult task of this degree much easier by their friendly relationship with me. I ask Allah to help them and make their future project of building a new mosque a success. My friends in the Saudi house showed me the true meaning of brotherhood. Our weakly meeting proved to be very important in bring us home far away from home. I would like also to thank my wife Nawal for her extreme support and patience with me and our five kids. Without her help, my goal could have been much more difficult to accomplish. My kids AbduRahman, Maha, Ahmad, Roba, and Arwa gave me the social stability which was also necessary to achieve this task. My good neighbor Sonya Brewster helped to proofread many chapters of my dissertation. Finally I would like to thank the Saudi government for their help and financial support of my research. vi TABLE OF CONTENTS Chapter Page 1. OBJECTIVE AND CONTRIBUTIONS ........................................................................1 1.1 Introduction .............................................................................................................. 1 1.2 Objective................................................................................................................... 1 1.3 Contributions ............................................................................................................ 1 1.4 Literature review ...................................................................................................... 2 1.5 Remaining chapters summary .................................................................................. 5 2. INTRODUCTION TO CHAOTIC SYSTEMS ..............................................................6 2.1 History of dynamical systems .................................................................................. 6 2.2 Dynamical definitions .............................................................................................. 7 2.2.1 Dynamical System ............................................................................................7 2.2.2 Equilibrium point ..............................................................................................8 2.2.3 Orbit of a dynamical system .............................................................................8 2.2.4 Attractor ............................................................................................................8 2.3 Characteristics of Chaotic Systems .......................................................................... 8 2.3.1 Sensitivity of chaotic systems to initial conditions ..........................................8 2.3.2 Chaotic signals look random but they are deterministic ...................................9 2.3.3 Chaotic systems have attractors with fractal dimension .................................11 2.3.3.1 The boxcounting dimension: .................................................................12 2.4 Examples of Chaotic systems................................................................................. 14 2.4.1 Weather system ...............................................................................................14 2.4.2 Biological models ...........................................................................................14 2.4.3 Laser signals ...................................................................................................15 2.4.4 Chaotic circuits ...............................................................................................15 2.4.5 Discrete chaotic systems: ................................................................................15 2.5 Chapter Summary................................................................................................... 15 vii 3. INTRODUCTION TO DYNAMICAL MODELING ..................................................17 3.1 Introduction ............................................................................................................ 17 3.2 Modeling................................................................................................................. 17 3.3 Applications of modeling ....................................................................................... 18 3.3.1 Detection of chaos ..........................................................................................18 3.3.2 Prediction ........................................................................................................18 3.3.3 Diagnosis ........................................................................................................19 3.4 Definitions .............................................................................................................. 19 3.4.1 Injection (1to1) function ..............................................................................19 3.4.1.1 Example: An injection ............................................................................19 3.4.1.2 Example: A function that is not injection ...............................................20 3.4.2 Immersion function..........................................................................................21 3.4.2.1 Example: An immersion .........................................................................21 3.4.2.2 Example: A function that is not immersion ............................................21 3.4.3 Embedding function ........................................................................................22 3.4.3.1 Example: An embedding ........................................................................22 3.4.3.2 Example: A nonembedding ...................................................................23 3.5 Modeling by Embedding ........................................................................................ 23 3.5.1 The delaycoordinate map ..............................................................................23 3.5.2 Example: The delaycoordinate map ..............................................................24 3.5.3 Example: Sufficient but not necessary condition for embedding ...................25 3.6 Chapter summary.................................................................................................... 26 4. INTRODUCTION TO THE MINIMUM EMBEDDING DIMENSION ESTIMATION 28 4.1 Introduction ............................................................................................................ 28 4.2 Modeling chaotic systems ...................................................................................... 28 4.3 Estimating the minimum embedding dimension.................................................... 30 4.3.1 The geometric technique .................................................................................30 4.3.1.1 The change of neighbors with dimension method (CND) ......................31 4.3.1.2 The change of distance with dimension (CDD) method .........................31 4.3.1.3 The change of distance with time method (CDT) ...................................32 4.3.2 The predictive technique: ...............................................................................34 4.4 Dynamic equivalence: ............................................................................................ 36 4.5 Chapter summary: .................................................................................................. 37 5. EXAMPLES OF THE MINIMUM EMBEDDING DIMENSION ESTIMATION ....38 5.1 Introduction ............................................................................................................ 38 5.2 The geometric technique ........................................................................................ 38 5.2.1 Using the change of neighbors with dimension method (CND) .....................39 5.2.2 Using the change of distance with dimension method (CDD) .......................43 viii 5.2.3 Using the change of distance with time method (CDT) .................................47 5.3 The Predictive technique ........................................................................................ 51 5.4 Chapter summary.................................................................................................... 52 6. ADVANCED ALGORITHMS FOR ESTIMATING THE MINIMUM EMBEDDING DIMENSION...............................................................................................................53 6.1 Introduction ............................................................................................................ 53 6.2 The delaytime (T).................................................................................................. 55 6.3 Two algorithms that use the local neighbor search ................................................ 57 6.3.1 The Algorithm ....................................................................................57 6.3.1.1 Example: Neighbor indices .....................................................................58 6.3.1.2 The vectorcoordinates projection method .............................................59 6.3.1.3 Example: Vectorcoordinates projection ................................................59 6.3.2 The Algorithm ....................................................................................61 6.3.2.1 The PCA projection method ...................................................................62 6.4 Limitations of the and the Algorithms ........................................... 66 6.5 Three new algorithms that use the global neighbor search method ....................... 68 6.5.1 The first new algorithm: The CND .................................................................68 6.5.2 The second new algorithm: The .........................................................71 6.5.3 The third new algorithm: The .............................................................74 6.6 The fourth new algorithm: The Predictive ............................................................. 77 6.7 Chapter summary.................................................................................................... 80 7. MINIMUM EMBEDDING DIMENSION RESULTS.................................................81 7.1 Introduction: ........................................................................................................... 81 7.2 Testing systems ...................................................................................................... 82 7.2.1 Noise free chaotic systems ..............................................................................82 7.2.2 Practical systems .............................................................................................86 7.3 Estimating for the noise free chaotic systems .................................................. 87 7.3.1 Using the algorithm ............................................................................87 7.3.2 Using the algorithm ............................................................................88 7.3.3 Using the CND algorithm ...............................................................................89 7.3.4 Using the algorithm ............................................................................90 7.3.5 Using the algorithm ............................................................................91 7.3.6 Using the predictive algorithm .......................................................................92 7.4 Estimating for the noisy systems ...................................................................... 94 7.4.1 Estimating for the noisy chaotic circuit (cc) .............................................94 CDDL CDTL CDDL CDTL CDDG CDTG dL CDDL CDTL CDDG CDTG dL dL ix 7.4.1.1 Using the algorithm ....................................................................94 7.4.1.2 Using the algorithm. ....................................................................95 7.4.1.3 Using the CND algorithm .......................................................................96 7.4.1.4 Using the algorithm ....................................................................96 7.4.1.5 Using the algorithm ....................................................................97 7.4.1.6 Using the predictive algorithm ...............................................................98 7.4.2 Estimating for the Santa Fe data sets ........................................................99 7.4.2.1 Using the algorithm ....................................................................99 7.4.2.2 Using the algorithm ...................................................................100 7.4.2.3 Using the CND algorithm .....................................................................101 7.4.2.4 Using the algorithm ..................................................................102 7.4.2.5 Using the algorithm ..................................................................103 7.4.2.6 Using the predictive algorithm .............................................................104 7.5 Testing the algorithms with random signals......................................................... 105 7.6 Tables and discussions ......................................................................................... 107 7.6.1 Tables ............................................................................................................107 7.6.2 Discussion of the Results ..............................................................................108 7.6.2.1 The effect of the original dimension of the system ..............................108 7.6.2.2 Local versus global neighbor search methods ......................................108 7.6.2.3 Estimation time .....................................................................................109 7.6.2.4 Sensitivity to the threshold value ..........................................................109 7.6.2.5 Noise detection .....................................................................................110 7.6.2.6 Dependence of the algorithms on the number of data points ...............110 7.6.3 General conclusions ......................................................................................110 7.7 Chapter Summary................................................................................................. 110 8. THEORY OF LYAPUNOV EXPONENTS...............................................................112 8.1 Introduction .......................................................................................................... 112 8.2 Sensitivity of some linear systems to initial conditions ....................................... 113 8.3 Lyapunov Exponent (LE) of a first order chaotic system .................................... 115 8.4 Lyapunov exponents for a multidimensional system........................................... 117 8.5 Invariant sets in modeling by embedding............................................................. 119 8.6 LEs from the Jacobian matrix............................................................................... 120 8.6.1 Oseledec theorem ..........................................................................................121 8.7 Linear Algebra definitions.................................................................................... 122 8.7.1 Definition: Inner product ..............................................................................122 8.7.2 Definition: Vector Norm ...............................................................................122 8.7.3 Definition: Matrix Norm ...............................................................................123 8.7.3.1 Some important properties of the matrix norm .....................................123 8.7.3.2 Lemma: Norm of a diagonal matrix .....................................................124 CDDL CDTL CDDG CDTG dL CDDL CDTL CDDG CDTG x 8.7.4 Singular value decomposition .......................................................................125 8.7.4.1 Important SVD properties .....................................................................125 8.7.5 Definition: Diagonalizable matrix ................................................................125 8.8 Multilinear algebra definitions ............................................................................. 126 8.8.1 Vector exterior products ...............................................................................126 8.8.2 Linear operator exterior power .....................................................................127 8.8.2.1 Definition: Adjoint of a linear operator ................................................128 8.8.2.2 Lemma: Adjoint of the wedge ..............................................................128 8.8.2.3 Linear operator properties .....................................................................129 8.8.2.4 Lemma: Eigen values of .........................................................130 8.9 The linear Oseledec matrix................................................................................... 131 8.9.1 Theorem: Eigen values of a linear Oseledec matrix .....................................131 8.9.2 For a symmetric matrix .................................................................................133 8.9.3 For a general matrix ......................................................................................134 8.9.3.1 Lemma: Eigen values of from singular values of .................134 8.9.3.2 Proposition: Limit of the root of ...................................135 8.10 Chapter summary................................................................................................ 141 9. ESTIMATING THE SET OF LYAPUNOV EXPONENTS......................................142 9.1 Introduction .......................................................................................................... 142 9.2 Estimating the LEs from Jacobian matrices ......................................................... 143 9.2.1 Estimating the LE by the QR decomposition ...............................................143 9.2.1.1 Example: Estimating LEs of the Henon map ........................................145 9.2.2 Spurious exponents .......................................................................................146 9.3 Estimating the LEs by the Geometric algorithms................................................. 147 9.3.1 Eckmann algorithm .......................................................................................147 9.3.2 Linear Oseledec algorithm ............................................................................149 9.4 Estimating the LEs by the predictive algorithm ................................................... 150 9.5 Pseudo codes of the LEs estimation algorithms ................................................... 154 9.5.1 Pseudo code of the Eckmann algorithm .......................................................155 9.5.2 Pseudo code of the linear Oseledec algorithm ..............................................157 9.5.3 Pseudo code of the predictive algorithm ......................................................158 9.6 Results of estimating the LEs by using the three algorithms ............................... 159 9.6.1 Discussion of the results ...............................................................................160 9.7 Chapter summary.................................................................................................. 161 10. SUMMARY, CONCLUSIONS, AND FUTURE RECOMMENDATIONS...........162 10.1 Summary............................................................................................................. 162 10.2 Conclusions ........................................................................................................ 164 10.3 Future recommendations .................................................................................... 165 (L*L)^q Ok (A)k kth σd (A)k ( ) xi LIST OF TABLES Table Page 5.1 The CND method (1/3) ...........................................................................................42 5.2 The CND method (2/3) ...........................................................................................42 5.3 The CND method (3/3) ...........................................................................................43 5.4 The CDD method tables for d = 1 ...........................................................................45 5.5 The CDD method tables for d = 2 ...........................................................................46 5.6 The CDD method tables for d = 3 ...........................................................................47 5.7 The CDT method for d = 1 ......................................................................................49 5.8 The CDT method for d = 2 ......................................................................................50 5.9 The CDT method for d = 3 ......................................................................................50 5.10 The neural network SSE as a function of d ............................................................52 7.1 Tabulation of the estimated for all the testing systems shown in Sections 7.3 through 7.5. The wrong estimates are circled. The abbreviation cc: chaotic circuit, and N. Able: not able ............................................................................................107 7.2 Six algorithms estimation time (in seconds) for of data set A of the Santa Fe competition, it shows the is the fastest and the predictive is the slowest. ......107 8.1 Estimating for . ..................................115 9.1 LEs estimation results using the three algorithms. The numbers inside the curly parenthesis are the estimated LEs and the number in bottom is the Lyapunov dimension ( ). ....................................................................................................................160 dL dL CDDG log(c) f(x(k)) = (3 ⁄ 4)(1 – 1 – 2x(k) ) dLyp xii LIST OF FIGURES Figure Page 2.1 Sensitivity of the chaotic tent map to initial conditions............................................ 9 2.2 The 60 Hz periodic wave and its frequency response............................................. 10 2.3 The random signal and its frequency response....................................................... 10 2.4 The chaotic signal and its frequency response........................................................ 11 2.5 The chaotic attractor of the Henon map.................................................................. 12 2.6 The slope of the line is the boxcounting dimension of the Henon map attractor. . 14 3.1 a) The attractor C in , b) the set of measurements, c) The reconstructed attractor by the delaycoordinate map................................................................................... 25 3.2 The circle in and its reconstruction in . ................................................ 26 4.1 Projection artifacts .................................................................................................. 30 4.2 The CDD method.................................................................................................... 32 4.3 The projection of a 3D curve (solid line) into 2D curve (dotted line) ................. 33 4.4 The CDT method: a) Distances between false neighbors increase with time, b) The intersection in the middle of the curve is a projection artifact................................ 33 4.5 The function approximation in the predictive technique.................................... 36 4.6 A neural network model with a TDL used to approximate the function , D is a one step delaytime........................................................................................................ 36 5.1 The nonlinear network used to estimate of the Henon map.............................. 52 5.2 LogPlot of the SSE versus d for the Henon map................................................... 52 6.1 The Lorenz model average mutual information ..................................................... 57 6.2 Pseudocode of the algorithm ...................................................................... 61 6.3 Pseudocode of the algorithm....................................................................... 65 6.4 Pseudocode of the CND algorithm ......................................................................... 70 6.5 Pseudocode of the algorithm...................................................................... 73 6.6 Pseudocode of the algorithm ...................................................................... 76 6.7 Pseudocode of the predictive algorithm.................................................................. 79 7.1 Lorenz chaotic attractor .......................................................................................... 83 7.2 The chaotic circuit diagram .................................................................................... 83 7.3 The chaotic circuit response ................................................................................... 84 7.4 The Rossler model attractor .................................................................................... 85 7.5 The estimated for the six noise free systems using the algorithm. The vertical axis is the percentage of the FNNs found from each dimension and the horizon ℜ3 S1 ℜ3 ℜ2 μ μ dL CDDL CDTL CDDG CDTG dL CDDL xiii tal axis is the dimension d, a) for the Lorenz model, , b) for the chaotic circuit, , c) for the Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, .. 88 7.6 The estimated for the six noise free systems using the algorithm. The vertical axis is the percentage of the FNNs and the horizontal axis is the dimension d, a) for Lorenz model, , b) for chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ......................................................................... 89 7.7 Estimating for the noise free systems using our first algorithm (CND). The vertical axis is the estimated and the horizontal axis is the number of neighbors used (w). a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, this algorithm does not give a stable answer.............. 90 7.8 Estimating for the noise free systems using our second algorithm ( ). The vertical axis is the percentage of the FNNs found and the horizontal axis is the dimension d. a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ................................................... 91 7.9 Estimating for the noise free systems using our third algorithm ( ). The vertical axis is the percentage of FNNs and the horizontal axis is the dimension d. a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, . ......................................................................... 92 7.10 The predictive algorithm results for the noise free systems, the vertical axis is the SSE of the prediction errors and the horizontal axis is the dimension d, a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ............................................................................................. 93 dL = 3 dL = 3 dL = 2 dL = 4 dL = 4 dL = 4 dL CDTL dL = 2 dL = 3 dL = 2 dL = 3 dL = 3 dL = 4 dL dL dL = 3 dL = 3 dL = 3 dL = 4 dL = 5 dL CDDG dL = 3 dL = 3 dL = 3 dL = 4 dL = 4 dL = 4 dL CDTG dL = 3 dL = 3 dL = 3 dL = 4 dL = 7 dL = 7 dL = 3 dL = 3 dL = 3 dL = 4 dL = 7 dL = 13 xiv 7.11 Estimating for the noisy cc. using Abarbanel et al first algorithm ( ). a) at 200 dB, , b) at 100 dB, , c) at 50 dB, , d) at 20 dB, the algorithm assumes the signal is noise............................................................................ 95 7.12 Estimating for the noisy cc. using algorithm. a) at 200 dB, , b) 100 dB, , c) at 50 dB, and d) at 20 dB the algorithm assumes the signal is noise . 95 7.13 Estimating for the noisy cc. using our first algorithm (the CND), a) at 200 dB, , b) at 100 dB, , c) at 50 dB, , d) at 20 dB, ......... 96 7.14 Estimating for the noisy cc. using our second algorithm ( ). It shows for SNR=200, 100, 50, and 20 dB ............................................................ 97 7.15 Estimated for the noisy cc. using the algorithm. For SNR=200, 100, and 50, . At SNR=20db, not stable................................................................... 98 7.16 Estimating using our fourth algorithm (Predictive), for 200, 100, and 50 dB, at 20 dB it assumes a random signal ................................................................ 98 7.17 Using the algorithm, a) for A data set, , b) for , , c) for , .................................................................................................................. 99 7.18 Using algorithm, a) for data set A, , b) for , , c) for , . .............................................................................................................. 100 7.19 Using the CND algorithm, a) for data set A, , b) for , , c) for , the result is not stable............................................................................................ 101 7.20 Using algorithm, a) for data set A, , b) for , , c) for , .............................................................................................................. 102 7.21 Using algorithm, a) for data set A, , b) for , , c) for , ................................................................................................................. 103 7.22 Using the predictive algorithm, a) for data set A, , b) for , , c) for , ....................................................................................................... 104 7.23 Testing Abarbanel et al two algorithms with a random signal, a) the algorithm fails to recognize the random signal, b) the algorithm succeeded in recognizing the random signal............................................................................................ 105 dL CDDL dL = 3 dL = 3 dL > 5 dL CDTL dL = 2 dL = 2 dL dL = 3 dL = 3 dL = 3 dL = 4 dL CDDG dL = 3 dL CDTG dL = 3 dL dL = 3 CDDL dL = 4 B1 dL = 4 D1 dL = 6 CDTL dL = 3 B1 dL = 6 D1 dL = 10 dL = 3 B1 dL = 5 D1 CDDG dL = 3 B1 dL = 5 D1 dL = 11 CDTG dL = 3 B1 dL = 5 D1 dL = 6 dL = 3 B1 dL = 4 D1 dL = 10 CDDL CDTL xv 7.24 Testing our four algorithms for estimating the dimension of the random signal, a) the CND algorithm recognizes the signal is random by not changing the estimated as w increases, b) for the algorithm, it did not recognize the random signal, c) for the algorithm, it was able to recognize the random signal, the same for the predictive algorithm in d). .................................................................................... 106 8.1 Linear systems can be sensitive to initial conditions ............................................ 113 8.2 Vectors a and b exterior product ( ) ............................................................... 126 9.1 The feed forward network used to estimate the LEs............................................. 152 9.2 The Eckmann algorithm Pseudo code .................................................................. 156 9.3 The linear Oseledec algorithm pseudo code ......................................................... 157 9.4 The predictive algorithm pseudo code.................................................................. 158 dL CDDG CDTG a^b xvi NOMENCLATURE Basic Concepts Scalars: small italic letters...a, b, Vectors: small bold nonitalic letters, or underbar symbol a, b, Matrices: capital Bold nonitalic letters ... A, B, C Language Vector means a column of numbers. Vector element Small italic with an index Matrix elements An element of matrix A: small italic showing the row then the column of the matrix: , or The column vector of a matrix I The row vector of a matrix I φ φ a(i) or ai a(i, j), aij A(i, j) kth ik kth i k xvii Norm Covariance matrix C Eigen value and eigen vector and (or bold small nonitalic letter) Set Empty set Subset of A, B, C Dimension Space dimension is superscript Box counting dimension Theoretical minimum embedding dimension Actual minimum embedding dimension Side length of a dcube (boxcounting dimension) a λ η ∅ ℜn ℜd dc dE dL σ xviii Interval in (boxcounting dimension) N Number of dcubes needed to cover a set (boxcounting dimension) Lyapunov dimension (Chapter 8) The number of Lyapunov exponents such that Power For a vector For a matrix Dynamics State vector in the original space Equilibrium point ℜ1 σ cσ d dLyp λj j = 1, 2, …, kL Σ > 0 kL (x)k (A)k x(n) x* xix The map for the state update in the original space Scalar Measurements Measurement function Delaytime T Delaycoordinate map (for embedding) Delayvector in the reconstructed attractor in The map for the state update in the reconstructed space Map for the output update in the reconstructed space The sampled measurement (predictive algorithm) , Equation (6.16) The sampled measurement with zero mean x(k + 1) = f(x(k)) y(m) y(m) = h(x(m)) yd(m) = Fd(h, x(n)) ℜd yd(m) y(m) y(m – T) … y(m – (d – 1)T) t = yd(m + 1) = φ(yd(m)) y(m + 1) = μ(yd(m)) ys(m) = y(1 + (m – 1)T) s(m) ys(m) 1 L  ys(k) k = 1 L = – Σ xx The delayvector created from Vectors distances matrix in the ddimensional space , see section 5.2.1 on page 39 Number of neighbors in used to search for the nearest neighbor of (for the CND method) w Indices matrix of the wnearest neighbors in the ddimensional space Nearest neighbor distances vector in ddimensional space Nearest neighbor distances in ddimensional space with distances computed in (d + 1)dimensional space Number of neighbors of Basis matrix for the PCA projection method ; are the eigen vectors of the covariance matrix C s(m) yd s (m) s(m) s(m – 1) … s(m – d + 1) t = Qd ℜd + 1 yd(m) Id rd ed + 1 yd(m) Nb Bd η1 η2 … ηd = ηi xxi The neighbors of the delayvector in : when the delaycoordinate projection method is used. Its element is : when the PCA projection method is used. It element is (see page 62) Index of the neighbor of Average vector in the space , the CDT pseudocode Average mutual information between two measurements I Threshold for significant FNN change (for the CND method) Threshold of distance increase (for the CDD method) Attractor mean value (for the CDT method) A fraction of attractor mean (for the CDT method) Nb yd m ( ) ℜd Yd Nb yd(m) Yd Nb, p yd p (m) kth yd(m) id k (m) ℜdE ydE (m) α ρ β ζ xxii Time step for the CDT method Distance between the nearest neighbors after time steps (CDT) Threshold for significant change in prediction error (Predictive technique) Signal variance Noise variance Jacobian matrix at or J Sampling time Delay time in MackeyGlass equation (Chapter 7) Lyapunov Exponent (global) Finite time Lyapunov Exponent (local) Δ Δ σΔ γ (σs)2 (σn)2 yd(m) Jyd(m) τs tf λ λt xxiii Infinitesimal distance Infinitesimal perturbation from to Abbreviations FNN: False Nearest Neighbor CND: Change of Neighbor with Dimension CDD: Change of Distance with Dimension CDT: Change of Distance with Time PCA: Principal Component Analysis SNR: Signal to Noise Ratio ECG: Electrocardiogram TDL: Tapped Delay Line SSE: Sum Squared Errors LE: Lyapunov Exponent SVD: Singular Value Decomposition Linear Algebra Singular value of a matrix or Multiplications of the n Jacobian matrices ε x(t) xε(t) δ(t) = xε(t) – x(t) σi σi(A) Jy n = Jy(n – 1)Jy(n – 2)…Jy(0) xxiv Perturbation vector from time m to time n Matrix of the perturbation vectors from Diagonal matrix D Diagonal matrix with elements magnitude of x ; Linear Oseledec matrix General Oseledec matrix Multilinear Algebra Vector exterior product The exterior power of a linear operator L Adjoint of a linear operator L Matrix representation of the wedge of the linear operator L ( ) δ(m) = x(n) – x(m) Nb yd(m) Em Da x ∈ Da Ok Λx ei^ej qth L^q L* L^q Bq 1 CHAPTER 1 OBJECTIVE AND CONTRIBUTIONS 1.1 Introduction In this chapter, we will show the objective of this research and give a brief list of our contributions. In Section 1.4, we give a literature review of some important publications in chaotic modeling. A summary of the remaining chapters is given in Section 1.5. 1.2 Objective The objective of this research is to model chaotic systems. We will build a chaotic model from a set of scalar measurements taken from the system. Evolution of the states in the model follows the evolution of the hidden states in the original system. To build a chaotic model, we start by estimating the model order. The next step is to estimate the set of Lyapunov exponents (the model parameters) that characterize the chaotic system. A good chaotic model depends on both the accuracy of the estimated model order and the estimate of Lyapunov exponent values. 1.3 Contributions Our contributions are: (1) Four new algorithms for estimating the model order. a) The Change of Neighbors with Dimension (CND). 2 b) The Change of Distance with Dimension using the Global neighbor search ( ). c) The Change of Distance with Time using the Global neighbor search ( ). d) The predictive. (2) The proof of a new theorem that relates the poles of a linear system to the set of Lyapunov exponents. (3) Implementation of the new theorem result to estimate the set of Lyapunov exponents for chaotic systems. (4) Development of an existing algorithm which uses a neural network to estimate the set of Lyapunov exponents. In total, we implemented four algorithms to estimate the model order and two algorithms to estimate the set of Lyapunov exponents for a chaotic system. The six algorithms were tested on different chaotic systems. The testing examples include noise free and noisy signals. In addition, the testing systems have different dimensions. 1.4 Literature review There are two primary areas of research in modeling chaotic systems. These areas, which will be discussed briefly, are: estimating the model order and estimating the set of Lyapunov exponents. In this section, we list some key references in these areas. We begin by listing literature dealing with estimating the model order. After that we list the literature that describes estimating the set of Lyapunov exponents. CDDG CDTG 3 We begin by reviewing publications that discuss estimating the model order. We found many papers by Abarbanel and his colleagues that talk about different methods for estimating the model order. Kennel, Brown, and Abarbanel [KBA92] introduced estimating the model order by the method of False Nearest Neighbors (FNNs). This method is a geometric one. It depends on measuring the distance between neighbors as the model order increases. The distance between false neighbors increases more than the distance between true neighbors as the model order increases. In 1993, Abarbanel and Kennel [AbKe93] developed the method of FNNs further and introduced a new method that uses time as well as distance to check for the existence of FNNs. In this proposal, we will present new algorithms that improve on the results of Abarbanel’s algorithms. In 1997, Cao [Cao97] suggested a new method for estimating the model order. His method is a geometric method and it is closely related to Abarbanel method [KBA92]. Instead of using a threshold value to determine false neighbors, he suggested the use of the mean of distance change as the model order increases. Abarbanel [Aba98] used a geometric method with a predictor to estimate the model order. Since the points inside the attractor have an evolution rule, their neighboring points will evolve according to a certain rule as well. To determine the model order, he used a predictor to estimate the evolution rule. When the model order is not large enough, prediction errors will increase. As the model order increases, the prediction errors drop. At certain point, further increase of the model order does not improve the prediction errors. At this point, the model order is found. In 2002, Kennel and Abarbanel [KeAb02] used the method of FNNs with a global neighbor search method. The projection method used in this paper is the delaycoordinate method. (Our second algorithm (CDDG ) used the method of 4 FNNs, with a global neighbor search, but the used the Principal Component Analysis projection method.) Now we review some papers that talk about estimating the set of Lyapunov exponents for a chaotic model. The most important paper in this field is the one by Oseledec [Ose68]. He proved that the set of Lyapunov exponents can be estimated from the product of an infinite number of Jacobian matrices along the attractor of the system. (In this proposal, we will present a new theorem that provides insight into the Oseledec theorem by relating Lyapunov exponents to the poles of a linear system.) In 1985, Eckmann and Ruelle [EcRu85] and Sano and Sawada [SnSd85] applied the Oseledec theorem to measurements taken from a chaotic system. Both papers use orthogonalization methods to overcome the numerical problem of multiplying a large number of matrices. Parlitz [Par92] introduced the identification of spurious Lyapunov exponents. He did that by reestimating the set of Lyapunov exponents from measurements backward in time. By doing that, he found that true exponents change their signs, while spurious ones change their values as well as signs. Darbyshire and Broomhead [DaBr96] estimated the set of Lyapunov exponents by the methods of least squares and total least squares. Their method uses the pseudoinverse to estimate the Jacobian matrices. After that, they applied the orthogonalization method proposed by Eckmann. Djamai and Coirault [DjCo02] introduced the use of neural networks to estimate the set of Jacobian matrices. After estimating the Jacobian matrices, they used the same method proposed by Eckmann to reorthogonalize the product of these matrices. CDDG 5 1.5 Remaining chapters summary In Chapter 2, we will present an introduction to chaotic dynamical systems. In Chapter 3, we introduce modeling chaotic systems and discuss modeling by embedding. Examples are given to illustrate the idea of embedding functions. In Chapter 4, we will present an introduction to four methods used to find the minimum dimension (model order) required for embedding. In Chapter 5, we implement the methods found in Chapter 4 to find the minimum embedding dimension of the Henon map. Six algorithms based on the four previous methods are presented in Chapter 6. (Four of the six algorithms represent our original work.) Full details of the algorithms are given, and more practical issues are discussed in Chapter 6. The results of implementing the six algorithms on nine chaotic systems are found in Chapter 7. The minimum embedding dimensions of the nine systems are listed in this chapter, and a comparison is made among the six algorithms. In Chapter 8, we will explore the theory of Lyapunov exponents (LEs) and prove a new theorem that relates the poles of a linear system to the set of LEs. In Chapter 9, we will discuss the estimation of the LE values and apply the result of the new theorem to estimate the LEs. In this chapter, we will also improve an existing algorithm that uses a neural network to estimate the LEs. We tested the two algorithms by estimating the LEs of different chaotic systems. We also compare the two algorithms to an existing algorithm using the test systems. In Chapter 10, we will give a conclusion of the dissertation and give some recommendations on possible future research. 6 CHAPTER 2 INTRODUCTION TO CHAOTIC SYSTEMS 2.1 History of dynamical systems Isaac Newton (16421727) introduced the idea of modeling the motion of objects by equations. Position, velocity, and acceleration were the fundamental parameters of his equations. From position, velocity, and acceleration, he could describe the state of a moving object at any given time. Later scientists modelled dynamical systems by using a set of differential equations. The solution of these equations describes the state of the system at any given time. If the solution of the set of differential equations remains in a bounded region, the sequence of states reduces to either i) a steady state, generally because of loss of energy by friction, or ii) a periodic or quasi periodic motion. An example of case (ii) is the motion of the moon around the earth and the motion of the earth around the sun. Case (i) and case (ii) remained the only two known bounded solutions until the use of modern computers made possible the numerical solution of sets of differential equations. In 1963 Edward Lorenz published a paper entitled “Deterministic non Periodic Flow” [Lo63]. He discovered a new bounded attractor that is not periodic but was filling a region in space. This was the first time the world knew about the third possible bounded attractor: iii) the chaotic attractor. The new solution (chaotic) is not simply periodic nor 7 quasi periodic with a large number of periods. Chaotic motion is possible with simple one dimensional systems. Although chaotic motion can be produced from simple systems, it is very complicated and becomes unpredictable after a short time. This happens because of the sensitivity of chaotic systems to changes in the initial conditions [ASY96]. In the next section, we define some terms related to dynamical systems. Three characteristics of chaotic systems are presented in Section 2.3. In Section 2.4, we present some examples of chaotic systems. Finally, Section 2.5 provides a summary of the chapter. 2.2 Dynamical definitions In the previous section, we gave a brief history of dynamical systems. In this section, we define some important dynamical terms that will be used in the remaining chapters. We start by defining the dynamical system. 2.2.1 Dynamical System The dynamical system is a system that consists of a sequence of states which are governed by a certain rule that determines the next state given the previous one. A dynamical system in can be described by either d first order ordinary differential equations (flow) or d difference equations (map). In the first case, the time is a continuous variable; , and the system is called a continuous dynamical system: . (2.1) In the second case, the time is a discrete variable; , and the system is called a discrete dynamical system: . (2.2) ℜd t ℜ1 ∈ d(x(t)) dt  = f(x(t)) n ∈ ℵ x(n + 1) = f(x(n)) 8 2.2.2 Equilibrium point A point in the state space is said to be an equilibrium point of a continuous dynamical system where . The equilibrium point of a discrete dynamical system, on the other hand, happens when . 2.2.3 Orbit of a dynamical system The sequence of points or that results from the solution of the set of equations representing the system is called the orbit (trajectory) of the dynamical system. The point or is called the initial condition of the system. 2.2.4 Attractor The attractor of a dynamical system is set containing the limits of all orbits that start sufficiently close to it. 2.3 Characteristics of Chaotic Systems In this section, we explore three characteristics of chaotic systems. These characteristics will help us to understand how to model these systems. 2.3.1 Sensitivity of chaotic systems to initial conditions Chaotic systems are known for their sensitivity to initial conditions. To illustrate this idea, let’s look into the tent map, which is a first order chaotic system: . (2.3) Figure 2.1 shows two curves (solid, and dotted with ‘x’) representing the responses of the system when two sightly different initial conditions are used. The solid curve was produced from the initial condition and the dotted curve with ‘x’ was produced from x(t) x* = dx(t) dt  x(t) x* = = 0 f x* ( ) x* = x(t) x(n) x0 x(0) x(n + 1) = f(x(n)) = (3 ⁄ 4)(1 – 1 – 2x(n) ) x1(0) = 0.23 9 . The dotted curve in the bottom of the figure is the difference between the two responses. We can see that the difference starts to grow after approximately 15 iterations. Notice that the responses of the system due to different initial conditions remain within a bounded region. Figure 2.1 Sensitivity of the chaotic tent map to initial conditions. 2.3.2 Chaotic signals look random but they are deterministic Before the discovery of chaotic systems, scientists thought of chaos as a random signal (noise). To illustrate this idea, let’s look into the frequency responses (Fourier spectrums [KaSc00]) of a periodic signal and a random signal, and then compare them to a chaotic signal. The left hand side of Figure 2.2 shows a 60 Hz periodic sinusoidal wave. The plot on the right hand side is its frequency response. As we can see, it has only one component at 60 Hz. x2(0) = 0.2301 0 5 10 15 20 25 30 35 40 45 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 The chaotic response of the tent map f(x(n))=x(n+1)=3/4(1−1−2x(n)) time n f x1(0) = 0.23 x2(0) = 0.2301 growth of the error 10 Figure 2.2 The 60 Hz periodic wave and its frequency response If we look into the frequency response of a random signal, shown in the right hand side of Figure 2.3, we can see that it has a component at every frequency from 0 through 120 Hz. Figure 2.3 The random signal and its frequency response We will now compare the frequency responses of the two signals above to the frequency response of a chaotic signal. The left hand side of Figure 2.4 shows the time domain plot of a chaotic signal (x component of the Lorenz model). Its frequency response is shown on the right hand side of the same figure. We can see that the chaotic signal has a component at every frequency from 0 through 120 Hz, which is similar to the random signal. 0 0.005 0.01 0.015 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 The periodic sinusoidal wave of period 60 Hz Time (sec) 20 40 60 80 100 0 50 100 150 Frequency content of the wave frequency (Hz) 0 200 400 600 800 1000 −5 −4 −3 −2 −1 0 1 2 3 4 5 The random signal (noise) Time 0 20 40 60 80 100 120 10−2 10−1 100 101 Frequency content of the signal frequency (Hz) 11 Figure 2.4 The chaotic signal and its frequency response Although the frequency response of the chaotic system looks random, we know that this system is deterministic. This is true since it is produced from a known set of differential equations. 2.3.3 Chaotic systems have attractors with fractal dimension To see that chaotic attractors have fractal dimensions, let’s look into the Henon map: , (2.4) , (2.5) where , , and the initial conditions are [Hen76]. In Figure 2.5, we can see that the attractor of the Henon map is not a simple line of dimension 1, nor it is a closed curve. It is also not a plane of dimension 2. It is an object that fills a region in the plane. 0 5 10 15 20 −30 −20 −10 0 10 20 30 The chaotic Lorenz signal Time (sec) 0 20 40 60 80 100 120 10−3 10−2 10−1 100 101 102 103 104 105 106 Frequency content of y frequency (Hz) x1(n + 1) 1 a(x1(n))2 = – + x2(n) x2(n + 1) = bx1(n) a = 1.4 b = 0.3 x1(0) = x2(0) = 0.5 12 Figure 2.5 The chaotic attractor of the Henon map. As a matter of fact, its dimension is not an integer, but rather is fractal. The dimension of the attractor of the Henon map can be found by the boxcounting dimension. 2.3.3.1 The boxcounting dimension: We can find the dimension of an interval , in , by dividing it into subintervals of length such that , where . Then the minimum number of subintervals needed to cover N is , which can be written as . Notice that the exponent in the expression is 1, which is also the dimension of N. For the case of a unit square in , , we can find the dimension of S by dividing it into small squares of side length (called squares), such that . Then the minimum number of squares needed to cover S is , which can be written as . Notice that the exponent in the expression is 2 which is also the dimension of S. The method used to find the −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 −1.5 −1 −0.5 0 0.5 1 1.5 The chaotic attractor of the Henon map x1 x2 N = [0, 1] ℜ1 Nj σ = 0.1 Nj ∩ Nj + 1 = ∅ j ∈ ℵ Nj cσ N = 10 cσ N (1 ⁄ σ)1 = ℜ2 S x y ℜ2 ∈ 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 ⎩ ⎭ ⎨ ⎬ ⎧ ⎫ = Sj σ = 0.1 σ Sj ∩ Sj + 1 = ∅ Sj cσ S 100 = cσ S (1 ⁄ σ)2 = 13 dimension of the sets in the previous two examples is called the boxcounting dimension. It is used to find the dimension of more complicated sets in , like fractal sets. The boxcounting dimension is denoted by which is the exponent in the relation: as , (2.6) where is the number of dcubes needed to cover the set in . Taking the limit of the log of Equation (2.6) as , we have . (2.7) To evaluate the boxcounting dimension ( ) in Equation (2.7), we draw a loglog plot of versus as . The value of is the slope of the resulting curve. Now we can apply the boxcounting dimension to the attractor of the Henon map. We start by choosing the side length of the squares to be some value . Then we count the minimum number of squares needed to cover the set of points in the attractor (which results from iterating equations (2.4) and (2.5) times). Next we decrease , and count the minimum number needed to cover the attractor for the new . Figure 2.6 shows the number of squares ( ) versus , for . ℜd dc cσ d (1 ⁄ σ)dc ≈ σ →0 cσ d σ ℜd σ → 0 dc cσ d ln ln(1 ⁄ σ)  σ →0 = lim dc cσ d 1 ⁄ σ σ → 0 dc σ σ < 1 105 σ σ σ cσ d σ σ = {0.1, 0.05, 0.01, 0.004, 0.002, 0.001} 14 Figure 2.6 The slope of the line is the boxcounting dimension of the Henon map attractor. The resulting boxcounting dimension of the Henon map was found to be , which is a fractal. After presenting three characteristics of the chaotic systems, in the next section we list some examples of chaotic systems. 2.4 Examples of Chaotic systems 2.4.1 Weather system We said in Section 2.1 that E. Lorenz was the first scientist to quantify chaos. He modelled the heat convection phenomena in fluids by using a set of three differential equations. His model represents earth’s atmosphere. He used it to forecast weather. 2.4.2 Biological models Many biological activities exhibit chaotic behaviors. One example is the Epileptic seizure. Epilepsy causes patients suffering from this disease to experience bouts of unconsciousness. Epileptic seizures result from an abnormal neuronal discharge from the central nervous system [Zyl01]. Another example of chaotic systems is the red blood cell production. Mackey and Glass [MG77] modelled red blood cell production and found that it ex 3 4 5 6 7 8 9 10 7 8 9 10 11 12 13 14 15 16 17 The box dimension of the Henon attractor σε = 0.1, 0.05, 0.01, 0.004, 0.002, 0.001 dc ≈ 1.189 15 hibits chaotic behavior at some parameter values of the model, as will be shown in Chapter 7. A third example of chaotic systems is the Electrocardiogram (ECG). In 1983, Glass et al [GGS83] experimented on the spontaneous beating of cells from embryonic chick hearts. They found that when these cells were stimulated by external periodic pulses, they showed chaotic motion (see also [Moo92]). 2.4.3 Laser signals Measurements taken from a laser output can be represented by a chaotic model. In Chapter 7, we will give more detail regarding the modeling of two data sets of laser outputs. 2.4.4 Chaotic circuits Chaos can be observed in electrical circuits as well. An example of this is the RLC circuit designed by Rulkov et al [RVRDV92]. 2.4.5 Discrete chaotic systems: In the previous section, we have shown two discrete chaotic systems: the tent and the Henon maps. Discrete chaotic systems are mainly used for analysis. In Chapter 7, we will show more discrete chaotic systems used to analyze chaos of different dimensions. 2.5 Chapter Summary In this chapter, we presented a historical background for dynamical systems. We defined some dynamical terms that will be used in subsequent chapters. Three characteristics of chaotic systems were illustrated with examples. We also explored a few examples of dynamical systems that show chaotic behaviors. In the next chapter, we will define dy 16 namical modeling, and show how to implement modeling of chaotic systems by using measurements taken from these systems. 17 CHAPTER 3 INTRODUCTION TO DYNAMICAL MODELING 3.1 Introduction In Chapter 2, we introduced dynamical chaotic systems. We have seen three characteristics of chaotic systems as well. We mentioned that these characteristics will help us to understand modeling of chaotic systems. In this chapter, we explore the modeling process. We also show some of its applications in real life. A theoretical background of the method used for modeling chaotic systems is illustrated with examples. In the next section, we present a brief introduction to modeling chaotic systems. Some applications of the modeling process are presented in Section 3.3. We do modeling by a technique called embedding. Some mathematical background on embedding is contained in Section 3.4. The use of the delaycoordinate map for modeling is illustrated in Section 3.5 with two examples. 3.2 Modeling We said in Chapter 2 that dynamical systems can be represented as a set of differential equations (for a continuous system), or a set of difference equations (for a discrete system). We also said that these equations determine the evolution of states, which converges to an attractor of the system. Knowing the evolution of the states is important for understanding the behavior of the system at a future time. Unfortunately these equations, 18 in most cases, are not known. All we can see from a chaotic system (like those listed in Chapter 2) is a set of scalar measurements. Modeling chaotic systems entails describing the hidden states of the system from these measurements only. From the set of measurements, the modeling process builds a new set of states, which are in some sense equivalent to the original hidden states. In Chapter 8, we will give more detail regarding the meaning of equivalent chaotic systems. 3.3 Applications of modeling 3.3.1 Detection of chaos In Section 2.3.2, we stated that, before the quantification of chaos, scientists thought of chaos as a random signal (noise). Given a set of scalar measurements, the modeling process enables us to distinguish between chaotic systems (deterministic) and random signals (noise). Not only that, but modeling can also extract the hidden deterministic part from a noisy signal. This application of the modeling process is important because, in most cases, the set of measurements is contaminated with noise from different sources. 3.3.2 Prediction The ability to predict the future has fascinated scientists for a long time. Modeling chaotic systems enables us to predict the hidden future states of the system. It does so using only the set of measurements. However, as we have said in Section 2.3.1, chaotic systems have sensitive dependence on the initial conditions. This problem limits the ability of chaotic modeling to make long term predictions of future states. In Chapter 8, we will give more detail. 19 3.3.3 Diagnosis Abarbanel [Aba98] conducted an experiment on the ECG of subjects undergoing a stress test for a specific pathology. He found that in the extreme pathology of ventricular fibrillation, characteristics of the model are different from those of a healthy person. This means the modeling process can be used to diagnose life threatening diseases (for more examples, see [Hol86 Chapters 9 and 11]). In the next two sections, we present some mathematical definitions, then we present the embedding method which is used to build models of chaotic systems. 3.4 Definitions Before we present the embedding method, let’s give a mathematical background of embedding functions. We begin by defining an injection function and an immersion function, then we define the embedding function. 3.4.1 Injection (1to1) function Let M and N be two sets, a function is an injection (1to1) if , (3.1) where . To understand the injection function, consider the next example. 3.4.1.1 Example: An injection Let the function be , (3.2) where . By applying the condition for an injection, we have g: M → N g x1 ( ) g x2 = ( ) x1 ⇒ x2 = x1 x2 , ∈ M g: ℜ2 ℜ2 → g(x) g1(x) g2(x) t x1 (x2)3 t = = x x1 x2 ℜ2 = ∈ 20 . (3.3) From the first element of the vectors on the right hand side we have , (3.4) while from the second element of the vectors we have . (3.5) In conclusion, we can see that , so the function g is an injection. On the other hand, let’s look at an example of a function which is not an injection. 3.4.1.2 Example: A function that is not injection Let the function be . (3.6) To test the condition for injection, we have . (3.7) From the first element of the vectors on the right hand side we have . (3.8) From the second element of the vectors we have . (3.9) Which means that does not have to equal , so the function g is not an injection. g x1 ( ) g x2 ( ) x1 1 x2 1 ( )3 t x1 2 x2 2 ( )3 t = = = x1 1 x1 2 = x2 1 ( )3 x2 2 ( )3 = x2 1 x2 ⇒ = 2 x1 x2 = g: ℜ2 ℜ2 → g(x) x1 (x2)2 t = g x1 ( ) g x2 ( ) x1 1 x2 1 ( )2 t x1 2 x2 2 ( )2 t = = = x1 1 x1 2 = x2 1 ( )2 x2 2 ( )2 = x2 1 ±x2 2 ⇒ = x1 x2 21 3.4.2 Immersion function Let M and N be two sets, a function is an immersion if its Jacobian matrix is an injection (full rank) [Nak90 p.149]. 3.4.2.1 Example: An immersion Let the function be , (3.10) where . The Jacobian matrix of g is . (3.11) The determinant of J is , therefore J is full rank. Hence, the function g is an immersion. 3.4.2.2 Example: A function that is not immersion Let the function be , (3.12) where . The Jacobian matrix of g is . (3.13) g: M → N g: ℜ2 ℜ2 → g(x) g1(x) g2(x) t x1 (x2)3 + x2 t = = x ℜ2 ∈ J 1 0 0 3(x2)2 + 1 x = x0 = 3(x2)2 + 1 ≠ 0 g: ℜ2 ℜ2 → g(x) x1 (x2)3 t = x ℜ2 ∈ J 1 0 0 3(x2)2 x = x0 = 22 The determinant of J is , which zero for . Therefore J is not full rank. Hence g is not an immersion. Now that we have defined the injection and the immersion functions, we will define the embedding function. 3.4.3 Embedding function Let M and N be two sets, a function is an embedding if it is an injection and an immersion [Nak90]. 3.4.3.1 Example: An embedding Let’s look at the function g in Example 3.4.2.1 ( ). In that example, we have seen that g is an immersion. To prove that it is an embedding, we need to show that it is also an injection. To do this, we start by assuming that . (3.14) If we look at the first element of the vectors on the right hand side of Equation (3.14), we can see that . While the second elements give . (3.15) By rearranging the above equation, we have . (3.16) Which can be written as . (3.17) 3(x2)2 x2 = 0 g: M → N g(x) x1 (x2)3 + x2 t = g x1 ( ) g x2 ( ) x1 1 x2 1 ( )3 x2 1 + t x1 2 x2 2 ( )3 x2 2 + t = = = x1 1 x1 2 = x2 1 ( )3 x2 1 + x2 2 ( )3 x2 2 = + x2 1 ( )3 x2 1 x2 2 ( )3 – x2 2 + – = 0 x2 1 x2 2 ( – ) x2 1 ( )2 x2 2 ( )2 x2 1x2 2 ( + + + 1) = 0 23 By looking at the second parenthesis of the left hand side of Equation (3.17), we can see that it can be simplified as follows . (3.18) So the only solution of Equation (3.17) is . As a result, we see that . Which means the function g is an injection. Since g is an injection and an immersion (Example 3.4.2.1), we see that g is an embedding. 3.4.3.2 Example: A nonembedding From Example 3.4.1.1, we have seen that is an injection. But from Example 3.4.2.2, we have seen that this function is not an immersion. So the function g is not an embedding. Notice that the main feature of the embedding is that it is an injection, and its linearization (Jacobian matrix) at every point along the attractor is also an injection. After presenting the embedding, let’s see now how we can use it for modeling. As we said above, we have only a set of scalar measurements from the system we want to model. To build a model for the system, we will use the delaycoordinate map, which is explained below. 3.5 Modeling by Embedding 3.5.1 The delaycoordinate map Let the state of the original system that we want to model be , and let the measurements taken from this system be governed by the map such that . (3.19) x2 1 1 2 x2 2 ⎝ + ⎠ ⎛ ⎞2 3 4  x2 2 ( ) 2 + + 1 ≠ 0 ∀x x2 1 x2 2 = x1 x2 = g(x) x1 (x2)3 t = x ℜk ∈ h: ℜk ℜ1 → y(n) = h(x(n)) 24 The delaycoordinate map [SYC91] is represented by , (3.20) where is the delaytime, and is the delayvector. The attractor built from the delayvectors is called the reconstructed attractor. According to the embedding theorem by Sauer at al [SYC91], the function is an embedding if , where is the boxcounting dimension (see Section 2.3.3.1) of the attractor of the original system. We will call the minimum dimension d that satisfies the embedding condition, the theoretical minimum embedding dimension: . The embedding map guarantees that evolution of the states in the original unknown attractor is equivalent to the evolution of the delayvectors in the reconstructed attractor (Chapter 8). In the next two examples, we illustrate the modeling process by using the delaycoordinate map. 3.5.2 Example: The delaycoordinate map Let the set represent an attractor in . Further, let the measurement function , for the purpose of explanation, be , where . The delaycoordinate map operating on C will be: ( ). The original attractor C is shown in Figure 3.1a, the measurement function h is shown in Figure 3.1b, and the reconstructed attractor using is shown in Figure 3.1c. Fd: ℜk ℜd → yd(n) Fd(h, x(n)) y(n) y(n – T) … y(n – (d – 1)T) t = = T ℵ ∈ yd n ( ) ℜd ∈ Fd d ≥ 2dc + 1 dc dE C y ℜ3 ∈ y1 cos t, y2 2t y3 sint cos2t t ℜ1 = { = = sin , = , ∈ } ℜ3 h:ℜ3 ℜ1 → y = h(x) = x1 + x2 + x3 x ℜ3 ∈ y3(n) F3(h, x(n)) y(n) y(n – 1) y(n – 2) t = = T = 1 F3 25 Figure 3.1 a) The attractor C in , b) the set of measurements, c) The reconstructed attractor by the delaycoordinate map. As we can see, this is a curve of dimension 1 in a three dimensional space. According to the embedding condition; . Which is the same as the dimension of the space required to see the curve C without ambiguity. 3.5.3 Example: Sufficient but not necessary condition for embedding Let the circle; represent an attractor in . By using the delaycoordinate map, we can reconstruct in . The attractor is shown in Figure 3.2 as a solid curve. The reconstructed attractor using the delay coordinate map ( ) is shown in the same figure as a dashed curve, the delaytime (T) is 1. −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 The image of the function Fd Operating on the Curve C in ℜ3 0 2 4 6 8 10 12 14 16 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 The set of observations from the function h:R3 → R1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 X The original attractor of the 1−D set C in ℜ3 Y Z h(x) a ) b ) c ) F3 x(n) y3(n) x2 x1 x3 y(n) y(n – 1) y(n – 2) ℜ3 dE = 2x1 + 1 = 3 S1 y ℜ3 ∈ y1 = cost, y2 = sint, y3 0.5 t ℜ1 = { = , ∈ } ℜ3 S1 ℜ3 S1 F3 26 Figure 3.2 The circle in and its reconstruction in . In Figure 3.2, we can see that is a curve of dimension 1 in a 3D space. Notice that we needed a 2D space to reconstruct this curve. This means the dimension required for embedding is 2 rather than 3, as suggested by the embedding condition ( ). In conclusion, we can see that the embedding condition gives a sufficient but not a necessary condition for embedding. In other words, it is possible that we can embed an attractor in a space of dimension less than . In the next chapter, the delaycoordinate map will be used to model chaotic systems. We will also present four methods used to find the actual, rather than the theoretical, minimum dimension required for embedding. The choice of T (the second parameter in the delaycoordinate map), will be shown in Chapter 6. 3.6 Chapter summary In this chapter, we presented an introduction to modeling chaotic systems. Some applications of modeling chaotic systems were given, mathematical definitions of the embedding functions were illustrated with examples, and we showed two examples of modeling −1 −1 −0.5 0 0.5 1 0 1 −0.5 0 0.5 1 Embedding of a circle by delay−coordinate map Original curve Reconstructed curve S1 S1 ℜ3 ℜ2 S1 d = 2 < dE = 3 dE 27 using the embedding method. We found that the embedding condition is a sufficient but not necessary condition. In the next chapter, we show different methods used to find the minimum dimension required for embedding. 28 CHAPTER 4 INTRODUCTION TO THE MINIMUM EMBEDDING DIMENSION ESTIMATION 4.1 Introduction In Chapter 3, we talked about modeling by the embedding method. Modeling is done from a time series of measurements taken from the system. We said that two parameters have to be found in order to model the system by embedding. Those parameters are the dimension of the delayvector (d) and the delaytime (T). They will be used to build the delayvectors. In this chapter we introduce four methods for estimating the value of d. Three of these methods are geometric and one method is predictive. This chapter will provide only a simple overview of the methods. Chapters 5 and 6 will provide more detail. In the next section, we present modeling of chaotic systems by using the delayvectors. In Section 4.3, we introduce two different approaches for estimating d: geometric and predictive. Then, in Section 4.4, we introduce equivalent chaotic systems. A chapter summary is given in Section 4.5. In Chapter 6, we will discuss the selection of the delaytime T. 4.2 Modeling chaotic systems The state evolution of a chaotic system in the space can be written in the form of the difference equation: ℜk 29 , (4.1) where is the state of the system, the time index , and N is the total number of points. In the steady state condition, the evolution of the state follows an attractor with a fractal dimension (chaotic attractor). Practically, the state and the function f of the system are invisible to us, and all we can see is a set of scalar measurements . We can write these measurements as (see Equations (3.19). To model a chaotic system, we need to build a system model which uses the delayvector . (4.2) The attractor of the system should be equivalent to the original attractor, and the state evolution from should follow that of the original attractor from [Hak98]. (A more careful definition of equivalence is contained in Chapter 8.) The key idea is to find d and T such that the two systems are equivalent. The embedding theorem guarantees this equivalence if . But in practice, we do not generally know the value of , so d has to be estimated. The projection from the original state to the delayvector is called the delaycoordinate map. As we said in the previous chapter, the embedding theorem gives a sufficient, but not a necessary, condition for embedding. In other words, it could be possible to find an embedding map at . Our goal is to estimate the minimum embedding dimension. We will label this dimension . In the next section, we will talk about estimating by using two techniques: a geometric and a predictive technique. In Chapter 6, we will discuss the delaytime T to complete the modeling process. x(m) = f(x(m – 1)) x m ( ) ℜk ∈ m = 1, 2, …, N x(m) dc x(m) y(m) y(m) = h(x(m)) yd(m) y(m) y(m – T) … y(m – (d – 1)T) t = yd(m) yd(m) → yd(m + 1) x(m) → x(m + 1) d ≥ 2dc + 1 dc x yd d ≤ 2dc dL dL 30 4.3 Estimating the minimum embedding dimension 4.3.1 The geometric technique To understand the geometric technique, let the circle in Figure 4.1 represent an attractor of a dynamical system in , and the line below it represent its projection into . Figure 4.1 Projection artifacts From the above figure, we can see that the points on the horizontal line (a, b, and c) do not occur in the same sequence as those on the circle ( ). That happened because the circle was projected into a space with insufficient dimension. That caused the distances between points to shrink, and the original order of points to change. We call this effect the projection artifact. We can see from this example that, if the dimension of the space is too small to represent the attractor of the system, the projection of the attractor to this space will have artifacts. Therefore, we need to increase the dimension of the space in order to get rid of these artifacts. One technique of estimating the minimum embedding dimension is the geometric technique. To estimate by using the geometric technique, we start from and find for every point its nearest neighbor. This neighbor has to be tested to see if it is a true neighbor or if it became a neighbor because of some projection artifacts. After that, d is increased and the previous steps are repeated. The value of d where all the neighbors are true neighbors, is the minimum embedding dimension . Under the geometric tech ℜ2 ℜ1 a′ b′ c′ a c b projection a′, b′, and c′ dL dL d = 1 yd(m) dL 31 nique, three methods can be used to detect the existence of projection artifacts. These methods are 1) the Change of Neighbors with Dimension method, or CND, 2) the Change of Distance with Dimension method, or CDD, and 3) the Change of Distance with Time method, or CDT. 4.3.1.1 The change of neighbors with dimension method (CND) One can perform a visual test using Figure 4.1 to determine that the horizontal line does not have enough dimension to represent the structure of the circle, and that we need to have a two dimensional space (meaning that ). However, we need to automate this test using an algorithm. To estimate using the first geometric method, the change of neighbors with dimension method (CND), the algorithm starts from the scalar measurements a, b, and c, and computes the distances between point a and the other two points. It will find that the nearest neighbor to a is c, while on the circle, the nearest neighbor to is . The algorithm concludes that the neighbors have changed, and that the points a and c became neighbors because of the projection artifact and not because they are true neighbors. That means the 1D space is not large enough to represent the structure of the circle since it has this projection artifact. By repeating the above steps on the 2D space, the algorithm would find that there are no projection artifacts left on the attractor in this space and concludes that . 4.3.1.2 The change of distance with dimension (CDD) method The second geometric method that can be used to estimate is the change of distance with dimension method (CDD). The CDD method is similar to the CND method except that it detects the existence of projection artifacts by comparing the distances between neighbors rather than comparing the neighbors themselves. To understand this method, let dL = 2 dL a′ b′ dL = 2 dL 32 the curve in Figure 4.2 represent an attractor of a dynamical system in , and the line below it represent its projection into . Figure 4.2 The CDD method We can see the points , , and on the curve in and their projection into the horizontal line a, b, and c respectively. As in the CND method, we want an algorithm to find that the 1D space is not large enough to represent the structure of the curve. From the horizontal line, the algorithm can find that b is the nearest neighbor of a and the distance between them in is . In , along the curve, the distance between and is . The algorithm can now compare the two distances to find that , therefore, b became a neighbor of a because of the projection artifact and not because they are true neighbors. As a result, the algorithm will find that the 1D space is not large enough to represent the structure of the curve since it has this projection artifact. By repeating the above steps on the 2D space, the algorithm would find that there are no projection artifacts left on the attractor in this space and concludes that . 4.3.1.3 The change of distance with time method (CDT) The third geometric method used to estimate is the change of distance with time method (CDT). This method is different from the previous two in that it tries to detect the existence of projection artifacts by moving neighbors forward in time and checking to see ℜ2 ℜ1 a' b′ c′ a b c ℜ2 ℜ1 projection v1 v2 a' b' c' ℜ2 ℜ1 v1 ℜ2 a′ b′ v2 v1 « v2 dL = 2 dL 33 if they will remain neighbors. To understand this method, let the nonintersecting solid curve in shown in Figure 4.3 represent an attractor of a dynamical system, and let the dotted curve with the “+” sign represent its projection into the XY plane ( space). Figure 4.3 The projection of a 3D curve (solid line) into 2D curve (dotted line) Figure 4.4 The CDT method: a) Distances between false neighbors increase with time, b) The intersection in the middle of the curve is a projection artifact The projection of the systems’s attractor shown in Figure 4.3 is redrawn in Figure 4.4 for further discussion. Notice that the middle of the curve in Figure 4.4a has a projection artifact, which can be found by looking into Figure 4.4b. If we start from point a, we can’t tell whether the next correct point in time is b or c. In other words, the curve intersection makes it impossible to determine the correct sequence of points. ℜ3 ℜ2 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 X 2−D projection of a curve in ℜ3 Y Z c a d b v1 v2 ℜ2 f a b f a ) b ) 34 By using the CDT method, we want an algorithm to be able to find that the 2D space is not large enough to represent the structure of the system’s attractor (which is in : ). This can be done by detecting the projection artifact in the middle of the curve. To do that, the algorithm should start by calculating the distance between point a and the other points on the curve. It will find that c is closest to a; let the distance between them be . Next it should move forward in time to find that a moves to b, and c moves to d. It should then compute the distance between b and d; let the distance be . Now it can compares the two distances and find that , which means that a and c are neighbors because of some projection artifacts and not because they are true neighbors. (See Figure 4.4a) Hence the algorithm concludes that the 2D space is not large enough to represent the structure of the system’s attractor since the curve has this projection artifact. By repeating this for the 3D space, the algorithm would find no projection artifacts left on the curve in this space and would conclude that . 4.3.2 The predictive technique: Now that we have talked about the three geometric methods, we will discuss the predictive technique. The predictive technique is also used to estimate the minimum embedding dimension such that the system of (using delay embedding) is equivalent to the system of (original system). For simplicity, we will set in the delayvector (see Equations (4.2)). By using the delaycoordinate map; , we can write (see Equations (3.20)). From the embedding definition (see Section 3.4.3), we know that if the map is an embedding, then it is an injection. Hence, the inverse map exists: . (4.3) ℜ3 dL = 3 v1 v2 v1 « v2 dL = 3 dL ydL (m) x(m) T = 1 yd(m) Fd yd(m) = Fd(x(m)) Fd Fd –1:ℜd ℜk → x(m) Fd –1 = (yd(m)) 35 We can substitute Equations (4.1) for the measurement function h at Equations (3.19) ( ) to write h as . If we let the composite function , we can write . (4.4) By substituting Equations (4.3) for Equations (4.4), we can write . (4.5) Now, if we let the composite function , we can write . (4.6) We conclude from Equations (4.6) that if the delayvector at time is known to us, we can predict the current measurement by approximating the unknown function . The predictive technique depends on the idea that when the system of is equivalent to the system of , we can use the delayvector to predict the current measurement . To be able to do that, we need to approximate the unknown function such that, . (4.7) To approximate we will use a neural network with a Tapped Delay Line (TDL) connected to its input. The TDL is used to produce the delayvector from the current measurement , as seen in Figure 4.6. Each tap of the TDL is a delayed version of and the total number of taps is d. The network is trained to predict when it is presented with d previous measurements of (coordinates of ). The resulting prediction error is recorded as a function of d. As the predictor order d increases, the prediction error will decrease. However, after a certain point, further increase of d results in only a very small decrease in the error. The minimum dimension where the error y(m) = h(x(m)) y(m) = h(f(x(m – 1))) h o f = q y(m) = q(x(m – 1)) y(m) q Fd –1 = ( (yd(m – 1))) q o Fd –1 = μ y(m) μ(yd(m – 1)) μ y(m – 1) y(m – 2) … y(m – d) t ⎝ ⎠ = = ⎛ ⎞ m – 1 y(m) μ yd(m) x(m) yd(m – 1) y(m) μ yˆ m ( ) μ ˆ = (yd(m – 1)) μ yd(m – 1) y(m) y(m) y(m) y(m) yd(m – 1) 36 does not improve any further is the minimum embedding dimension . At this point, the approximation of is accurate and the systems of and are equivalent to each other. Figures 4.5 and 4.6 show the block diagram of the function approximation and the neural network model used to create the approximation. Figure 4.5 The function approximation in the predictive technique Figure 4.6 A neural network model with a TDL used to approximate the function , D is a one step delaytime 4.4 Dynamic equivalence The goal of modeling by embedding is to find a chaotic model that is equivalent to the original unknown system. (The notion of equivalent chaotic systems will be covered in detail in later chapters.) To understand this idea, we know that in the original system the hidden state evolution is governed by the map f , (4.8) dL μ yd(m) x(m) μ + y(m) yˆ (m)  + μ ˆ (.) e(m) (prediction error) TDL μ(.) yd(m – 1) μ D y(m) yˆ y(m – 1) (m) y(m – 2) y(m – d) D D Neural Network Model μ ˆ (.) yd(m – 1) μ x(m + 1) = f(x(m)) 37 (see Equations (4.1)), while all we can see from the original system is a set of measurements which are governed by the map h: , (4.9) (see Equations (3.19)). On the other hand, the states in the reconstructed space are the delayvectors (see Equations (4.2)), while the output update in this space is governed by the map : , (4.10) (see Equations (4.6)). So the delayvector can be written as , (4.11) ( ). If the reconstructed model given by Equations (4.10) and (4.11) is equivalent to the original model given by Equations (4.8) and (4.9), we can use the evolution of in place of the evolution of to gain an understanding of the system characteristics. In other words, if the two systems are equivalent, we will be able to use to estimate the original system dimension and parameters. In the coming chapters, we will find the delayvector dimension d and the delaytime T which guarantee this dynamic equivalence. 4.5 Chapter summary In this chapter, we provided an introduction to modeling chaotic systems. We introduced two techniques (geometric and predictive) that can be used to estimate the minimum embedding dimension that is required for chaotic systems modeling. We also gave a brief introduction to equivalent chaotic systems. In Chapter 5, we will see the application of both techniques to the chaotic Henon map. In Chapter 6, we will provide complete and detailed algorithms for estimating d. We will also discuss practical considerations for implementing the algorithms, including the choice of the delaytime (T). y(m) = h(x(m)) yd(m) μ y(m) = μ(yd(m – 1)) yd(m) μ(yd(m – 1)) y(m – 1) … y(m – (d – 1)) t = T = 1 yd(m) x(m) yd(m) 38 CHAPTER 5 Examples of the minimum embedding dimension estimation 5.1 Introduction In Chapter 3, we introduced modeling of chaotic systems. We used delayvectors created from measurements to model these systems. We also introduced four methods to estimate the minimum embedding dimension: three geometric methods one predictive method. In this chapter, we will use the four methods to estimate the minimum embedding dimension of the Henon map. The Henon map is a chaotic system created from a set of two difference equations (see Section 2.3.3). By using the embedding method (see Section 3.5), the theoretical minimum embedding dimension is . But we suspect that the actual minimum embedding dimension is since the system’s dynamics are generated from a set of two difference equations. In Sections 5.2 and 5.3, we show the estimated minimum embedding dimension of the Henon map by using the geometric and the predictive techniques, respectively. In Section 5.4, we provide a chapter summary. 5.2 The geometric technique To show how we can estimate by using the geometric technique, we take 15 points from the coordinate of the Henon map to represent the measurements from this dE = 2dc = 2x1.189 = 3 dL = 2 dL X1 39 system. That means the measurements where . In the next step, we use to construct which is the delayvector (where ). We said in Chapter 4 that the idea of estimating by using the geometric technique depends on detecting the existence of the projection artifacts on the system’s attractor. We can do that by using any of the three methods that we mentioned before: the CND, the CDD, or the CDT method. We need first to compute the distance between the reference vector; and every other vector in the space . The nearest neighbor of is the vector with the shortest distance. 5.2.1 Using the change of neighbors with dimension method (CND) Before we get into the details of the CND method, we need to introduce some notation. First, the nearest neighbor of will be denoted . The nearest neighbor of will be denoted . (This means there are vectors that are closer to .) We will indicate the time index of the nearest neighbor as . (5.1) For example, if the nearest neighbor of is , then the nearest neighbor index is . To estimate by using the CND method, we need to check if the nearest neighbor of will remain a neighbor as the dimension d grows to . In other words, we need to check if the nearest neighbor of will appear as the first, second, , or neighbor of . To do that, we need to build the matrix which has the elements where and . If , we label as y(m) = x1(m) m = 1, 2, …, 15 y(m) yd(m) y(m) y(m – 1) … y(m – (d – 1)) t = T = 1 dL yd m ( ) ℜd yd(m) yd(m) y d 1 m ( ) ) kth yd(m) y d k (m) ) k – 1 yd(m) kth id k (m) index y d k = [ (m)] ) yd(5) y d 1 (5) = yd(9) ) id 1 (5) index y d 1 = [ (5)] = index[yd(9)] = 9 ) dL yd(m) d + 1 yd(m) … wth yd + 1(m) Qd qd(i, j) = yd(i) – yd(j) i, j = 1, 2, …, M i ≠ j i = j qd(i, j) 40 not a number (NaN) since we are not interested in the distance between a vector and itself. We need also to find the vector whose element is the time index of the nearest neighbor of . Further, we need to compute the through the neighbors of each and save their indices in the row of the matrix . In other words, , and are defined as follows: , (5.2) , (5.3) where , (5.4) and . (5.5) Each row of the matrix will be labeled by . So, the maid 1 mth id 1(m) yd(m) 1st wth yd + 1(m) mth Id + 1 Qd id k Id Qd qd(1, 1) qd(1, 2) … qd(1, M) qd(2, 1) qd(2, 2) … qd(2, M) : qd(M, 1) qd(M, 2) … qd(M, M) = Qd NaN yd(1) – yd(2) … yd(1) – yd(M) yd(2) – yd(1) NaN … yd(2) – yd(M) : : yd(M) – yd(1) yd(M) – yd(2) … NaN = id k id k (1) id k (2) : id k (M) index y d k [ (1)] index y d k [ (2)] : index y d k [ (M)] = = ) ) ) k = 1, 2, …, w Id id 1 id 2 … id w id 1 (1) id 2 (1) … id w (1) id 1 (2) id 2 (2) … id w (2) : : id 1 (M) id 2 (M) … id w (M) = = Id i m d id 1 (m) id 2 (m) … id w = (m) 41 trix can be written as: . In Tables 5.1a and 5.1c and Tables 5.2a and 5.2c, we can see listed for through 4. In Tables 5.1d, 5.2b and 5.2d, we can see through showing the indices of the three neighbors ( ) of each . Tables 5.3a, 5.3b, and 5.3c summarize the neighbors time indices found for through 4. For example, in Table 5.3a we can see that the nearest neighbor of is and it fails to appear as the first, second, or third neighbor of , so we label as a FNN. On the other hand, we can see that the nearest neighbor of is and when the dimension increases to 2, it appears as the second neighbor of , so we label as a true neighbor. As a conclusion, we can see in Table 5.3d that the total number of FNNs drops from 2 at to 0 at and remains 0 at . That means the minimum embedding dimension is , where there are no projection artifacts left on the attractor of the system. Notice that in this method we could have chosen , so that only the nearest neighbor is considered. We have found through experimentation that using w greater than 1 provides a more robust algorithm. This is a new result. In Chapter 6, we will provide more detail on the practical implementation aspects of the various algorithms. Id Id i 1 d i 2 d : i M d = Qd d = 1 I2 I4 w = 3 yd + 1(m) d = 1 y1(1) y1(8) y2(1) y1(8) y1(8) y1(9) y2(8) y1(9) d = 1 d = 2 d = 3 dL = 2 w = 1 42 Table 5.1 The CND method (1/3) Table 5.2 The CND method (2/3) Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 NaN 0.65 1.2015 0.1561 0.313 0.6479 1.2886 0.0262 0.0508 0.3597 0.4 0.7439 1.6362 0.9341 0.1046 8 2 0.65 NaN 1.8515 0.4939 0.963 0.0021 1.9386 0.6762 0.7008 0.2903 1.05 0.0939 2.2862 1.5841 0.7546 6 3 1.2015 1.8515 NaN 1.3576 0.8885 1.8494 0.0871 1.1753 1.1507 1.5612 0.8015 1.9454 0.4347 0.2674 1.0969 7 4 0.1561 0.4939 1.3576 NaN 0.4691 0.4918 1.4446 0.1823 0.2069 0.2036 0.556 0.5878 1.7922 1.0901 0.2607 1 5 0.313 0.963 0.8885 0.4691 NaN 0.9609 0.9755 0.2868 0.2622 0.6727 0.0869 1.0569 1.3231 0.621 0.2084 11 6 0.6479 0.0021 1.8494 0.4918 0.9609 NaN 1.9364 0.6741 0.6987 0.2882 1.0478 0.096 2.284 1.5819 0.7525 2 7 1.2886 1.9386 0.0871 1.4446 0.9755 1.9364 NaN 1.2624 1.2377 1.6483 0.8886 2.0325 0.3476 0.3545 1.1839 3 8 0.0262 0.6762 1.1753 0.1823 0.2868 0.6741 1.2624 NaN 0.0246 0.3859 0.3737 0.7701 1.61 0.9078 0.0784 9 9 0.0508 0.7008 1.1507 0.2069 0.2622 0.6987 1.2377 0.0246 NaN 0.4105 0.3491 0.7947 1.5853 0.8832 0.0538 8 10 0.3597 0.2903 1.5612 0.2036 0.6727 0.2882 1.6483 0.3859 0.4105 NaN 0.7596 0.3842 1.9959 1.2937 0.4643 4 11 0.4 1.05 0.8015 0.556 0.0869 1.0478 0.8886 0.3737 0.3491 0.7596 NaN 1.1438 1.2362 0.5341 0.2953 5 12 0.7439 0.0939 1.9454 0.5878 1.0569 0.096 2.0325 0.7701 0.7947 0.3842 1.1438 NaN 2.3801 1.6779 0.8485 2 13 1.6362 2.2862 0.4347 1.7922 1.3231 2.284 0.3476 1.61 1.5853 1.9959 1.2362 2.3801 NaN 0.7021 1.5316 7 14 0.9341 1.5841 0.2674 1.0901 0.621 1.5819 0.3545 0.9078 0.8832 1.2937 0.5341 1.6779 0.7021 NaN 0.8294 3 15 0.1046 0.7546 1.0969 0.2607 0.2084 0.7525 1.1839 0.0784 0.0538 0.4643 0.2953 0.8485 1.5316 0.8294 NaN 9 a) b) Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 NaN 1.9623 1.2991 0.9756 0.313 2.044 1.4552 0.7013 0.2947 1.1099 0.4108 2.4041 2.2773 1.2008 9 5 11 2 1.9623 NaN 2.2959 1.0165 2.0851 0.0871 2.267 1.3346 1.7113 0.8525 2.2106 0.4447 2.3018 1.9268 6 12 10 3 1.2991 2.2959 NaN 1.4363 1.0155 2.3467 0.202 1.1934 1.1685 1.6572 0.994 2.6451 1.1736 0.3735 7 14 11 4 0.9756 1.0165 1.4363 NaN 1.0693 1.0925 1.4728 0.3193 0.7038 0.2214 1.1942 1.4478 1.8968 1.1099 10 8 9 5 0.313 2.0851 1.0155 1.0693 NaN 2.1617 1.1858 0.7553 0.3896 1.2452 0.1295 2.5167 2.0623 0.9757 11 1 9 6 2.044 0.0871 2.3467 1.0925 2.1617 NaN 2.3116 1.4094 1.7902 0.9342 2.2867 0.3606 2.3114 1.9759 2 12 10 7 1.4552 2.267 0.202 1.4728 1.1858 2.3116 NaN 1.2626 1.2965 1.6901 1.1759 2.5928 0.9721 0.3631 3 14 13 8 0.7013 1.3346 1.1934 0.3193 0.7553 1.4094 1.2626 NaN 0.4113 0.5204 0.8782 1.7625 1.8363 0.9094 4 9 10 9 0.2947 1.7113 1.1685 0.7038 0.3896 1.7902 1.2965 0.4113 NaN 0.8635 0.5191 2.1483 2.0462 0.9978 1 5 8 10 1.1099 0.8525 1.6572 0.2214 1.2452 0.9342 1.6901 0.5204 0.8635 NaN 1.3731 1.2945 2.0661 1.327 4 8 2 11 0.4108 2.2106 0.994 1.1942 0.1295 2.2867 1.1759 0.8782 0.5191 1.3731 NaN 2.6407 2.0842 1.0026 5 1 9 12 2.4041 0.4447 2.6451 1.4478 2.5167 0.3606 2.5928 1.7625 2.1483 1.2945 2.6407 NaN 2.4815 2.2718 6 2 10 13 2.2773 2.3018 1.1736 1.8968 2.0623 2.3114 0.9721 1.8363 2.0462 2.0661 2.0842 2.4815 NaN 1.0867 7 14 3 14 1.2008 1.9268 0.3735 1.1099 0.9757 1.9759 0.3631 0.9094 0.9978 1.327 1.0026 2.2718 1.0867 NaN 7 3 8 c) d) i1 1 I2 a ) b ) c ) d ) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 11 9 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 10 4 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 9 13 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 8 2 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 11 9 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 12 10 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 13 8 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 10 9 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 8 5 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 8 2 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 1 9 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 2 10 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 3 8 a) b) Q4 1 2 3 4 5 6 7 8 9 10 11 12 1 NaN 2.4318 1.6489 2.5414 0.3725 2.3668 1.8663 1.7995 1.0368 2.8685 1.2434 2.433 5 9 11 2 2.4318 NaN 2.5327 1.4923 2.5528 0.331 2.3738 1.3529 2.0868 1.6802 2.9129 1.1956 6 12 8 3 1.6489 2.5327 NaN 2.5954 1.5612 2.4653 0.4389 1.7247 1.1757 3.0134 2.2894 2.8193 7 9 5 4 2.5414 1.4923 2.5954 NaN 2.5469 1.7832 2.3182 0.9872 2.3925 0.4232 2.6017 2.4496 10 8 2 5 0.3725 2.5528 1.5612 2.5469 NaN 2.5035 1.757 1.8512 1.2387 2.8763 0.9807 2.5428 1 11 9 6 2.3668 0.331 2.4653 1.7832 2.5035 NaN 2.3479 1.5024 1.9941 1.9947 2.9327 0.9783 2 12 8 7 1.8663 2.3738 0.4389 2.3182 1.757 2.3479 NaN 1.5296 1.3966 2.7334 2.36 2.7782 3 9 8 8 1.7995 1.3529 1.7247 0.9872 1.8512 1.5024 1.5296 NaN 1.4334 1.3952 2.245 2.2062 4 2 10 9 1.0368 2.0868 1.1757 2.3925 1.2387 1.9941 1.3966 1.4334 NaN 2.7782 2.1478 2.3707 1 3 5 10 2.8685 1.6802 3.0134 0.4232 2.8763 1.9947 2.7334 1.3952 2.7782 NaN 2.836 2.6148 4 8 2 11 1.2434 2.9129 2.2894 2.6017 0.9807 2.9327 2.36 2.245 2.1478 2.836 NaN 2.8555 5 1 9 12 2.433 1.1956 2.8193 2.4496 2.5428 0.9783 2.7782 2.2062 2.3707 2.6148 2.8555 NaN 6 2 8 c) d) I3 I4 a ) c ) d ) b ) 43 Table 5.3 The CND method (3/3) 5.2.2 Using the change of distance with dimension method (CDD) As we said in Chapter 4, the CDD method can also be used to estimate . This method depends on the idea that false nearest neighbor distances increase significantly as the dimension of the space increases. To estimate of the Henon map using the CDD method, we need to compute the matrix and find the vector as in the CND method. In addition, we need to find the vector whose element is the distance between and its nearest neighbor : . We also need to find out how much the distances to the nearest neighbors grow as the dimension increases. For these new distances we define the vector whose elements represent the distance between and . (Notice that is not the same as . It is with one component added.) For example, when in Table 5.4a, the nearest neighbor of is and the distance between them is . At in Table 5.4d, the distance between and is . If is a true nearest neighbor of , will be close to . To apply this idea, we need to see how d d d d FNN time 1 2 time 2 3 time 3 4 1 2 1 8 9 5 11 FNN 1 9 5 11 9 1 5 5 9 11 2 0 2 6 6 12 10 2 6 6 10 4 2 6 6 12 8 3 0 3 7 7 14 11 3 7 7 9 13 3 7 7 9 5 d) 4 1 10 8 9 FNN 4 10 10 8 2 4 10 10 8 2 5 11 11 1 9 5 11 1 11 9 5 1 1 11 9 6 2 2 12 10 6 2 2 12 10 6 2 2 12 8 7 3 3 14 13 7 3 3 13 8 7 3 3 9 8 8 9 4 9 10 8 4 4 10 9 8 4 4 2 10 9 8 1 5 8 9 1 1 8 5 9 1 1 3 5 10 4 4 8 2 10 4 4 8 2 10 4 4 8 2 11 5 5 1 9 11 5 5 1 9 11 5 5 1 9 12 2 6 2 10 12 6 6 2 10 12 6 6 2 8 13 7 7 14 3 13 7 7 3 8 13 7 14 3 7 3 8 14 7 c) 15 9 aa)) b )b) c ) d ) dL dL Qd id 1 rd rd(m) yd(m) y d 1 (m) yd id 1 = ( (m)) ) rd(m) yd(m) y d 1 m ( ) – = ) ed + 1 ed + 1(m) yd + 1(m) yd + 1 id 1 ( (m)) yd + 1 id 1 ( (m)) y d + 1 1 (m) ) y d 1 (m) ) d = 1 y1(1) y1(8) r1(1) = 0.026 d = 2 y2(1) y2(8) e2(1) = 0.701 yd id 1 ( (m)) yd(m) yd + 1 id 1 ( (m)) yd + 1(m) 44 much the distance between the nearest neighbors grows as d increases to . We can do that by forming the vector which has the elements . (5.6) If where is some predefined threshold, we label the nearest neighbor of as a false nearest neighbor (FNN). For instance, let . We can see from our example that as seen at the second column of Table 5.4g. Since , we label as a FNN. The results in Tables 5.4g, 5.5g, and 5.6g show for through 3. In Table 5.6h, we summarize the results of the previous tables. In it, we can see that at , the number of FNNs is 5, while at , the number of FNNs drops to 0 and remains 0 at . Therefore, the minimum embedding dimension is . d + 1 cd cd(m) ed + 1(m) – rd(m) rd(m) =  cd(m) > ρ ρ yd(m) ρ = 10 c1(1) = 25.9731 c1(1) > ρ = 10 y1(8) cd d = 1 d = 1 d = 2 d = 3 dL = 2 45 Table 5.4 The CDD method tables for d = 1 Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 i1 r1 1 NaN 0.650 1.202 0.156 0.313 0.648 1.289 0.026 0.051 0.360 0.400 0.744 1.636 0.934 0.105 8 0.026 2 0.650 NaN 1.852 0.494 0.963 0.002 1.939 0.676 0.701 0.290 1.050 0.094 2.286 1.584 0.755 6 0.002 3 1.202 1.852 NaN 1.358 0.889 1.849 0.087 1.175 1.151 1.561 0.802 1.945 0.435 0.267 1.097 7 0.087 4 0.156 0.494 1.358 NaN 0.469 0.492 1.445 0.182 0.207 0.204 0.556 0.588 1.792 1.090 0.261 1 0.156 5 0.313 0.963 0.889 0.469 NaN 0.961 0.976 0.287 0.262 0.673 0.087 1.057 1.323 0.621 0.208 11 0.087 6 0.648 0.002 1.849 0.492 0.961 NaN 1.936 0.674 0.699 0.288 1.048 0.096 2.284 1.582 0.753 2 0.002 7 1.289 1.939 0.087 1.445 0.976 1.936 NaN 1.262 1.238 1.648 0.889 2.033 0.348 0.355 1.184 3 0.087 8 0.026 0.676 1.175 0.182 0.287 0.674 1.262 NaN 0.025 0.386 0.374 0.770 1.610 0.908 0.078 9 0.025 9 0.051 0.701 1.151 0.207 0.262 0.699 1.238 0.025 NaN 0.411 0.349 0.795 1.585 0.883 0.054 8 0.025 10 0.360 0.290 1.561 0.204 0.673 0.288 1.648 0.386 0.411 NaN 0.760 0.384 1.996 1.294 0.464 4 0.204 11 0.400 1.050 0.802 0.556 0.087 1.048 0.889 0.374 0.349 0.760 NaN 1.144 1.236 0.534 0.295 5 0.087 12 0.744 0.094 1.945 0.588 1.057 0.096 2.033 0.770 0.795 0.384 1.144 NaN 2.380 1.678 0.849 2 0.094 13 1.636 2.286 0.435 1.792 1.323 2.284 0.348 1.610 1.585 1.996 1.236 2.380 NaN 0.702 1.532 7 0.348 14 0.934 1.584 0.267 1.090 0.621 1.582 0.355 0.908 0.883 1.294 0.534 1.678 0.702 NaN 0.829 3 0.267 15 0.105 0.755 1.097 0.261 0.208 0.753 1.184 0.078 0.054 0.464 0.295 0.849 1.532 0.829 NaN 9 0.054 a) b) c) Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 i1 e2 1 NaN 1.962 1.299 0.976 0.313 2.044 1.455 0.701 0.295 1.110 0.411 2.404 2.277 1.201 8 0.701 2 1.962 NaN 2.296 1.017 2.085 0.087 2.267 1.335 1.711 0.853 2.211 0.445 2.302 1.927 6 0.087 3 1.299 2.296 NaN 1.436 1.016 2.347 0.202 1.193 1.169 1.657 0.994 2.645 1.174 0.374 7 0.202 4 0.976 1.017 1.436 NaN 1.069 1.093 1.473 0.319 0.704 0.221 1.194 1.448 1.897 1.110 1 0.976 5 0.313 2.085 1.016 1.069 NaN 2.162 1.186 0.755 0.390 1.245 0.130 2.517 2.062 0.976 11 0.130 6 2.044 0.087 2.347 1.093 2.162 NaN 2.312 1.409 1.790 0.934 2.287 0.361 2.311 1.976 2 0.087 7 1.455 2.267 0.202 1.473 1.186 2.312 NaN 1.263 1.297 1.690 1.176 2.593 0.972 0.363 3 0.202 8 0.701 1.335 1.193 0.319 0.755 1.409 1.263 NaN 0.411 0.520 0.878 1.763 1.836 0.909 9 0.411 9 0.295 1.711 1.169 0.704 0.390 1.790 1.297 0.411 NaN 0.864 0.519 2.148 2.046 0.998 8 0.411 10 1.110 0.853 1.657 0.221 1.245 0.934 1.690 0.520 0.864 NaN 1.373 1.295 2.066 1.327 4 0.221 11 0.411 2.211 0.994 1.194 0.130 2.287 1.176 0.878 0.519 1.373 NaN 2.641 2.084 1.003 5 0.130 12 2.404 0.445 2.645 1.448 2.517 0.361 2.593 1.763 2.148 1.295 2.641 NaN 2.482 2.272 2 0.445 13 2.277 2.302 1.174 1.897 2.062 2.311 0.972 1.836 2.046 2.066 2.084 2.482 NaN 1.087 7 0.972 14 1.201 1.927 0.374 1.110 0.976 1.976 0.363 0.909 0.998 1.327 1.003 2.272 1.087 NaN 3 0.374 d) e) f) Time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 r1 0.026 0.002 0.087 0.156 0.087 0.002 0.087 0.025 0.025 0.204 0.087 0.094 0.348 0.267 e2 0.701 0.087 0.202 0.976 0.130 0.087 0.202 0.411 0.411 0.221 0.130 0.445 0.972 0.374 c1 25.9731 42.55 1.32184 5.25385 0.48851 42.55 1.32184 15.452 15.452 0.08529 0.48851 3.73085 1.79339 0.39888 FNN FNN FNN FNN FNN g) i1 1 e2 r1 i1 1 a ) d ) b ) c ) e ) f ) g ) 46 Table 5.5 The CDD method tables for d = 2 Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 i2 r2 1 NaN 1.962 1.299 0.976 0.313 2.044 1.455 0.701 0.295 1.110 0.411 2.404 2.277 1.201 9 0.295 2 1.962 NaN 2.296 1.017 2.085 0.087 2.267 1.335 1.711 0.853 2.211 0.445 2.302 1.927 6 0.087 3 1.299 2.296 NaN 1.436 1.016 2.347 0.202 1.193 1.169 1.657 0.994 2.645 1.174 0.374 7 0.202 4 0.976 1.017 1.436 NaN 1.069 1.093 1.473 0.319 0.704 0.221 1.194 1.448 1.897 1.110 10 0.221 5 0.313 2.085 1.016 1.069 NaN 2.162 1.186 0.755 0.390 1.245 0.130 2.517 2.062 0.976 11 0.13 6 2.044 0.087 2.347 1.093 2.162 NaN 2.312 1.409 1.790 0.934 2.287 0.361 2.311 1.976 2 0.087 7 1.455 2.267 0.202 1.473 1.186 2.312 NaN 1.263 1.297 1.690 1.176 2.593 0.972 0.363 3 0.202 8 0.701 1.335 1.193 0.319 0.755 1.409 1.263 NaN 0.411 0.520 0.878 1.763 1.836 0.909 4 0.319 9 0.295 1.711 1.169 0.704 0.390 1.790 1.297 0.411 NaN 0.864 0.519 2.148 2.046 0.998 1 0.295 10 1.110 0.853 1.657 0.221 1.245 0.934 1.690 0.520 0.864 NaN 1.373 1.295 2.066 1.327 4 0.221 11 0.411 2.211 0.994 1.194 0.130 2.287 1.176 0.878 0.519 1.373 NaN 2.641 2.084 1.003 5 0.13 12 2.404 0.445 2.645 1.448 2.517 0.361 2.593 1.763 2.148 1.295 2.641 NaN 2.482 2.272 6 0.361 13 2.277 2.302 1.174 1.897 2.062 2.311 0.972 1.836 2.046 2.066 2.084 2.482 NaN 1.087 7 0.972 14 1.201 1.927 0.374 1.110 0.976 1.976 0.363 0.909 0.998 1.327 1.003 2.272 1.087 NaN 7 0.363 a) b) c) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 i2 e3 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 9 0.854 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 0.202 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 0.331 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 0.241 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 0.371 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 0.202 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 0.331 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 0.43 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 0.854 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 0.241 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 0.371 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 0.977 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 0.974 d) e) f) Time 1 2 3 4 5 6 7 8 9 10 11 12 13 r2 0.295 0.087 0.202 0.221 0.13 0.087 0.202 0.319 0.295 0.221 0.13 0.361 0.972 e3 0.854 0.202 0.331 0.2413 0.3709 0.202 0.331 0.4301 0.854 0.2413 0.3709 0.9768 0.9736 c2 1.89492 1.32184 0.63861 0.09186 1.85308 1.32184 0.63861 0.34828 1.89492 0.09186 1.85308 1.70582 0.00165 g) i2 1 e3 i2 1 r2 a ) b ) c ) d ) e ) f ) g ) 47 Table 5.6 The CDD method tables for d = 3 5.2.3 Using the change of distance with time method (CDT) As we said in Chapter 4, the CDT method can also be used to estimate . This method depends on the idea that false nearest neighbor distances grow significantly as time increases. As in the previous methods, we begin by computing the nearest neighbors for . Then we track how the distance between the original nearest neighbors increases with time. This method can be applied by computing the distance in dimension , rather than in dimension . This has the advantage that the measured distance won’t be affected by projection artifacts. Next, we need to have a method for determining whether or not the distances between two vectors can be considered large. For our purpose, we will use the average vector Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 i3 r3 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 0.325 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 0.202 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 0.331 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 0.241 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 0.325 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 0.202 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 0.331 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 0.43 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 0.854 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 0.241 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 0.371 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 0.977 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 0.97 a) b) c) Q4 1 2 3 4 5 6 7 8 9 10 11 12 i3 e4 1 NaN 2.4318 1.6489 2.5414 0.3725 2.3668 1.8663 1.7995 1.0368 2.8685 1.2434 2.433 5 0.3725 2 2.4318 NaN 2.5327 1.4923 2.5528 0.331 2.3738 1.3529 2.0868 1.6802 2.9129 1.1956 6 0.331 3 1.6489 2.5327 NaN 2.5954 1.5612 2.4653 0.4389 1.7247 1.1757 3.0134 2.2894 2.8193 7 0.4389 4 2.5414 1.4923 2.5954 NaN 2.5469 1.7832 2.3182 0.9872 2.3925 0.4232 2.6017 2.4496 10 0.4232 5 0.3725 2.5528 1.5612 2.5469 NaN 2.5035 1.757 1.8512 1.2387 2.8763 0.9807 2.5428 1 0.3725 6 2.3668 0.331 2.4653 1.7832 2.5035 NaN 2.3479 1.5024 1.9941 1.9947 2.9327 0.9783 2 0.331 7 1.8663 2.3738 0.4389 2.3182 1.757 2.3479 NaN 1.5296 1.3966 2.7334 2.36 2.7782 3 0.4389 8 1.7995 1.3529 1.7247 0.9872 1.8512 1.5024 1.5296 NaN 1.4334 1.3952 2.245 2.2062 4 0.9872 9 1.0368 2.0868 1.1757 2.3925 1.2387 1.9941 1.3966 1.4334 NaN 2.7782 2.1478 2.3707 1 1.0368 10 2.8685 1.6802 3.0134 0.4232 2.8763 1.9947 2.7334 1.3952 2.7782 NaN 2.836 2.6148 4 0.4232 11 1.2434 2.9129 2.2894 2.6017 0.9807 2.9327 2.36 2.245 2.1478 2.836 NaN 2.8555 5 0.9807 12 2.433 1.1956 2.8193 2.4496 2.5428 0.9783 2.7782 2.2062 2.3707 2.6148 2.8555 NaN 6 0.9783 d) e) f) Time d FNN 1 2 3 4 5 6 7 8 9 10 11 12 1 5 r3 0.3249 0.202 0.331 0.2413 0.3249 0.202 0.331 0.4301 0.854 0.2413 0.3709 0.9768 2 0 e4 0.3725 0.331 0.4389 0.4232 0.3725 0.331 0.4389 0.9872 1.0368 0.4232 0.9807 0.9783 3 0 c2 0.14651 0.63861 0.32598 0.75383 0.14651 0.63861 0.32598 1.29528 0.21405 0.75383 1.64411 0.00154 h) g) i3 1 e4 i3 r3 1 a ) b ) c ) d ) e ) f ) g ) h ) dL d = 1 dE d 48 length as a bench mark. If the distance between two vectors is larger than the average vector length, we will consider the distance to be large. The average vector length is defined as . (5.7) If we let be the nearest neighbor of , we can check if is a FNN by measuring the distance between and as time increases ( ). For example, we can see from the row of Table 5.7a that the nearest neighbor of is . If we look at Table 5.7b, we can see that the distance between and is (see the circled value). But, after one time step ahead, we can see that the distance between and is which is greater than as indicated by the arrow. We conclude that is a FNN since the distance between and has grown more than as time increases. A second example of the FNNs can be seen from the first row of Table 5.7a. We can see that the nearest neighbor of is , while in Table 5.7b, the distance between and is that means is a FNN. On the other hand, the nearest neighbor of is , as seen in the second row of the last column in Table 5.7a, while if we look at Table 5.7b, we can see that the distance between and is . If we move one step ahead, we can see that the distance between and is which is still less than . A conclusion that can be reached is that is a true nearest neighbor of . Notice here that we choose one as the number of forward time steps that is used to check for FNN. Table 5.9b summarizes the number of FNNs found for through 3. It shows that at , the number of FNNs found is 4, while at , the number of FNNs is 0 β 1 M  yd(k) k = 1 M = Σ yd(n) yd(m) yd(n) y3(m) y3(n) dE = 3 8th y1(8) y1(9) y3(8) y3(9) q3(8, 9) = 0.8638 < β = 1.3404 y3(9) y3(10) q3(9, 10) = 1.4332 β = 1.3404 y1(9) y3(8) y3(9) β y1(1) y1(8) y3(1) y3(8) 1.7115 > β y1(8) y1(2) y1(6) y3(2) y3(6) 0.202 < β y3(3) y3(7) 0.331 β y1(6) y1(2) d = 1 d = 1 d = 2 49 and remains 0 at . That means the minimum embedding dimension is . Finally, notice that the nearest neighbor of shown in the last row of the Tables 5.7b, 5.8b, and 5.9a are labeled as not decidable (nd). The reason for this is that at the first instance of time, the distance between the two vectors is less than , but we do not have enough data to check if the distance between the two vectors after one time step ahead is greater than . So, we have to discard this neighbor from the count of the FNNs and label it as not decidable (nd). Table 5.7 The CDT method for d = 1 d = 3 dL = 2 yd(13) β β The threshold is Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 NaN 0.65 1.2015 0.1561 0.313 0.6479 1.2886 0.0262 0.0508 0.3597 0.4 0.7439 1.6362 0.9341 0.1046 8 2 0.65 NaN 1.8515 0.4939 0.963 0.0021 1.9386 0.6762 0.7008 0.2903 1.05 0.0939 2.2862 1.5841 0.7546 6 3 1.2015 1.8515 NaN 1.3576 0.8885 1.8494 0.0871 1.1753 1.1507 1.5612 0.8015 1.9454 0.4347 0.2674 1.0969 7 4 0.1561 0.4939 1.3576 NaN 0.4691 0.4918 1.4446 0.1823 0.2069 0.2036 0.556 0.5878 1.7922 1.0901 0.2607 1 5 0.313 0.963 0.8885 0.4691 NaN 0.9609 0.9755 0.2868 0.2622 0.6727 0.0869 1.0569 1.3231 0.621 0.2084 11 6 0.6479 0.0021 1.8494 0.4918 0.9609 NaN 1.9364 0.6741 0.6987 0.2882 1.0478 0.096 2.284 1.5819 0.7525 2 7 1.2886 1.9386 0.0871 1.4446 0.9755 1.9364 NaN 1.2624 1.2377 1.6483 0.8886 2.0325 0.3476 0.3545 1.1839 3 8 0.0262 0.6762 1.1753 0.1823 0.2868 0.6741 1.2624 NaN 0.0246 0.3859 0.3737 0.7701 1.61 0.9078 0.0784 9 9 0.0508 0.7008 1.1507 0.2069 0.2622 0.6987 1.2377 0.0246 NaN 0.4105 0.3491 0.7947 1.5853 0.8832 0.0538 8 10 0.3597 0.2903 1.5612 0.2036 0.6727 0.2882 1.6483 0.3859 0.4105 NaN 0.7596 0.3842 1.9959 1.2937 0.4643 4 11 0.4 1.05 0.8015 0.556 0.0869 1.0478 0.8886 0.3737 0.3491 0.7596 NaN 1.1438 1.2362 0.5341 0.2953 5 12 0.7439 0.0939 1.9454 0.5878 1.0569 0.096 2.0325 0.7701 0.7947 0.3842 1.1438 NaN 2.3801 1.6779 0.8485 2 13 1.6362 2.2862 0.4347 1.7922 1.3231 2.284 0.3476 1.61 1.5853 1.9959 1.2362 2.3801 NaN 0.7021 1.5316 7 14 0.9341 1.5841 0.2674 1.0901 0.621 1.5819 0.3545 0.9078 0.8832 1.2937 0.5341 1.6779 0.7021 NaN 0.8294 3 15 0.1046 0.7546 1.0969 0.2607 0.2084 0.7525 1.1839 0.0784 0.0538 0.4643 0.2953 0.8485 1.5316 0.8294 NaN 9 a) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 8 FNN 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 1 FNN 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 9 FNN 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 8 FNN 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 2 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd b) At d =1, the number of FNN = 4 nd: not decidable β = 1.3404 i1 1 i1 1 a ) b ) 50 Table 5.8 The CDT method for d = 2 Table 5.9 The CDT method for d = 3 The circled distances in the above tables are the distances between the reference points and their nearest neighbors at the first instance of time in . The arrows represent the direction where the reference points and their nearest neighbors move after one time step. Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 NaN 1.9623 1.2991 0.9756 0.313 2.044 1.4552 0.7013 0.2947 1.1099 0.4108 2.4041 2.2773 1.2008 9 2 1.9623 NaN 2.2959 1.0165 2.0851 0.0871 2.267 1.3346 1.7113 0.8525 2.2106 0.4447 2.3018 1.9268 6 3 1.2991 2.2959 NaN 1.4363 1.0155 2.3467 0.202 1.1934 1.1685 1.6572 0.994 2.6451 1.1736 0.3735 7 4 0.9756 1.0165 1.4363 NaN 1.0693 1.0925 1.4728 0.3193 0.7038 0.2214 1.1942 1.4478 1.8968 1.1099 10 5 0.313 2.0851 1.0155 1.0693 NaN 2.1617 1.1858 0.7553 0.3896 1.2452 0.1295 2.5167 2.0623 0.9757 11 6 2.044 0.0871 2.3467 1.0925 2.1617 NaN 2.3116 1.4094 1.7902 0.9342 2.2867 0.3606 2.3114 1.9759 2 7 1.4552 2.267 0.202 1.4728 1.1858 2.3116 NaN 1.2626 1.2965 1.6901 1.1759 2.5928 0.9721 0.3631 3 8 0.7013 1.3346 1.1934 0.3193 0.7553 1.4094 1.2626 NaN 0.4113 0.5204 0.8782 1.7625 1.8363 0.9094 4 9 0.2947 1.7113 1.1685 0.7038 0.3896 1.7902 1.2965 0.4113 NaN 0.8635 0.5191 2.1483 2.0462 0.9978 1 10 1.1099 0.8525 1.6572 0.2214 1.2452 0.9342 1.6901 0.5204 0.8635 NaN 1.3731 1.2945 2.0661 1.327 4 11 0.4108 2.2106 0.994 1.1942 0.1295 2.2867 1.1759 0.8782 0.5191 1.3731 NaN 2.6407 2.0842 1.0026 5 12 2.4041 0.4447 2.6451 1.4478 2.5167 0.3606 2.5928 1.7625 2.1483 1.2945 2.6407 NaN 2.4815 2.2718 6 13 2.2773 2.3018 1.1736 1.8968 2.0623 2.3114 0.9721 1.8363 2.0462 2.0661 2.0842 2.4815 NaN 1.0867 7 14 1.2008 1.9268 0.3735 1.1099 0.9757 1.9759 0.3631 0.9094 0.9978 1.327 1.0026 2.2718 1.0867 NaN 7 a) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 9 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd b) At d =2, the number of FNN is 0 i2 1 i2 1 a ) b ) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd a) At d =3, the number of FNNs is 0 also. d FNNs 1 4 2 0 3 0 b) i3 1 a ) b ) ℜ3 51 5.3 The Predictive technique Now that we have discussed the three geometric methods (CND, CDD, and CDT), we will now talk about the predictive technique. We said in Chapter 4 that the predictive technique can also be used to estimate by approximating the function in Equation (4.6). We also said that we will approximate by using a neural network. Since chaotic systems are nonlinear, we need to use a nonlinear network to make the approximation. To estimate of the Henon map by using the predictive technique, we will use a nonlinear neural network which consists of two layers. The first layer has 2 neurons with sigmoid transfer functions and a Tapped Delay Line (TDL) connected to it. The second layer has one neuron with a linear transfer function. The TDL is fed with the measurements to produce the delayvectors . The network structure is shown in Figure 5.1. While the number of taps in the TDL (d) is changed from 1 to 2 to 3, the network is trained to predict the current measurement . After the end of the training process, the sum of the squares of the prediction errors (SSE) is recorded as a function of d. We applied this technique to 100 points from the coordinate of the Henon map. As we can see from Table 5.10, the SSE changes significantly as d increases from 1 to 2 and remains almost the same at . We conclude that the network has accurately approximated and that . The logplot of the SSE versus d is shown in Figure 5.2. dL μ μ dL y(m) yd(m – 1) y(m) X1 d = 3 μ dL = 2 52 Figure 5.1 The nonlinear network used to estimate of the Henon map Table 5.10 The neural network SSE as a function of d Figure 5.2 LogPlot of the SSE versus d for the Henon map. 5.4 Chapter summary In this chapter, we demonstrated the estimation of the minimum embedding dimension of the Henon map by using the geometric and the predictive techniques. The purpose of this chapter was to provide some insight into the operation of the algorithms. In the next chapter, we will present the complete algorithms in full detail and discuss practical issues in using the algorithms on more complex systems than the Henon map. dx1 1 1x1 y(m) n1 (m) a1 (m) 1 b1 2x1 2xd b2 LW2, 1 n2 (m) TDL tansig purelin IW1, 1 1x2 a2 (m) dL d SSE(d ) 1 4.0321 2 4.7117x107 3 4.0824x107 1 2 3 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 SSE of the prediction error 53 CHAPTER 6 ADVANCED ALGORITHMS FOR ESTIMATING THE MINIMUM EMBEDDING DIMENSION 6.1 Introduction In Chapter 3, we explained that two parameters are needed in order to apply the embedding theorem to model a chaotic system. These two parameters are the dimension (d) of the delayvector ( ) and the delaytime (T) between the delayvector coordinates. We presented in Chapter 4 three geometric methods (CND, CDD, and CDT) and one predictive method through which we estimate the minimum dimension ( ) required to embed the system’s attractor. More details were given in Chapter 5 regarding the application of the four methods to estimate of the Henon map. The purpose of this chapter is to give full detail of six algorithms which use the four methods mentioned above to estimate . We will also show a method used to find the delaytime T. Before that, we first give a summary of the four methods. The CND, CDD, and the CDT geometric methods depend on the idea that if the dimension of the space is not large enough to represent the attractor of the system, projection artifacts will appear in the projected attractor. These artifacts cause points on the attractor to be falsely projected close to each other and produce False Nearest Neighbors (FNN) (see Section 4.3 for more detail). The dimension can be estimated as the minimum dimenyd dL dL dL dL 54 sion where the percentage of the FNNs does not change significantly with further increase in dimension. The CND method detects the existence of FNNs by checking to see if the nearest neighbors in the space of dimension d remain neighbors in dimension . On the other hand, the CDD method detects the existence of FNNs by checking to see if the distance between the nearest neighbors in dimension d will increase significantly as the dimension increases to . For the case of the CDT method, detection of the existence of FNNs is done by checking to see if the distance between the nearest neighbors in dimension d will change significantly as time increases. On the other hand, in the predictive method, the estimation of is done by approximating the function that operates on the reconstructed attractor. is approximated by using a neural network with a Tapped Delay Line (TDL) connected to its input. As the number of taps in the TDL (d) increases, the prediction error decreases. At one point, further increase of d does not improve the prediction error. At this point, is found. In the remaining parts of this chapter, we present in Section 6.2 a method used to find the delaytime T. In Section 6.3, we present two different algorithms based on the CDD and the CDT methods which were proposed by Abarbanel et al to estimated . In Section 6.4, we discuss some limitations of the previous two algorithms and suggest four new algorithms to overcome these limitations. In Section 6.5, we present three of the four algorithms, which are based on the geometric technique. The fourth algorithm, which is based on the predictive techniques, is presented in Section 6.6. At the end, Section 6.7 concludes with a chapter summary. d + 1 d + 1 dL μ: ℜd ℜ1→ μ dL dL 55 6.2 The delaytime (T) We explained in Section 4.2 that we wanted to find the values of the parameters d and T such that the system we developed
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Modeling Chaotic Systems 
Date  20050701 
Author  AlMughadhawi, Khaled Marzoug 
Department  Electrical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  The purpose of this research is to model chaotic systems. From a set of scalar measurements taken from a chaotic system, we want to describe the hidden states of the system. The modeling process has two parts: determining the model order, and determining the model parameters. The major characteristic of chaotic systems is their sensitivity to changes in initial conditions. This characteristic limits the ability to make a long term prediction in these systems. To be able to make an accurate model, we need to accurately estimates of the model parameters, which requires an accurate estimate of the model order. Findings and conclusions./ In this research we found four new algorithms to estimate the model order. These algorithms were tested on different chaotic systems. Three of these algorithms are geometric and one is predictive. We checked on the ability of these algorithms to make accurate estimates of the model order. We also proved a new theorem that relates the Lyapunov exponents (model parameters) to the linear system poles. We introduced a new algorithm that implement the result of this theorem to estimate the model parameters of a chaotic model. Beside this algorithm, we improved an existing algorithm that uses a neural network to estimate the model parameters. Both algorithms were tested on different chaotic systems. Our finding with respect to the algorithms used to estimate the model order can be summarized in the following: (i)�the predictive algorithm gave the best estimate of model order as long as the signal to noise ratio is not too low, (ii)�the global neighbor search algorithms proposed in this research gave better estimate of model order than the local neighbor search algorithms, (iii)�the use of more than one algorithm to estimate model order increases our confidence in the estimated value. With respect to the algorithms used to estimate the model parameters, the main conclusions are: (i)�the predictive algorithm gives good estimate of the model parameters for signals with high signal to noise ratio, (ii)�the algorithm that implement the new theorem result gives good estimate of the model parameters in some cases, and (iii)�when using more than one algorithm to estimate the model parameters, one can have confidence in the estimate values if they are in agreement with each other. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  MODELING CHAOTIC SYSTEMS By Khaled Marzoug AlMughadhawi Bachelor of Science in Electrical Engineering King Fahad University of Petroleum and Minerals Dhahran, Saudi Arabia 1989 Master of Science in Electrical Engineering Southern Methodist University Dallas, Texas 1999 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY In Electrical Engineering July, 2005 ii COPYRIGHT By Khaled Marzoug AlMughadhawi July, 2005 iii MODELING CHAOTIC SYSTEMS Thesis Approved: Dr. Martin T. Hagan Thesis Adviser Dr. Gary E. Young Dr. Rao. K. Yarlagadda Dr. Carl D. Latino Dr. A. Gordon Emslie Dean of the Graduate College iv ACKNOWLEDGMENTS First and foremost, I would like to thank Allah (the God) for helping me and making it possible for me to reach this point of knowledge. Peace and blessing of Allah are on his messenger Muhammad, his last prophet and the seal of all messengers. My thanks extend to my parents: Marzoug, and Mohelah for their support, patience, and extreme love. My father grew up as an orphan in a small town known as Abar Ali (south of Madina, Saudi Arabia). When I was born, in 1966, he had to stop his education in the Madina teaching institute. He did that in order to support his new born child. He use to say to me “Khaled, I want you to compensate my loss of education by bringing me the highest degree in your major”. “I hope I did that Dad.” I would like to thank my advisor Dr. Martin Hagan for his support, patience, and help. He use to arrange weakly seminars for me, often lasting for three hours. He listened very carefully to my presentations and corrected me, suggesting new ways to solve the research problems. My thanks extend to all his graduate students who attended these seminars. I would like also to thank Dr. Yong Hu, and Dr. Tagi Menhaj for their valuable suggestions and comments. My thanks extend also to the dissertation committee members for their comments and suggestions. v Three professors in the Mathematics department gave me a much needed help and support. Dr. Saleh Mehdi gave me three lectures in differential geometry which helped me in understanding the embedding theorem. He and I corresponded often regarding the mathematical challenges I faced. Dr. Roger Zierau was very helpful to me. We spent long hours for many weeks discussing the linear Oseledec theorem. His support was very valuable to me. Dr. Leticia Barchini’s suggestions and comments made the proof of the theorem possible as well. My friends in the Islamic society of Stillwater made the difficult task of this degree much easier by their friendly relationship with me. I ask Allah to help them and make their future project of building a new mosque a success. My friends in the Saudi house showed me the true meaning of brotherhood. Our weakly meeting proved to be very important in bring us home far away from home. I would like also to thank my wife Nawal for her extreme support and patience with me and our five kids. Without her help, my goal could have been much more difficult to accomplish. My kids AbduRahman, Maha, Ahmad, Roba, and Arwa gave me the social stability which was also necessary to achieve this task. My good neighbor Sonya Brewster helped to proofread many chapters of my dissertation. Finally I would like to thank the Saudi government for their help and financial support of my research. vi TABLE OF CONTENTS Chapter Page 1. OBJECTIVE AND CONTRIBUTIONS ........................................................................1 1.1 Introduction .............................................................................................................. 1 1.2 Objective................................................................................................................... 1 1.3 Contributions ............................................................................................................ 1 1.4 Literature review ...................................................................................................... 2 1.5 Remaining chapters summary .................................................................................. 5 2. INTRODUCTION TO CHAOTIC SYSTEMS ..............................................................6 2.1 History of dynamical systems .................................................................................. 6 2.2 Dynamical definitions .............................................................................................. 7 2.2.1 Dynamical System ............................................................................................7 2.2.2 Equilibrium point ..............................................................................................8 2.2.3 Orbit of a dynamical system .............................................................................8 2.2.4 Attractor ............................................................................................................8 2.3 Characteristics of Chaotic Systems .......................................................................... 8 2.3.1 Sensitivity of chaotic systems to initial conditions ..........................................8 2.3.2 Chaotic signals look random but they are deterministic ...................................9 2.3.3 Chaotic systems have attractors with fractal dimension .................................11 2.3.3.1 The boxcounting dimension: .................................................................12 2.4 Examples of Chaotic systems................................................................................. 14 2.4.1 Weather system ...............................................................................................14 2.4.2 Biological models ...........................................................................................14 2.4.3 Laser signals ...................................................................................................15 2.4.4 Chaotic circuits ...............................................................................................15 2.4.5 Discrete chaotic systems: ................................................................................15 2.5 Chapter Summary................................................................................................... 15 vii 3. INTRODUCTION TO DYNAMICAL MODELING ..................................................17 3.1 Introduction ............................................................................................................ 17 3.2 Modeling................................................................................................................. 17 3.3 Applications of modeling ....................................................................................... 18 3.3.1 Detection of chaos ..........................................................................................18 3.3.2 Prediction ........................................................................................................18 3.3.3 Diagnosis ........................................................................................................19 3.4 Definitions .............................................................................................................. 19 3.4.1 Injection (1to1) function ..............................................................................19 3.4.1.1 Example: An injection ............................................................................19 3.4.1.2 Example: A function that is not injection ...............................................20 3.4.2 Immersion function..........................................................................................21 3.4.2.1 Example: An immersion .........................................................................21 3.4.2.2 Example: A function that is not immersion ............................................21 3.4.3 Embedding function ........................................................................................22 3.4.3.1 Example: An embedding ........................................................................22 3.4.3.2 Example: A nonembedding ...................................................................23 3.5 Modeling by Embedding ........................................................................................ 23 3.5.1 The delaycoordinate map ..............................................................................23 3.5.2 Example: The delaycoordinate map ..............................................................24 3.5.3 Example: Sufficient but not necessary condition for embedding ...................25 3.6 Chapter summary.................................................................................................... 26 4. INTRODUCTION TO THE MINIMUM EMBEDDING DIMENSION ESTIMATION 28 4.1 Introduction ............................................................................................................ 28 4.2 Modeling chaotic systems ...................................................................................... 28 4.3 Estimating the minimum embedding dimension.................................................... 30 4.3.1 The geometric technique .................................................................................30 4.3.1.1 The change of neighbors with dimension method (CND) ......................31 4.3.1.2 The change of distance with dimension (CDD) method .........................31 4.3.1.3 The change of distance with time method (CDT) ...................................32 4.3.2 The predictive technique: ...............................................................................34 4.4 Dynamic equivalence: ............................................................................................ 36 4.5 Chapter summary: .................................................................................................. 37 5. EXAMPLES OF THE MINIMUM EMBEDDING DIMENSION ESTIMATION ....38 5.1 Introduction ............................................................................................................ 38 5.2 The geometric technique ........................................................................................ 38 5.2.1 Using the change of neighbors with dimension method (CND) .....................39 5.2.2 Using the change of distance with dimension method (CDD) .......................43 viii 5.2.3 Using the change of distance with time method (CDT) .................................47 5.3 The Predictive technique ........................................................................................ 51 5.4 Chapter summary.................................................................................................... 52 6. ADVANCED ALGORITHMS FOR ESTIMATING THE MINIMUM EMBEDDING DIMENSION...............................................................................................................53 6.1 Introduction ............................................................................................................ 53 6.2 The delaytime (T).................................................................................................. 55 6.3 Two algorithms that use the local neighbor search ................................................ 57 6.3.1 The Algorithm ....................................................................................57 6.3.1.1 Example: Neighbor indices .....................................................................58 6.3.1.2 The vectorcoordinates projection method .............................................59 6.3.1.3 Example: Vectorcoordinates projection ................................................59 6.3.2 The Algorithm ....................................................................................61 6.3.2.1 The PCA projection method ...................................................................62 6.4 Limitations of the and the Algorithms ........................................... 66 6.5 Three new algorithms that use the global neighbor search method ....................... 68 6.5.1 The first new algorithm: The CND .................................................................68 6.5.2 The second new algorithm: The .........................................................71 6.5.3 The third new algorithm: The .............................................................74 6.6 The fourth new algorithm: The Predictive ............................................................. 77 6.7 Chapter summary.................................................................................................... 80 7. MINIMUM EMBEDDING DIMENSION RESULTS.................................................81 7.1 Introduction: ........................................................................................................... 81 7.2 Testing systems ...................................................................................................... 82 7.2.1 Noise free chaotic systems ..............................................................................82 7.2.2 Practical systems .............................................................................................86 7.3 Estimating for the noise free chaotic systems .................................................. 87 7.3.1 Using the algorithm ............................................................................87 7.3.2 Using the algorithm ............................................................................88 7.3.3 Using the CND algorithm ...............................................................................89 7.3.4 Using the algorithm ............................................................................90 7.3.5 Using the algorithm ............................................................................91 7.3.6 Using the predictive algorithm .......................................................................92 7.4 Estimating for the noisy systems ...................................................................... 94 7.4.1 Estimating for the noisy chaotic circuit (cc) .............................................94 CDDL CDTL CDDL CDTL CDDG CDTG dL CDDL CDTL CDDG CDTG dL dL ix 7.4.1.1 Using the algorithm ....................................................................94 7.4.1.2 Using the algorithm. ....................................................................95 7.4.1.3 Using the CND algorithm .......................................................................96 7.4.1.4 Using the algorithm ....................................................................96 7.4.1.5 Using the algorithm ....................................................................97 7.4.1.6 Using the predictive algorithm ...............................................................98 7.4.2 Estimating for the Santa Fe data sets ........................................................99 7.4.2.1 Using the algorithm ....................................................................99 7.4.2.2 Using the algorithm ...................................................................100 7.4.2.3 Using the CND algorithm .....................................................................101 7.4.2.4 Using the algorithm ..................................................................102 7.4.2.5 Using the algorithm ..................................................................103 7.4.2.6 Using the predictive algorithm .............................................................104 7.5 Testing the algorithms with random signals......................................................... 105 7.6 Tables and discussions ......................................................................................... 107 7.6.1 Tables ............................................................................................................107 7.6.2 Discussion of the Results ..............................................................................108 7.6.2.1 The effect of the original dimension of the system ..............................108 7.6.2.2 Local versus global neighbor search methods ......................................108 7.6.2.3 Estimation time .....................................................................................109 7.6.2.4 Sensitivity to the threshold value ..........................................................109 7.6.2.5 Noise detection .....................................................................................110 7.6.2.6 Dependence of the algorithms on the number of data points ...............110 7.6.3 General conclusions ......................................................................................110 7.7 Chapter Summary................................................................................................. 110 8. THEORY OF LYAPUNOV EXPONENTS...............................................................112 8.1 Introduction .......................................................................................................... 112 8.2 Sensitivity of some linear systems to initial conditions ....................................... 113 8.3 Lyapunov Exponent (LE) of a first order chaotic system .................................... 115 8.4 Lyapunov exponents for a multidimensional system........................................... 117 8.5 Invariant sets in modeling by embedding............................................................. 119 8.6 LEs from the Jacobian matrix............................................................................... 120 8.6.1 Oseledec theorem ..........................................................................................121 8.7 Linear Algebra definitions.................................................................................... 122 8.7.1 Definition: Inner product ..............................................................................122 8.7.2 Definition: Vector Norm ...............................................................................122 8.7.3 Definition: Matrix Norm ...............................................................................123 8.7.3.1 Some important properties of the matrix norm .....................................123 8.7.3.2 Lemma: Norm of a diagonal matrix .....................................................124 CDDL CDTL CDDG CDTG dL CDDL CDTL CDDG CDTG x 8.7.4 Singular value decomposition .......................................................................125 8.7.4.1 Important SVD properties .....................................................................125 8.7.5 Definition: Diagonalizable matrix ................................................................125 8.8 Multilinear algebra definitions ............................................................................. 126 8.8.1 Vector exterior products ...............................................................................126 8.8.2 Linear operator exterior power .....................................................................127 8.8.2.1 Definition: Adjoint of a linear operator ................................................128 8.8.2.2 Lemma: Adjoint of the wedge ..............................................................128 8.8.2.3 Linear operator properties .....................................................................129 8.8.2.4 Lemma: Eigen values of .........................................................130 8.9 The linear Oseledec matrix................................................................................... 131 8.9.1 Theorem: Eigen values of a linear Oseledec matrix .....................................131 8.9.2 For a symmetric matrix .................................................................................133 8.9.3 For a general matrix ......................................................................................134 8.9.3.1 Lemma: Eigen values of from singular values of .................134 8.9.3.2 Proposition: Limit of the root of ...................................135 8.10 Chapter summary................................................................................................ 141 9. ESTIMATING THE SET OF LYAPUNOV EXPONENTS......................................142 9.1 Introduction .......................................................................................................... 142 9.2 Estimating the LEs from Jacobian matrices ......................................................... 143 9.2.1 Estimating the LE by the QR decomposition ...............................................143 9.2.1.1 Example: Estimating LEs of the Henon map ........................................145 9.2.2 Spurious exponents .......................................................................................146 9.3 Estimating the LEs by the Geometric algorithms................................................. 147 9.3.1 Eckmann algorithm .......................................................................................147 9.3.2 Linear Oseledec algorithm ............................................................................149 9.4 Estimating the LEs by the predictive algorithm ................................................... 150 9.5 Pseudo codes of the LEs estimation algorithms ................................................... 154 9.5.1 Pseudo code of the Eckmann algorithm .......................................................155 9.5.2 Pseudo code of the linear Oseledec algorithm ..............................................157 9.5.3 Pseudo code of the predictive algorithm ......................................................158 9.6 Results of estimating the LEs by using the three algorithms ............................... 159 9.6.1 Discussion of the results ...............................................................................160 9.7 Chapter summary.................................................................................................. 161 10. SUMMARY, CONCLUSIONS, AND FUTURE RECOMMENDATIONS...........162 10.1 Summary............................................................................................................. 162 10.2 Conclusions ........................................................................................................ 164 10.3 Future recommendations .................................................................................... 165 (L*L)^q Ok (A)k kth σd (A)k ( ) xi LIST OF TABLES Table Page 5.1 The CND method (1/3) ...........................................................................................42 5.2 The CND method (2/3) ...........................................................................................42 5.3 The CND method (3/3) ...........................................................................................43 5.4 The CDD method tables for d = 1 ...........................................................................45 5.5 The CDD method tables for d = 2 ...........................................................................46 5.6 The CDD method tables for d = 3 ...........................................................................47 5.7 The CDT method for d = 1 ......................................................................................49 5.8 The CDT method for d = 2 ......................................................................................50 5.9 The CDT method for d = 3 ......................................................................................50 5.10 The neural network SSE as a function of d ............................................................52 7.1 Tabulation of the estimated for all the testing systems shown in Sections 7.3 through 7.5. The wrong estimates are circled. The abbreviation cc: chaotic circuit, and N. Able: not able ............................................................................................107 7.2 Six algorithms estimation time (in seconds) for of data set A of the Santa Fe competition, it shows the is the fastest and the predictive is the slowest. ......107 8.1 Estimating for . ..................................115 9.1 LEs estimation results using the three algorithms. The numbers inside the curly parenthesis are the estimated LEs and the number in bottom is the Lyapunov dimension ( ). ....................................................................................................................160 dL dL CDDG log(c) f(x(k)) = (3 ⁄ 4)(1 – 1 – 2x(k) ) dLyp xii LIST OF FIGURES Figure Page 2.1 Sensitivity of the chaotic tent map to initial conditions............................................ 9 2.2 The 60 Hz periodic wave and its frequency response............................................. 10 2.3 The random signal and its frequency response....................................................... 10 2.4 The chaotic signal and its frequency response........................................................ 11 2.5 The chaotic attractor of the Henon map.................................................................. 12 2.6 The slope of the line is the boxcounting dimension of the Henon map attractor. . 14 3.1 a) The attractor C in , b) the set of measurements, c) The reconstructed attractor by the delaycoordinate map................................................................................... 25 3.2 The circle in and its reconstruction in . ................................................ 26 4.1 Projection artifacts .................................................................................................. 30 4.2 The CDD method.................................................................................................... 32 4.3 The projection of a 3D curve (solid line) into 2D curve (dotted line) ................. 33 4.4 The CDT method: a) Distances between false neighbors increase with time, b) The intersection in the middle of the curve is a projection artifact................................ 33 4.5 The function approximation in the predictive technique.................................... 36 4.6 A neural network model with a TDL used to approximate the function , D is a one step delaytime........................................................................................................ 36 5.1 The nonlinear network used to estimate of the Henon map.............................. 52 5.2 LogPlot of the SSE versus d for the Henon map................................................... 52 6.1 The Lorenz model average mutual information ..................................................... 57 6.2 Pseudocode of the algorithm ...................................................................... 61 6.3 Pseudocode of the algorithm....................................................................... 65 6.4 Pseudocode of the CND algorithm ......................................................................... 70 6.5 Pseudocode of the algorithm...................................................................... 73 6.6 Pseudocode of the algorithm ...................................................................... 76 6.7 Pseudocode of the predictive algorithm.................................................................. 79 7.1 Lorenz chaotic attractor .......................................................................................... 83 7.2 The chaotic circuit diagram .................................................................................... 83 7.3 The chaotic circuit response ................................................................................... 84 7.4 The Rossler model attractor .................................................................................... 85 7.5 The estimated for the six noise free systems using the algorithm. The vertical axis is the percentage of the FNNs found from each dimension and the horizon ℜ3 S1 ℜ3 ℜ2 μ μ dL CDDL CDTL CDDG CDTG dL CDDL xiii tal axis is the dimension d, a) for the Lorenz model, , b) for the chaotic circuit, , c) for the Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, .. 88 7.6 The estimated for the six noise free systems using the algorithm. The vertical axis is the percentage of the FNNs and the horizontal axis is the dimension d, a) for Lorenz model, , b) for chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ......................................................................... 89 7.7 Estimating for the noise free systems using our first algorithm (CND). The vertical axis is the estimated and the horizontal axis is the number of neighbors used (w). a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, this algorithm does not give a stable answer.............. 90 7.8 Estimating for the noise free systems using our second algorithm ( ). The vertical axis is the percentage of the FNNs found and the horizontal axis is the dimension d. a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ................................................... 91 7.9 Estimating for the noise free systems using our third algorithm ( ). The vertical axis is the percentage of FNNs and the horizontal axis is the dimension d. a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, . ......................................................................... 92 7.10 The predictive algorithm results for the noise free systems, the vertical axis is the SSE of the prediction errors and the horizontal axis is the dimension d, a) for Lorenz model, , b) for the chaotic circuit, , c) for Rossler model, , d) for MG of dimension 4, , e) for MG of dimension 7, , f) for MG of dimension 13, ............................................................................................. 93 dL = 3 dL = 3 dL = 2 dL = 4 dL = 4 dL = 4 dL CDTL dL = 2 dL = 3 dL = 2 dL = 3 dL = 3 dL = 4 dL dL dL = 3 dL = 3 dL = 3 dL = 4 dL = 5 dL CDDG dL = 3 dL = 3 dL = 3 dL = 4 dL = 4 dL = 4 dL CDTG dL = 3 dL = 3 dL = 3 dL = 4 dL = 7 dL = 7 dL = 3 dL = 3 dL = 3 dL = 4 dL = 7 dL = 13 xiv 7.11 Estimating for the noisy cc. using Abarbanel et al first algorithm ( ). a) at 200 dB, , b) at 100 dB, , c) at 50 dB, , d) at 20 dB, the algorithm assumes the signal is noise............................................................................ 95 7.12 Estimating for the noisy cc. using algorithm. a) at 200 dB, , b) 100 dB, , c) at 50 dB, and d) at 20 dB the algorithm assumes the signal is noise . 95 7.13 Estimating for the noisy cc. using our first algorithm (the CND), a) at 200 dB, , b) at 100 dB, , c) at 50 dB, , d) at 20 dB, ......... 96 7.14 Estimating for the noisy cc. using our second algorithm ( ). It shows for SNR=200, 100, 50, and 20 dB ............................................................ 97 7.15 Estimated for the noisy cc. using the algorithm. For SNR=200, 100, and 50, . At SNR=20db, not stable................................................................... 98 7.16 Estimating using our fourth algorithm (Predictive), for 200, 100, and 50 dB, at 20 dB it assumes a random signal ................................................................ 98 7.17 Using the algorithm, a) for A data set, , b) for , , c) for , .................................................................................................................. 99 7.18 Using algorithm, a) for data set A, , b) for , , c) for , . .............................................................................................................. 100 7.19 Using the CND algorithm, a) for data set A, , b) for , , c) for , the result is not stable............................................................................................ 101 7.20 Using algorithm, a) for data set A, , b) for , , c) for , .............................................................................................................. 102 7.21 Using algorithm, a) for data set A, , b) for , , c) for , ................................................................................................................. 103 7.22 Using the predictive algorithm, a) for data set A, , b) for , , c) for , ....................................................................................................... 104 7.23 Testing Abarbanel et al two algorithms with a random signal, a) the algorithm fails to recognize the random signal, b) the algorithm succeeded in recognizing the random signal............................................................................................ 105 dL CDDL dL = 3 dL = 3 dL > 5 dL CDTL dL = 2 dL = 2 dL dL = 3 dL = 3 dL = 3 dL = 4 dL CDDG dL = 3 dL CDTG dL = 3 dL dL = 3 CDDL dL = 4 B1 dL = 4 D1 dL = 6 CDTL dL = 3 B1 dL = 6 D1 dL = 10 dL = 3 B1 dL = 5 D1 CDDG dL = 3 B1 dL = 5 D1 dL = 11 CDTG dL = 3 B1 dL = 5 D1 dL = 6 dL = 3 B1 dL = 4 D1 dL = 10 CDDL CDTL xv 7.24 Testing our four algorithms for estimating the dimension of the random signal, a) the CND algorithm recognizes the signal is random by not changing the estimated as w increases, b) for the algorithm, it did not recognize the random signal, c) for the algorithm, it was able to recognize the random signal, the same for the predictive algorithm in d). .................................................................................... 106 8.1 Linear systems can be sensitive to initial conditions ............................................ 113 8.2 Vectors a and b exterior product ( ) ............................................................... 126 9.1 The feed forward network used to estimate the LEs............................................. 152 9.2 The Eckmann algorithm Pseudo code .................................................................. 156 9.3 The linear Oseledec algorithm pseudo code ......................................................... 157 9.4 The predictive algorithm pseudo code.................................................................. 158 dL CDDG CDTG a^b xvi NOMENCLATURE Basic Concepts Scalars: small italic letters...a, b, Vectors: small bold nonitalic letters, or underbar symbol a, b, Matrices: capital Bold nonitalic letters ... A, B, C Language Vector means a column of numbers. Vector element Small italic with an index Matrix elements An element of matrix A: small italic showing the row then the column of the matrix: , or The column vector of a matrix I The row vector of a matrix I φ φ a(i) or ai a(i, j), aij A(i, j) kth ik kth i k xvii Norm Covariance matrix C Eigen value and eigen vector and (or bold small nonitalic letter) Set Empty set Subset of A, B, C Dimension Space dimension is superscript Box counting dimension Theoretical minimum embedding dimension Actual minimum embedding dimension Side length of a dcube (boxcounting dimension) a λ η ∅ ℜn ℜd dc dE dL σ xviii Interval in (boxcounting dimension) N Number of dcubes needed to cover a set (boxcounting dimension) Lyapunov dimension (Chapter 8) The number of Lyapunov exponents such that Power For a vector For a matrix Dynamics State vector in the original space Equilibrium point ℜ1 σ cσ d dLyp λj j = 1, 2, …, kL Σ > 0 kL (x)k (A)k x(n) x* xix The map for the state update in the original space Scalar Measurements Measurement function Delaytime T Delaycoordinate map (for embedding) Delayvector in the reconstructed attractor in The map for the state update in the reconstructed space Map for the output update in the reconstructed space The sampled measurement (predictive algorithm) , Equation (6.16) The sampled measurement with zero mean x(k + 1) = f(x(k)) y(m) y(m) = h(x(m)) yd(m) = Fd(h, x(n)) ℜd yd(m) y(m) y(m – T) … y(m – (d – 1)T) t = yd(m + 1) = φ(yd(m)) y(m + 1) = μ(yd(m)) ys(m) = y(1 + (m – 1)T) s(m) ys(m) 1 L  ys(k) k = 1 L = – Σ xx The delayvector created from Vectors distances matrix in the ddimensional space , see section 5.2.1 on page 39 Number of neighbors in used to search for the nearest neighbor of (for the CND method) w Indices matrix of the wnearest neighbors in the ddimensional space Nearest neighbor distances vector in ddimensional space Nearest neighbor distances in ddimensional space with distances computed in (d + 1)dimensional space Number of neighbors of Basis matrix for the PCA projection method ; are the eigen vectors of the covariance matrix C s(m) yd s (m) s(m) s(m – 1) … s(m – d + 1) t = Qd ℜd + 1 yd(m) Id rd ed + 1 yd(m) Nb Bd η1 η2 … ηd = ηi xxi The neighbors of the delayvector in : when the delaycoordinate projection method is used. Its element is : when the PCA projection method is used. It element is (see page 62) Index of the neighbor of Average vector in the space , the CDT pseudocode Average mutual information between two measurements I Threshold for significant FNN change (for the CND method) Threshold of distance increase (for the CDD method) Attractor mean value (for the CDT method) A fraction of attractor mean (for the CDT method) Nb yd m ( ) ℜd Yd Nb yd(m) Yd Nb, p yd p (m) kth yd(m) id k (m) ℜdE ydE (m) α ρ β ζ xxii Time step for the CDT method Distance between the nearest neighbors after time steps (CDT) Threshold for significant change in prediction error (Predictive technique) Signal variance Noise variance Jacobian matrix at or J Sampling time Delay time in MackeyGlass equation (Chapter 7) Lyapunov Exponent (global) Finite time Lyapunov Exponent (local) Δ Δ σΔ γ (σs)2 (σn)2 yd(m) Jyd(m) τs tf λ λt xxiii Infinitesimal distance Infinitesimal perturbation from to Abbreviations FNN: False Nearest Neighbor CND: Change of Neighbor with Dimension CDD: Change of Distance with Dimension CDT: Change of Distance with Time PCA: Principal Component Analysis SNR: Signal to Noise Ratio ECG: Electrocardiogram TDL: Tapped Delay Line SSE: Sum Squared Errors LE: Lyapunov Exponent SVD: Singular Value Decomposition Linear Algebra Singular value of a matrix or Multiplications of the n Jacobian matrices ε x(t) xε(t) δ(t) = xε(t) – x(t) σi σi(A) Jy n = Jy(n – 1)Jy(n – 2)…Jy(0) xxiv Perturbation vector from time m to time n Matrix of the perturbation vectors from Diagonal matrix D Diagonal matrix with elements magnitude of x ; Linear Oseledec matrix General Oseledec matrix Multilinear Algebra Vector exterior product The exterior power of a linear operator L Adjoint of a linear operator L Matrix representation of the wedge of the linear operator L ( ) δ(m) = x(n) – x(m) Nb yd(m) Em Da x ∈ Da Ok Λx ei^ej qth L^q L* L^q Bq 1 CHAPTER 1 OBJECTIVE AND CONTRIBUTIONS 1.1 Introduction In this chapter, we will show the objective of this research and give a brief list of our contributions. In Section 1.4, we give a literature review of some important publications in chaotic modeling. A summary of the remaining chapters is given in Section 1.5. 1.2 Objective The objective of this research is to model chaotic systems. We will build a chaotic model from a set of scalar measurements taken from the system. Evolution of the states in the model follows the evolution of the hidden states in the original system. To build a chaotic model, we start by estimating the model order. The next step is to estimate the set of Lyapunov exponents (the model parameters) that characterize the chaotic system. A good chaotic model depends on both the accuracy of the estimated model order and the estimate of Lyapunov exponent values. 1.3 Contributions Our contributions are: (1) Four new algorithms for estimating the model order. a) The Change of Neighbors with Dimension (CND). 2 b) The Change of Distance with Dimension using the Global neighbor search ( ). c) The Change of Distance with Time using the Global neighbor search ( ). d) The predictive. (2) The proof of a new theorem that relates the poles of a linear system to the set of Lyapunov exponents. (3) Implementation of the new theorem result to estimate the set of Lyapunov exponents for chaotic systems. (4) Development of an existing algorithm which uses a neural network to estimate the set of Lyapunov exponents. In total, we implemented four algorithms to estimate the model order and two algorithms to estimate the set of Lyapunov exponents for a chaotic system. The six algorithms were tested on different chaotic systems. The testing examples include noise free and noisy signals. In addition, the testing systems have different dimensions. 1.4 Literature review There are two primary areas of research in modeling chaotic systems. These areas, which will be discussed briefly, are: estimating the model order and estimating the set of Lyapunov exponents. In this section, we list some key references in these areas. We begin by listing literature dealing with estimating the model order. After that we list the literature that describes estimating the set of Lyapunov exponents. CDDG CDTG 3 We begin by reviewing publications that discuss estimating the model order. We found many papers by Abarbanel and his colleagues that talk about different methods for estimating the model order. Kennel, Brown, and Abarbanel [KBA92] introduced estimating the model order by the method of False Nearest Neighbors (FNNs). This method is a geometric one. It depends on measuring the distance between neighbors as the model order increases. The distance between false neighbors increases more than the distance between true neighbors as the model order increases. In 1993, Abarbanel and Kennel [AbKe93] developed the method of FNNs further and introduced a new method that uses time as well as distance to check for the existence of FNNs. In this proposal, we will present new algorithms that improve on the results of Abarbanel’s algorithms. In 1997, Cao [Cao97] suggested a new method for estimating the model order. His method is a geometric method and it is closely related to Abarbanel method [KBA92]. Instead of using a threshold value to determine false neighbors, he suggested the use of the mean of distance change as the model order increases. Abarbanel [Aba98] used a geometric method with a predictor to estimate the model order. Since the points inside the attractor have an evolution rule, their neighboring points will evolve according to a certain rule as well. To determine the model order, he used a predictor to estimate the evolution rule. When the model order is not large enough, prediction errors will increase. As the model order increases, the prediction errors drop. At certain point, further increase of the model order does not improve the prediction errors. At this point, the model order is found. In 2002, Kennel and Abarbanel [KeAb02] used the method of FNNs with a global neighbor search method. The projection method used in this paper is the delaycoordinate method. (Our second algorithm (CDDG ) used the method of 4 FNNs, with a global neighbor search, but the used the Principal Component Analysis projection method.) Now we review some papers that talk about estimating the set of Lyapunov exponents for a chaotic model. The most important paper in this field is the one by Oseledec [Ose68]. He proved that the set of Lyapunov exponents can be estimated from the product of an infinite number of Jacobian matrices along the attractor of the system. (In this proposal, we will present a new theorem that provides insight into the Oseledec theorem by relating Lyapunov exponents to the poles of a linear system.) In 1985, Eckmann and Ruelle [EcRu85] and Sano and Sawada [SnSd85] applied the Oseledec theorem to measurements taken from a chaotic system. Both papers use orthogonalization methods to overcome the numerical problem of multiplying a large number of matrices. Parlitz [Par92] introduced the identification of spurious Lyapunov exponents. He did that by reestimating the set of Lyapunov exponents from measurements backward in time. By doing that, he found that true exponents change their signs, while spurious ones change their values as well as signs. Darbyshire and Broomhead [DaBr96] estimated the set of Lyapunov exponents by the methods of least squares and total least squares. Their method uses the pseudoinverse to estimate the Jacobian matrices. After that, they applied the orthogonalization method proposed by Eckmann. Djamai and Coirault [DjCo02] introduced the use of neural networks to estimate the set of Jacobian matrices. After estimating the Jacobian matrices, they used the same method proposed by Eckmann to reorthogonalize the product of these matrices. CDDG 5 1.5 Remaining chapters summary In Chapter 2, we will present an introduction to chaotic dynamical systems. In Chapter 3, we introduce modeling chaotic systems and discuss modeling by embedding. Examples are given to illustrate the idea of embedding functions. In Chapter 4, we will present an introduction to four methods used to find the minimum dimension (model order) required for embedding. In Chapter 5, we implement the methods found in Chapter 4 to find the minimum embedding dimension of the Henon map. Six algorithms based on the four previous methods are presented in Chapter 6. (Four of the six algorithms represent our original work.) Full details of the algorithms are given, and more practical issues are discussed in Chapter 6. The results of implementing the six algorithms on nine chaotic systems are found in Chapter 7. The minimum embedding dimensions of the nine systems are listed in this chapter, and a comparison is made among the six algorithms. In Chapter 8, we will explore the theory of Lyapunov exponents (LEs) and prove a new theorem that relates the poles of a linear system to the set of LEs. In Chapter 9, we will discuss the estimation of the LE values and apply the result of the new theorem to estimate the LEs. In this chapter, we will also improve an existing algorithm that uses a neural network to estimate the LEs. We tested the two algorithms by estimating the LEs of different chaotic systems. We also compare the two algorithms to an existing algorithm using the test systems. In Chapter 10, we will give a conclusion of the dissertation and give some recommendations on possible future research. 6 CHAPTER 2 INTRODUCTION TO CHAOTIC SYSTEMS 2.1 History of dynamical systems Isaac Newton (16421727) introduced the idea of modeling the motion of objects by equations. Position, velocity, and acceleration were the fundamental parameters of his equations. From position, velocity, and acceleration, he could describe the state of a moving object at any given time. Later scientists modelled dynamical systems by using a set of differential equations. The solution of these equations describes the state of the system at any given time. If the solution of the set of differential equations remains in a bounded region, the sequence of states reduces to either i) a steady state, generally because of loss of energy by friction, or ii) a periodic or quasi periodic motion. An example of case (ii) is the motion of the moon around the earth and the motion of the earth around the sun. Case (i) and case (ii) remained the only two known bounded solutions until the use of modern computers made possible the numerical solution of sets of differential equations. In 1963 Edward Lorenz published a paper entitled “Deterministic non Periodic Flow” [Lo63]. He discovered a new bounded attractor that is not periodic but was filling a region in space. This was the first time the world knew about the third possible bounded attractor: iii) the chaotic attractor. The new solution (chaotic) is not simply periodic nor 7 quasi periodic with a large number of periods. Chaotic motion is possible with simple one dimensional systems. Although chaotic motion can be produced from simple systems, it is very complicated and becomes unpredictable after a short time. This happens because of the sensitivity of chaotic systems to changes in the initial conditions [ASY96]. In the next section, we define some terms related to dynamical systems. Three characteristics of chaotic systems are presented in Section 2.3. In Section 2.4, we present some examples of chaotic systems. Finally, Section 2.5 provides a summary of the chapter. 2.2 Dynamical definitions In the previous section, we gave a brief history of dynamical systems. In this section, we define some important dynamical terms that will be used in the remaining chapters. We start by defining the dynamical system. 2.2.1 Dynamical System The dynamical system is a system that consists of a sequence of states which are governed by a certain rule that determines the next state given the previous one. A dynamical system in can be described by either d first order ordinary differential equations (flow) or d difference equations (map). In the first case, the time is a continuous variable; , and the system is called a continuous dynamical system: . (2.1) In the second case, the time is a discrete variable; , and the system is called a discrete dynamical system: . (2.2) ℜd t ℜ1 ∈ d(x(t)) dt  = f(x(t)) n ∈ ℵ x(n + 1) = f(x(n)) 8 2.2.2 Equilibrium point A point in the state space is said to be an equilibrium point of a continuous dynamical system where . The equilibrium point of a discrete dynamical system, on the other hand, happens when . 2.2.3 Orbit of a dynamical system The sequence of points or that results from the solution of the set of equations representing the system is called the orbit (trajectory) of the dynamical system. The point or is called the initial condition of the system. 2.2.4 Attractor The attractor of a dynamical system is set containing the limits of all orbits that start sufficiently close to it. 2.3 Characteristics of Chaotic Systems In this section, we explore three characteristics of chaotic systems. These characteristics will help us to understand how to model these systems. 2.3.1 Sensitivity of chaotic systems to initial conditions Chaotic systems are known for their sensitivity to initial conditions. To illustrate this idea, let’s look into the tent map, which is a first order chaotic system: . (2.3) Figure 2.1 shows two curves (solid, and dotted with ‘x’) representing the responses of the system when two sightly different initial conditions are used. The solid curve was produced from the initial condition and the dotted curve with ‘x’ was produced from x(t) x* = dx(t) dt  x(t) x* = = 0 f x* ( ) x* = x(t) x(n) x0 x(0) x(n + 1) = f(x(n)) = (3 ⁄ 4)(1 – 1 – 2x(n) ) x1(0) = 0.23 9 . The dotted curve in the bottom of the figure is the difference between the two responses. We can see that the difference starts to grow after approximately 15 iterations. Notice that the responses of the system due to different initial conditions remain within a bounded region. Figure 2.1 Sensitivity of the chaotic tent map to initial conditions. 2.3.2 Chaotic signals look random but they are deterministic Before the discovery of chaotic systems, scientists thought of chaos as a random signal (noise). To illustrate this idea, let’s look into the frequency responses (Fourier spectrums [KaSc00]) of a periodic signal and a random signal, and then compare them to a chaotic signal. The left hand side of Figure 2.2 shows a 60 Hz periodic sinusoidal wave. The plot on the right hand side is its frequency response. As we can see, it has only one component at 60 Hz. x2(0) = 0.2301 0 5 10 15 20 25 30 35 40 45 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 The chaotic response of the tent map f(x(n))=x(n+1)=3/4(1−1−2x(n)) time n f x1(0) = 0.23 x2(0) = 0.2301 growth of the error 10 Figure 2.2 The 60 Hz periodic wave and its frequency response If we look into the frequency response of a random signal, shown in the right hand side of Figure 2.3, we can see that it has a component at every frequency from 0 through 120 Hz. Figure 2.3 The random signal and its frequency response We will now compare the frequency responses of the two signals above to the frequency response of a chaotic signal. The left hand side of Figure 2.4 shows the time domain plot of a chaotic signal (x component of the Lorenz model). Its frequency response is shown on the right hand side of the same figure. We can see that the chaotic signal has a component at every frequency from 0 through 120 Hz, which is similar to the random signal. 0 0.005 0.01 0.015 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 The periodic sinusoidal wave of period 60 Hz Time (sec) 20 40 60 80 100 0 50 100 150 Frequency content of the wave frequency (Hz) 0 200 400 600 800 1000 −5 −4 −3 −2 −1 0 1 2 3 4 5 The random signal (noise) Time 0 20 40 60 80 100 120 10−2 10−1 100 101 Frequency content of the signal frequency (Hz) 11 Figure 2.4 The chaotic signal and its frequency response Although the frequency response of the chaotic system looks random, we know that this system is deterministic. This is true since it is produced from a known set of differential equations. 2.3.3 Chaotic systems have attractors with fractal dimension To see that chaotic attractors have fractal dimensions, let’s look into the Henon map: , (2.4) , (2.5) where , , and the initial conditions are [Hen76]. In Figure 2.5, we can see that the attractor of the Henon map is not a simple line of dimension 1, nor it is a closed curve. It is also not a plane of dimension 2. It is an object that fills a region in the plane. 0 5 10 15 20 −30 −20 −10 0 10 20 30 The chaotic Lorenz signal Time (sec) 0 20 40 60 80 100 120 10−3 10−2 10−1 100 101 102 103 104 105 106 Frequency content of y frequency (Hz) x1(n + 1) 1 a(x1(n))2 = – + x2(n) x2(n + 1) = bx1(n) a = 1.4 b = 0.3 x1(0) = x2(0) = 0.5 12 Figure 2.5 The chaotic attractor of the Henon map. As a matter of fact, its dimension is not an integer, but rather is fractal. The dimension of the attractor of the Henon map can be found by the boxcounting dimension. 2.3.3.1 The boxcounting dimension: We can find the dimension of an interval , in , by dividing it into subintervals of length such that , where . Then the minimum number of subintervals needed to cover N is , which can be written as . Notice that the exponent in the expression is 1, which is also the dimension of N. For the case of a unit square in , , we can find the dimension of S by dividing it into small squares of side length (called squares), such that . Then the minimum number of squares needed to cover S is , which can be written as . Notice that the exponent in the expression is 2 which is also the dimension of S. The method used to find the −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 −1.5 −1 −0.5 0 0.5 1 1.5 The chaotic attractor of the Henon map x1 x2 N = [0, 1] ℜ1 Nj σ = 0.1 Nj ∩ Nj + 1 = ∅ j ∈ ℵ Nj cσ N = 10 cσ N (1 ⁄ σ)1 = ℜ2 S x y ℜ2 ∈ 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 ⎩ ⎭ ⎨ ⎬ ⎧ ⎫ = Sj σ = 0.1 σ Sj ∩ Sj + 1 = ∅ Sj cσ S 100 = cσ S (1 ⁄ σ)2 = 13 dimension of the sets in the previous two examples is called the boxcounting dimension. It is used to find the dimension of more complicated sets in , like fractal sets. The boxcounting dimension is denoted by which is the exponent in the relation: as , (2.6) where is the number of dcubes needed to cover the set in . Taking the limit of the log of Equation (2.6) as , we have . (2.7) To evaluate the boxcounting dimension ( ) in Equation (2.7), we draw a loglog plot of versus as . The value of is the slope of the resulting curve. Now we can apply the boxcounting dimension to the attractor of the Henon map. We start by choosing the side length of the squares to be some value . Then we count the minimum number of squares needed to cover the set of points in the attractor (which results from iterating equations (2.4) and (2.5) times). Next we decrease , and count the minimum number needed to cover the attractor for the new . Figure 2.6 shows the number of squares ( ) versus , for . ℜd dc cσ d (1 ⁄ σ)dc ≈ σ →0 cσ d σ ℜd σ → 0 dc cσ d ln ln(1 ⁄ σ)  σ →0 = lim dc cσ d 1 ⁄ σ σ → 0 dc σ σ < 1 105 σ σ σ cσ d σ σ = {0.1, 0.05, 0.01, 0.004, 0.002, 0.001} 14 Figure 2.6 The slope of the line is the boxcounting dimension of the Henon map attractor. The resulting boxcounting dimension of the Henon map was found to be , which is a fractal. After presenting three characteristics of the chaotic systems, in the next section we list some examples of chaotic systems. 2.4 Examples of Chaotic systems 2.4.1 Weather system We said in Section 2.1 that E. Lorenz was the first scientist to quantify chaos. He modelled the heat convection phenomena in fluids by using a set of three differential equations. His model represents earth’s atmosphere. He used it to forecast weather. 2.4.2 Biological models Many biological activities exhibit chaotic behaviors. One example is the Epileptic seizure. Epilepsy causes patients suffering from this disease to experience bouts of unconsciousness. Epileptic seizures result from an abnormal neuronal discharge from the central nervous system [Zyl01]. Another example of chaotic systems is the red blood cell production. Mackey and Glass [MG77] modelled red blood cell production and found that it ex 3 4 5 6 7 8 9 10 7 8 9 10 11 12 13 14 15 16 17 The box dimension of the Henon attractor σε = 0.1, 0.05, 0.01, 0.004, 0.002, 0.001 dc ≈ 1.189 15 hibits chaotic behavior at some parameter values of the model, as will be shown in Chapter 7. A third example of chaotic systems is the Electrocardiogram (ECG). In 1983, Glass et al [GGS83] experimented on the spontaneous beating of cells from embryonic chick hearts. They found that when these cells were stimulated by external periodic pulses, they showed chaotic motion (see also [Moo92]). 2.4.3 Laser signals Measurements taken from a laser output can be represented by a chaotic model. In Chapter 7, we will give more detail regarding the modeling of two data sets of laser outputs. 2.4.4 Chaotic circuits Chaos can be observed in electrical circuits as well. An example of this is the RLC circuit designed by Rulkov et al [RVRDV92]. 2.4.5 Discrete chaotic systems: In the previous section, we have shown two discrete chaotic systems: the tent and the Henon maps. Discrete chaotic systems are mainly used for analysis. In Chapter 7, we will show more discrete chaotic systems used to analyze chaos of different dimensions. 2.5 Chapter Summary In this chapter, we presented a historical background for dynamical systems. We defined some dynamical terms that will be used in subsequent chapters. Three characteristics of chaotic systems were illustrated with examples. We also explored a few examples of dynamical systems that show chaotic behaviors. In the next chapter, we will define dy 16 namical modeling, and show how to implement modeling of chaotic systems by using measurements taken from these systems. 17 CHAPTER 3 INTRODUCTION TO DYNAMICAL MODELING 3.1 Introduction In Chapter 2, we introduced dynamical chaotic systems. We have seen three characteristics of chaotic systems as well. We mentioned that these characteristics will help us to understand modeling of chaotic systems. In this chapter, we explore the modeling process. We also show some of its applications in real life. A theoretical background of the method used for modeling chaotic systems is illustrated with examples. In the next section, we present a brief introduction to modeling chaotic systems. Some applications of the modeling process are presented in Section 3.3. We do modeling by a technique called embedding. Some mathematical background on embedding is contained in Section 3.4. The use of the delaycoordinate map for modeling is illustrated in Section 3.5 with two examples. 3.2 Modeling We said in Chapter 2 that dynamical systems can be represented as a set of differential equations (for a continuous system), or a set of difference equations (for a discrete system). We also said that these equations determine the evolution of states, which converges to an attractor of the system. Knowing the evolution of the states is important for understanding the behavior of the system at a future time. Unfortunately these equations, 18 in most cases, are not known. All we can see from a chaotic system (like those listed in Chapter 2) is a set of scalar measurements. Modeling chaotic systems entails describing the hidden states of the system from these measurements only. From the set of measurements, the modeling process builds a new set of states, which are in some sense equivalent to the original hidden states. In Chapter 8, we will give more detail regarding the meaning of equivalent chaotic systems. 3.3 Applications of modeling 3.3.1 Detection of chaos In Section 2.3.2, we stated that, before the quantification of chaos, scientists thought of chaos as a random signal (noise). Given a set of scalar measurements, the modeling process enables us to distinguish between chaotic systems (deterministic) and random signals (noise). Not only that, but modeling can also extract the hidden deterministic part from a noisy signal. This application of the modeling process is important because, in most cases, the set of measurements is contaminated with noise from different sources. 3.3.2 Prediction The ability to predict the future has fascinated scientists for a long time. Modeling chaotic systems enables us to predict the hidden future states of the system. It does so using only the set of measurements. However, as we have said in Section 2.3.1, chaotic systems have sensitive dependence on the initial conditions. This problem limits the ability of chaotic modeling to make long term predictions of future states. In Chapter 8, we will give more detail. 19 3.3.3 Diagnosis Abarbanel [Aba98] conducted an experiment on the ECG of subjects undergoing a stress test for a specific pathology. He found that in the extreme pathology of ventricular fibrillation, characteristics of the model are different from those of a healthy person. This means the modeling process can be used to diagnose life threatening diseases (for more examples, see [Hol86 Chapters 9 and 11]). In the next two sections, we present some mathematical definitions, then we present the embedding method which is used to build models of chaotic systems. 3.4 Definitions Before we present the embedding method, let’s give a mathematical background of embedding functions. We begin by defining an injection function and an immersion function, then we define the embedding function. 3.4.1 Injection (1to1) function Let M and N be two sets, a function is an injection (1to1) if , (3.1) where . To understand the injection function, consider the next example. 3.4.1.1 Example: An injection Let the function be , (3.2) where . By applying the condition for an injection, we have g: M → N g x1 ( ) g x2 = ( ) x1 ⇒ x2 = x1 x2 , ∈ M g: ℜ2 ℜ2 → g(x) g1(x) g2(x) t x1 (x2)3 t = = x x1 x2 ℜ2 = ∈ 20 . (3.3) From the first element of the vectors on the right hand side we have , (3.4) while from the second element of the vectors we have . (3.5) In conclusion, we can see that , so the function g is an injection. On the other hand, let’s look at an example of a function which is not an injection. 3.4.1.2 Example: A function that is not injection Let the function be . (3.6) To test the condition for injection, we have . (3.7) From the first element of the vectors on the right hand side we have . (3.8) From the second element of the vectors we have . (3.9) Which means that does not have to equal , so the function g is not an injection. g x1 ( ) g x2 ( ) x1 1 x2 1 ( )3 t x1 2 x2 2 ( )3 t = = = x1 1 x1 2 = x2 1 ( )3 x2 2 ( )3 = x2 1 x2 ⇒ = 2 x1 x2 = g: ℜ2 ℜ2 → g(x) x1 (x2)2 t = g x1 ( ) g x2 ( ) x1 1 x2 1 ( )2 t x1 2 x2 2 ( )2 t = = = x1 1 x1 2 = x2 1 ( )2 x2 2 ( )2 = x2 1 ±x2 2 ⇒ = x1 x2 21 3.4.2 Immersion function Let M and N be two sets, a function is an immersion if its Jacobian matrix is an injection (full rank) [Nak90 p.149]. 3.4.2.1 Example: An immersion Let the function be , (3.10) where . The Jacobian matrix of g is . (3.11) The determinant of J is , therefore J is full rank. Hence, the function g is an immersion. 3.4.2.2 Example: A function that is not immersion Let the function be , (3.12) where . The Jacobian matrix of g is . (3.13) g: M → N g: ℜ2 ℜ2 → g(x) g1(x) g2(x) t x1 (x2)3 + x2 t = = x ℜ2 ∈ J 1 0 0 3(x2)2 + 1 x = x0 = 3(x2)2 + 1 ≠ 0 g: ℜ2 ℜ2 → g(x) x1 (x2)3 t = x ℜ2 ∈ J 1 0 0 3(x2)2 x = x0 = 22 The determinant of J is , which zero for . Therefore J is not full rank. Hence g is not an immersion. Now that we have defined the injection and the immersion functions, we will define the embedding function. 3.4.3 Embedding function Let M and N be two sets, a function is an embedding if it is an injection and an immersion [Nak90]. 3.4.3.1 Example: An embedding Let’s look at the function g in Example 3.4.2.1 ( ). In that example, we have seen that g is an immersion. To prove that it is an embedding, we need to show that it is also an injection. To do this, we start by assuming that . (3.14) If we look at the first element of the vectors on the right hand side of Equation (3.14), we can see that . While the second elements give . (3.15) By rearranging the above equation, we have . (3.16) Which can be written as . (3.17) 3(x2)2 x2 = 0 g: M → N g(x) x1 (x2)3 + x2 t = g x1 ( ) g x2 ( ) x1 1 x2 1 ( )3 x2 1 + t x1 2 x2 2 ( )3 x2 2 + t = = = x1 1 x1 2 = x2 1 ( )3 x2 1 + x2 2 ( )3 x2 2 = + x2 1 ( )3 x2 1 x2 2 ( )3 – x2 2 + – = 0 x2 1 x2 2 ( – ) x2 1 ( )2 x2 2 ( )2 x2 1x2 2 ( + + + 1) = 0 23 By looking at the second parenthesis of the left hand side of Equation (3.17), we can see that it can be simplified as follows . (3.18) So the only solution of Equation (3.17) is . As a result, we see that . Which means the function g is an injection. Since g is an injection and an immersion (Example 3.4.2.1), we see that g is an embedding. 3.4.3.2 Example: A nonembedding From Example 3.4.1.1, we have seen that is an injection. But from Example 3.4.2.2, we have seen that this function is not an immersion. So the function g is not an embedding. Notice that the main feature of the embedding is that it is an injection, and its linearization (Jacobian matrix) at every point along the attractor is also an injection. After presenting the embedding, let’s see now how we can use it for modeling. As we said above, we have only a set of scalar measurements from the system we want to model. To build a model for the system, we will use the delaycoordinate map, which is explained below. 3.5 Modeling by Embedding 3.5.1 The delaycoordinate map Let the state of the original system that we want to model be , and let the measurements taken from this system be governed by the map such that . (3.19) x2 1 1 2 x2 2 ⎝ + ⎠ ⎛ ⎞2 3 4  x2 2 ( ) 2 + + 1 ≠ 0 ∀x x2 1 x2 2 = x1 x2 = g(x) x1 (x2)3 t = x ℜk ∈ h: ℜk ℜ1 → y(n) = h(x(n)) 24 The delaycoordinate map [SYC91] is represented by , (3.20) where is the delaytime, and is the delayvector. The attractor built from the delayvectors is called the reconstructed attractor. According to the embedding theorem by Sauer at al [SYC91], the function is an embedding if , where is the boxcounting dimension (see Section 2.3.3.1) of the attractor of the original system. We will call the minimum dimension d that satisfies the embedding condition, the theoretical minimum embedding dimension: . The embedding map guarantees that evolution of the states in the original unknown attractor is equivalent to the evolution of the delayvectors in the reconstructed attractor (Chapter 8). In the next two examples, we illustrate the modeling process by using the delaycoordinate map. 3.5.2 Example: The delaycoordinate map Let the set represent an attractor in . Further, let the measurement function , for the purpose of explanation, be , where . The delaycoordinate map operating on C will be: ( ). The original attractor C is shown in Figure 3.1a, the measurement function h is shown in Figure 3.1b, and the reconstructed attractor using is shown in Figure 3.1c. Fd: ℜk ℜd → yd(n) Fd(h, x(n)) y(n) y(n – T) … y(n – (d – 1)T) t = = T ℵ ∈ yd n ( ) ℜd ∈ Fd d ≥ 2dc + 1 dc dE C y ℜ3 ∈ y1 cos t, y2 2t y3 sint cos2t t ℜ1 = { = = sin , = , ∈ } ℜ3 h:ℜ3 ℜ1 → y = h(x) = x1 + x2 + x3 x ℜ3 ∈ y3(n) F3(h, x(n)) y(n) y(n – 1) y(n – 2) t = = T = 1 F3 25 Figure 3.1 a) The attractor C in , b) the set of measurements, c) The reconstructed attractor by the delaycoordinate map. As we can see, this is a curve of dimension 1 in a three dimensional space. According to the embedding condition; . Which is the same as the dimension of the space required to see the curve C without ambiguity. 3.5.3 Example: Sufficient but not necessary condition for embedding Let the circle; represent an attractor in . By using the delaycoordinate map, we can reconstruct in . The attractor is shown in Figure 3.2 as a solid curve. The reconstructed attractor using the delay coordinate map ( ) is shown in the same figure as a dashed curve, the delaytime (T) is 1. −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 The image of the function Fd Operating on the Curve C in ℜ3 0 2 4 6 8 10 12 14 16 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 The set of observations from the function h:R3 → R1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 X The original attractor of the 1−D set C in ℜ3 Y Z h(x) a ) b ) c ) F3 x(n) y3(n) x2 x1 x3 y(n) y(n – 1) y(n – 2) ℜ3 dE = 2x1 + 1 = 3 S1 y ℜ3 ∈ y1 = cost, y2 = sint, y3 0.5 t ℜ1 = { = , ∈ } ℜ3 S1 ℜ3 S1 F3 26 Figure 3.2 The circle in and its reconstruction in . In Figure 3.2, we can see that is a curve of dimension 1 in a 3D space. Notice that we needed a 2D space to reconstruct this curve. This means the dimension required for embedding is 2 rather than 3, as suggested by the embedding condition ( ). In conclusion, we can see that the embedding condition gives a sufficient but not a necessary condition for embedding. In other words, it is possible that we can embed an attractor in a space of dimension less than . In the next chapter, the delaycoordinate map will be used to model chaotic systems. We will also present four methods used to find the actual, rather than the theoretical, minimum dimension required for embedding. The choice of T (the second parameter in the delaycoordinate map), will be shown in Chapter 6. 3.6 Chapter summary In this chapter, we presented an introduction to modeling chaotic systems. Some applications of modeling chaotic systems were given, mathematical definitions of the embedding functions were illustrated with examples, and we showed two examples of modeling −1 −1 −0.5 0 0.5 1 0 1 −0.5 0 0.5 1 Embedding of a circle by delay−coordinate map Original curve Reconstructed curve S1 S1 ℜ3 ℜ2 S1 d = 2 < dE = 3 dE 27 using the embedding method. We found that the embedding condition is a sufficient but not necessary condition. In the next chapter, we show different methods used to find the minimum dimension required for embedding. 28 CHAPTER 4 INTRODUCTION TO THE MINIMUM EMBEDDING DIMENSION ESTIMATION 4.1 Introduction In Chapter 3, we talked about modeling by the embedding method. Modeling is done from a time series of measurements taken from the system. We said that two parameters have to be found in order to model the system by embedding. Those parameters are the dimension of the delayvector (d) and the delaytime (T). They will be used to build the delayvectors. In this chapter we introduce four methods for estimating the value of d. Three of these methods are geometric and one method is predictive. This chapter will provide only a simple overview of the methods. Chapters 5 and 6 will provide more detail. In the next section, we present modeling of chaotic systems by using the delayvectors. In Section 4.3, we introduce two different approaches for estimating d: geometric and predictive. Then, in Section 4.4, we introduce equivalent chaotic systems. A chapter summary is given in Section 4.5. In Chapter 6, we will discuss the selection of the delaytime T. 4.2 Modeling chaotic systems The state evolution of a chaotic system in the space can be written in the form of the difference equation: ℜk 29 , (4.1) where is the state of the system, the time index , and N is the total number of points. In the steady state condition, the evolution of the state follows an attractor with a fractal dimension (chaotic attractor). Practically, the state and the function f of the system are invisible to us, and all we can see is a set of scalar measurements . We can write these measurements as (see Equations (3.19). To model a chaotic system, we need to build a system model which uses the delayvector . (4.2) The attractor of the system should be equivalent to the original attractor, and the state evolution from should follow that of the original attractor from [Hak98]. (A more careful definition of equivalence is contained in Chapter 8.) The key idea is to find d and T such that the two systems are equivalent. The embedding theorem guarantees this equivalence if . But in practice, we do not generally know the value of , so d has to be estimated. The projection from the original state to the delayvector is called the delaycoordinate map. As we said in the previous chapter, the embedding theorem gives a sufficient, but not a necessary, condition for embedding. In other words, it could be possible to find an embedding map at . Our goal is to estimate the minimum embedding dimension. We will label this dimension . In the next section, we will talk about estimating by using two techniques: a geometric and a predictive technique. In Chapter 6, we will discuss the delaytime T to complete the modeling process. x(m) = f(x(m – 1)) x m ( ) ℜk ∈ m = 1, 2, …, N x(m) dc x(m) y(m) y(m) = h(x(m)) yd(m) y(m) y(m – T) … y(m – (d – 1)T) t = yd(m) yd(m) → yd(m + 1) x(m) → x(m + 1) d ≥ 2dc + 1 dc x yd d ≤ 2dc dL dL 30 4.3 Estimating the minimum embedding dimension 4.3.1 The geometric technique To understand the geometric technique, let the circle in Figure 4.1 represent an attractor of a dynamical system in , and the line below it represent its projection into . Figure 4.1 Projection artifacts From the above figure, we can see that the points on the horizontal line (a, b, and c) do not occur in the same sequence as those on the circle ( ). That happened because the circle was projected into a space with insufficient dimension. That caused the distances between points to shrink, and the original order of points to change. We call this effect the projection artifact. We can see from this example that, if the dimension of the space is too small to represent the attractor of the system, the projection of the attractor to this space will have artifacts. Therefore, we need to increase the dimension of the space in order to get rid of these artifacts. One technique of estimating the minimum embedding dimension is the geometric technique. To estimate by using the geometric technique, we start from and find for every point its nearest neighbor. This neighbor has to be tested to see if it is a true neighbor or if it became a neighbor because of some projection artifacts. After that, d is increased and the previous steps are repeated. The value of d where all the neighbors are true neighbors, is the minimum embedding dimension . Under the geometric tech ℜ2 ℜ1 a′ b′ c′ a c b projection a′, b′, and c′ dL dL d = 1 yd(m) dL 31 nique, three methods can be used to detect the existence of projection artifacts. These methods are 1) the Change of Neighbors with Dimension method, or CND, 2) the Change of Distance with Dimension method, or CDD, and 3) the Change of Distance with Time method, or CDT. 4.3.1.1 The change of neighbors with dimension method (CND) One can perform a visual test using Figure 4.1 to determine that the horizontal line does not have enough dimension to represent the structure of the circle, and that we need to have a two dimensional space (meaning that ). However, we need to automate this test using an algorithm. To estimate using the first geometric method, the change of neighbors with dimension method (CND), the algorithm starts from the scalar measurements a, b, and c, and computes the distances between point a and the other two points. It will find that the nearest neighbor to a is c, while on the circle, the nearest neighbor to is . The algorithm concludes that the neighbors have changed, and that the points a and c became neighbors because of the projection artifact and not because they are true neighbors. That means the 1D space is not large enough to represent the structure of the circle since it has this projection artifact. By repeating the above steps on the 2D space, the algorithm would find that there are no projection artifacts left on the attractor in this space and concludes that . 4.3.1.2 The change of distance with dimension (CDD) method The second geometric method that can be used to estimate is the change of distance with dimension method (CDD). The CDD method is similar to the CND method except that it detects the existence of projection artifacts by comparing the distances between neighbors rather than comparing the neighbors themselves. To understand this method, let dL = 2 dL a′ b′ dL = 2 dL 32 the curve in Figure 4.2 represent an attractor of a dynamical system in , and the line below it represent its projection into . Figure 4.2 The CDD method We can see the points , , and on the curve in and their projection into the horizontal line a, b, and c respectively. As in the CND method, we want an algorithm to find that the 1D space is not large enough to represent the structure of the curve. From the horizontal line, the algorithm can find that b is the nearest neighbor of a and the distance between them in is . In , along the curve, the distance between and is . The algorithm can now compare the two distances to find that , therefore, b became a neighbor of a because of the projection artifact and not because they are true neighbors. As a result, the algorithm will find that the 1D space is not large enough to represent the structure of the curve since it has this projection artifact. By repeating the above steps on the 2D space, the algorithm would find that there are no projection artifacts left on the attractor in this space and concludes that . 4.3.1.3 The change of distance with time method (CDT) The third geometric method used to estimate is the change of distance with time method (CDT). This method is different from the previous two in that it tries to detect the existence of projection artifacts by moving neighbors forward in time and checking to see ℜ2 ℜ1 a' b′ c′ a b c ℜ2 ℜ1 projection v1 v2 a' b' c' ℜ2 ℜ1 v1 ℜ2 a′ b′ v2 v1 « v2 dL = 2 dL 33 if they will remain neighbors. To understand this method, let the nonintersecting solid curve in shown in Figure 4.3 represent an attractor of a dynamical system, and let the dotted curve with the “+” sign represent its projection into the XY plane ( space). Figure 4.3 The projection of a 3D curve (solid line) into 2D curve (dotted line) Figure 4.4 The CDT method: a) Distances between false neighbors increase with time, b) The intersection in the middle of the curve is a projection artifact The projection of the systems’s attractor shown in Figure 4.3 is redrawn in Figure 4.4 for further discussion. Notice that the middle of the curve in Figure 4.4a has a projection artifact, which can be found by looking into Figure 4.4b. If we start from point a, we can’t tell whether the next correct point in time is b or c. In other words, the curve intersection makes it impossible to determine the correct sequence of points. ℜ3 ℜ2 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 X 2−D projection of a curve in ℜ3 Y Z c a d b v1 v2 ℜ2 f a b f a ) b ) 34 By using the CDT method, we want an algorithm to be able to find that the 2D space is not large enough to represent the structure of the system’s attractor (which is in : ). This can be done by detecting the projection artifact in the middle of the curve. To do that, the algorithm should start by calculating the distance between point a and the other points on the curve. It will find that c is closest to a; let the distance between them be . Next it should move forward in time to find that a moves to b, and c moves to d. It should then compute the distance between b and d; let the distance be . Now it can compares the two distances and find that , which means that a and c are neighbors because of some projection artifacts and not because they are true neighbors. (See Figure 4.4a) Hence the algorithm concludes that the 2D space is not large enough to represent the structure of the system’s attractor since the curve has this projection artifact. By repeating this for the 3D space, the algorithm would find no projection artifacts left on the curve in this space and would conclude that . 4.3.2 The predictive technique: Now that we have talked about the three geometric methods, we will discuss the predictive technique. The predictive technique is also used to estimate the minimum embedding dimension such that the system of (using delay embedding) is equivalent to the system of (original system). For simplicity, we will set in the delayvector (see Equations (4.2)). By using the delaycoordinate map; , we can write (see Equations (3.20)). From the embedding definition (see Section 3.4.3), we know that if the map is an embedding, then it is an injection. Hence, the inverse map exists: . (4.3) ℜ3 dL = 3 v1 v2 v1 « v2 dL = 3 dL ydL (m) x(m) T = 1 yd(m) Fd yd(m) = Fd(x(m)) Fd Fd –1:ℜd ℜk → x(m) Fd –1 = (yd(m)) 35 We can substitute Equations (4.1) for the measurement function h at Equations (3.19) ( ) to write h as . If we let the composite function , we can write . (4.4) By substituting Equations (4.3) for Equations (4.4), we can write . (4.5) Now, if we let the composite function , we can write . (4.6) We conclude from Equations (4.6) that if the delayvector at time is known to us, we can predict the current measurement by approximating the unknown function . The predictive technique depends on the idea that when the system of is equivalent to the system of , we can use the delayvector to predict the current measurement . To be able to do that, we need to approximate the unknown function such that, . (4.7) To approximate we will use a neural network with a Tapped Delay Line (TDL) connected to its input. The TDL is used to produce the delayvector from the current measurement , as seen in Figure 4.6. Each tap of the TDL is a delayed version of and the total number of taps is d. The network is trained to predict when it is presented with d previous measurements of (coordinates of ). The resulting prediction error is recorded as a function of d. As the predictor order d increases, the prediction error will decrease. However, after a certain point, further increase of d results in only a very small decrease in the error. The minimum dimension where the error y(m) = h(x(m)) y(m) = h(f(x(m – 1))) h o f = q y(m) = q(x(m – 1)) y(m) q Fd –1 = ( (yd(m – 1))) q o Fd –1 = μ y(m) μ(yd(m – 1)) μ y(m – 1) y(m – 2) … y(m – d) t ⎝ ⎠ = = ⎛ ⎞ m – 1 y(m) μ yd(m) x(m) yd(m – 1) y(m) μ yˆ m ( ) μ ˆ = (yd(m – 1)) μ yd(m – 1) y(m) y(m) y(m) y(m) yd(m – 1) 36 does not improve any further is the minimum embedding dimension . At this point, the approximation of is accurate and the systems of and are equivalent to each other. Figures 4.5 and 4.6 show the block diagram of the function approximation and the neural network model used to create the approximation. Figure 4.5 The function approximation in the predictive technique Figure 4.6 A neural network model with a TDL used to approximate the function , D is a one step delaytime 4.4 Dynamic equivalence The goal of modeling by embedding is to find a chaotic model that is equivalent to the original unknown system. (The notion of equivalent chaotic systems will be covered in detail in later chapters.) To understand this idea, we know that in the original system the hidden state evolution is governed by the map f , (4.8) dL μ yd(m) x(m) μ + y(m) yˆ (m)  + μ ˆ (.) e(m) (prediction error) TDL μ(.) yd(m – 1) μ D y(m) yˆ y(m – 1) (m) y(m – 2) y(m – d) D D Neural Network Model μ ˆ (.) yd(m – 1) μ x(m + 1) = f(x(m)) 37 (see Equations (4.1)), while all we can see from the original system is a set of measurements which are governed by the map h: , (4.9) (see Equations (3.19)). On the other hand, the states in the reconstructed space are the delayvectors (see Equations (4.2)), while the output update in this space is governed by the map : , (4.10) (see Equations (4.6)). So the delayvector can be written as , (4.11) ( ). If the reconstructed model given by Equations (4.10) and (4.11) is equivalent to the original model given by Equations (4.8) and (4.9), we can use the evolution of in place of the evolution of to gain an understanding of the system characteristics. In other words, if the two systems are equivalent, we will be able to use to estimate the original system dimension and parameters. In the coming chapters, we will find the delayvector dimension d and the delaytime T which guarantee this dynamic equivalence. 4.5 Chapter summary In this chapter, we provided an introduction to modeling chaotic systems. We introduced two techniques (geometric and predictive) that can be used to estimate the minimum embedding dimension that is required for chaotic systems modeling. We also gave a brief introduction to equivalent chaotic systems. In Chapter 5, we will see the application of both techniques to the chaotic Henon map. In Chapter 6, we will provide complete and detailed algorithms for estimating d. We will also discuss practical considerations for implementing the algorithms, including the choice of the delaytime (T). y(m) = h(x(m)) yd(m) μ y(m) = μ(yd(m – 1)) yd(m) μ(yd(m – 1)) y(m – 1) … y(m – (d – 1)) t = T = 1 yd(m) x(m) yd(m) 38 CHAPTER 5 Examples of the minimum embedding dimension estimation 5.1 Introduction In Chapter 3, we introduced modeling of chaotic systems. We used delayvectors created from measurements to model these systems. We also introduced four methods to estimate the minimum embedding dimension: three geometric methods one predictive method. In this chapter, we will use the four methods to estimate the minimum embedding dimension of the Henon map. The Henon map is a chaotic system created from a set of two difference equations (see Section 2.3.3). By using the embedding method (see Section 3.5), the theoretical minimum embedding dimension is . But we suspect that the actual minimum embedding dimension is since the system’s dynamics are generated from a set of two difference equations. In Sections 5.2 and 5.3, we show the estimated minimum embedding dimension of the Henon map by using the geometric and the predictive techniques, respectively. In Section 5.4, we provide a chapter summary. 5.2 The geometric technique To show how we can estimate by using the geometric technique, we take 15 points from the coordinate of the Henon map to represent the measurements from this dE = 2dc = 2x1.189 = 3 dL = 2 dL X1 39 system. That means the measurements where . In the next step, we use to construct which is the delayvector (where ). We said in Chapter 4 that the idea of estimating by using the geometric technique depends on detecting the existence of the projection artifacts on the system’s attractor. We can do that by using any of the three methods that we mentioned before: the CND, the CDD, or the CDT method. We need first to compute the distance between the reference vector; and every other vector in the space . The nearest neighbor of is the vector with the shortest distance. 5.2.1 Using the change of neighbors with dimension method (CND) Before we get into the details of the CND method, we need to introduce some notation. First, the nearest neighbor of will be denoted . The nearest neighbor of will be denoted . (This means there are vectors that are closer to .) We will indicate the time index of the nearest neighbor as . (5.1) For example, if the nearest neighbor of is , then the nearest neighbor index is . To estimate by using the CND method, we need to check if the nearest neighbor of will remain a neighbor as the dimension d grows to . In other words, we need to check if the nearest neighbor of will appear as the first, second, , or neighbor of . To do that, we need to build the matrix which has the elements where and . If , we label as y(m) = x1(m) m = 1, 2, …, 15 y(m) yd(m) y(m) y(m – 1) … y(m – (d – 1)) t = T = 1 dL yd m ( ) ℜd yd(m) yd(m) y d 1 m ( ) ) kth yd(m) y d k (m) ) k – 1 yd(m) kth id k (m) index y d k = [ (m)] ) yd(5) y d 1 (5) = yd(9) ) id 1 (5) index y d 1 = [ (5)] = index[yd(9)] = 9 ) dL yd(m) d + 1 yd(m) … wth yd + 1(m) Qd qd(i, j) = yd(i) – yd(j) i, j = 1, 2, …, M i ≠ j i = j qd(i, j) 40 not a number (NaN) since we are not interested in the distance between a vector and itself. We need also to find the vector whose element is the time index of the nearest neighbor of . Further, we need to compute the through the neighbors of each and save their indices in the row of the matrix . In other words, , and are defined as follows: , (5.2) , (5.3) where , (5.4) and . (5.5) Each row of the matrix will be labeled by . So, the maid 1 mth id 1(m) yd(m) 1st wth yd + 1(m) mth Id + 1 Qd id k Id Qd qd(1, 1) qd(1, 2) … qd(1, M) qd(2, 1) qd(2, 2) … qd(2, M) : qd(M, 1) qd(M, 2) … qd(M, M) = Qd NaN yd(1) – yd(2) … yd(1) – yd(M) yd(2) – yd(1) NaN … yd(2) – yd(M) : : yd(M) – yd(1) yd(M) – yd(2) … NaN = id k id k (1) id k (2) : id k (M) index y d k [ (1)] index y d k [ (2)] : index y d k [ (M)] = = ) ) ) k = 1, 2, …, w Id id 1 id 2 … id w id 1 (1) id 2 (1) … id w (1) id 1 (2) id 2 (2) … id w (2) : : id 1 (M) id 2 (M) … id w (M) = = Id i m d id 1 (m) id 2 (m) … id w = (m) 41 trix can be written as: . In Tables 5.1a and 5.1c and Tables 5.2a and 5.2c, we can see listed for through 4. In Tables 5.1d, 5.2b and 5.2d, we can see through showing the indices of the three neighbors ( ) of each . Tables 5.3a, 5.3b, and 5.3c summarize the neighbors time indices found for through 4. For example, in Table 5.3a we can see that the nearest neighbor of is and it fails to appear as the first, second, or third neighbor of , so we label as a FNN. On the other hand, we can see that the nearest neighbor of is and when the dimension increases to 2, it appears as the second neighbor of , so we label as a true neighbor. As a conclusion, we can see in Table 5.3d that the total number of FNNs drops from 2 at to 0 at and remains 0 at . That means the minimum embedding dimension is , where there are no projection artifacts left on the attractor of the system. Notice that in this method we could have chosen , so that only the nearest neighbor is considered. We have found through experimentation that using w greater than 1 provides a more robust algorithm. This is a new result. In Chapter 6, we will provide more detail on the practical implementation aspects of the various algorithms. Id Id i 1 d i 2 d : i M d = Qd d = 1 I2 I4 w = 3 yd + 1(m) d = 1 y1(1) y1(8) y2(1) y1(8) y1(8) y1(9) y2(8) y1(9) d = 1 d = 2 d = 3 dL = 2 w = 1 42 Table 5.1 The CND method (1/3) Table 5.2 The CND method (2/3) Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 NaN 0.65 1.2015 0.1561 0.313 0.6479 1.2886 0.0262 0.0508 0.3597 0.4 0.7439 1.6362 0.9341 0.1046 8 2 0.65 NaN 1.8515 0.4939 0.963 0.0021 1.9386 0.6762 0.7008 0.2903 1.05 0.0939 2.2862 1.5841 0.7546 6 3 1.2015 1.8515 NaN 1.3576 0.8885 1.8494 0.0871 1.1753 1.1507 1.5612 0.8015 1.9454 0.4347 0.2674 1.0969 7 4 0.1561 0.4939 1.3576 NaN 0.4691 0.4918 1.4446 0.1823 0.2069 0.2036 0.556 0.5878 1.7922 1.0901 0.2607 1 5 0.313 0.963 0.8885 0.4691 NaN 0.9609 0.9755 0.2868 0.2622 0.6727 0.0869 1.0569 1.3231 0.621 0.2084 11 6 0.6479 0.0021 1.8494 0.4918 0.9609 NaN 1.9364 0.6741 0.6987 0.2882 1.0478 0.096 2.284 1.5819 0.7525 2 7 1.2886 1.9386 0.0871 1.4446 0.9755 1.9364 NaN 1.2624 1.2377 1.6483 0.8886 2.0325 0.3476 0.3545 1.1839 3 8 0.0262 0.6762 1.1753 0.1823 0.2868 0.6741 1.2624 NaN 0.0246 0.3859 0.3737 0.7701 1.61 0.9078 0.0784 9 9 0.0508 0.7008 1.1507 0.2069 0.2622 0.6987 1.2377 0.0246 NaN 0.4105 0.3491 0.7947 1.5853 0.8832 0.0538 8 10 0.3597 0.2903 1.5612 0.2036 0.6727 0.2882 1.6483 0.3859 0.4105 NaN 0.7596 0.3842 1.9959 1.2937 0.4643 4 11 0.4 1.05 0.8015 0.556 0.0869 1.0478 0.8886 0.3737 0.3491 0.7596 NaN 1.1438 1.2362 0.5341 0.2953 5 12 0.7439 0.0939 1.9454 0.5878 1.0569 0.096 2.0325 0.7701 0.7947 0.3842 1.1438 NaN 2.3801 1.6779 0.8485 2 13 1.6362 2.2862 0.4347 1.7922 1.3231 2.284 0.3476 1.61 1.5853 1.9959 1.2362 2.3801 NaN 0.7021 1.5316 7 14 0.9341 1.5841 0.2674 1.0901 0.621 1.5819 0.3545 0.9078 0.8832 1.2937 0.5341 1.6779 0.7021 NaN 0.8294 3 15 0.1046 0.7546 1.0969 0.2607 0.2084 0.7525 1.1839 0.0784 0.0538 0.4643 0.2953 0.8485 1.5316 0.8294 NaN 9 a) b) Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 NaN 1.9623 1.2991 0.9756 0.313 2.044 1.4552 0.7013 0.2947 1.1099 0.4108 2.4041 2.2773 1.2008 9 5 11 2 1.9623 NaN 2.2959 1.0165 2.0851 0.0871 2.267 1.3346 1.7113 0.8525 2.2106 0.4447 2.3018 1.9268 6 12 10 3 1.2991 2.2959 NaN 1.4363 1.0155 2.3467 0.202 1.1934 1.1685 1.6572 0.994 2.6451 1.1736 0.3735 7 14 11 4 0.9756 1.0165 1.4363 NaN 1.0693 1.0925 1.4728 0.3193 0.7038 0.2214 1.1942 1.4478 1.8968 1.1099 10 8 9 5 0.313 2.0851 1.0155 1.0693 NaN 2.1617 1.1858 0.7553 0.3896 1.2452 0.1295 2.5167 2.0623 0.9757 11 1 9 6 2.044 0.0871 2.3467 1.0925 2.1617 NaN 2.3116 1.4094 1.7902 0.9342 2.2867 0.3606 2.3114 1.9759 2 12 10 7 1.4552 2.267 0.202 1.4728 1.1858 2.3116 NaN 1.2626 1.2965 1.6901 1.1759 2.5928 0.9721 0.3631 3 14 13 8 0.7013 1.3346 1.1934 0.3193 0.7553 1.4094 1.2626 NaN 0.4113 0.5204 0.8782 1.7625 1.8363 0.9094 4 9 10 9 0.2947 1.7113 1.1685 0.7038 0.3896 1.7902 1.2965 0.4113 NaN 0.8635 0.5191 2.1483 2.0462 0.9978 1 5 8 10 1.1099 0.8525 1.6572 0.2214 1.2452 0.9342 1.6901 0.5204 0.8635 NaN 1.3731 1.2945 2.0661 1.327 4 8 2 11 0.4108 2.2106 0.994 1.1942 0.1295 2.2867 1.1759 0.8782 0.5191 1.3731 NaN 2.6407 2.0842 1.0026 5 1 9 12 2.4041 0.4447 2.6451 1.4478 2.5167 0.3606 2.5928 1.7625 2.1483 1.2945 2.6407 NaN 2.4815 2.2718 6 2 10 13 2.2773 2.3018 1.1736 1.8968 2.0623 2.3114 0.9721 1.8363 2.0462 2.0661 2.0842 2.4815 NaN 1.0867 7 14 3 14 1.2008 1.9268 0.3735 1.1099 0.9757 1.9759 0.3631 0.9094 0.9978 1.327 1.0026 2.2718 1.0867 NaN 7 3 8 c) d) i1 1 I2 a ) b ) c ) d ) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 11 9 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 10 4 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 9 13 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 8 2 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 11 9 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 12 10 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 13 8 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 10 9 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 8 5 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 8 2 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 1 9 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 2 10 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 3 8 a) b) Q4 1 2 3 4 5 6 7 8 9 10 11 12 1 NaN 2.4318 1.6489 2.5414 0.3725 2.3668 1.8663 1.7995 1.0368 2.8685 1.2434 2.433 5 9 11 2 2.4318 NaN 2.5327 1.4923 2.5528 0.331 2.3738 1.3529 2.0868 1.6802 2.9129 1.1956 6 12 8 3 1.6489 2.5327 NaN 2.5954 1.5612 2.4653 0.4389 1.7247 1.1757 3.0134 2.2894 2.8193 7 9 5 4 2.5414 1.4923 2.5954 NaN 2.5469 1.7832 2.3182 0.9872 2.3925 0.4232 2.6017 2.4496 10 8 2 5 0.3725 2.5528 1.5612 2.5469 NaN 2.5035 1.757 1.8512 1.2387 2.8763 0.9807 2.5428 1 11 9 6 2.3668 0.331 2.4653 1.7832 2.5035 NaN 2.3479 1.5024 1.9941 1.9947 2.9327 0.9783 2 12 8 7 1.8663 2.3738 0.4389 2.3182 1.757 2.3479 NaN 1.5296 1.3966 2.7334 2.36 2.7782 3 9 8 8 1.7995 1.3529 1.7247 0.9872 1.8512 1.5024 1.5296 NaN 1.4334 1.3952 2.245 2.2062 4 2 10 9 1.0368 2.0868 1.1757 2.3925 1.2387 1.9941 1.3966 1.4334 NaN 2.7782 2.1478 2.3707 1 3 5 10 2.8685 1.6802 3.0134 0.4232 2.8763 1.9947 2.7334 1.3952 2.7782 NaN 2.836 2.6148 4 8 2 11 1.2434 2.9129 2.2894 2.6017 0.9807 2.9327 2.36 2.245 2.1478 2.836 NaN 2.8555 5 1 9 12 2.433 1.1956 2.8193 2.4496 2.5428 0.9783 2.7782 2.2062 2.3707 2.6148 2.8555 NaN 6 2 8 c) d) I3 I4 a ) c ) d ) b ) 43 Table 5.3 The CND method (3/3) 5.2.2 Using the change of distance with dimension method (CDD) As we said in Chapter 4, the CDD method can also be used to estimate . This method depends on the idea that false nearest neighbor distances increase significantly as the dimension of the space increases. To estimate of the Henon map using the CDD method, we need to compute the matrix and find the vector as in the CND method. In addition, we need to find the vector whose element is the distance between and its nearest neighbor : . We also need to find out how much the distances to the nearest neighbors grow as the dimension increases. For these new distances we define the vector whose elements represent the distance between and . (Notice that is not the same as . It is with one component added.) For example, when in Table 5.4a, the nearest neighbor of is and the distance between them is . At in Table 5.4d, the distance between and is . If is a true nearest neighbor of , will be close to . To apply this idea, we need to see how d d d d FNN time 1 2 time 2 3 time 3 4 1 2 1 8 9 5 11 FNN 1 9 5 11 9 1 5 5 9 11 2 0 2 6 6 12 10 2 6 6 10 4 2 6 6 12 8 3 0 3 7 7 14 11 3 7 7 9 13 3 7 7 9 5 d) 4 1 10 8 9 FNN 4 10 10 8 2 4 10 10 8 2 5 11 11 1 9 5 11 1 11 9 5 1 1 11 9 6 2 2 12 10 6 2 2 12 10 6 2 2 12 8 7 3 3 14 13 7 3 3 13 8 7 3 3 9 8 8 9 4 9 10 8 4 4 10 9 8 4 4 2 10 9 8 1 5 8 9 1 1 8 5 9 1 1 3 5 10 4 4 8 2 10 4 4 8 2 10 4 4 8 2 11 5 5 1 9 11 5 5 1 9 11 5 5 1 9 12 2 6 2 10 12 6 6 2 10 12 6 6 2 8 13 7 7 14 3 13 7 7 3 8 13 7 14 3 7 3 8 14 7 c) 15 9 aa)) b )b) c ) d ) dL dL Qd id 1 rd rd(m) yd(m) y d 1 (m) yd id 1 = ( (m)) ) rd(m) yd(m) y d 1 m ( ) – = ) ed + 1 ed + 1(m) yd + 1(m) yd + 1 id 1 ( (m)) yd + 1 id 1 ( (m)) y d + 1 1 (m) ) y d 1 (m) ) d = 1 y1(1) y1(8) r1(1) = 0.026 d = 2 y2(1) y2(8) e2(1) = 0.701 yd id 1 ( (m)) yd(m) yd + 1 id 1 ( (m)) yd + 1(m) 44 much the distance between the nearest neighbors grows as d increases to . We can do that by forming the vector which has the elements . (5.6) If where is some predefined threshold, we label the nearest neighbor of as a false nearest neighbor (FNN). For instance, let . We can see from our example that as seen at the second column of Table 5.4g. Since , we label as a FNN. The results in Tables 5.4g, 5.5g, and 5.6g show for through 3. In Table 5.6h, we summarize the results of the previous tables. In it, we can see that at , the number of FNNs is 5, while at , the number of FNNs drops to 0 and remains 0 at . Therefore, the minimum embedding dimension is . d + 1 cd cd(m) ed + 1(m) – rd(m) rd(m) =  cd(m) > ρ ρ yd(m) ρ = 10 c1(1) = 25.9731 c1(1) > ρ = 10 y1(8) cd d = 1 d = 1 d = 2 d = 3 dL = 2 45 Table 5.4 The CDD method tables for d = 1 Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 i1 r1 1 NaN 0.650 1.202 0.156 0.313 0.648 1.289 0.026 0.051 0.360 0.400 0.744 1.636 0.934 0.105 8 0.026 2 0.650 NaN 1.852 0.494 0.963 0.002 1.939 0.676 0.701 0.290 1.050 0.094 2.286 1.584 0.755 6 0.002 3 1.202 1.852 NaN 1.358 0.889 1.849 0.087 1.175 1.151 1.561 0.802 1.945 0.435 0.267 1.097 7 0.087 4 0.156 0.494 1.358 NaN 0.469 0.492 1.445 0.182 0.207 0.204 0.556 0.588 1.792 1.090 0.261 1 0.156 5 0.313 0.963 0.889 0.469 NaN 0.961 0.976 0.287 0.262 0.673 0.087 1.057 1.323 0.621 0.208 11 0.087 6 0.648 0.002 1.849 0.492 0.961 NaN 1.936 0.674 0.699 0.288 1.048 0.096 2.284 1.582 0.753 2 0.002 7 1.289 1.939 0.087 1.445 0.976 1.936 NaN 1.262 1.238 1.648 0.889 2.033 0.348 0.355 1.184 3 0.087 8 0.026 0.676 1.175 0.182 0.287 0.674 1.262 NaN 0.025 0.386 0.374 0.770 1.610 0.908 0.078 9 0.025 9 0.051 0.701 1.151 0.207 0.262 0.699 1.238 0.025 NaN 0.411 0.349 0.795 1.585 0.883 0.054 8 0.025 10 0.360 0.290 1.561 0.204 0.673 0.288 1.648 0.386 0.411 NaN 0.760 0.384 1.996 1.294 0.464 4 0.204 11 0.400 1.050 0.802 0.556 0.087 1.048 0.889 0.374 0.349 0.760 NaN 1.144 1.236 0.534 0.295 5 0.087 12 0.744 0.094 1.945 0.588 1.057 0.096 2.033 0.770 0.795 0.384 1.144 NaN 2.380 1.678 0.849 2 0.094 13 1.636 2.286 0.435 1.792 1.323 2.284 0.348 1.610 1.585 1.996 1.236 2.380 NaN 0.702 1.532 7 0.348 14 0.934 1.584 0.267 1.090 0.621 1.582 0.355 0.908 0.883 1.294 0.534 1.678 0.702 NaN 0.829 3 0.267 15 0.105 0.755 1.097 0.261 0.208 0.753 1.184 0.078 0.054 0.464 0.295 0.849 1.532 0.829 NaN 9 0.054 a) b) c) Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 i1 e2 1 NaN 1.962 1.299 0.976 0.313 2.044 1.455 0.701 0.295 1.110 0.411 2.404 2.277 1.201 8 0.701 2 1.962 NaN 2.296 1.017 2.085 0.087 2.267 1.335 1.711 0.853 2.211 0.445 2.302 1.927 6 0.087 3 1.299 2.296 NaN 1.436 1.016 2.347 0.202 1.193 1.169 1.657 0.994 2.645 1.174 0.374 7 0.202 4 0.976 1.017 1.436 NaN 1.069 1.093 1.473 0.319 0.704 0.221 1.194 1.448 1.897 1.110 1 0.976 5 0.313 2.085 1.016 1.069 NaN 2.162 1.186 0.755 0.390 1.245 0.130 2.517 2.062 0.976 11 0.130 6 2.044 0.087 2.347 1.093 2.162 NaN 2.312 1.409 1.790 0.934 2.287 0.361 2.311 1.976 2 0.087 7 1.455 2.267 0.202 1.473 1.186 2.312 NaN 1.263 1.297 1.690 1.176 2.593 0.972 0.363 3 0.202 8 0.701 1.335 1.193 0.319 0.755 1.409 1.263 NaN 0.411 0.520 0.878 1.763 1.836 0.909 9 0.411 9 0.295 1.711 1.169 0.704 0.390 1.790 1.297 0.411 NaN 0.864 0.519 2.148 2.046 0.998 8 0.411 10 1.110 0.853 1.657 0.221 1.245 0.934 1.690 0.520 0.864 NaN 1.373 1.295 2.066 1.327 4 0.221 11 0.411 2.211 0.994 1.194 0.130 2.287 1.176 0.878 0.519 1.373 NaN 2.641 2.084 1.003 5 0.130 12 2.404 0.445 2.645 1.448 2.517 0.361 2.593 1.763 2.148 1.295 2.641 NaN 2.482 2.272 2 0.445 13 2.277 2.302 1.174 1.897 2.062 2.311 0.972 1.836 2.046 2.066 2.084 2.482 NaN 1.087 7 0.972 14 1.201 1.927 0.374 1.110 0.976 1.976 0.363 0.909 0.998 1.327 1.003 2.272 1.087 NaN 3 0.374 d) e) f) Time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 r1 0.026 0.002 0.087 0.156 0.087 0.002 0.087 0.025 0.025 0.204 0.087 0.094 0.348 0.267 e2 0.701 0.087 0.202 0.976 0.130 0.087 0.202 0.411 0.411 0.221 0.130 0.445 0.972 0.374 c1 25.9731 42.55 1.32184 5.25385 0.48851 42.55 1.32184 15.452 15.452 0.08529 0.48851 3.73085 1.79339 0.39888 FNN FNN FNN FNN FNN g) i1 1 e2 r1 i1 1 a ) d ) b ) c ) e ) f ) g ) 46 Table 5.5 The CDD method tables for d = 2 Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 i2 r2 1 NaN 1.962 1.299 0.976 0.313 2.044 1.455 0.701 0.295 1.110 0.411 2.404 2.277 1.201 9 0.295 2 1.962 NaN 2.296 1.017 2.085 0.087 2.267 1.335 1.711 0.853 2.211 0.445 2.302 1.927 6 0.087 3 1.299 2.296 NaN 1.436 1.016 2.347 0.202 1.193 1.169 1.657 0.994 2.645 1.174 0.374 7 0.202 4 0.976 1.017 1.436 NaN 1.069 1.093 1.473 0.319 0.704 0.221 1.194 1.448 1.897 1.110 10 0.221 5 0.313 2.085 1.016 1.069 NaN 2.162 1.186 0.755 0.390 1.245 0.130 2.517 2.062 0.976 11 0.13 6 2.044 0.087 2.347 1.093 2.162 NaN 2.312 1.409 1.790 0.934 2.287 0.361 2.311 1.976 2 0.087 7 1.455 2.267 0.202 1.473 1.186 2.312 NaN 1.263 1.297 1.690 1.176 2.593 0.972 0.363 3 0.202 8 0.701 1.335 1.193 0.319 0.755 1.409 1.263 NaN 0.411 0.520 0.878 1.763 1.836 0.909 4 0.319 9 0.295 1.711 1.169 0.704 0.390 1.790 1.297 0.411 NaN 0.864 0.519 2.148 2.046 0.998 1 0.295 10 1.110 0.853 1.657 0.221 1.245 0.934 1.690 0.520 0.864 NaN 1.373 1.295 2.066 1.327 4 0.221 11 0.411 2.211 0.994 1.194 0.130 2.287 1.176 0.878 0.519 1.373 NaN 2.641 2.084 1.003 5 0.13 12 2.404 0.445 2.645 1.448 2.517 0.361 2.593 1.763 2.148 1.295 2.641 NaN 2.482 2.272 6 0.361 13 2.277 2.302 1.174 1.897 2.062 2.311 0.972 1.836 2.046 2.066 2.084 2.482 NaN 1.087 7 0.972 14 1.201 1.927 0.374 1.110 0.976 1.976 0.363 0.909 0.998 1.327 1.003 2.272 1.087 NaN 7 0.363 a) b) c) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 i2 e3 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 9 0.854 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 0.202 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 0.331 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 0.241 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 0.371 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 0.202 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 0.331 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 0.43 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 0.854 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 0.241 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 0.371 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 0.977 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 0.974 d) e) f) Time 1 2 3 4 5 6 7 8 9 10 11 12 13 r2 0.295 0.087 0.202 0.221 0.13 0.087 0.202 0.319 0.295 0.221 0.13 0.361 0.972 e3 0.854 0.202 0.331 0.2413 0.3709 0.202 0.331 0.4301 0.854 0.2413 0.3709 0.9768 0.9736 c2 1.89492 1.32184 0.63861 0.09186 1.85308 1.32184 0.63861 0.34828 1.89492 0.09186 1.85308 1.70582 0.00165 g) i2 1 e3 i2 1 r2 a ) b ) c ) d ) e ) f ) g ) 47 Table 5.6 The CDD method tables for d = 3 5.2.3 Using the change of distance with time method (CDT) As we said in Chapter 4, the CDT method can also be used to estimate . This method depends on the idea that false nearest neighbor distances grow significantly as time increases. As in the previous methods, we begin by computing the nearest neighbors for . Then we track how the distance between the original nearest neighbors increases with time. This method can be applied by computing the distance in dimension , rather than in dimension . This has the advantage that the measured distance won’t be affected by projection artifacts. Next, we need to have a method for determining whether or not the distances between two vectors can be considered large. For our purpose, we will use the average vector Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 i3 r3 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 0.325 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 0.202 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 0.331 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 0.241 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 0.325 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 0.202 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 0.331 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 0.43 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 0.854 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 0.241 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 0.371 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 0.977 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 0.97 a) b) c) Q4 1 2 3 4 5 6 7 8 9 10 11 12 i3 e4 1 NaN 2.4318 1.6489 2.5414 0.3725 2.3668 1.8663 1.7995 1.0368 2.8685 1.2434 2.433 5 0.3725 2 2.4318 NaN 2.5327 1.4923 2.5528 0.331 2.3738 1.3529 2.0868 1.6802 2.9129 1.1956 6 0.331 3 1.6489 2.5327 NaN 2.5954 1.5612 2.4653 0.4389 1.7247 1.1757 3.0134 2.2894 2.8193 7 0.4389 4 2.5414 1.4923 2.5954 NaN 2.5469 1.7832 2.3182 0.9872 2.3925 0.4232 2.6017 2.4496 10 0.4232 5 0.3725 2.5528 1.5612 2.5469 NaN 2.5035 1.757 1.8512 1.2387 2.8763 0.9807 2.5428 1 0.3725 6 2.3668 0.331 2.4653 1.7832 2.5035 NaN 2.3479 1.5024 1.9941 1.9947 2.9327 0.9783 2 0.331 7 1.8663 2.3738 0.4389 2.3182 1.757 2.3479 NaN 1.5296 1.3966 2.7334 2.36 2.7782 3 0.4389 8 1.7995 1.3529 1.7247 0.9872 1.8512 1.5024 1.5296 NaN 1.4334 1.3952 2.245 2.2062 4 0.9872 9 1.0368 2.0868 1.1757 2.3925 1.2387 1.9941 1.3966 1.4334 NaN 2.7782 2.1478 2.3707 1 1.0368 10 2.8685 1.6802 3.0134 0.4232 2.8763 1.9947 2.7334 1.3952 2.7782 NaN 2.836 2.6148 4 0.4232 11 1.2434 2.9129 2.2894 2.6017 0.9807 2.9327 2.36 2.245 2.1478 2.836 NaN 2.8555 5 0.9807 12 2.433 1.1956 2.8193 2.4496 2.5428 0.9783 2.7782 2.2062 2.3707 2.6148 2.8555 NaN 6 0.9783 d) e) f) Time d FNN 1 2 3 4 5 6 7 8 9 10 11 12 1 5 r3 0.3249 0.202 0.331 0.2413 0.3249 0.202 0.331 0.4301 0.854 0.2413 0.3709 0.9768 2 0 e4 0.3725 0.331 0.4389 0.4232 0.3725 0.331 0.4389 0.9872 1.0368 0.4232 0.9807 0.9783 3 0 c2 0.14651 0.63861 0.32598 0.75383 0.14651 0.63861 0.32598 1.29528 0.21405 0.75383 1.64411 0.00154 h) g) i3 1 e4 i3 r3 1 a ) b ) c ) d ) e ) f ) g ) h ) dL d = 1 dE d 48 length as a bench mark. If the distance between two vectors is larger than the average vector length, we will consider the distance to be large. The average vector length is defined as . (5.7) If we let be the nearest neighbor of , we can check if is a FNN by measuring the distance between and as time increases ( ). For example, we can see from the row of Table 5.7a that the nearest neighbor of is . If we look at Table 5.7b, we can see that the distance between and is (see the circled value). But, after one time step ahead, we can see that the distance between and is which is greater than as indicated by the arrow. We conclude that is a FNN since the distance between and has grown more than as time increases. A second example of the FNNs can be seen from the first row of Table 5.7a. We can see that the nearest neighbor of is , while in Table 5.7b, the distance between and is that means is a FNN. On the other hand, the nearest neighbor of is , as seen in the second row of the last column in Table 5.7a, while if we look at Table 5.7b, we can see that the distance between and is . If we move one step ahead, we can see that the distance between and is which is still less than . A conclusion that can be reached is that is a true nearest neighbor of . Notice here that we choose one as the number of forward time steps that is used to check for FNN. Table 5.9b summarizes the number of FNNs found for through 3. It shows that at , the number of FNNs found is 4, while at , the number of FNNs is 0 β 1 M  yd(k) k = 1 M = Σ yd(n) yd(m) yd(n) y3(m) y3(n) dE = 3 8th y1(8) y1(9) y3(8) y3(9) q3(8, 9) = 0.8638 < β = 1.3404 y3(9) y3(10) q3(9, 10) = 1.4332 β = 1.3404 y1(9) y3(8) y3(9) β y1(1) y1(8) y3(1) y3(8) 1.7115 > β y1(8) y1(2) y1(6) y3(2) y3(6) 0.202 < β y3(3) y3(7) 0.331 β y1(6) y1(2) d = 1 d = 1 d = 2 49 and remains 0 at . That means the minimum embedding dimension is . Finally, notice that the nearest neighbor of shown in the last row of the Tables 5.7b, 5.8b, and 5.9a are labeled as not decidable (nd). The reason for this is that at the first instance of time, the distance between the two vectors is less than , but we do not have enough data to check if the distance between the two vectors after one time step ahead is greater than . So, we have to discard this neighbor from the count of the FNNs and label it as not decidable (nd). Table 5.7 The CDT method for d = 1 d = 3 dL = 2 yd(13) β β The threshold is Q1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 NaN 0.65 1.2015 0.1561 0.313 0.6479 1.2886 0.0262 0.0508 0.3597 0.4 0.7439 1.6362 0.9341 0.1046 8 2 0.65 NaN 1.8515 0.4939 0.963 0.0021 1.9386 0.6762 0.7008 0.2903 1.05 0.0939 2.2862 1.5841 0.7546 6 3 1.2015 1.8515 NaN 1.3576 0.8885 1.8494 0.0871 1.1753 1.1507 1.5612 0.8015 1.9454 0.4347 0.2674 1.0969 7 4 0.1561 0.4939 1.3576 NaN 0.4691 0.4918 1.4446 0.1823 0.2069 0.2036 0.556 0.5878 1.7922 1.0901 0.2607 1 5 0.313 0.963 0.8885 0.4691 NaN 0.9609 0.9755 0.2868 0.2622 0.6727 0.0869 1.0569 1.3231 0.621 0.2084 11 6 0.6479 0.0021 1.8494 0.4918 0.9609 NaN 1.9364 0.6741 0.6987 0.2882 1.0478 0.096 2.284 1.5819 0.7525 2 7 1.2886 1.9386 0.0871 1.4446 0.9755 1.9364 NaN 1.2624 1.2377 1.6483 0.8886 2.0325 0.3476 0.3545 1.1839 3 8 0.0262 0.6762 1.1753 0.1823 0.2868 0.6741 1.2624 NaN 0.0246 0.3859 0.3737 0.7701 1.61 0.9078 0.0784 9 9 0.0508 0.7008 1.1507 0.2069 0.2622 0.6987 1.2377 0.0246 NaN 0.4105 0.3491 0.7947 1.5853 0.8832 0.0538 8 10 0.3597 0.2903 1.5612 0.2036 0.6727 0.2882 1.6483 0.3859 0.4105 NaN 0.7596 0.3842 1.9959 1.2937 0.4643 4 11 0.4 1.05 0.8015 0.556 0.0869 1.0478 0.8886 0.3737 0.3491 0.7596 NaN 1.1438 1.2362 0.5341 0.2953 5 12 0.7439 0.0939 1.9454 0.5878 1.0569 0.096 2.0325 0.7701 0.7947 0.3842 1.1438 NaN 2.3801 1.6779 0.8485 2 13 1.6362 2.2862 0.4347 1.7922 1.3231 2.284 0.3476 1.61 1.5853 1.9959 1.2362 2.3801 NaN 0.7021 1.5316 7 14 0.9341 1.5841 0.2674 1.0901 0.621 1.5819 0.3545 0.9078 0.8832 1.2937 0.5341 1.6779 0.7021 NaN 0.8294 3 15 0.1046 0.7546 1.0969 0.2607 0.2084 0.7525 1.1839 0.0784 0.0538 0.4643 0.2953 0.8485 1.5316 0.8294 NaN 9 a) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 8 FNN 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 1 FNN 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 9 FNN 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 8 FNN 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 2 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd b) At d =1, the number of FNN = 4 nd: not decidable β = 1.3404 i1 1 i1 1 a ) b ) 50 Table 5.8 The CDT method for d = 2 Table 5.9 The CDT method for d = 3 The circled distances in the above tables are the distances between the reference points and their nearest neighbors at the first instance of time in . The arrows represent the direction where the reference points and their nearest neighbors move after one time step. Q2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 NaN 1.9623 1.2991 0.9756 0.313 2.044 1.4552 0.7013 0.2947 1.1099 0.4108 2.4041 2.2773 1.2008 9 2 1.9623 NaN 2.2959 1.0165 2.0851 0.0871 2.267 1.3346 1.7113 0.8525 2.2106 0.4447 2.3018 1.9268 6 3 1.2991 2.2959 NaN 1.4363 1.0155 2.3467 0.202 1.1934 1.1685 1.6572 0.994 2.6451 1.1736 0.3735 7 4 0.9756 1.0165 1.4363 NaN 1.0693 1.0925 1.4728 0.3193 0.7038 0.2214 1.1942 1.4478 1.8968 1.1099 10 5 0.313 2.0851 1.0155 1.0693 NaN 2.1617 1.1858 0.7553 0.3896 1.2452 0.1295 2.5167 2.0623 0.9757 11 6 2.044 0.0871 2.3467 1.0925 2.1617 NaN 2.3116 1.4094 1.7902 0.9342 2.2867 0.3606 2.3114 1.9759 2 7 1.4552 2.267 0.202 1.4728 1.1858 2.3116 NaN 1.2626 1.2965 1.6901 1.1759 2.5928 0.9721 0.3631 3 8 0.7013 1.3346 1.1934 0.3193 0.7553 1.4094 1.2626 NaN 0.4113 0.5204 0.8782 1.7625 1.8363 0.9094 4 9 0.2947 1.7113 1.1685 0.7038 0.3896 1.7902 1.2965 0.4113 NaN 0.8635 0.5191 2.1483 2.0462 0.9978 1 10 1.1099 0.8525 1.6572 0.2214 1.2452 0.9342 1.6901 0.5204 0.8635 NaN 1.3731 1.2945 2.0661 1.327 4 11 0.4108 2.2106 0.994 1.1942 0.1295 2.2867 1.1759 0.8782 0.5191 1.3731 NaN 2.6407 2.0842 1.0026 5 12 2.4041 0.4447 2.6451 1.4478 2.5167 0.3606 2.5928 1.7625 2.1483 1.2945 2.6407 NaN 2.4815 2.2718 6 13 2.2773 2.3018 1.1736 1.8968 2.0623 2.3114 0.9721 1.8363 2.0462 2.0661 2.0842 2.4815 NaN 1.0867 7 14 1.2008 1.9268 0.3735 1.1099 0.9757 1.9759 0.3631 0.9094 0.9978 1.327 1.0026 2.2718 1.0867 NaN 7 a) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 9 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 11 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd b) At d =2, the number of FNN is 0 i2 1 i2 1 a ) b ) Q3 1 2 3 4 5 6 7 8 9 10 11 12 13 1 NaN 2.3861 1.5738 2.0909 0.3249 2.3578 1.8552 1.7115 0.854 2.2397 0.5981 2.419 2.5277 5 2 2.3861 NaN 2.3433 1.1293 2.5366 0.202 2.2764 1.3501 1.7993 1.0355 2.8459 1.1773 2.3165 6 3 1.5738 2.3433 NaN 1.7281 1.4082 2.3642 0.331 1.3699 1.1718 1.9656 1.6549 2.717 1.1919 7 4 2.0909 1.1293 1.7281 NaN 2.212 1.2837 1.6302 0.4301 1.2623 0.2413 2.5774 2.1445 2.0406 10 5 0.3249 2.5366 1.4082 2.212 NaN 2.5033 1.7141 1.8131 0.9703 2.3836 0.3709 2.5416 2.378 1 6 2.3578 0.202 2.3642 1.2837 2.5033 NaN 2.3117 1.4613 1.8288 1.2107 2.7966 0.9768 2.3127 2 7 1.8552 2.2764 0.331 1.6302 1.7141 2.3117 NaN 1.3277 1.3427 1.8676 1.9738 2.7391 0.9736 3 8 1.7115 1.3501 1.3699 0.4301 1.8131 1.4613 1.3277 NaN 0.8638 0.6468 2.1805 2.1863 1.8941 4 9 0.854 1.7993 1.1718 1.2623 0.9703 1.8288 1.3427 0.8638 NaN 1.4332 1.3408 2.2137 2.0674 1 10 2.2397 1.0355 1.9656 0.2413 2.3836 1.2107 1.8676 0.6468 1.4332 NaN 2.7477 2.1193 2.2335 4 11 0.5981 2.8459 1.6549 2.5774 0.3709 2.7966 1.9738 2.1805 1.3408 2.7477 NaN 2.7324 2.5864 5 12 2.419 1.1773 2.717 2.1445 2.5416 0.9768 2.7391 2.1863 2.2137 2.1193 2.7324 NaN 2.6164 6 13 2.5277 2.3165 1.1919 2.0406 2.378 2.3127 0.9736 1.8941 2.0674 2.2335 2.5864 2.6164 NaN 7 nd a) At d =3, the number of FNNs is 0 also. d FNNs 1 4 2 0 3 0 b) i3 1 a ) b ) ℜ3 51 5.3 The Predictive technique Now that we have discussed the three geometric methods (CND, CDD, and CDT), we will now talk about the predictive technique. We said in Chapter 4 that the predictive technique can also be used to estimate by approximating the function in Equation (4.6). We also said that we will approximate by using a neural network. Since chaotic systems are nonlinear, we need to use a nonlinear network to make the approximation. To estimate of the Henon map by using the predictive technique, we will use a nonlinear neural network which consists of two layers. The first layer has 2 neurons with sigmoid transfer functions and a Tapped Delay Line (TDL) connected to it. The second layer has one neuron with a linear transfer function. The TDL is fed with the measurements to produce the delayvectors . The network structure is shown in Figure 5.1. While the number of taps in the TDL (d) is changed from 1 to 2 to 3, the network is trained to predict the current measurement . After the end of the training process, the sum of the squares of the prediction errors (SSE) is recorded as a function of d. We applied this technique to 100 points from the coordinate of the Henon map. As we can see from Table 5.10, the SSE changes significantly as d increases from 1 to 2 and remains almost the same at . We conclude that the network has accurately approximated and that . The logplot of the SSE versus d is shown in Figure 5.2. dL μ μ dL y(m) yd(m – 1) y(m) X1 d = 3 μ dL = 2 52 Figure 5.1 The nonlinear network used to estimate of the Henon map Table 5.10 The neural network SSE as a function of d Figure 5.2 LogPlot of the SSE versus d for the Henon map. 5.4 Chapter summary In this chapter, we demonstrated the estimation of the minimum embedding dimension of the Henon map by using the geometric and the predictive techniques. The purpose of this chapter was to provide some insight into the operation of the algorithms. In the next chapter, we will present the complete algorithms in full detail and discuss practical issues in using the algorithms on more complex systems than the Henon map. dx1 1 1x1 y(m) n1 (m) a1 (m) 1 b1 2x1 2xd b2 LW2, 1 n2 (m) TDL tansig purelin IW1, 1 1x2 a2 (m) dL d SSE(d ) 1 4.0321 2 4.7117x107 3 4.0824x107 1 2 3 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 SSE of the prediction error 53 CHAPTER 6 ADVANCED ALGORITHMS FOR ESTIMATING THE MINIMUM EMBEDDING DIMENSION 6.1 Introduction In Chapter 3, we explained that two parameters are needed in order to apply the embedding theorem to model a chaotic system. These two parameters are the dimension (d) of the delayvector ( ) and the delaytime (T) between the delayvector coordinates. We presented in Chapter 4 three geometric methods (CND, CDD, and CDT) and one predictive method through which we estimate the minimum dimension ( ) required to embed the system’s attractor. More details were given in Chapter 5 regarding the application of the four methods to estimate of the Henon map. The purpose of this chapter is to give full detail of six algorithms which use the four methods mentioned above to estimate . We will also show a method used to find the delaytime T. Before that, we first give a summary of the four methods. The CND, CDD, and the CDT geometric methods depend on the idea that if the dimension of the space is not large enough to represent the attractor of the system, projection artifacts will appear in the projected attractor. These artifacts cause points on the attractor to be falsely projected close to each other and produce False Nearest Neighbors (FNN) (see Section 4.3 for more detail). The dimension can be estimated as the minimum dimenyd dL dL dL dL 54 sion where the percentage of the FNNs does not change significantly with further increase in dimension. The CND method detects the existence of FNNs by checking to see if the nearest neighbors in the space of dimension d remain neighbors in dimension . On the other hand, the CDD method detects the existence of FNNs by checking to see if the distance between the nearest neighbors in dimension d will increase significantly as the dimension increases to . For the case of the CDT method, detection of the existence of FNNs is done by checking to see if the distance between the nearest neighbors in dimension d will change significantly as time increases. On the other hand, in the predictive method, the estimation of is done by approximating the function that operates on the reconstructed attractor. is approximated by using a neural network with a Tapped Delay Line (TDL) connected to its input. As the number of taps in the TDL (d) increases, the prediction error decreases. At one point, further increase of d does not improve the prediction error. At this point, is found. In the remaining parts of this chapter, we present in Section 6.2 a method used to find the delaytime T. In Section 6.3, we present two different algorithms based on the CDD and the CDT methods which were proposed by Abarbanel et al to estimated . In Section 6.4, we discuss some limitations of the previous two algorithms and suggest four new algorithms to overcome these limitations. In Section 6.5, we present three of the four algorithms, which are based on the geometric technique. The fourth algorithm, which is based on the predictive techniques, is presented in Section 6.6. At the end, Section 6.7 concludes with a chapter summary. d + 1 d + 1 dL μ: ℜd ℜ1→ μ dL dL 55 6.2 The delaytime (T) We explained in Section 4.2 that we wanted to find the values of the parameters d and T such that the system we developed 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


