VERIFYING C O M P U T E R SOLUTIONS TO DISCRETE-TIME D Y N A M I C A L SYSTEMS by DAVID JOHN U R M I N S K Y B.Sc. McMaster University, 1997 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required staadard T H E UNIVERSITY OF BRITISH C O L U M B I A April 2000 © David John Urminsky, 2000 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Mathematics The University of British Columbia Vancouver, Canada Abstract Chaotic dynamical systems exhibit sensitive dependence on ini t ia l conditions. Round-off errors introduced in computer simulations may cause a computed orbit to quickly diverge from the true orbit. One may therefore question the validity of these computations. Shadowing provides a means for studying the validity of a computation. If we are able to show that a true solution with a different ini t ial condition, called a shadowing orbit, stays close to a computed solution, which we call a pseudo-orbit, then perhaps the computed solution has some validity. We wi l l investigate two methods for proving the existence of a shadowing orbit. The first method comes from two theorems by Chow and Palmer. These theorems provide us with sufficient conditions for when a pseudo-orbit is shadowed by a true orbit and provide us with a bound on the shadowing distance of a true orbit. The second method is by Grebogi, Hammel, Yorke and Sauer. Their method contains a pseudo-orbit in a carefully constructed sequence of parallelograms which helps to prove the existence of a shadowing orbit. These two methods wi l l then be used to prove the existence of a shadowing orbit for examples in one and two dimensions. F ina l ly we wi l l discuss a refinement technique which takes a 'noisy' orbit and produces a less 'noisy' orbit which shadows the original orbit. This technique w i l l then be applied to some examples to find a numerical shadow for a pseudo-orbit. i i Table of Contents Abst rac t i i Table of Contents i i i Acknowledgements iv Chapter 1. In t roduct ion 1 1.1 Numerical solutions and shadowing 1 1.2 The basic idea 2 1.3 Is shadowing always possible? 6 1.4 Suggestions for measures of output error 7 1.5 Why shadowing? 8 1.6 Summary of this thesis 9 Chapter 2. Shadowing 11 2.1 Introduction 11 2.2 The one-dimensional case 12 2.2.1 Existence of a shadowing orbit 12 2.2.2 The computation 17 2.2.3 Existence of a shadowing orbit: another method 23 2.2.4 Discussion of one-dimensional shadowing 26 2.3 The two-dimensional case 28 2.3.1 Two-dimensional shadowing: Chow and Palmer's method 29 2.3.2 The computation 36 2.3.3 Two-dimensional shadowing: the GHYS method 43 2.3.4 The computation 44 2.3.5 Discussion of two-dimensional shadowing 47 2.4 Higher-dimensional shadowing 48 Chapter 3. Refinement 50 3.1 Introduction 50 3.2 The G H Y S refinement algorithm 52 3.2.1 One dimensional case 52 3.2.2 The two-dimensional case 55 3.2.3 Generalizing the GHYS algorithm for Hamiltonian systems 58 3.2.4 Discussion of the GHYS algorithm 60 Chapter 4. Conclusions and future work 62 Bib l iography 65 ii i Acknowledgements First and foremost, I would like to thank my supervisor, Dr. Wayne Nagata, for all of his help in putting this thesis together. He has shown me how to write in the language of mathematics. Secondly, I would like to thank Dr. Brian Wetton for his suggestions and comments and for getting my thesis back to me quickly. Special thanks also goes out to the Institute of Applied Mathematics, for the use of their computer facilities. Especially after several unfortunate computer crashes some of my programs caused. Thanks also goes to Bridget Gilbride and Valerie Dion for their many insightful conversations. Last of all, I would like to thank Jonathan Patrick and Marc Mailhot (my house mates). When my thesis got the better of me, they were there to wrestle me back. iv Chapter 1 Introduction There is something that is nothing, but it has a name. It joins our walks; it joins our talks; it plays in every game. What is it? A shadow. -Traditional 1 . 1 N u m e r i c a l s o l u t i o n s a n d s h a d o w i n g Numerical simulations are often used to investigate the long term behavior of solutions of ordinary or partial differential and discrete time dynamical systems. A basic requirement needed to interpret the result of a simulation is that the behavior of a solution determined by the computation is the same as the behavior of some true solution of the system. For many dynamical systems, a numerical solution afflicted with small truncation and roundoff errors will diverge away from the true solution of the system over a small interval of time. If there is a difference in behavior between a computed solution and the actual solution, then the investigator must question the validity of his solution. The problem of computing "accurate" solutions to a dynamical system is made even harder when the system is chaotic. In this case, orbits exhibit sensitive dependence on initial conditions: two orbits with initial conditions which are close together will tend to diverge away from each other exponentially. Because of this behavior, any numerical truncation or rounding error made at any step in the computation will be greatly mag-1 nified in the future evolution of the system. In view of this, it is natural to ask: "under which conditions will a computed orbit be close to a true orbit of the model? Perhaps a better question to ask is: can we trust the computed orbit to represent faithfully the dynamics of the system? There has been much work done to answer this question. Several computational techniques for verifying computer-generated orbits have been developed for low dimen-sional dynamical systems. Verifying a computer generated solution means we produce a computer-assisted proof of the existence of a true orbit of the system which stays close to the computed solution. This true orbit we call a shadowing orbit. The distance between the computer generated orbit and the shadowing orbit remains small over a long period of time. Now we can give an answer to the question asked above. It will be the view of this thesis that we can trust a computed orbit to represent the true dynamics of a system if we can show that there exists a true orbit that shadows the computed orbit for a long period of time. 1 . 2 T h e b a s i c i d e a The study of discrete-time dynamical systems has flourished over the last half a century. During the same time frame, we have also seen many advances in computer technology. As a result, computers have become a main tool in studying dynamical systems. Given the map xi+l = f{%i) with an initial condition Xo, we ask the computer to calculate the orbit {xz}^=0. In general, the computer is not able to calculate the orbit exactly. The computer will only keep so many digits of accuracy and makes a round-off error with every calculation. This rounding error the computer makes is referred to as noise. Instead of calculating {xi}^=0, the computer actually calculates {pi}^ which we call a pseudo-orbit. 2 Let us consider the quadratic map, xi+i = f(xi) = \Xi(l - x^, (1.1) for 0 ^ i ^ N = 38, A = 3.8 and x0 = 0.6. The maximum error the computer makes when it performs any of the operations of addition, multiplication, division and subtracting is 2~d. The exponent d is the number of binary digits in the floating point precision. Following the ideas by Wilkinson in [10], when the computer tries to calculate x\ we find that the computer actually calculates Pi = Xx0(l - x0)(l + e u ) ( l + € 2 i ) ( l + e3i)- (1-2) where |e,;i| ^ 2~d. The factors (1 + en),( l + e 2 i) and (1 + tzi) arise from the error when calculating \xQ, 1 — x0 and (Ax 0 )( l — x0). Let pi = (1 + ei;)(l + £ 2 i ) ( l + £31) be the floating point error from calculating from pi. The maximum value .x,; can obtain for A = 3.8 is So large errors resulting from the subtraction of two numbers which are close together is avoided because the maximum value xi can take is bounded well away from 1. For A values close to 4, we would obtain large errors calculating (1 — Xj) when Xi is close to 1. Let us examine how the error accumulates at each time step. We can write the error at the % + 1 time step as, ei+l = — Pi+1 = Xl+l - Pif{pi). Expanding f(x) about Xi using Taylor's Theorem, we can show that el+i = xi+i - Pi[f(xi) + f'(xi)(pi - xz) + f ^l\pi - x^2 + 0((pi - x^3) . = [xi+1 - P;.[{:>:,)] + Pif'{x^e, + 0(e2) = ( l - p i ) / ( x i ) + p i/'(a; i)e i+.0(e?). (1.3) The term is expected to be small, so 0(ej) is generally small compared to f'(xi)ei. Also, |(1 — Pi)f(xi)\ is bounded by (1 + 2 _ d ) 3 — 1, and |/'(.x)j is on average greater than 3 one o n the i n t e rva l [0,1]. So \pif'(x)ei\ can grow large as i increases, a n d leads to an error w h i c h t y p i c a l l y grows exponent ia l ly , at least i n i t i a l l y . F i g u r e 1.1 represents the d is tance between two c o m p u t e d orbi ts for equa t ion (1.1) where A = 3.8. B o t h o rb i t s have the same i n i t i a l c o n d i t i o n 0.6, bu t one o rb i t is ca l cu la t ed us ing doub le p rec i s ion a n d the o ther us ing s ingle prec is ion . N o t i c e how the d is tance between the orb i t s is inc reas ing a p p r o x i m a t e l y l ike an exponen t i a l f unc t i on w h i c h agrees our resul ts f r o m (1.3). A f t e r 37 t i m e steps, the d is tance between the two c o m p u t e d orb i t s becomes as large as the a c t u a l values of the orb i t s . Distance between a double and single precision orbit F i g u r e 1.1: Th i s graph represents the distance between two orbits for (1.1), one is computed in single precision and the other is computed i n double precision D e s p i t e the p o s s i b i l i t y tha t our n u m e r i c a l so lu t ion m a y diverge e x p o n e n t i a l l y f rom the a c t u a l so lu t i on , we can s t i l l pu t the n u m e r i c a l i n f o r m a t i o n to g o o d use. C o n s i d e r a s o l u t i o n Xi of equa t ion (1.1) where we take xN = pN. O u r s o l u t i o n Xi is c a l cu l a t ed us ing a h igher p rec i s ion t h a n W e i terate (1.1) backwards , to o b t a i n a new o rb i t Xi for 0 ^ i ^ N. O u r new orb i t satisfies (1.1) except w i t h a different i n i t i a l c o n d i t i o n 4 x0 = 0.60000002552983. This new computed orbit satisfies \xi-pi\ ^ 2.9680201e-07, for 0 ^ ^ 37. Using the same method we can find a computed orbit which stays within 9.75213521e-05 for 0 ^ i ^ 100, 000. These orbits stay close to the original orbit because they were calculated backwards in time. In general, the orbits for (1.1) with A = 3.8 diverge exponentially forwards in time. If we iterate (1.1) backwards in time we can keep the method stable and orbits tend to^stay close together. Let us consider the map xl+x — kxi, (1.4) for k > 0, 0 ^ i ^ JV, and x0 = 1. It can be shown using a similar technique as above that the error at time step i + 1 is bounded by | e l + 1 K 2" d|A;x| ^ | ( 1 + 2~d)k\^j , where x = max Xi. If 1(1 + 2 d)k\ < 1, then the sum above is a geometric series which is bounded by 1_\r1+2-d)k\ a n i ^ x w ^ be bounded by x0. So the error will be bounded and we can compute the maximum shadowing distance. For |(1 + 2~d)/c| > 1, the series diverges and x will be large for large values of N. So for |(1 + 2~d)k\ > 1, the error between the true orbit and the computed orbit may grow exponentially. Now let us look at the system of equations given by u i + i = Dm, O ^ i ^ N , (1.5) where D = diag(/c 1 ; k2), = [xi,yi], \ki\ < 1 and \k2\ > 1. Errors to this system will be damped in the x direction and will grow in the y direction. But we can find a more accurate solution by prescribing x0 and yN and solving for the resulting pseudo-orbit. For more general linear systems u i + 1 = Aui, wi th A being a d x d matr ix wi th no eigenvalues on the imaginary axis, the normal modes are either exponentially growing or exponentially decreasing. Shadowing of a numerical solution is possible by prescribing the ini t ia l value of the decreasing modes and the final value of the growing modes. 1 . 3 I s s h a d o w i n g a l w a y s p o s s i b l e ? In Section 1.2 we demonstrated shadowing for a simple one-dimensional map. B y spec-ifying the finial condition for the shadowing orbit we produced a less noisy orbit which shadowed the original orbit. Can we safely say that shadowing wi l l always work? The answer is no. This is especially true in chaotic systems. After some time it may not be possible to find a true orbit that remains close to a computed noisy orbit. Let us consider the quadratic map (1.1) with, A = 4 and x0 = 0.4999, using single precision floating point arithmetic where the computer only uses 8 digits of accuracy. A t the first time step, the computer maps the pseudo-orbit outside the interval [0,1]. Once this happens, the pseudo-orbit cannot be shadowed since any true orbit should always be contained in the interval [0,1]. The occurrence of such a phenomenon is called a glitch. Glitches can cause numerical errors which can introduce new behavior into the system. York and Sauer (1991) suggest that a similar phenomenon occurs in higher-dimensional systems because of the folds created by the tangencies of homoclinic orbits of stable and unstable manifolds. Shadowing was originally discussed for hyperbolic dynamical systems. Systems clas-sified as hyperbolic have the property that each point in an orbit has a stable manifold and an unstable manifold. Infinitesimal displacements in the unstable direction grow exponentially when followed forward in time and infinitesimal displacement in the stable direction grow exponentially as time goes backwards. The idea is similar to what we discussed for equation (1.5). Here, the unstable direction is given by the vector [0,1] 6 and the stable direction is given by [1,0]. To be hyperbolic, it is required that the angle between the stable and unstable direction is bounded away from zero. If these conditions are met, then it is possible to show (Anosov, 1967, Bowen, 1975) that a true orbit can be found near the noisy numerical orbit for arbitrarily long times. Unfortunately many dynamical systems are not hyperbolic so the theorem does not apply. Surprisingly, we are still able to show that non-hyperbolic systems can be shadowed for long periods of times. We will see cases of this later on in this thesis. 1.4 Suggestions for measures of output error It is a good idea to distinguish between the desired properties of a computation and the deviations that computations make from those properties, i.e., the output errors they make. 1. We could demand that the computation precisely follow the exact evolution of the true solution to the system. This would be possible for discrete-time dynamical systems but not for ODE's. For iterated maps we can use programs such as Maple that can use exact arithmetic to calculate the orbits. However, this is usually infeasible since each iteration typically requires more memory than the previous one. For ODE's there is no numerical integration routine which will give exact solutions to arbitrary ODE's. 2. We could demand that the phase space distance between the computed solution and the true solution with the same initial condition be bounded by a small constant. This condition is currently infeasible. O D E integrations have truncation and round-off errors. For many dynamical systems the magnitude of these errors will be magnified on a short time scale. Similarly, for iterated maps round-off errors can be magnified making our calculations questionable. 7 3. Somewhere in between 1 and 2 is that we could demand, with some small error, that the simulation follow an the exact phase space solution with a slightly different initial condition. The idea of shadowing is essentially this condition, and most of the work of this thesis will be to study this situation. 1 . 5 Why shadowing? Now that we have an introduction to the idea of shadowing, we can now look at the idea from a wider perspective. In particular I will answer some questions that people have asked me about the subject. 1. Q Why should we care about finding an exact solution for a specific initial condi-tion? A We actually don't care about finding an exact solution. We just want to know that a true solution follows close to our computed solution. Finding the true solu-tion is usually not possible for most dynamical systems. Once we know that a true solution exists, we can show that a cheaply calculated (low-precision) solution has a more expensively calculated (high-precision) solution close by. If we cannot prove such an expensive solution exists, then we may want to suspect the reliability of the cheap solution. 2. Q Why not just use higher precision arithmetic to calculate the solution to the system? A This seems like a good idea. However, if the system exhibits chaotic behavior eventually any round-off error we make will cause the solution to diverge exponen-tially from the true solution. Also, if we use too many digits of accuracy, we may run into problems storing the numbers in the computer. 8 3. Q Computers almost always create round-off errors in their calculations. Does this mean that we cannot use computer to study these types of systems? A This question wi l l hopefully be answered to some degree in this thesis. C o m -puters are the best tools which we have to solve these types of systems quickly and efficiently. Al though they may not be able to provide us with an exact solution to a dynamical system, they can sti l l provide us with many insights to the system. So i f we are able to show that there is some true solution which stays close to our computed solution, then perhaps our computer solution has some merit. 4. Q What is the goal of shadowing? A Firs t of al l we are answering question 3 above. We want to verify the reliability of our computer simulations. As discussed earlier, places where our computer solutions may be suspicious are at glitches. We want to be able to identify glitches in our simulations, learn about them, and perhaps try to eliminate them. Also, we wish to show that we can use low precision arithmetic for our calculations and our solutions represent some kind of true behavior of the system. 1.6 Summary of this thesis In Chapter 2 we wi l l define what a shadowing orbit is and define what it means for a pseudo-orbit to be shadowed by a true orbit. In particular we wi l l investigate two theorems by Chow and Palmer given in [3] and [4]. These theorems provide us with sufficient conditions for when a pseudo-orbit is shadowed by a true orbit and provide us wi th a bound on the shadowing distance of a true orbit. We wi l l also look at algorithms for these theorems which carefully calculate the necessary quantities for proving the two theorems. Also in Chapter 2, we wi l l study methods of containing a pseudo-orbit presented in [6] and [7] by Grebogi, Ffammel, Yorke and Sauer. These methods carefully contain a given pseudo-orbit in intervals or parallelograms. From these intervals and 9 parallelograms we can determine bounds on the shadowing distance of a true orbit. Being able to determine whether a true orbit shadows a pseudo-orbit is important for computer simulations of chaotic systems. So how can we compute an orbit that is closer to this true orbit? In chapter 3 we discuss a noise reduction technique by Grebogi, Hammel, Yorke and Sauer in [7] which reduces the dynamical noise of the system. Appl ica t ion of the technique is applied to an one- and a two-dimensional system. Final ly, in Chapter 4 we discuss some open questions and further avenues of research. 10 Chapter 2 Shadowing "Time flies over us, but leaves its shadow behind." - N . Hawthorne 2.1 Introduction W i t h a lit t le understanding of the problems faced when running a computer simulation of a dynamical system, we now make some quantitative observations about the solutions we compute. Most importantly we are interested in the existence of a true orbit which stays close to the computed orbit. Once the existence of such an orbit can be established, the maximal distance between the true and computed orbit can be determined. Several theorems have been proven which measure certain quantities of the computed orbit and give bounds on the maximal shadowing distance. I wi l l outline some of these theorems and provide some proofs. A most useful method of proving a true orbit exists and stays close to a computed orbit is called containment. Containment is a rigorous method to prove the existence of a true orbit. In this chapter I wi l l outline a method by Yorke, Hammel, Sauer and Grebogi [7] in which a special sequence of parallelograms is constructed which contain both the true orbit and the computed orbit. This method can be extended to systems of three or more dimensions. Other methods for showing that a true orbit shadows a pseudo-orbit 11 have been developed in [3] and [4] by Chow and Palmer. These methods involve rigorous calculations of the expansiveness of the pseudo-orbit. We now introduce some definitions. Some of the following definitions are taken, with some modifications, from Grebogi, Hammel, Yorke, and Sauer [7] (GHYS) Def in i t ion A true orbit {.Xj}-^ 0 for f satisfies xi+i — f{xj) for 0 ^ i < N. Defini t ion {pi}^L0 is 5/-pseudo-orbit (also called a noisy orbit) for f, if\pi+i — f(pi)\ < 5j for 0 ^ i < N, where 5f is the noise amplitude. Defin i t ion For 0 ^ % < N, the one-step error made between step i and step i + 1 of the pseudo-orbit {Pi}"=Q is ei+l = pi+i - f(p{). Thus, a true orbit is a sequence whose one-step errors are identically zero. Def in i t ion The true orbit {xi}^L0 8x-shadows {pi}$L0 on 0 ^ i < N if \xi — ^ 5X for 0 ^ i < N. 2.2 The one-dimensional case Some interesting results can be easily obtained from simple one-dimensional maps. By analyzing one-step errors, we can determine the existence of a shadowing orbit and ac-tually find an upper bound for the shadowing distance. The methods outlined in this section can be extended to higher-dimensional systems. 2.2.1 Existence of a shadowing orbit The first method we will discuss is outlined in Chow and Palmer [3]. Their method is as follows: Let / : [0,1] [0,1] be a C2 function and {pi}^=0 be a pseudo-orbit of this map. 12 We want to find a true orbit {xi}fL0 such that the one step errors \xi — pi\ are small for o < ?; < N - 1. We define the following quantities: N-l a= sup '£\Df(pi)-1Df(pi+1)-1---Df(pj)-1\ 0<i<N-l . 3=i (2-1) and T = SUp 0<i<N-l N-l Y^DfiPiY'Dfip^)-1 • ••i)f(P:i) ; J V l - f(Pj; (2.2) These quantities estimate how large the one-step errors are at each time step along the pseudo-orbit. We can think of these quantities as measures of expansiveness since they will give bounds on how far the true orbit (if it exists) and the pseudo-orbit are from each other. Theorem 2.1 (Chow and Palmer , 1991) Let f : [0,1] -> [0,1] be a C2 function with M = sup{\D2f(x)\ : 0 ^ x ^ 1}. Let {pi}iL0 be a pseudo-orbit of f such that 2Mor < 1. Then there is an true orbit {xi}fL0 with -l ^ sup \xi - pi| ^ 2 r ( l + Vl - 2Mar)~1 2(1 + V l - 2MCTT). In particular, the true orbit 5x-shadows the pseudo-orbit with 5X ^ 2r(l + v / l — 2Mar)~l. Proof: This proof is similar to the one given in [3] with a few corrections and clarifica-tions. The true orbit {xi}fLQ has to satisfy xi+i = f{xi), for 0 ^ i ^ T V - 1. 13 We expect a true orbit, if one exists, to be close to the the pseudo-orbit. So we set Xi = Pi + and find a bound on ej. We note ej satisfies <2j+l — xi+i — pi+i = f{Xi)-Pi+l = f(pi + ei)+pi+1. (2.3) If we let &(e) = f(Pi + e) - ffa) - Dfipje, then equation (2.3) becomes el+l = Df(pi)ei + /fa) +pi+i + cji(ei). (2.4) Now we must show that there exists a sequence {e;}^ 0 which satisfies equation (2.4) such that 2r , , e , K ^ , (2.5) for 0 ^ ?; ^ N. Define S to be the set of all sequences e = {ej}^ 0 such that |ej| ^ e for 0 ^ i ^ N. Now we will define the mapping T by N-l (Te)i = -J^Dfip^Dfip^y1 • ••Df(pj)-1\f(pj)-pi+1+gj(ej)], (2.6) j—i for 0 ^ i ^ N — 1. [The reason for choosing T in this way is to measure by how much the backwards errors are shrinking. For example, suppose we wanted to find a true orbit for which eN = 0 at time step N. Then we iterate (2.4) backwards and obtain (2.6).] In fact, T is a continuous mapping of S into itself. To show this, observe that if we approximate /(j/i + ej) using a linear approximation we can use Taylor's theorem to show 14 that \gi(e)\ ^ \M\e\2 for all i, and e. Now we show that Te is bounded for e e S. \N-1 \(Te)i\ < J2Df(Pir1Df(Pl+1)-1 • • • Df(Pj)-l[f(Pj) - pj+1] N-l •YilDfipi^Df^)-1 • • • Df{p3)~l\ • l-M\e3 3=1 ^ T + -Mae2 = e, 2 (2.7) for 0 ^ i ^ N. Hence T maps 5 into itself. So from Brouwer's fixed point theorem we know that there must be a fixed point, The sequence {e{}?L0 is a solution to (2.4). By iterating (2.4) backwards we can show that, N-l et = -J^Dfip^Dfip^)-1 • • • Df(p3)-l[,f(Pj) - P j + 1 + 9]{e3)l j=i for i = 0,1, . . . , iV — 1. Using this fact, it follows that N-l Y;Df(pi)-1Df(Pi+l)-1 • ••Df(Pj)-l{f(P])-Pj+1} j=i N-l = -e{ - J2Df(Pi)-lDf(Pl+l)-1 • ••Df(Pj)-1gJ(e]). (2.8) 3=1 Now we define sup \e[ \ ^ e, 0<i<N-l (2.9) 15 and by using (2.8) and (2.9) we can bound r as follows: T ^ \e{\ + N-l 3=1 ^ He'll + ^Ma\\ef\\2 ^ \\ef\\ + ]-Mae\\ef\\ * ( 1 + ^ ) l l e / H 1 + (since 2MOT < 1) 2(l + v / l - 2 M a r ) J We solve for ||e'|| to get: e'll > 1 - i 1 2(1 + V I - 2MOT)\ T. (2.10) Combining (2.5) and (2.10) we get, - l 1 I 1 + ^ sup \xi -pi\ < 2r 1 + y/\ - 2Mar) 0<i<N-l V / 2 ( l + v l - 2M(TT) This completes the proof of the theorem. • To be able to use this theorem we have to calculate a and r. In practice this is a time consuming process. If we are able to get bounds on these quantities then we can still bound the quantity 2Mar. In [3], Chow and Palmer show that a and r can be bounded by, a ^ (1 - / J , s ) "V S , T ^ (l-Ps) -1, (2.11) (2.12) where Ps = sup \Df(pl)-1---Df(Pl+s)-1\, O^i^N-s min(i+s,N) S U P , E \Df(pl)-1Df(Pi+l)-1 j=i min(i+s,N) ] T Df(PlylDf{Pl+l)-1 0<i<N-l sup 0<i<N-l •Df(Pj)-'{pJ+1-f(Pj) 16 and s is a positive integer in the interval [0, N] chosen so that fis < 1. Now, the problem is reduced to the calculation of fj,s, as, and TS. 2.2.2 T h e c o m p u t a t i o n Given a map f(x) and a pseudo-orbit {pi}^L0, we wish to apply Theorem 2.1 to determine i f a true orbit shadows the pseudo-orbit and also how closely the pseudo-orbit is shadowed by the true orbit. The following steps outline a procedure to do this. S t e p 1 First we wish to calculate upper bounds for the quantities: M = sup{\D2f{x)\:Q^x^l} 5 = sup{|./'(.x) -f(x)\ : 0 x ^ 1} A = sup{|D/(.x) - Df(x)\ : 0 ^ x ^ 1} where f(x) and Df(x) are f(x) and Df(x) calculated in double precision. The second two quantities are used when we correct for round-off errors in Step 3. S t e p 2 Compute the pseudo-orbit {pi}f=0 in both double and single precision while simulta-neously computing the quantities sup \Df{pi 0<i<N-s Df(Pl+s) - n 5 = min{i+s,N} S U P , E \Df(Pl)-1Df(pl+l)-1 j=i min{i+s,N} 0 < i < i V - l sup sup \pi+i-f(Pi)\ 0<i<N-l Df{p3)-l\ Df{p3yl[p3+l-f{p3) 17 The computer actually works with the numbers f(x) and Df(x) and tries to calculate the numbers 7JJ, aj, and 5. By doing this, the computer makes round-off errors and actually calculates the quantities u.'s, a'sl r's, and 6'. At this point we will assume that a suitable value for s.has been chosen so that 717 < 1. The number s is found by trial and error and can be a time consuming process. Step 3 Now we need estimates of a and r. In order to do this we first must place a bound on the round off error the computer makes while trying to calculate ~JTS. Following the methods in [10], we see that when the computer multiplies two floating point numbers x\ and X2, the computer actually calculates xix2(l + e) where 1 - 2~d ^ 1 + e 1 + 2'd Similarly, when the computer calculates x~l the computer actually calculates x~l(l + e), where |e| < 2~d. So in the calculation of JTS, the computer tries to calculate Df(Pl)-l-..Df(Pl+s)-1 but instead calculates Df(Pl)-l---Df(pl+s)-l(l + es), where (1 _ 2 - ^ + 1 ^ 1 + 6 s <; (1 + 2-d)2s+l. Hence we can bound the calculation of ]TS by JTs^P's(l-2-d)-2s-\ (2.13) We can bound the error the computer makes while calculating o~s in the same way. When the computer adds two floating point numbers xx, x2 it actually computes (xi + 18 x2)(l + ei) where |ei| ^ 2 d. So when the computer tries to calculate X\ + x2 + • • • + xn, it actually calculates .Ti(l + 77!) + X2(l + 7)2) + h + ??n) where, (l-2-d)n-1 ^ l + ?7i ^ (l + 2 - d ) n - 1 (1 - d~d)n+l-r ^ 1 + V r ^ (1 + 2-d)n+l-r for 2 ^ r ^ n. Hence, instead of calculating as, the computer actually calculates min(i+.s,./V) ]T \Df(pi)-1Df(Pi+1)-1 • • • Df{P])-l\{l + € i i ) ( l + tTij) j=i where 1 + arises from the roundoff error of the products, 1 + 77^ arises from the round off error when adding. We can bound each of these errors by (1 _ 2-rf)2(J-0+i <c 1 + 6 z j <c (1 _ 2 -*)2(J-i)+i and (1 _ 2-dy+2-(J-i-l) ^ 1 + m j ^ (1 + 2-dy+2-(J-i-l) ^ and combining them we get ' (1 _ 2-dy+U-i)+4 ^ ( 1 + £ . + ^ < f l + 2-dy+(j^)+\ Since (j — i) ^ p, we get (1 - 2-d)2s+i <: (1 + ei3)(l + 77^ ) ^ (1 + 2 " d ) 2 s + 4 , and the number a's, which the computer obtains when it tries to calculate a's satisfies a 7 ^ ( l - 2 - r f ) 2 s + 4 . (2.14) 19 Using similar methods as described above, it can be shown that the number T'S, which the computer obtains when it tries to calculate T~S satisfies T-S^T'S + a;5'{(l + 2-d)2s+5 - 1 } (2.15) where 5' = sup \Pl+l-f(Pl)\ (2.16) Equations (2.13), (2.14) and (2.15) are the numbers obtained by the computer to bound the numbers that we should obtain if there were no rounding error. Equations (2.11) and (2.12) give us a relationship between CTJ and 77 to a and T. How are we able to relate a and r to a and r? It is shown in [3] that a and a are norms of linear operators, and that ff*(dkj ( 2 1 7 ) and r T [(1 - a A ) - : a 2 A + a(2d - l ) " 1 ] 5' + 06 (2.18) Step 4 Now calculate 2MOT. If the value of this is less than 1, then Theorem 2.1 says there exists a shadowing orbit {xi}*LQ such that 2r \Xj — Vi ^ =, lor 0 < i <^ N. 1 1 1 + V I -2Mor Example : The logistic map The method outlined above will now be applied to the logistic map f(x) = Aa;(l - x), where 1 ^ A ^ 4. (2.19) 20 For this specific example I will take A = 3.8 and p0 — 0.6. The maximum round-off error the computer makes while computing the pseudo-orbit is 2~ 2 7 . Also, it has been determined that s = 30 is sufficient to satisfy /J S < 1. A l l calculations were done in using a Gnu C compiler on a Intel Celeron, 466 MHz processor. Step 1 To-estimate the quantity M, note that D2f (x) = -2A, so for our example M = 2A = 7.6. To estimate 8, let x represent a number stored in the computer where 0 ^ x ^ 1. Then m = xx(i-x)(i+o, where 1 + £ arises from the round off error caused by multiplying x and (1 — x). Also, note that ( l - 2 - r f ) 2 ^ l + £ ^ ( l + 2" f 7) 2. Thus, | 7 R - / ( x ) | < \Xx(l-x)(\ $ ^((1 + 2 - ^ - 1 ) . So for our example, we can take * = ^ ((1 + 2" d ) 2 - 1) . We can find an upper bound for A in a similar way to get A = 3.8(3 + 2- d + 1 )2~ r f . 21 Step 2 Now we obtain estimates for a and r. The pseudo-orbit was calculated in single precision and the all other calculations were done in double precision for which d = 54. These are the same values used by Chow and Palmer in [3]. For N = 426, 000 I obtained A4O = 1-372777636932113e-02 a'30 = 1.315562128027402e+01 T£0 = 4.579665661388762e-08 Step 3 The final estimates for a and r, were found to be a ^ 1.333873241593024e+01 r < 5.966314241498407e-07 Step 4 With the values of a and r from step 3, we obtain IMor ^ 0.000010157748749. Since this quantity is less than 1, we can conclude from Theorem 2.1 that there exists a shadowing orbit {xi}^L0 such that \pi - Xi\ ^ 5.966327901 le-07. for i = 0 . . . 426, 000. Thus the pseudo-orbit is £ x-shadowed by a true orbit with Sx = 6 x 10~7. This results is slightly better than that obtained by Chow and Palmer in [3], where 5X = 9 x 1 0 - 6 for a pseudo-orbit with the same initial condition. This difference could result from the fact that different computers will cause different round-off errors so the resulting pseudo-orbit could be different on different computers. 22 2.2.3 Existence of a shadowing orbit: another method. Another method for proving the existence of a shadowing orbit for one-dimensional maps is outlined in Grebogi et al. [6]. The idea of this method stems from the fact that for chaotic attractors, distances between two nearby orbits, grow approximately exponen-tially with every iterate. If we consider the inverse of such a map, the distance between two nearby orbits shrinks exponentially as we iterate backwards. By taking advantage of this fact, we construct a sequence of intervals which contain the pseudo-orbit. These intervals are constructed in such a way that by iterating an interval forward using the mapping, the mapped interval contains the next interval in the sequence. If we are able to construct a sequence of intervals which satisfies the description above, we can bound the shadowing distance of the pseudo-orbit. A procedure for constructing the sequence of intervals is described in detail in the proof of Theorem 2.2. Using a Gnu C compiler on a Intel Celeron 466 MHz processor, we compute a pseudo-orbit for the logistic map f(x) = Xx(l-x). with A = 3.8 starting with initial point p0 — 0.4. Let P i denote the ith point of the pseudo-orbit. The following theorem can be proven: Theorem 2.2 For N=1,000,000, the pseudo-orbit {Pi}f=0 with A = 3.8 and p0 = 0.4 is 5x-shadowed by a true orbit {xi}f=0 where 8X = 10~7. M e t h o d of proof: The pseudo-orbit {pi}f=Q is calculated in single precision. A l l other calculations are done using double precision. A true orbit { . T j } ^ 0 , is found by finding a sequence of intervals {Ii}f=0, such that .Xj e f, for 0 ^ i ^ N. Without actually calculating X j we can bound the location of each Xi and get a bound on the shadowing distance for the pseudo-orbit. 23 The intervals {Ii}^L0 are chosen such that f{h-i) 5 Ii, for 1 ^ i <i N, (2.20) and are found by first defining 1^, and then working backwards to find IN-I, I N - 2 , and so on. If the intervals {Ii}fLQ can be successfully determined, then given .T; e there exists Xi-i G such that X i = / ( X i _ x ) and {xi}^LQ is a true orbit with x^ e 7 n for a l H = 0 , . . . , N. • In order to prove the theorem above, we must be able to construct the sequence of intervals {Ii}^=0- If U-i is chosen sufficiently large, then condition (2.20) can be verified numerically. However, if I^i is chosen too large, then it may be impossible to verify condition (2.20). So the objective is to chose 7i_i just big enough so that (2.20) will hold. We will start by defining the interval I% — [s~,sf]. The endpoint interval is defined to be IN = [PN — £N>PN + £«]• Given an interval Ii, we must show how to construct 7j_x such that we satisfy condition (2.20). This is done taking the inverse of the endpoints of a "thickened" version of The thickened version of f is defined by Ii = [s,~ — Si, sf + We have freedom to choose E\ > 0, for 0 ^ i ^ N. For the proof of this particular theorem, Ei is chosen to be constant Si = 2" 4 8 , for l ^ i ^ N . The reason for choosing this value for the e^s is that we cannot expect to find a shadowing distance that is smaller than the precision we do our double precision calculations in. The double precision floating point error the computer makes is 2 - 5 4 . By making £j = 2~ 4 8 for all i, we can numerically check that the intervals are in fact thicker. Perhaps a method which would vary et depending on the size of the one-step error would be better. If the 24 one-step error at time step i is e^ , then we would create f such that > 2 1 | . Future research into this would be a good idea. The inverse of the map, is defined in such a way that at time step i — 1, f~1(pi) and f~l(Ii) lie on the same side of the critical point 0.5. With this definition for the inverse map we calculate the previous interval in the sequence by ii-1 = r1(ii). This is a first approximation for If (2.21) then we are in good shape. If we cannot verify (2.21), then we further thicken 77;_i until we can satisfy (2.21). As before, the computer does not actually calculate / and Computer calculations will cause some round off error and the computer is in fact calculating f(x) and f~l(x). An upper bound for the difference between f(Ii-\) and f{Ii-i) can be found. For 0 < x < 1 and 1 < A < 4 we get f(x) - 2 % < J{x) < f(x) + 2 5eM, where e ; i = 2~ 5 4 is the double precision machine epsilon. The extra factor of 2 5 is added to ensure that, we can place bounds on the double precision error. Once we verified (2.21) for an interval — [s~_l, s^i^, we correct for the double precision error. We can do this with the formula / i = [ S l - - 2 - 4 9 , ^ + 2- 4 9 ] . To get a upper bound on 5X we first define the quantity 5x(i) = m a x O , - s~\, \Pi - s+\}. (2.22) 25 The above equation calculates the maximum distance from P i to the furthest end point of the interval I%. If we take the maximum of 5x(i) for 1 ^ i ^ JV we get, (L ^ max \5T(i)}, which gives a bound for the shadowing distance of the true orbit. • Note that the pseudo-orbit was calculated in single precision where the noise level is 2~ 2 7 , and all other calculations were done in double precision where the maximum round-off error is 2 - 5 4 . Using these values we were able to successfully contain the pseudo-orbit. 2.2.4 Discussion of one-dimensional shadowing Sections 2.2.1 and 2.2.3 demonstrate two different methods for shadowing one-dimensional maps. Perhaps the biggest difference between the two methods is that the first method works forward in time and the other works backwards in time from the end of the orbit. The procedure by Chow and Palmer in Section 2.2.1 calculates the pseudo-orbit forwards in time while at the same time calculating the necessary quantities needed for the proof of the Theorem 2.1. At any time step N we can decide whether our theorem applies, and if it does, we can calculate the shadowing error. Another way of looking at this is that we can continually check whether a true orbit is shadowing the pseudo-orbit. As soon as the theorem fails we stop the program. The procedure by Grebogi et al. in Section 2.2.3 first calculates the entire orbit and then works backwards in time to obtain the shadowing error. So this procedure is not as dynamic as the first because we restrict ourselves to orbits of a fixed length N. Another important difference between the two methods is that one makes explicit use of the inverse of the map while the other does not. This is a result of one method working forward and the other backward in time. The method in [6] makes uses the inverse of the map to iterate intervals backwards so as to contain the pseudo-orbit. The 26 problem is that the inverse map may not always be well defined, and special care must be taken so that intervals are mapped backwards properly. The method by Chow and Palmer does not use the inverse of the map in any of the calculations, thus avoiding the the difficulties discussed above. The versatility that the Grebogi et al. method lacks by running backwards in time is sometimes gained back by being easily applicable to different pseudo-orbits. The method does not need to be changed significantly in order to be applied to other pseudo-orbits. A l l we need to do is change the initial condition. Thus many simulations can be run soon after each other. The method by Chow and Palmer lacks this feature. The cause of this is in finding a suitable choice of s. For many pseudo-orbits choosing a suitable s can be a time consuming process. In fact, it was found that as A increases past the period doubling cascade, it became increasingly difficult to find a suitable s to ensure that ns < 1 for the logistic map. A difference between Theorem 2.1 and 2.2 is that the first is true for any C2 function on the unit interval [0,1] whereas the second is specific to the logistic map, although the method of proof for Theorem 2.2 can be applied to other maps to obtain similar results. For the Baker map, where po — 0.6, the following theorem can be proven using this method. Theorem 2.3 For N = 1,000,000, the pseudo-orbit {Pl}f=0 with P o = 0.6, is shadowed by a true orbit {xi}^LQ where 8X ^ I O - 7 . It is my belief that the method of proof for Theorem 2.2 is applicable to any C2 map of the unit interval / : [0,1] —> [0,1], provided a suitable / _ 1 can be chosen. 2 ( 1 - x ) , if 0 ^ x ^ 1/2 if 1/2 < x ^ 1 27 As mentioned earlier, if / is not hyperbolic then it is not necessarily true that ev-ery pseudo-orbit is shadowable. Tests were done using the method in Theorem 2.2 for different parameter values. It was found that for po = 0.6 we can show that for A = 3.6, 3.73, 3.81, 3.87, and 3.9, a true orbit shadows the pseudo-orbit for N = 330, 000. However, for p0 = 0.6 and A = 3.84, 3.83, we could not successfully construct a sequence of intervals which contain the pseudo-orbit for N = 1,000,000. This may suggest that there does not exist a true orbit which shadows the pseudo-orbit for these parameter val-ues. However, when we apply the Chow and Palmer method to these parameter values, we find that a true orbit does indeed exist! So, where one method failed to show that a true orbit exists, the other was able to provide us with the necessary information. In fact, it was found that the method of Chow and Palmer worked for many simulations where the other method failed. For a rigorous proof of the existence of a shadowing true-orbit for a pseudo-orbit, it is recommended that the method of Chow and Palmer be used. This seemed to be a more reliable method for obtaining shadowing errors. The method takes extra care in determining roundoff errors for every step of the calculation. Also, this method is nice since we do not have to calculate the inverse of the map. The problem with this method is sometimes it can run significantly slower than the second. If quick results are required, then the second method is recommended. Just beware that if that method fails, it does not necessarily mean that a shadowing orbit does not exist. If this happens, try the first method to obtain a better result. 2.3 The two-dimensional case The methods discussed in-Section 2.2 can be easily extended for two-dimensional maps. In fact, Theorem 2.1 can be changed for C2 maps, / , where / : 5Rfc —> Also the method of proof for Theorem 2.2 can be changed for the two-dimensional well. 28 Instead of mapping intervals backwards in time, we map parallelograms forwards in time to contain the pseudo-orbit. In this section we will discuss these methods and provide some computer aided proofs. 2.3.1 Two-dimensional shadowing: Chow and Palmer's method The first method we will discuss, is a continuation of the study of one-dimensional shad-owing by Chow and Palmer. Let / : -> 5Rfc be a C2 function and let {pl}iLo be a pseudo-orbit of / with noise level 5. We want to show that there is a true orbit {xi}f_0 such that, \xi -pi\ ^ e, for n = 0 , . . . , N — 1. To do this we use the following theorem by Chow and Palmer [4]. Theorem 2.4 Let f : 3?fc -» 3ftfe be a C2 function with M = sup{|£> 2/(z)l : x E 3?*}. Let {pj}^ 0 be a 5-pseudo-orbit of f with 2 M | | L - 1 | | 5 ^ 1, where L~l is a right inverse of L which is defined in the proof of this theorem. Then there is an true orbit {xi}f=0 of f such that \Xi - P i| ^ 25IIL"1!! ( l + V/1-2MHL-1H25)" 1, for i = 0,1,...,N. Proof: The proof of this theorem is similar to the proof of Theorem 2.1. The following proof is presented by Chow and Palmer in [4]. 29 The true orbit satisfies xi+i — f(xi) for 0 ^ i ^ N - 1. If we set xl — pt + eiy we can show that &i satisfies el+1 = Df{Pi)ei + gi(ei), (2.23) where 9i(e) = fipi) - Pi+i + f(p* + e) - f(pi) - Df(Pi)e. We want to show that there exists a sequence {ei}fL0 such that |e,| ^ e = 25\\L-1\\ (l + ^1 - 2 M H L " 1 1 | 2 5 ) , for i = 0 ,1 , . . . ,N . Define S to be the set of all sequences e = {ei}f=Q such that max lej ^ e. Note that S is a compact convex subset of $tk(N+lK Rewrite equation (2.23) as Le = g(e) where (Le)i = ei+1 - Df(pi)ei. Now we define the mapping T on S by Te = L-lg(e), where L~l is a right inverse of L. The choice of L - 1 is important and so is its existence. For an explanation of the choice of L~l refer the discussion in the following subsection. 30 The mapping T is a continuous mapping of S into $lk(N l \ Also, T maps S into itself. To show this, note that g is a continuous function and \9i(e)\ ^ \f{pi)-Pi+i\ + \fiPi + e)-f{pi)-Df(pi)e\ ^ S + ^M\e\2. Also, if e € S then, \(Te)i\^ \\L-l\\ (S + \ME^ =e. So, T maps S into itself. So, from Brouwer's fixed point theorem, T has a fixed point e = {e;}^_0. Since LL~X = I, e satisfies Le = g(e). This means that e; is a solution to (2.23) such that j | ^ e for i = 0 , 1 , . . . , N, and Xi = Pi + e*. • The factor plays the same role that a and r, had in Theorem 2.1. We can think of H-t' - 11| as a measure of expansiveness. Also, if 8 is the maximum error made by the map at each time step, then | | L - 1 | | 5 is approximately the distance away from the true orbit. So, if we are able to choose L~v in such a way to make | | £ - 1 | | as small as possible, we can use Theorem 2.4 to obtain a bound on the shadowing orbit for the pseudo-orbit. Notice that Theorem 2.4 holds for any /c-dimensional maps. For the remainder of this section we will restrict ourselves to two-dimensional maps. Choice of L Our problem is to determine a suitable L - 1 such that 2M||L - 1 | |<5 ^ 1. One assumption we will make is that the map / is somewhat hyperbolic although it may not be uniformly hyperbolic. Then we look at the linear map, ul+1 = Df(pi)ui, (2.24) for i = 0 , 1 , . . . , TV — 1. Given any sequence {hi}^1 in RkN, the difference equation ui+1 = Df(pi)ul + hi, 31 has many solutions for i = 0,1, 2 , . . . , N — 1. So the linear operator L : $lk(N+1) —> ynkN defined by (Lu)i = ui+i - Df(pi)ui, for i = 0,1, 2 , . . . , N — 1 is onto and has a right inverse. The subspace generated by (2.24) converges to the unstable subspace for the pseudo-orbit. So we set Pi Ui and choose /3t to be orthogonal to rji such that the matrix Si = [Pi Vi has determinant equal to 1. If we make the transformation Ui — Si'Vi, (2.25) then equation (2.24) becomes Vi+l (2.26) where a, = \Df{pl)pi\, bi operator I by o-i bi PT+iDf(Pi)Vi, and c, = vf+iDf(Pi)Vi- Define the linear (lv)i = vi+i -^ bi 0 a Vi. (2.27) With I defined in this way, we can see that for L defined as in the proof the the theorem that {Lu)i = Si+i(lv)i. 32 So, if / 1 is a right inverse of /, we define the right inverse of L by where h is in K 2 / v Since Si is orthogonal we can show from (2.28) that Hi/- 1! ; - H (2.28) (2.29) Now the problem is choosing a suitable / 1 . To do this, we need to solve Vi+l = for Vi so that v.[ is as small as possible. If Vi = [ v l i \ v ^ ] T , then we can show from (2.30) that, o,i bi 0 d vz + hi, (2) f n—l Ylci-\---cihf iii ^i^N, 3=0 0 if i = 0. and N-l - £ a7l • • • a f ib3vf] + hf\ i f 0 ^ i ^ N - 1 3=1 0 if i = N. Equations (2.30), (2.31) and (2.32) define a right inverse ,(2) 0 < i < N. (2.30) (2.31) (2.32) (2.33) Our goal is compute \\l To do this, define the quantities i-l SUp E l c i -1 • • • Cj | if N ^ 1 , a{N,2) = { i ^ ^ = o if AT = 0, (2.34) 33 and sup NEK1---ajl\{\bJ\a{j,2) + l) if TV > 1, a{N,l) = { °$K*-ij=i (2.35) 0 if N = 0. It is easy to see from (2.34) and (2.35) that |wj 2 )| ^ a{N, 2) sup \hj\ for Q^i^N, "3\ and ,0) vf'\ ^ a{N, 1) sup \hj\ for 0 < i < TV. Hence, ^ ^ ( / V , 1)2 + <7(/V, 2)2||AI||, (2.36) and from (2.29) we get HL"1!! ^ y/a(N,l)* + o-{N, 2)2. (2.37) Now we have reduced the problem of determining to the calculation of a(N, 1) and a (TV, 2). The calculation of a(N, 1), and a(N,2) may be difficult if TV is too large. We can truncate the sums in (2.34) and (2.35) in much the same way as we did in the one-dimensional case. We calculate instead { min(i+s,N-l) max £ | a r i . . . a T i | ( | 6 i | f f , ( j , 2 ) + l ) if JV > 1, 0 if TV = 0 and crs(TV, 2) i - l max \ci-i---Cj\ if TV ^ 1, 1^^ iVJ=max(0,i-s-l) (2.39) 0 if TV = 0, 34 where s is a positive integer such that 0 ^ s ^ TV. We can use these truncated sums to approximate a(N, 1) and a(N,2) in a similar way as we did for the one-dimensional case. Define { sup \a~l... a~+s\ if TV ^ s + 1, 0 if 0 ^ TV ^ s, and sup \ d . . . C i + s \ i f ' T V ^ s + l , s<i<N-l 0 if 0 ^ TV s, then it can be shown that cr(TV, 1) and a(N, 2) satisfy a(N,l) ^ ( l - / i s ( / V , l ) ) - 1 ( l - / i s ( i V , 2 ) ) - V , ( / V , l ) , (2.40) a(TV,2) ^ ( l - / i s ( 7 V , 2 ) ) - V s ( / V , 2 ) . (2.41) These truncated sums are much faster to compute. We can further bound a, (TV, 1) by noticing that m'm(i+s,N-l) • as(N,l) = max ••• aj11(^^,2) + 1) j=i min(i+s,N-l) ^ o l t E \a;l...ajl\(\bJ\as(N,2) + l) (2.42) j=i = < ( i V , l ) (2.43) for TV > s + 1. We improve the computing time for calculating as(N,l) by making this approximation since the calculation of the term as(j, 1) requires the most processing time as j becomes large. The values of as(j, 1) for 0 ^ j^ TV, could be stored in a vector of length TV or even a file, but if TV was too large we could run into memory problems or it would slow down our process to read off a file. It was found that for all pseudo-orbits where a shadowing orbit exists, that the bound for e did not differ by a significant amount by making this approximation. 35 2.3.2 The computation Given a map / : 5R2 —> 5ft2 and an initial value p0, we wish to use the results in Section 2.3.1 to determine how closely the pseudo-orbit is shadowed by a true orbit. To do this, we go through the following steps: Step 1 First calculate the upper bounds for the quantities: 5 = sup|/(pi) - f{pi)\, A = sup\Dr&)-Df(pi)\, B = sup |D/(p n ) | , M = sup | / J 2 / ( P n ) | . These quantities will be used for bounding round-off error effects in Step 4. Step 2 Compute the pseudo-orbit pi, the unit vectors /?,; and 'qi, and compute the quantities a„i, bu and ct for 0 ^ % ^ N. Remember that the computer actually calculates rf^ a';, b't, and c[ because of round-off errors. Step 3 Compute the quantities p's(N,l), p's(N,2), a's(N,l) and a's(N,2). We must choose s in such a way so that p's(N,l) < 1 and p's(N,2) < 1. Since we are using f(x), and Df(x) instead of f(x) and Df(x), we eventually want to obtain estimates for cr(JV, 1) and a(N,2). We can proceed like we did in Step 3 in Section 2.2.2 to obtain estimates 36 for fj,s(N, 1), fj,s(N, 2), a-s(N, 1), and a(N, 2). It can be shown that /is(N,l) < [j,'s(N,l)(l-2-d)-2s-1 (2.44) /X(iV72) ^ / / s ( i V , 2 ) ( l - 2 - r f ) - s (2.45) and a s(7V,l) <; a's(N,l)(l-2-d)-*°-\ (2.46) ^ ( iV,2) ^ a's(N,2)(l-2-d)-2s. (2.47) Equations (2.44)-(2.47) estimate the round-off error effect. After estimating the round-off error effect we then estimate the truncation error. From equations (2.40) and (2.41) we get, a(N,l) ^ (l-tisiN,!))-1^-^^))-1*^,!) (2.48) a(N,2) <: ( l - / x s (7V,2 ) ) -V s ( iV ,2 ) . ' (2.49) Step 4 Compute an upper bound for \\L taking into account all the calculation errors. We can identify three sources of error in our calculations: 1. The round-off error and truncation error involved in the calculation of || / 2. The error due to the fact that we are actually estimating whereas what we want to estimate is ||/ 3. The error involved in estimating from \\l To estimate the first source of error, note that is estimated by |rx|| <; y/o-s(N,i)* + o-s{N,2y. 37 Since the computer works with f(x) and Df(x), we can estimate by II / - 1 ! ! ^ yJa'N,l) + a(N,2). It is shown in Appendix 1 of [4] that 7-1 - z - 1 ! ! ^ V 3 ( ^? + ( 2 - 5 ° ) where 0i = (1 + h? h = (1 + A, 2 -d+2 = ( 1 - 2~d) I^_^_+A+B[(1+2-,r_1]) + 2 A + 2 1 ( 1 + 2 ~ r - 1 1 ( s + A ) ) + B ( 1 + w < f e It is also shown in [4] that ||Z l\\ satisfies l i r ^ K a - i i z - R i i i i R i i ) - 1 ! ! ^ ! . (2.51) Using (2.51) we can estimate the second source of error for the calculation of | |Z | | _ 1 . The third source of error arises because the matrix Si may not be exactly orthogonal for 0 ^ 1 ^ N. The operators and are related by ( £ - 1 h ) i = Sfr'S-'h), for 0 ^ i ^ N — 1. So, since Si may not be exactly orthogonal it is clear that llZr 1 !! ^ max ISA- max IS1, -1! • IIZ x II. (2.52) O ^ i V - l Now we try to bound for 0 ^ i ^ N — 1. Matrix Si = [/?,;, 77^] was chosen so that and rji are orthogonal and of the same length. It follows that, \Si\ ^ IAI 38 It is shown in [4] that we can estimate |/34| by Hence, and finally, 1 + 3 ( II 7 — 1 I 1 (2.53) Equation (2.53) will give us our final estimate for \\L Step 5 Calculate 2M| |Z/ _ 1 ! |5 . If this is not greater than 1, then Theorem 2.4 says that there is a true orbit {xi}^L0 such that ^ 25HL-1!! ( l + v"! - 2 M | | L - 1 | | j ) - l for 0 < i < N. Example : The Henon map The method outlined above will be applied to the Henon map f(xi,x2) 1 — ax\ + x2 —bx\ (2.54) for a = 1.4, b = —0.3, and p0 = [0;0]. For this example, all calculations were done using double precision arithmetic where the maximum error the computer makes when adding, subtracting, multiplying and dividing is at most 2 ~ D . Also, a quick calculation has already determined that for s — 40 the quantities /J,s(N, 1) and LIs(N, 2) are less than 1. This choice of s is the same used by Chow and Palmer in [4] 39 Step 1 To estimate 5, let xx and x2 be machine numbers, and the numbers 1 - ax\ + x2 and -bxi denote the numbers the computer obtains when it tries to calculate 1 — a.x\ + x2 and —bxi. Now, following the methods in [10] we can show that, 1 - ax\ + x2 = {[l- ax\{\ + ei)( l + e2)] (1 + e3) + x2}(l + eA) where |e,| ^ 2~d. So, |1 — ax\ + x2 — (1 — ax\ + x2)\ = |e 4 (l - ax\ + x2) + [e3(l + ax\) -axl(et + e2 + eie2)(l + e3)](l + £4)! ^ 2'd(l + \a\\x\2 + |x|) + [2" r f(l + |a||x|2) + |a||x|2(2-t/+1 + 2"M)(1 + 2-d)](l + 2~d) = 2'd(2 + 2'd + \x\ + [4 + 6 • 2-d + 4 • Td + 2 _ 3 d ] |a | |x | 2 ) Similarly, we can show that —bx\ = —bxi(l '+ e5), and \-bxi - ^ \ - bxi{l - e5) - (-bxi)\ ^ \bxi I e5 ^ \b\\x\2-d. So, if ^ R for 0 ^ 1 ^ N, we have 5 ^ 2 " V ( 2 + 2" ( / + + [4 + 6 • 2-d + 4 • 2~ d + 2-M]|a||i?|2)2 + ft2^2. (2.55) We can estimate A and B in a similar way to show that A ^ 2- d(4 + 2" d- 1)|a| JR, (2.56) B ^ v /4a 2 JR 2 + l + 62. (2.57) 40 Also, since only one second partial derivative is non-zero we have, M = 2\a\. (2.58) Step 2 The pseudo-orbit pi, the unit vectors /3j and ry,;, and the quantities a*, 6,: and cz were calculated for 0 ^ i ^ N, where JV = 330, 000. A l l calculations were done in double precision using a Gnu C compiler on an I B M compatible computer with an Intel Celeron 466 MHz processor. The maximum error the computer makes in double precision is 2~ 5 4 . Step 3 It was found that fi's(N,l) = 2.082322390632631e-09, p'S{N,2) = 4.725229628391131e-29, CT'S{N,1) = 8.688439701823398e+02, <r'S{N,2) = 1.368108060968983e+01. So, from inequalities (2.44), (2.45), (2.46) and (2.47) we get fia(N,l) ^ 2.082322390632669e-09, HS{N, 2) 4.725229628391173e-29, aS(N,l) <C 8.688439701823476e+02, as(N, 2) 1.368108060969008e+01. Finally, equations (2.48) and (2.49) give a(N,l) ^ 8.688439719915608e+02, a(N, 2) sC 1.368108060969008e+01. 41 Step 4 In order to calculate we first note that A ^ 1.666372867159099e-15, M = 2.8, 8 ^ 2.976135796265091e-15, B ^ 3.844481929370897e+00. Using the results from Step 3, we can show HF1!! ^ 8.689516785372596e+02. Using the above quantities, it follows from (2.50) and (2.53) that — ^ 2.299604736526933e-14. \\L"l\\ ^ 8.689516785546249e+02. Step 5 From the information obtained in Step 4 we get, 2M\\L~l\\5 ^ 1.258438578407314e-08. Since this quantity is less than 1, the theorem says there is a true orbit, x;, such that \xi -Pi\ ^ 5.172236358997404e-ll for 0 ^ i K. N. The results from Step 4 and Step 5 show that our pseudo-orbit, which is calculated with an error of 5 ^ 3e-15, has a true orbit which lies within e < 6e- l l . So after N = 330, 000 time steps, we only lose about 4 digits of accuracy. 42 2.3.3 Two-dimensional shadowing: the GHYS method In Grebogi et al. 1990 [7], the authors take the ideas explained in Section 2.2.3 and adapt them for two-dimensional dynamical systems. They call their method containment and it is described below. The containment of a true orbit {xi}fL0, which shadows a pseudo-orbit {pi}^=0, re-quires the construction of a sequence of parallelograms {Mj}£f_0. The parallelogram M j is centered at p,L for 0 < i ^ Af. Also, these parallelograms are constructed such that / ( M j ) lies across M i + 1 in a particular way. First, we designate two parallel sides of each M j as the expanding sides and the other two parallel sides as the contracting sides. The image of the expanding sides of M j must intersect the contracting sides of M i + 1 and cannot intersect the expanding sides of M i + 1 , see Figure 2.1. The expanding sides for each M j are chosen to be parallel to the unstable unit vector at the point pi along the pseudo-orbit. When we are able to construct such a sequence of parallelograms we say that it is a containing sequence of parallelograms. Now we will argue that if we have a containing sequence of parallelograms, { M , } ^ , there must be a true orbit {xi}fL0 contained within { M j } ^ 0 where Xj is contained in M j for 0 <.i < N. Let 70 be a curve that lies in M 0 . Now, / (70 ) contains a curve 71 which lies totally in M x . The curve ji that runs from one contracting side of M.i to the other contracting side. In fact, we can find a sequence of curves, {ji}^=0, where yt is contained in fiji-i) and lies wholly in M j for 0 < i < N. Now choose a point xN which lies on the final curve 77V. We can find the previous point £JV_I by the inverse map j f - 1 ( . X j v ) and continuing in this way we can define Xi to be f~l(xi+i) for 0 < i < N. Thus we have shown that the sequence of parallelograms { M j } ^ 0 contains the true orbit. Notice there in fact exist infinitely many true trajectories since the final point xN was arbitrary and could have been any point along the curve 7 ^ . 43 Image of Expanding f (Bp) Figure 2.1: Containment of a true orbit. The shadowing distance at time step i is bounded by the distance from the point pi of the original pseudo-orbit to the furthest point on the i'th parallelogram M j . To find the maximal shadowing distance we take the maximum or the these distance over the whole orbit. For dynamics that which chaotic and not uniformly hyperbolic, containment may not continue forever. There can be some point where we won't be able to construct a parallelogram M w + i such that f(MN) will intersect it in the manner described above. This usually happens if the angle between the expanding and contracting sides of the parallelogram becomes small in comparison to the noise level of the orbit. So in essence, we lose a dimension of the parallelogram and hence containment of any orbit is unreliable. 2.3.4 T h e computation Step 1 Calculate the pseudo-orbit pi, the unstable unit direction uu and the stable unit direction Sj, for 0 ^ i ^ N in double precision (maximum floating point error is |e| < 2~d). 44 The unstable and stable unit directions can be approximated by Ui+i = 7 ^ 7 , (2.59) «i = 7 7 ^ 1 > ( 2 - 6 ° ) where V,Q and SJV are arbitrary unit vectors and Jj = Df(Pi). The calculations could be improved if we calculated u0 and exactly. Step 2 Construct the parallelogram M 0 centered around p0 so that two parallel sides are in the same direction as UQ and the other two parallel sides are in the same direction as So- Make sure that the dimensions of the parallelograms are small. For our purpose, we construct M 0 so that each side has length 1 0 _ i l . The reason for this choice was arbitrary and it seemed to work well. The only limitation to the size of the i ' th parallelogram is that it should contain the circle centered at pi with a radius slightly larger than the noise level S of the pseudo-orbit. Step 3 Now map the parallelogram using / . This gives us a new quadrilateral Q = / (Mo) , hopefully containing the point p\. If not," we must enlarge M 0 so that its image under / contains px. The object Q is probably not a true quadrilateral. Instead it is something which resembles a quadrilateral but with sides which are curved. However, the sides are on such a small scale that the curvature is negligible. So I do not think much is lost by assuming the object is a quadrilateral. The parallelogram M 0 under the map / is stretched in the unstable direction and contracted in the stable direction. So if we construct a new parallelogram M i , centered around pi whose dimensions are similar to that of M 0 . The quadrilaterals Q and M t should intersect in the way shown in Figure 2.1. If they do not intersect in this way, then 45 change M j so that they do. This is done by stretching the stable sides and shrinking the unstable sides of M x until the proper orientation is achieved. Step 4 Repeat Step 3 for the entire length of the orbit. To obtain an estimate for the shadowing distance, first measure the distance between P i to the furthest point cfo on the of the zth parallelogram by Sf(i) = \pi-qi\. (2.61) As before, we have to estimate for the effect of round-off errors since the computer is actually computing f(x). So if 5 = sup \f(pi) - f(pi)\, then we can estimate the shadowing distance by, 6f = max5f(i) + 25. (2.62) Example : The Henon map We wish to apply the method outlined above to the Henon map (2.54), where a =- 1.4, b = -0.3, p0 = [0; 0], and TV = 330, 000. Step 1 The pseudo-orbit {Pi}f=0 and the unit vectors, ut and Sj, were calculated for 0 ^ i ^ N. A l l calculations were done in double precision using MatLab on an I B M compatible computer with an Intel Celeron 466 MHz chip. The maximum error the computer makes in double precision is 2 - 5 4 . The vectors u0 and were chosen to be [1,0] and [0,1] respectively. Similar calculations were done using different choices for u0 and S/v but it led to no significant difference on the bound for the shadowing distance. 46 Step 2 The parallelogram M 0 was constructed such that each side has length I O - 1 1 . Step 3 - Step 4 The sequence of parallelograms, { M i } ^ 0 was successfully constructed for N = 330, 000. It was also found that there has to be a true orbit, {xi}fL0, such that |:ri-#1 ^ 2.1e-10 (2.63) for 0 i ^ 330, 000. 2.3.5 Discussion of two-dimensional shadowing Sections 2.3.1 and 2.3.3 demonstrate two different methods for determining if a true orbit shadows a pseudo-orbit. For both methods we must first calculate the pseudo-orbit. This is different from their corresponding one-dimensional cases where one method worked forward in time and we could do our calculations as we calculate the pseudo-orbit, and the other needed to run backwards in time and required the entire pseudo-orbit before we applied the technique. The reason for this is that both methods require the knowledge of the stable and unstable manifolds at every time step for the calculations, and both used the same technique for determining them. If we could use a method which calculated the stable and unstable manifolds without determining the entire pseudo-orbit, then both methods could be much more efficient. Perhaps the main problem with the GHYS method is that the image Q of the paral-lelogram M j under the map / may not necessarily be a quadrilateral since / is usually non-linear. So there is some difficulty in making sure that the next parallelogram in the sequence intersects Q in the way shown in Figure 2.1. This does not seem to present too many difficulties if we make sure the parallelograms are small enough. If the sides of the 47 parallelograms are small and our map / is 'nice', then the curvature over this distance may not be too great. Hence, it would not be too difficult to bound this curvature and construct the parallelogram M i + i . Besides the fact that the GHYS method is much more intuitive than the C P method, it does not seem to be any better than the C P method. In fact, it is my opinion that for iterative maps, the C P method is superior. The algorithms seem to run much faster, even if the algorithms are not written to be very efficient and especially if we use the approximation (2.42). Also, the C P method is more rigorous than the G H Y S method. In the G H Y S method, round-off errors may cause some points to be mapped out of the parallelograms and thus we cannot safely say that a shadowing orbit exists. The CP method measure quantities associated with the pseudo-orbit and uses them to determine if a true-orbit shadows the given pseudo-orbit. In Grebogi et. al. [7], the authors make the following conjecture: For a typical two-dimensional Hamiltonian ma,p yielding chaotic trajectories with a small noise amplitude 5j- > 0, we expect to find Sx ^ ^lJ5~f for a trajectory of length N = l/y/Sj . A l l of our results for the C P and GHYS method are consistent with this conjecture. In fact, the C P method gave significantly better results. 2.4 Higher-dimensional shadowing There are many interesting examples of chaotic maps whose dimensions are three or greater. It would be nice if we could apply the methods we used in the previous sections to these systems. Although we haven't looked into these problems specifically, note that Theorem 2.4 does not specify the dimension of the map / . So we should be able to apply this theorem for maps whose dimensions are greater than two. The method has to be altered slightly. Again, the first step in doing this is to triangularize the right inverse for L for the map / : Rk —> Rk. We do this by defining an orthonormal basis ew, e2o, • • •, ek0, 48 then we define So = [eio62o • • • e/co]-Then define a sequence Si of such orthogonal matrices recursively by choosing an orthog-onal matrix Si+X and an upper triangular matrix R4 such that Df{pi)Si = Si+1Ri. If we make the transformation Ui = SiVi, then the linear map (2.24) becomes vi+i - RiVi. Then our right inverse is found by solving the difference system vi+i = RiVi + hi for 0 ^ 7 ^ N. Following the methods in Section 2.3.1 we can obtain quantities similar to a(N, 1) and a(N, 2). Then we proceed in a similar fashion to bound all the errors that arise in our calculations. To extend the GHYS method for a system with dimension three may not be too diffi-cult. If we could construct a sequence of cubes that intersect in a similar way as in Figure 2.1 we could contain the pseudo-orbit and say that there is a true orbit which shadows the pseudo-obit. We would have to determine all the unstable and stable directions for the pseudo-orbit at each time step. To actually construct the the rc-dimension parallel-ogram for k > 3 may be a little more difficult than for the two- or three-dimensional case. Also, to determine whether the parallelograms intersect in a way similar to the two-dimensional case may be more challenging. Certainly this would be an interesting further area of research. 49 Chapter 3 Refinement "Beware lest you lose the substance by grasping at the shadow." -Saint-Exupery 3.1 Introduction In the previous chapter we discussed the problem of determining whether a pseudo-trajectory can be shadowed by a true trajectory. What if it was important to know a more accurate solution? If we were interested in finding a more accurate solution we could always run the simulation again but this time we could use a higher precision method to get a more accurate solution. The problem with this is that our new solution may be quite different than the solution calculated with lower precision, especially if the dynamics are chaotic. Now we would have two solutions which are inconsistent with each other. It would be more interesting if we could find a more accurate solution which shadows the less accurate solution. It might be wise at this moment to make a distinction between dynamical noise and observational noise. Observational noise does not affect the future evolution of the system. Errors made in laboratory measurements are usually of this type. Another example is computer output that prints fewer digits than are represented internally in the computer. On the other hand, dynamical noise does effect the future evolution of the 50 system. The round off errors introduced by numerical solutions of a system is an example of dynamical noise. Quinlan and Tremaine [8] have studied some existing noise reduction methods to refine noisy trajectories, and found that none worked as well as the method presented in [7]. The reason why this method works so well is that the G H Y S procedure reduces the dynamical noise instead of the observational noise. It is the dynamical noise which we will be interested in reducing. The method of reducing dynamical noise is called refinement. Refinement is an iterative process that perturbs each point of the pseudo-trajectory in an attempt to produce a nearby trajectory with less noise. We say a refinement iteration is successful if the trajectory before the refinement has noise of magnitude Sj and after the refinement it has noise Sj where Sj < and \i G (0,1). Let p = {Pi}iLo be a pseudo-trajectory with N steps which has noise level 5 > ne^ where is the machine precision, and n is some constant significantly greater than 1. The constant 77 should be chosen so that there can be room to make improvements to the trajectory p towards machine precision. Next, we let el+x = — f(pi) be the one-step error at step i +1, where |ej+i| < 6 for 0 ^ i < N. The set of one-step errors, e = {ei}^Lx, is estimated by a numerical integration technique that has a higher accuracy than that used to compute the noisy trajectory p. A l l of this describes a function g which takes the trajectory p as input and whose output is the set of one-step errors e, i.e. g(p) = e. Since the one-step errors are assumed to be small, we would expect p to be close to a zero of g if one exists. If g(p) = 0, then p would represent a trajectory with one-step errors which are all zero, i.e. p would represent a true orbit. Usually it is next to impossible to find a true orbit for a dynamical system even after we refine the orbit because we can only hope to reduce the noise level to eM after refinement. So it makes sense to introduce a definition for a numerical shadow for p. 51 Defini t ion A trajectory p' is a numerical shadow for a pseudo-trajectory p if 1. p' has a noise level 8', where 8' < 8 and 8 is the noise level of the pseudo-trajectory. 2. The distance between p' and p is bounded by some constant 8pi. The main refinement procedure that will be studied in the thesis is one first presented in Grebogi et al [7]. 3.2 The GHYS refinement algorithm 3.2.1 One dimensional case This section presents the G H Y S refinement algorithm for a one-dimensional map. This algorithm reduces the dynamical noise of the pseudo-trajectory and produces a less noisy trajectory. The method is as follows: Given a pseudo-trajectory p = {pi}fL0 we want to find a less noisy pseudo-trajectory p = {pj^g which remains uniformly near p. The one-step errors are given by e»+i =Pi+\ ~ f{pi)- (3-1) We will construct the refined orbit by Pi =Pi + Ci (3.2) where c,; is a correction term for step i. Combining equations (3.1) and (3.2) we solve for Cj + i to obtain Ci+l = Pi+l - Pi+l = / ( p l ) - ( e m + / ( K ) ) , (3.3) where pi+\ ~ f(pi)- Assuming Cj to be small, we expand f{pi) about pi to get f(Pi) « f(Pi) + Ji<k, (3-4) 52 where J; = Df(pi). Using equation (3.4) we can write the final approximation for the corrections as ci+l ~ JiCi — (3.5) Since Jj tends to magnify errors as the map moves forward in time, to counter this effect solve for Q in equation (3.5) to get By iterating (3.6) backwards we can avoid the problem of the linear map magnifying the one-step errors. The choice of cN is arbitrary, and we will take its value to be = 0. The a lgor i thm Step 1 Calculate the pseudo-orbit {pi}^=0 using single precision arithmetic. A l l other calcu-lations will be done in double precision. Step 2 Calculate ejv and J/v-i and set C/v = 0. Calculate the correction term C / v - i using (3.6) and calculate PN-I using (3.2). Step 3 Repeat Step 2 for i = N - 1, N - 2 , . . . , 1, 0. If at any time step i we have J; = 0, equation (3.6) is undefined and we are unable to use this method find a less noisy trajectory. Once we find a less noisy orbit, we can apply the method again to reduce the noise level even more. In fact the method described above is super-convergent: the number [ci+l + Ci+l, • (3.6) 53 of significant digits typically doubles on each iteration of the process. Using computer software such as Maple, we can try to reduce the one-step errors to even smaller values. Example: The quadratic map revisited In Section 2.2.2 we proved the existence of a shadowing orbit for the quadratic map f(x) = \x(l-x), for A = 3.8 with an initial value of p0 = 0.6 and the noise level of the pseudo-orbit is 2~ 2 7 . It was also shown that the the pseudo-orbit was shadowed by a true orbit with a shadowing distance of about 6e-07. To actually find the true orbit would be impossible, but can we find a numerical shadow for the pseudo-orbit? Using the method described above, we found a numerical shadow for the pseudo-orbit computed in Section 2.2.2 whose shadowing distance is 8' ^ 5.842828e-06. A l l calcula-tions were done in double precision using a Gnu C compiler except for the generation of the pseudo-orbit which was done in single precision. The maximum errors the com-puter makes in single and double precision are 2~ 2 7 and 2~ 5 4 respectively. The process was iterated only 3 times before the method converged and the method was not able to obtain a better numerical shadow. The method did not find a numerical shadow whose shadowing distance was comparable to the shadowing distance of the true orbit. In Sec-tion 2.2.2 we showed that a true orbit exists where 8X < 6e-07, but the numerical shadow has a shadowing distance 8' < 5.9e-06. Although the method did not quite reach the shadowing distance set in Section 2.2.2, it was still a very good numerical shadow for the particular pseudo-orbit. It was found that when the orbit was near the critical point x = | , that the numerical orbit was furthest away from the pseudo-orbit. This is because at points close to the critical point, || J " 1 ! ! will be very large. Hence equation (3.6) has a factor which may be very large. If J; becomes too large, the correction may cause the numerical shadow to 54 lie on the opposite side of the critical point. If this happens the two orbits will quickly diverge. Sometimes computer errors may cause the numerical shadow to lie outside the interval [0,1]. If this happens the numerical shadow will blow up and would not make a very good shadow for the original pseudo-orbit. There is another method which can be used to find a less noisy orbit which shadows the pseudo-orbit for this example. In Section 2.2.3 a sequence of intervals {Ii}^L0 was constructed to contain the orbit and provide us with a shadowing distance for a true orbit. We can construct a numerical shadow by choosing a end point, p'N, for the numerical shadow first. We choose p'N so that it is contained in the last interval, IN, of the sequence of intervals. Now, we use the inverse map f~l to map p'N backwards in time. The numerical shadow is obtain by iterating the equation Pi = r\p'i+l), forz = ( i V - l ) , . . . , 0 , backwards. Note, all calculations were done in double precision except for the original pseudo-orbit. At each time step i of the orbit, we must ensure that the point p\ lies in the interval If we can find and orbit {p^}^ 0, we can call it an numerical shadow of the pseudo-orbit {pi}^=0. Using this method found a numerical shadow with noise level 8' < 2~ 5 4 which shadows the pseudo-orbit in Section 2.2.3 to within 6e-08. 3.2.2 The two-dimensional case We are able to apply the method discussed in the previous section to dynamical systems in two dimensions. The method is very similar right up until equation (3.5) with the only major difference being that Jj now represents the Jacobian matrix of the map f(x). If the system f(x) is not chaotic, the correction terms cz could be computed directly using (3.5). If the dynamics are chaotic then J; will amplify any errors which occur in the unstable direction and after a few iterations equation (3.5) will produce nothing but 55 noise. To solve this problem, GHYS suggest splitting the error and correction terms into components in the stable (s.j) and unstable (UJ) directions at each time step, to get, &i cUiUi -f~ e S i Sj (3.7) Ci = cUiUi + eSiSi. (3.8) By splitting the error in this way we can better control the corrections made in the unstable direction at each time step. The unstable unit directions can be found by iterating the linearized map ui+x = J j i t j , ui+i = T~~~~7 i (3.9) K + i l forwards in time. We need to choose an initial unit vector u0 for the above equation. This vector can be arbitrary because after a few iterations the unstable vectors will point roughly in the actual unstable direction (unless we were unlucky and actually chose uo to lie exactly in the stable direction at step 0). Similarly, we can calculate the stable unit directions Sj by iterating the map Si = Jilsi+l, sl = T ^ T , (3.10) N backwards in time choosing an arbitrary vector sN. The inverse of the linear map, J " " 1 can be calculated by inverting Jl. Substituting equations (3.7) and (3.8) into equation (3.5) gives Cui+lui+i + cSi+1si+1 ^ Ji{cUiUi + cS is,) - (eUi+lui+l + eSi+lsl+x) (3.11) For the same reasons that J j magnifies errors in the unstable directions, it diminishes errors in the stable directions. Likewise, J~l diminishes errors in the unstable directions and magnifies errors in the stable directions. So by taking components of equation (3.11) in the unstable directions at time step i + we can find the unstable components of the 56 correction coefficients by iterating backwards on ( c u i + i + e U i + 1 ) / | JiUi (3.12) Similarly, by taking components in the stable directions, we find the stable components of the correction coefficients by iterating forwards on The initial choices for cSQ and cUN are arbitrary as long as they are both small, i.e. smaller than the maximum shadowing distance. GHYS chooses them both to be 0. This seems as good as any since if there exists one shadow then there exist infinitely many shadowing trajectories so we would just end up finding a different shadowing trajectory for different choices of cso and cUN. An alternative way of looking at these end point conditions is that we want to 'twinge' the growing components of the end point of the pseudo-orbit and the backwards growing end point of the initial point. This means that we only introduce changes to the unstable directions of the initial conditions for the shadow and noisy trajectories. If we were to introduce any changes to the initial condition in the stable direction, the perturbation would die out leaving the new orbit to follow the true orbit that passes through our initial conditions. In other words, the shadow orbit will follow the orbit which is already known to diverge exponentially. Example: The Henon map The above procedure was applied to the Henon map (2.54) where a = 1.4, b = —0.3, and Po = [0; 0] for i = 0 , . . . , 330000. The pseudo-orbit was calculated in single precision where the maximum round-off error the computer makes is 2~ 2 7 . A l l other calculations were done in double precision where the maximum round-off error is 2~ 5 4 . A Gnu C compiler was used on a Intel Celeron 466 MHz processor. It was found that the GHYS algorithm converged to a numerical shadow whose noise level was 2~ 5 4 in only 3 iterations. The c s i + i — \JiSi\CSi (3.13) 57 numerical shadow was within 10 5 of the first (noisy) pseudo-orbit for i = 1,..., 33000. This suggests we only lost about 3 digits of accuracy. 3.2.3 Generalizing the GHYS algorithm for Hamiltonian sys-tems This section describes a generalization of the G H Y S algorithm for 2D-dimensional Hamil-tonian systems. The reason for restricting the procedure to Hamiltonian systems is that for this class of systems, the linear map J has D expanding directions and D contracting directions. For more general dynamical systems, the number of expanding and contract-ing sides may not be the same and the procedure described below would have to be modified. The following procedure is a modification of the GHYS algorithm taken from [8]: At time step i, let {u°i}^=-i be the D unstable unit vectors, and let {s{}f=1 be the D stable unit vectors. The vectors are evolved exactly as before except we use Gram-Schmidt orthonormalization to produce two sets of D-orthonormal vectors at each time step. Since we do not know which directions are stable and unstable at each time step, we choose an arbitrary orthonormal basis uJ0 at time zero, and an arbitrary orthonormal basis sJN at time step N. The unstable vectors at step i 4- 1 can be found as follows: 1 - W ui+l = [I -u\+lu\+l) • Jiu: (i - uil+lu\+l) • jlUfy (1 - u\+lu}+l ttj+iMJ+i) • JjUJ 1 + 1 1 ' l - ui+lui+l ui+lui+l) • Ji'ul I (j = 3 , . . . , D ) , 58 and the stable directions at time step i are given by s J " V \Ji S\+l\ "1.2 Si 2 _ 0--Sisl)-J;~Si+1 (1 S i S } ) ' Ji Si+l\ j _ (1 ~ s}s} - • • • - s\ s{ ) • Jj sJi+L _ i d - ^ j — s r l s n - j r i s u [J "•" j-After a few time steps i forward from i = 0, we would expect to find that u\ is nearly aligned with the most unstable direction, uf is aligned with the second most unstable and so on. Also we would expect to find that, after stepping backward in time from i = N, s} is nearly aligned with the most stable direction, s 2 is aligned with the second most stable direction, and so on. The higher dimensional equivalent of the error and correction vectors are given by D Yp ,-,r? -+-i = i D i = i The equation for the correction coefficients for the unstable subspace at time i + 1 is D D E ^ + 1 + e u i > i a = E c ^ 3 = 1 3=1 If we project out the component along uk+l we find that D cuki+1 + e u k i + l — EC«i U W ' ^iUi • j=k To find the unstable correction coefficients, we set CUJ = 0 for j = 1,..., D, and compute the coefficients backwards using °ui ,„k "i+l . J k ( °uk+l + e u k + 1 - EC«H+1 ' JlUi ) 1 1 \ j>k / 59 We first solve for {cup}fL0 which does not require knowledge of the other cu coefficients, then we solve for {c D-i}fLQ, and so on. 'i We can find the stable correction coefficients in a similar way by iterating forwards on 6 z Ji bz+l \ j > k J We solve for the cs coefficients the same way that we solved for the cu coefficients except now we start at time step 0 by setting CSJ = 0 for j = 1,... ,D. Then we compute { c s p } ^ 0 , which does not require the knowledge of the other cs coefficients, then we compute {csp-i}iL0, and so on. 3.2.4 Discussion of the GHYS algorithm It is important to note that there is no guarantee that refinement converges towards a true orbit. If there was, then all noisy orbits would be shadowable. In fact refinement alone does not prove that a true shadow exists, it only shows that a numerical shadow exists. In [8], it was shown for simple chaotic maps if the GHYS refinement procedure refined the trajectory to 1-step errors of about 1 0 - 1 5 , then successful refinements could be continued down to I O - 1 0 0 . It is reasonable to assume from this that refinement would continue to decrease the noise, converging in the limit to a true trajectory. To rigorously prove that a true shadow exists we must use methods like the shadowing procedures discussed in Chapter 2. There is also no guarantee that the GHYS refinement procedure, if it does converge, will converge to a reasonable shadow of our pseudo-trajectory. It may converge to a true orbit that is far from our original pseudo-trajectory. This could happen if the 1-step errors become so large that the linear map becomes invalid for computing corrections. Inaccuracies in the computation of the corrections can quickly amplify the noise instead 60 of decreasing it. 61 Chapter 4 Conclusions and future work "One sees chasing after shadows more fools than one can count." - La Fontaine, Fables. (1668) In Chapter 1 we asked the question: under what conditions will a computed orbit be close to a true orbit of the model? We then said that we can trust a computed solution if we can show that there exists a true orbit that shadows the computed orbit for a long period of time. Chapter 2 was dedicated to methods on how we can determine if a true orbit shadows a pseudo-orbit. We showed that despite errors growing exponentially in chaotic systems, many computer simulations have true orbits which stay close for long periods of time. Chapter 3 gave a noise reduction technique which reduces the dynamical noise of a system. If the noise reduction is able to reduce the noise level to machine precision, it suggests that perhaps a true solution does shadow a given pseudo-orbit. The observations in the previous two chapters suggest to us that computers are good tools for studying systems in which errors grow exponentially with time. Although all of the work in this thesis has been on discrete-time dynamical systems, similar ideas can be applied to ordinary differential equations. Both the G H Y S method and the C P method can be extended to ordinary differential equations. One difference in the G H Y S method is the the Jacobian J{ becomes the Jacobian of the integral of the ODE. A method similar to the CP method is given by Van Vleck [9] for ordinary 62 differential equations. The discretization error caused by the integration technique is much larger than the errors caused by floating point arithmetic. In [9], Van Vleck studies the discretization error and how it accumulates as time progresses. He then proves a shadowing theorem which involves the numerical method used to solve the ODE. The proof of the shadowing theorem in [9] requires a bound on the norm of an inverse of a linear operator. This linear operator has the same form as the linear operator in Sections 2.2.1 and 2.3.1. Although shadowing seems to work well for many pseudo-orbits, there are still many for which we cannot find a true orbit. Despite our efforts to prove that a shadow exists we still cannot show when a shadow does not exist. It is just as important to be able to determine when a shadow does not exist. For the logistic map we have shown how a round-off error introduces new behavior into the system when the pseudo-orbit is close to the critical point. For simple one-dimensional examples glitches are not too difficult to spot, but for higher-dimensional maps it may not be easy to determine where we have a glitch and why. How can we recognize glitches? How can we deal with them? When do glitches occur? These are many questions which have not been looked into in this thesis. Some work has been done in [5] to try to answer these questions. In this paper the authors show that the unshadowability of trajectories can be caused by Lyapunov exponents of the system fluctuating about zero. Much work has been done on the application of the GHYS refinement algorithm to the N-body problem. In [8], Quinlan and Tremaine show that the refinement procedure fails when a close encounter of two bodies occurs. This suggests that the pseudo-trajectory cannot be shadowed. If indeed this is the case, perhaps we can show that glitches which occur in other systems arise from different physical properties of the system such as energy, momentum, etc. 63 For the two-dimensional examples studied in this thesis, the stable and unstable manifolds for the pseudo-orbit were required to calculate the shadowing distance. The method used to calculate these vectors is similar to the power method for determining the eigenvector which corresponds to the largest eigenvalue. In order to determine the stable directions we need to iterate the method backwards which requires the knowledge of the entire pseudo-orbit. If we could determine the stable and unstable unit directions while we calculate the pseudo-orbit perhaps more efficient algorithms could be written to determine shadowing distances. Perhaps an eigenvector decomposition of the Jacobian matrix would be a better method to try. The G H Y S Refinement procedure has many similarities as a Newton's method pro-cedure for finding roots to an equation. If you recall the GHYS refinement procedure is like a function g which takes the pseudo-orbit and produces the one-step errors, ie. g{p) = e. Since the one-step error are assumed to be small, |e| is small. Hence p is close to a zero of g. And a zero of the function would represent a true orbit. Could it be that the G H Y S procedure is a Newton's method in disguise? They have many similarities: The resolvent is a Jacobian; the method is iterative; there is geometric convergence. If the G H Y S procedure could be arranged to look exactly like a Newton's method, then the massive amount of literature and code devoted to Newton's method may be brought to bear on the problem. 64 Bibliography [1] D. Anosov, Geodesic flows and closed Riemannian manifolds with negative curvature. Proceedings of the Steklov Institute of Mathematics, 90 (1967), pp. 209-235. [2] R. Bowen, co-limit sets for axiom A diffeomorphisms. Journal of Differential Equa-tions, 18 (1975), pp. 333-339. [3] S. Chow, K . J . Palmer, On the numerical computation of orbits of dynamical systems: the one dimensional case. Journal of Dynamics and Differential equations, 3 (1991), pp. 1527-1530. [4] S. Chow, K . J . Palmer, On the numerical computation of orbits of dynamical systems: the higher dimensional case, Journal of Complexity, 8 (1992), pp. 398-423. [5] S. Dawson, C. Grebogi, T. Sauer, J.A. Yorke, Obstructions to shadowing when a Lyapunov exponent fluctuates about zero. Physical Review Letters, 73(14) (1994), pp. 1927-1930. [6] C. Grebogi, S.M. Hammel, J.A. Yorke, Do numerical orbits of chaotic dynamical processes represent true orbits? Journal of Complexity, 3 (1987), pp. 136-145. [7] C. Grebogi, S.M. Hammel, J.A. Yorke, T. Sauer. Shadowing of physical trajectories in chaotic dynamics: containment and refinement. Physical Review Letters. 65(13) (1990), pp. 1527-1530. [8] G.D. Quinlan and S. Tremaine. On the Reliability of Gravitational N-Body Integra-tions. Monthly Notices of the Royal Astronomical Society, 259 (1992), pp. 505-5180. [9] E.S. Van Vleck, Numerical shadowing near hyperbolic trajectories. S IAM Journal on Scientific Computations, 16(5) (1995), pp. 1177-1189. [10] J .H. Wilkinson, Rounding Errors in Algebraic Processes. Prentice-Hall, Englewood Cliffs, NJ . (1963). 65
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Verifying computer solutions to discrete-time dynamical...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Verifying computer solutions to discrete-time dynamical systems Urminsky, David John 2000
pdf
Page Metadata
Item Metadata
Title | Verifying computer solutions to discrete-time dynamical systems |
Creator |
Urminsky, David John |
Date Issued | 2000 |
Description | Chaotic dynamical systems exhibit sensitive dependence on initial conditions. Round-off errors introduced in computer simulations may cause a computed orbit to quickly diverge from the true orbit. One may therefore question the validity of these computations. Shadowing provides a means for studying the validity of a computation. If we are able to show that a true solution with a different initial condition, called a shadowing orbit, stays close to a computed solution, which we call a pseudo-orbit, then perhaps the computed solution has some validity. We will investigate two methods for proving the existence of a shadowing orbit. The first method comes from two theorems by Chow and Palmer. These theorems provide us with sufficient conditions for when a pseudo-orbit is shadowed by a true orbit and provide us with a bound on the shadowing distance of a true orbit. The second method is by Grebogi, Hammel, Yorke and Sauer. Their method contains a pseudo-orbit in a carefully constructed sequence of parallelograms which helps to prove the existence of a shadowing orbit. These two methods will then be used to prove the existence of a shadowing orbit for examples in one and two dimensions. Finally we will discuss a refinement technique which takes a 'noisy' orbit and produces a less 'noisy' orbit which shadows the original orbit. This technique will then be applied to some examples to find a numerical shadow for a pseudo-orbit. |
Extent | 2693264 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0080030 |
URI | http://hdl.handle.net/2429/10521 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0301.pdf [ 2.57MB ]
- Metadata
- JSON: 831-1.0080030.json
- JSON-LD: 831-1.0080030-ld.json
- RDF/XML (Pretty): 831-1.0080030-rdf.xml
- RDF/JSON: 831-1.0080030-rdf.json
- Turtle: 831-1.0080030-turtle.txt
- N-Triples: 831-1.0080030-rdf-ntriples.txt
- Original Record: 831-1.0080030-source.json
- Full Text
- 831-1.0080030-fulltext.txt
- Citation
- 831-1.0080030.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080030/manifest