%%% Formated for "World Scientific Publishing Co" \magnification=\magstep1 %*************** TO GET SMALLER FONT FAMILIES ***************** \newskip\ttglue % ********** EIGHT POINT ************** \def\eightpoint{\def\rm{\fam0\eightrm} \textfont0=\eightrm \scriptfont0=\sixrm \scriptscriptfont0=\fiverm \textfont1=\eighti \scriptfont1=\sixi \scriptscriptfont1=\fivei \textfont2=\eightsy \scriptfont2=\sixsy \scriptscriptfont2=\fivesy \textfont3=\tenex \scriptfont3=\tenex \scriptscriptfont3=\tenex \textfont\itfam=\eightit \def\it{\fam\itfam\eightit} \textfont\slfam=\eightsl \def\sl{\fam\slfam\eightsl} \textfont\ttfam=\eighttt \def\tt{\fam\ttfam\eighttt} \textfont\bffam=\eightbf \scriptfont\bffam=\sixbf \scriptscriptfont\bffam=\fivebf \def\bf{\fam\bffam\eightbf} \tt \ttglue=.5em plus.25em minus.15em \normalbaselineskip=9pt \setbox\strutbox=\hbox{\vrule height7pt depth2pt width0pt} \let\sc=\sixrm \let\big=\eightbig \normalbaselines\rm} \font\eightrm=cmr8 \font\sixrm=cmr6 \font\fiverm=cmr5 \font\eighti=cmmi8 \font\sixi=cmmi6 \font\fivei=cmmi5 \font\eightsy=cmsy8 \font\sixsy=cmsy6 \font\fivesy=cmsy5 \font\eightit=cmti8 \font\eightsl=cmsl8 \font\eighttt=cmtt8 \font\eightbf=cmbx8 \font\sixbf=cmbx6 \font\fivebf=cmbx5 \def\eightbig#1{{\hbox{$\textfont0=\ninerm\textfont2=\ninesy \left#1\vbox to6.5pt{}\right.\enspace$}}} %************** NINE POINT ***************** \def\ninepoint{\def\rm{\fam0\ninerm} \textfont0=\ninerm \scriptfont0=\sixrm \scriptscriptfont0=\fiverm \textfont1=\ninei \scriptfont1=\sixi \scriptscriptfont1=\fivei \textfont2=\ninesy \scriptfont2=\sixsy \scriptscriptfont2=\fivesy \textfont3=\tenex \scriptfont3=\tenex \scriptscriptfont3=\tenex \textfont\itfam=\nineit \def\it{\fam\itfam\nineit} \textfont\slfam=\ninesl \def\sl{\fam\slfam\ninesl} \textfont\ttfam=\ninett \def\tt{\fam\ttfam\ninett} \textfont\bffam=\ninebf \scriptfont\bffam=\sixbf \scriptscriptfont\bffam=\fivebf \def\bf{\fam\bffam\ninebf} \tt \ttglue=.5em plus.25em minus.15em \normalbaselineskip=11pt \setbox\strutbox=\hbox{\vrule height8pt depth3pt width0pt} \let\sc=\sevenrm \let\big=\ninebig \normalbaselines\rm} \font\ninerm=cmr9 \font\sixrm=cmr6 \font\fiverm=cmr5 \font\ninei=cmmi9 \font\sixi=cmmi6 \font\fivei=cmmi5 \font\ninesy=cmsy9 \font\sixsy=cmsy6 \font\fivesy=cmsy5 \font\nineit=cmti9 \font\ninesl=cmsl9 \font\ninett=cmtt9 \font\ninebf=cmbx9 \font\sixbf=cmbx6 \font\fivebf=cmbx5 \def\ninebig#1{{\hbox{$\textfont0=\tenrm\textfont2=\tensy \left#1\vbox to7.25pt{}\right.$}}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \overfullrule=0pt \hsize=6.0truein \vsize=8.5truein \parindent=3truepc %%%% NO PAGE NUMBERS ON FINAL COPY !!!!!!!!!!!!!!!!!!!! %\nopagenumbers %%% Added macros (combs) % for Remarks, Proofs, Definitions, etc. \def\demo#1.{\medskip \noindent {\it #1.}\enspace} \outer\def\exercise#1#2{\medskip\noindent{\it Exercise #1.} \enspace{\rm#2} \medskip} \def\section#1{\vskip.6truecm\noindent{\bf#1}\vskip.4truecm} \def\qed{~\hfill\eop\medskip} \def\eop{\hbox{\vrule width6pt height7pt depth1pt}} \def\cite#1{$^{#1}$} \def\myskip{\noalign{\vskip6pt}} % extra space multi-line eqns. \def\frac#1#2{\textstyle{#1\over#2}} \def\bigmid{{\ \Big|\ }} \def\un#1{{\underline{#1}}} \def\iitem{\itemitem} % \def\balpha{{\underline{\alpha}}} \def\balpha{\mathop{\alpha\kern-.6em{\alpha}}\nolimits} % poormans' bold \def\bpsi{\mathop{\psi\kern-.73em{\psi}}\nolimits} % poormans' bold \def\bvarphi{\mathop{\varphi\kern-.71em{\varphi}}\nolimits}% poormans' bold \def\bomega{\mathop{\omega\kern-.59em{\omega}}\nolimits} \def\Ball{\mathop{\rm Ball}\nolimits} \def\diam{\mathop{\rm diam}\nolimits} \def\dist{\mathop{\rm dist}\nolimits} \def\Id{\mathop{\rm Id}\nolimits} \def\Im{\mathop{\rm Im}\nolimits} \def\mod{\mathop{\rm mod}\nolimits} \def\Vol{\mathop{\rm Vol}\nolimits} \def\triv{{|\!|\!|}} % ||| \def\varep{\varepsilon} \def\nat{{\bf N}} % natural \def\real{{\bf R}} % real \def\zed{{\bf Z}} % integer \def\bzero{{\bf 0}} \def\bA{{\bf A}} \def\bG{{\bf G}} \def\bY{{\bf Y}} \def\ba{{\bf a}} \def\bb{{\bf b}} \def\bk{{\bf k}} \def\bp{{\bf p}} \def\bq{{\bf q}} \def\bx{{\bf x}} \def\bz{{\bf z}} \def\A{{\cal A}} \def\B{{\cal B}} \def\C{{\cal C}} \def\F{{\cal F}} \def\G{{\cal G}} \def\H{{\cal H}} \def\L{{\cal L}} \def\R{{\cal R}} \def\S{{\cal S}} %%% End of added macros %\headline={\ifnum\pageno=1\firstheadline\else %\ifodd\pageno\rightheadline \else\leftheadline\fi\fi} %\def\firstheadline{\hfil} %\def\rightheadline{\hfil} %\def\leftheadline{\hfil} % \footline={\ifnum\pageno=1\firstfootline\else\otherfootline\fi} %\def\firstfootline{\rm\hss\folio\hss} %\def\otherfootline{\hfil} %\font\tenbf=cmbx10 %\font\tenrm=cmr10 %\font\tenit=cmti10 %\font\elevenbf=cmbx10 scaled\magstep 1 %\font\elevenrm=cmr10 scaled\magstep 1 %\font\elevenit=cmti10 scaled\magstep 1 %\font\ninebf=cmbx9 %\font\ninerm=cmr9 %\font\nineit=cmti9 %\font\eightbf=cmbx8 %\font\eightrm=cmr8 %\font\eightit=cmti8 %\font\sevenrm=cmr7 %\TagsOnRight %%%% NO PAGE NUMBERS ON FINAL COPY !!!!!!!!!!!!!!!!!!!! %\nopagenumbers %\line{\hfil } %\vglue 1cm \topinsert\vskip1truecm\endinsert \centerline{\ninebf INTRODUCTION TO THE THEORY} \vskip.1in \centerline{\ninebf OF COMPUTATIONAL DYNAMICS} \vskip 1.0truecm \centerline{\ninerm RAFAEL DE LA LLAVE} \vskip 1pt \centerline{\nineit Mathematics Department, The University of Texas at Austin} \centerline{\nineit Austin, Texas 78712-1082, USA} \vskip 0.8truecm %\centerline{\ninerm ABSTRACT} %\vskip 0.3truecm %{\narrower\noindent % \ninerm %%%%% put text here if there is an abstract % \noindent %\vskip 0.8truecm } \baselineskip=14pt \section{1. Introduction} It is a basic philosophical principle that for many physical systems, the state at a certain time determines the state one unit of time later. If the laws governing the system are independent of time we can write a function $f$, so that, if a system is in an state $x$ at time zero, $f^n (x) = \underbrace{f(f\cdots f}_{n\hbox{\sevenrm -times}} (x))$ will be the state occupied at $n$ times the basic unit of time. The sequence of $f^n(x)$ is usually called the orbit of $x$ and represents the states that the system will occupy at successive times. The map $f$ is called a discrete dynamical system in the Mathematical literature. If such a system were given to us in a practical problem (e.g., the solar system, whose laws of evolution are well known) there are some natural questions to ask. For example, we may inquire if there are certain types of orbits (e.g., starting near the Earth and going close to all the planets, remaining close to the Earth but not drifting away due to the periodic pull of other planets and moon) or, in a more philosophical mood we may inquire whether all the orbits are going to keep their shape forever except for small oscillations. For the solar system, the later is a non-trivial question, the planets do eject and capture small bodies such as meterorites all the time. Could it happen that the Earth captures the moon due to the influences of the other planets? Which is the time scale of such phenomena? This type of questions, which even now are not completely answered, were among the first posed to Mathematical Physics. To the previous questions, we could also add the worry of what would happen if the laws of evolution were slightly modified. Of course, in a physical system there is always the possibility that we have neglected some effect. More importantly, when we introduce a system in a computer or treat it theoretically, we make modifications to it to accommodate the limitations of the computer or the methods. Does this invalidate the conclusions?. Notice for example, that even if a harmonic oscillator has periodic solutions, the introduction of even the smallest friction drives it to equilibrium. Let us start by discussing away two naive approaches. The first is to try to find a solution in ``closed'' form and the second is to observe that since $f$ can be written in a computer we can easily iterate any initial condition we like. To the first we just answer that finding such closed form solutions turns out to be impossible for many dynamical systems (see e.g., \cite{1} for a discussion of impossibility of finding closed form solution of $\ddot y +y-y^2 +f(t)=0$). Even if we found those solutions, given that many of the questions we are interested in are qualitative in nature, it is not clear that we would profit much from such solutions as long formulas. For the second approach we observe that when the number of parameters is even as small as 6 --- 3 positions, 3 velocities --- a very crude examination of 1\% precision in each coordinate (nowhere near enough the precision required for an orbit of a satellite) would require $10^{12}$ orbits which is well beyond the reach of today's fastest computers. Even if by some remarkable progress in computer technology --- it is a very parallelizable problem --- we got to the level where enough orbits could be generated, we would still have the problem of selecting those that perform the motions we want or more difficult one of realizing that some of motions which we are computing may be useful for some goal. Neither of these two naive methods offers any clue on how to start considering the effects caused by studying modifications of the system. The fact that simple methods do not work is not a coincidence. Many people heuristically identify deterministic with predictable and predictable with simple. This is merely an abuse of language. It is becoming increasingly recognized that even very simple systems can have very complicated behaviour. Naive iteration of many initial conditions can make one realize very quickly how complicated things can be. The best answer to the type of questions that we have asked seems to be a program that could be attributed to Poincar\'e. It roughly, proceeds along the following steps \item{$\ast$} Find landmarks that organize the long term behaviour. For example: \iitem{$\bullet$} periodic orbits \iitem{$\bullet$} invariant tori \iitem{$\bullet$} attractors \iitem{$\bullet$} bifurcations \item{$\ast$} Perform a local analysis around those landmarks so that we control a non-negligible part of phase space. \iitem{$\bullet$} normal forms \iitem{$\bullet$} local stable manifolds and foliations \item{$\ast$} Join all the local information together into global structures \iitem{$\bullet$} index theory \iitem{$\bullet$} global bifurcation \iitem{$\bullet$} basins of attraction Part of the beauty of this program is that in it, one can make progress using at the same time theoretical analysis and numerical exploration. If the observer has the right theoretical training, even the naive direct iteration of random points may suggest the presence of features or landmarks which can then be computed precisely with more detailed numerical methods. After some more work, one can use some global topological argument to discover that the features already found imply the existence of another one, which once it is found and its details elucidated, implies some others, etc.. After a few iterations of the process, one may hope to arrive to an understanding sufficient for the applications. In this lecture, we would like to give some of the flavor of these techniques, emphasizing as corresponds to the title of the school the computational techniques that can be used to find the objects we describe theoretically. Given the limitations of time and space, we can only hope to discuss the most elementary phenomena and be just a guide to the literature. In particular, we will not discuss much about Global Behaviour. To a large extent, the global phenomena is not understood systematically yet, hence it is one of the most interesting areas of research and it is certainly one of the areas in which computer explorations is having a larger effect. Good introductions to some global problems from the point of view of applications exist in the literature \cite{2} \cite{3}. It is also instructive to see how this methods can be brought to bear on a concrete example \cite{4}. With the availability of powerful computers at affordable prices it has become possible to develop programs that display the graphical information in real time and allow the user to perform some of the standard diagnostics. Two such programs that are particularly appropriate for the level of these lectures are DSTOOL and DYNAMICS which are available with the source code. Reading the code of these programs can be a very instructive process that we recommend. \demo Exercise 1. Take one of the interactive packages mentioned in the text and run at least three examples observing the behaviour of different orbits. \demo Remark. As a historical note we may add that topology --- under the name of ``Analysis situs'' --- was developed by Poincar\'e as part of this program. This included simplicial homology and many other index theories. Even now topology is defined sometimes as the science of passing local information to global one. \demo Remark. It goes without saying that this program we have outlined is not the only one that is useful to study dynamics. (Not even if one wants to study the systems computationally.) We just mention the ``ergodic program'' or the Smale's generic program. In ergodic theory, we try to study the statistical properties of orbits. In the topological classification program one tries to identify properties that are ``typical'' and obtain a topological classification that applies to all systems. There are already good books and review papers \cite{5} \cite{6} \cite{7} \cite{8} about these subjects. In particular \cite{7} \cite{8} contain a very modern point of view that uses ergodic theory and the geoetric program at the same time. \demo Remark. Very often the laws of dynamics appears given as a differential equation in an space of finite dimensions. Standard theory of ordinary differential equations\cite{9} shows that under very mild regularity assumptions, the two formulations are equivalent. There is a well developed branch of numerical analysis devoted to the computation of evolution maps out of the differential equations. In this lecture we will emphasize maps. The issue of how to approximate differential equations by maps -- the only thing that is easy to implement on the computer -- is a large field of numerical analysis. Up to date reviews are the books \cite{10} \cite{11}. \section{2. Reliability of computer exploration.} Unfortunatley, if one is working with computers, all computer operations involve approximations. Computers work with a finite numbner of digits and, hence, all the operations are approximate. The most naive answer to this problem is just to increase the precision till we can obtain sufficiently many significant digits. This is not a satisfactory answer on two counts. First of all, it is very difficult to analyze theoretically what is the precision lost at every iteration. Second, it is true that in many cases the number of figures lost is proportional to the number of iterations. When one has to perform operations with numbers which are not supported directly by the hardware the efficiency drops very quickly. There are indeed methods to work with computers in such a way that one can get definitive answers to numerical computations. More important for dynamical systems is the fact that many systems satisfy the so-called shadowing property. That is: \demo Definition. We say that a system $f$ has the shadowing property whenever given a sequence of points $\{x_n\}$ such that $$|f(x_n) - x_{n+1}|\le\varep \quad (\varep \le\varep_0)$$ then, there exist a point $y$ such that $$|f^n (y)-x_n| \le\delta (\varep)\quad \bigl(\delta (\varep)\to0 \ \hbox{ when }\ \varep\to0\bigr)$$ A sequence $x_n$ such that $|f(x_n)-x_{n+1}|\le \varep$ is usually called an $\varep$ pseudo-orbit. Notice that a computer is perfectly capable of producing $\varep$ pseudo-orbits with $\varep$ about the round-off error. This result is quite satisfactory from the point of view of an analysis based --- as Poincar\'e program --- in the qualitative properties of some orbits. Nevertheless, it is important not to read too much from it. For example, it is not clear at all if properties considered ``typical'' on a system will remain typical in the system of pseudo-orbits produced by introducing approximations. Hence, this result is somewhat useless from the point of view of the ergodic program. For the so-called Anosov systems --- which we will define later --- in two dimensions, it is possible to show that if two of them assign positive probabilities to the same observables both of them differ just by a change of variables as smooth as they are. \demo Example. $x\to 2x$ mod 1. This system satisfies the shadowing property. It is possible to show that a typical orbit of this system is equally distributed over $[0,1)$. Nevertheless, on a binary machine, all computed orbits end at zero. \exercise{1} {Compute reliable 10 figures of $f^{100} (0.3)$ where $f(x) = 1-2x^2$.} \exercise{2} {The magnification of the error in one data after $N$ iterations can be estimated by $|{d\over dx} f^N(x)|$. For one dimensional maps this can be computed by $$\ln \left| {d\over dx} f^n(x)\right| = \ln \prod_{i=0}^{N-1} |f'(f^i (x))| = \sum_{i=0}^{N-1} \ln |f' (f^i (x))|\ .$$ \vskip.1in Compute this estimate of the error for several functions of the form $$f_\lambda = 1-\lambda x^2 \qquad \lambda \in [0,2]\qquad \hbox{(Take $N=500$, $x=0.3$)}$$} \exercise{3} {How does the previous result change if you take different $x$? What happens if you take a larger $N$?} \exercise{4} {Consider $f(x) = 2x$ mod 1. Show that if $f(x_n) - x_{n+1}=\varep_n$, then $y=x_0 +\sum (\varep_n/2^n)$ is a shadowing orbit.} If rather than starting with initial conditions $x$, we start with initial conditions $x+\varep e$, after $N$ iterations $f^N (x+\varep e)\approx f^N (x) +\varep Df^N (x) v$ so that $$\lim_{N\to\infty} {1\over N} \ln \bigl( \|Df^N(x) v\|/\|v\|\bigr) = \Lambda (x,v)$$ is an adequate measure of the effect of an error in the initial conditions. Even if it is not true that these limits always exist, a deep theorem of Osledec's --- see \cite{Ru} for a more modern proof --- tells us that for a ``typical'' initial condition $x$, these limits exist for all vectors. Moreover, if we define $E_\lambda = \{v\mid \Lambda (x,v) \le \lambda\}$ we can see it is a linear space $E_\lambda \subseteq E_{\lambda+\delta}$, $\delta>0$, so that we can except to have numbers $\lambda_1,\ldots,\lambda_k$ ($k\le$ dimension) that appear. Such numbers are called the Lyapunov exponents. Notice also that $\Lambda (f(x),Df(x)v) = \Lambda (x,v)$ so that the set --- and multiplicity --- of Lyapunov exponents is invariant under the map. In many cases this implies that the exponents are constant over large region of phase space. (This is related to the ``ergodic'' hypothesis for this region.) Intuitively, positive Lyapunov exponents mean that small causes will have large effect -- growing exponentially in time -- and, hence, we expect very widely varying types of long term behaviour intermingled. This intuition can be made precise and quantitative and there are theorems that relate the Lyapunov exponents to the ``entropy'' --- a measure of the complication of the set of orbits --- and to the dimension of the attractors --- a measure of the complication of the geometry of the set described by the orbits\cite{13}. Making all this more precise is a very active area of research now. For a typical orbit, we can find a basis $v_1,\ldots,v_n$ in such a way that $\|Df^N (x)v_i\| \approx e^{\lambda_i N} C$ both for $N\to \pm\infty$ (several $\lambda_i$'s here could be the same). This suggest that we can compute the largest Lyapunov exponent by a method very similar to the ``power method'' to compute the modulus of the largest eigenvalue of a matrix. We pick a vector $u_0$, $|u_0|=1$ and recursively set $u_{n+1} = Df(f^n(x))u_n$. Unless $v_0$ does not have any component on the direction of the leading eigenvalue, we will have $\lambda_{\max} \approx {1\over N} \ln \|u_N\|$. In practice, even the tiniest round-off errors will take us to this situation. Computation of other Lyapunov exponents is much harder. A good algorithm --- used in \cite{14} to give an alternate proof of Osledec's theorem but discussed in detail in discussed in \cite{7} --- is: set $Df(x) = Q_0R_0$ where $a_0$ is orthogonal $R_0$ is upper triangular (this is the well known $QR$ factorization of linear algebra, implemented reliably in many good numerical packages such as EISPACK). Then recursively set $Df(f^i) Q_{i-1}= Q_iR_i$. It is clear that then we have $Df (f^{N-1}x) \cdots Df^n(x)= Q_N R_{N-1} \cdots R_0$. For low dimensions, there are shortcuts. For example, in dimension 2, just notice that $\det Df^N(x) = \det f(f^{N-1}x) \cdots \det f(x)$ and that the rate of growth of the determinant should be the sum of the two Lyapunov exponents. We also point out that the concept of Lyapunov exponent can be generalized for approximate orbits. One can show that for approximate orbits with non-zero Lyapunov exponents can be shadowed (we will discuss this in more detail for periodic orbits). If a system has abundance of true orbits with non-zero Lyapunov exponents, these conditions will be satisfied. This seems paradoxical because we cannot compute well individual trajectories, we should believe qualitatively what we compute! Nevertheless (see exercise~4) what we are saying is that for these chaotic systems, anything can happen, in particular what we computed! \section{3. Periodic orbits} Periodic orbits are the simplest landmarks to compute. On the other hand, many features of dynamical systems can be approximated by periodic orbits. This was already predicted by Poincar\'e \cite{15} Vol I p. 82. \vskip 2 em {\it \narrower Il y a m\^eme plus: voici un fait que je n'ai pu demontrer rigorousment, mais qui me parait pourtant tr\`es vrisemblable. \narrower \`Etant donn\'ees des \`equations de la forme definie dans le ${\rm n}^{\rm o}$ 13 [ Hamilton eq. ] et unne solution particuli\`ere quelconque de ces equations, on peut toujours trouver une solution p\'eriodique (dont la p\'eriode peut, il est vrai,\^etre tr\`es longue), tel que la diff\'erence entre les deux solutions soit aussi petite qu'on le veut, pendant un temps, assi long qu'on le veut. D'allieurs, ce qui nous rende ces solutions p\'eriodiques si precieuses, c'est qu'elles sont, pour ansi dire, la seule br\`eche par o\`u nous puissions essayer de p\'en\'etrer dans une place jusqu'ici r\'eput\'ee inabordable. } The simplest method to compute periodic orbits is to search for solutions of $f^N (x)=x$ using, e.g., a Newton method or any other non-linear equation solver. As mentioned before, the numerical calculation of the derivative becomes delicate since numerical errors tend to line up all the vectors into the dominant Lyapunov exponents. A more reliable method is based on the study of the operator acting on sequences $\bx = \{x_i\}_{i=1}^N$ $$\C (\bx) = \bigl( f(x_N),f(x_1),\ldots,f(x_{N-1})\bigr)$$ We just point out that $|\C (\bx) - \bx | \le \varep$ (with the supremum norm) is the same as saying that $\bx$ is an $\varep$ pseudo-orbit. If one could show that $D\C (x)- \Id$ is invertible and $\varep$ is sufficiently small, then, a Newton method would converge to a true solution. This is a version of the shadowing theorem mentioned before. Notice that $D\C (x)-\Id$ is a very sparse matrix. If one picks an adequate solver for sparse matrices, it could be that the method is better conditioned than the naive one. As it turns out $D\C(x)-\Id$ is invertible when and only when the pseudo-orbit has no zero Lyapunov exponents. Hence, if there are recurrent orbits and non-zero Lyapunov exponents, there should be periodic orbits close by. That is, periodic orbits appear as tracers of chaotic behaviour. (A precise version of this are the ``ergodic'' closing lemmas\cite{16}\cite{17}\cite{18}.) For many systems that have special structures, one can use more effective methods. For example, the symplectic integrators considered in Professor Jauslin's lectures when applied to Hamiltonians $H= {p^2\over2} + V(q)$ lead to maps of the form $$\eqalign{ p_{n+1} & = p_n - h\nabla V(q_n)\cr q_{n+1} & = q_n + p_{n+1}\cr} \eqno(1)$$ One can see that these equations are variational equations for: $$S(\bq) = \sum_{i=1}^N (q_{n+1} - q_n)^2 \frac12 + hV(q_n)$$ (where $q_{N+1} = q_1$). This variational principle is rather well behaved and using a minimization algorithm it is possible to obtain good solutions even for fairly long periods. Of course, the existence of this variational principle is a discrete version of Hamilton's stationary action principle. For maps $f$ that have the property that $f^{-1} = I\circ f\circ I$ with $I^2=\Id$, usually called reversible --- very often in Physics this arises when $I$ is just reversing the velocity and keeping the positions--- one can find a very effective method to compute periodic orbits Notice that $f= I(If)$ and that $(If) (If)=\Id$ so that we can write $f=IJ$ where $I^2 = J^2 = \Id$. \demo Exercise. Compute explicitly the maps $I,J$ for (1). \vskip 1 em The following algorithm can be found in \cite{19}. \proclaim Algorithm. Let $\Sigma$ be the set of points $x$ such that $J(x) = x$. If $x\in\Sigma$ and $f^N (x)\in\Sigma$ then $f^{2N} (x)=x$. \demo Proof. $$\eqalign{f^N(x) & = (IJ)^{N-1} IJ(x) =\cr &= J(IJ)^{N-1} I(x) =\cr &= f^{-N} (x)\cr}$$ where in the first equality we have used $J(x) =x$, $Jf^N(x)=f^N(x)$ and in the second we have used $J^{-1} =J$, $I^{-1}=I$. Notice that $\Sigma$ is a manifold of half the dimension of the space and the equation $f^N(x)\in\Sigma$ amounts to a system of equations of half the dimension. In particular, for 2-dimensional reversible systems, this leads to solving 1-equation with one unknown which is quite easy since in one dimension finding zeros of a function can be done very systematically exploiting the fact that a change of sign implies the existence of a zero in between. For a complete discussion of algorithms for periodic orbits for symplectic and reversible mappings we refer to \cite{20}. The computation of periodic orbits for differential equations is essentially more difficult than for discrete systems since the periodic also has to be computed. Moreover, all the points in the periodic orbit are periodic points with the same period so that the solutions are never isolated, i.e., a Newton method will not work. As an added complication, computing the derivatives of the time $T$ map --- needed for a Newton method --- requires to integrate the so-called variation equation (see \cite{9}). A method which is quite analogous to the fixed point method we discussed for maps is the so-called ``multiple shooting method.'' One tries to compute a vector of points in which we impose that each one is in the image of the next. A convenient way of fixing the period is to impose that some point belongs to a transversal section of the flow and the last step is not a fixed time, but rather a projection on the transversal surface. (See \cite{21} for a discussion of the multiple shooting method for boundary value problems.) A good numerical integration that includes study of the variational equations is ODESSA. For Hamiltonian systems one still Lagrange's variational principle for periodic orbits: $$\delta \int_0^T (\dot q,q)\,dt = 0$$ among all functions which satisfy $q(t+T) = q(t)$. For Hamiltonian systems, one also has that the energy is a conserved quantity and often the orbits appear in families for which the period depends on the value of the energy. That is, we can prescribe the period of any preassigned value. If we discretize the space of periodic functions --- e.g., using Fourier series --- $$q(t) = \sum \hat q_k e^{2\pi ikt/T}$$ we can transform $$\int_0^T L(\dot q,q)\,dt = \S (\hat q_k)$$ which when truncated to a finite number of modes becomes a function of finitely many variables whose critical points can be computed by conventional methods. See \cite{22,23} for a more complete discussion of these variational methods for Hamiltonian systems. Notice also that, once we know that we can prescribe the period, studying a discrete version of the system (1) may produce a solution which is approximate up to the truncation error of the integration (related to the step $h$) --- one needs precautions against the problem of multiple solutions. If we require a solution with a prescribed energy, e.g., if we are computing several families of orbits to study their interactions, we may want to refer them to the same energy we still need to study the function $E(T)$ and invert it. \demo Exercise. For equations of the form $y'' = f(y)$, denote by $y_n = y(nh)$ and analogously for other functions notice that $$\eqalign{ &y_{n+1} + y_{n-1} - 2y_n = h^2 y''_n + {h^4 \over 12} y_n^{(iv} + o(h^6)\cr &y''_{n+1} + y''_{n-1} -2y''_n = h^2 y_n^{(iv} + o(h^4)\cr}$$ Show that $$\left( y_{n+1} - {h^2\over 12} f(y_{n+1})\right) + \left( y_{n-1} - {h^2\over 12} f(y_n)\right) = 2\left( y_n -{h^2\over 12} f(y_n)\right) = f(y_n)$$ is equivalent to the original equation up to errors of order $h^6$ and that if we consider the $y_n$'s as unknowns, this gives a system with sparse coefficients that can be solved efficiently. (This is a non-linear version of the so-called Numerov method for Schr\"odinger equations.) \vskip 2 em Before computers were available, one of the few methods to compute periodic orbits were perturbative methods. The most natural one for is the so-called Lindstedt-Poincar\'e method. We will discuss it in detail for the case of differential equations, where one has to guess not only the solution but also the period. The basic idea is that we search for a solution of the differential equation $ y_\varep(y) = u_\varep (T_\varep t)$ where $u_\varep = u_0 +\varep u_1 + \varep^2 u_2 + \cdots$, \ \ $u_\varep (t+2\pi) = u_\varep (t)$, \ \ $T_\varep = T_0 +\varep T_1 +\varep^2T_2 +\cdots$. We substitute in the original equations and require that the equations have solutions at all orders, the coefficients in the expansion of $T_\varep$ will appear as compatibility conditions in the expansion of $u$. For example, if we consider $$y'' = y+\varep N(y)$$ (a non-linear oscillator with friction). The Lindstedt-Poincar\'e equations become: $$(T_\varep)^2 u''_\varep (t) = -u_\varep (t) +\varep N\bigl(u(t)\bigr)$$ At order zero we have $$(T_0)^2 u''_0 = -u_0$$ This equation says $u_0 = A\sin (T_0 t+\varphi)$. We will set $\varphi =0$ since it physically corresponds to choosing the origin of time. We also obtain $T_0=1$ from the periodicity conditions. At order 1, we have $$2T_1 u_0 +u''_1 = -u_1 +N(u_0)$$ If we consider this as a (linear!) equation for $u_1$ we see that it is forced with the term $N(u_0) - 2Tu_0$. If this forcing term contained a term of the form $\sin (t)$ or $\cos (t)$ the theory of linear systems with constant coefficients tells us that the solution should contain a term of the form $Ct\sin (t)$ or $Ct \cos (t)$. This would violate the condition that the solutions are periodic. If we impose that these resonant forcings do not appear, this amounts to conditions on $A,T_1$. If these equations can be solved, then, we proceed to next order. It is by no means obvious that these equations can be solved. Indeed for $$y'' = -y - \varep (y')^3$$ they cannot. Physically, this is a system with friction and, indeed, as soon as friction appears and is not balanced by an external force, periodic motion becomes impossible. Nevertheless, it is quite frequently the case that, after a few order, the recursive equations develop some structure that makes it clear that they can be solved for all higher orders. This sometimes can even be proved using more elaborate theory. We just mention that there are symbolic programs that can perform these calculations \cite{24}. Usual numerical programs do not work very well since it is necessary to perform Taylor expansions, gather terms, etc. \demo Exercise. Compute by hand the Lindstedt-Poincar\'e series of to order 2 for the Van-der Pol oscillator $$y''=y +\varep (1-y^2) \dot y$$ \demo Exercise. Obtain one of the symbolic manipulation programs mentioned above and compute the same series. An complication on periodic orbits is quasi-periodic orbits. Recall that a quasi-periodic orbit is one that has an expansion: $$y(t) = \sum e^{2\pi i(k_1\omega_1+\cdots + k_n\omega_n)t} \hat y_{k_1\ldots k_n}$$ where $\omega_1,\ldots,\omega_n$ are independent frequencies. Unfortunately, the sense in which this expansion converges is quite delicate. We will just mention that the numerical computation of such orbits seems to be best done by transforming the equation for these orbits into an equation for the coefficients $\hat y_{k_1\ldots y_n}$. When $n>2$, they are practically intractable. Nevertheless, given the importance of quasi-periodic orbits, there is considerable effort done in their calculation. \section{4. Local theory near a fixed point} If we have a fixed point --- we will assume for simplicity it is the origin --- our goal is to obtain a detailed description of the dynamics of the map in a neighborhood of the origin. Since $f(x) = Df (0) x + o(|x|^2)$, it is natural to start by comparing the dynamics of the full map to that of the linear part. The dynamics of the linear part is very easy to study. Since $Df(0)=V^{-1}JV$ where $J$ is a (real) Jordan form we see $Df(0)^n = V^{-1}J^nV$ so that the long term behaviours in the dynamics of the linear approximations can be studied by studying those of Jordan blocks. The blocks corresponding to eigenvalues of modulus bigger than 1 will blow up exponentially fast, those corresponding to eigenvalues of modulus smaller than 1 will shrink exponentially fast and those corresponding to modulus 1 will remain bounded or blow up (polynomially) both in the future and in the past. Sample examples such as $f(x) = x+x^{2n+1}$, $f(x) = x-x^{2n+1}$ show that no matter how high order the perturbation, the dynamics of maps with eigenvalue 1 can be radically changed. We define a map to be hyperbolic if no eigenvalue has modulus 1. A very effective way to compare the dynamics of two maps is to find a change of variables that sends one into the other. The situation for hyperbolic maps is very clean thanks to the following theorem of Hartman (differential equation) and Grobman (maps). A modern proof can be found in \cite{25}. \proclaim Theorem. Let $f$ be a $C^1$ map on a hyperbolic fixed point. Then, there exists a continuous change of variables such that $$h^{-1} fh= Df (0)$$ in a neighborhood of the origin. The map produced the the Hartman-Grobman theorem is, unfortunately, not very computable and it is, in general not very smooth. To uderstand why it is instructive to try to compute it formally. This computation will also give hints of some other phenomena that we will need to study in the chapter on bifurcations. For simplicity, we will assume that the matrix $Df(0)$ is diagonal and will denote the components by subindices. Here we will have $$f^i (x) = \lambda^i x^i + f_{jk}^i x^i x^k + f_{jk\ell}^i x^j x^k x^\ell +\cdots$$ (No convention summation for $\lambda^i x^i$, but we use the convention for all the other indices.) If we set $h^i (x) = x^i + h_j^i x^j$ and try to solve $hf=Df(0)h$ up to third order we obtain $$(h_{jk}^i \lambda^i \lambda^k + f_{jk}^i - \lambda^i h_{jk})=0$$ so that we can solve the equation provided $\lambda^i - \lambda^j \lambda^i \ne 0$. When equalities such as this happen, they are usually called ``resonances'' and it is clear that, if we choose examples for which $f_{jk}^i \ne0$, it is impossible to obtain $h$ which are $C^2$ and which reduce the map to its linear part. (In \cite{26} one finds obstructions for $C^1$ maps!) The map $\tilde f = hfh^{-1}$ satisfies $$D\tilde f (0)= Df (0)\quad ,\quad \tilde f = Df(0) + o(|x|^3)$$ We will work by induction and show that if $$f^i (x) = \lambda^i x^i + f_{j_1\ldots j_n}^i x^{j_1}\cdots x^{j_n} + o(|x|^{n+1})$$ and if $\lambda^{j_1} \cdots \lambda^{j_n}\ne \lambda^i$ we can find $h^i (x) = x^i + h_{j_1\ldots j_n}^i x^{j_1}\cdots x^{j_n}$ in such a way that $$hf = Df(0) h + o(|x|^{n+1})$$ Hence $hfh^{-1}=\lambda^i x^i + \tilde f_{j_1\ldots j_{n+1}}^i x^{j_1}\cdots x^{j_{n+1}} + o(|x|^{n+2})$. The composition of all the changes of variables obtained by the induction process will achieve the desired reduction. If we expand $hf-Df(0)h$ we obtain that the coefficient of order $n$ is $$0 = h_{j_1\ldots j_n}^i \lambda^{j_1}\cdots \lambda^{j_n} + f_{j_1\ldots j_n}^i - \lambda^i h_{j_1\ldots j_n}^i$$ so that, we can solve provided that $\lambda^{j_1\cdots j_n} \ne \lambda^i$ and, if not, it is possible to find examples. A more complete treatment can be found in \cite{27}, \cite{28} among other places, but we recommend the reader that he/she tries the following exercises. \demo Exercise. Study the case where the matrix has a non-trivial Jordan block. Show that the equation always has solution provided that the product of eigenvalues is not an eigenvalue (i.e., the same as in the case discussed). \demo Exercise. Discuss the case of complex eigenvalues and show that even in that case, the conclusions are not altered. \demo Exercise. Show that in the case that we have one resonant term, we can still eliminate all others provided that they are not resonant. All this leads to the Poincar\'e-Dulac theorem. \proclaim Theorem. Let $f$ be a formal power expansion. Then, there exists a power expansion $h$ such that $$h\circ f\circ h^{-1} = Df(0) + R$$ where $R$ is a poser series expansion that only contains resonant terms (i.e., terms $R_{j_1\ldots j_n}^i$ where $\lambda^i = \lambda^{j_1}\cdots \lambda^{j_n}$). The fact that these formal power series have a real --- not purely formal meaning is given by the celebrated Sternberg theorem. \proclaim Theorem. Assume that $f$ as before is non-resonant, and hyperbolic. Then, one can find a $C^\infty$ change of variables reducing it to the linear part. There are many variants of the Sternberg theorem which cover cases in which there are resonances or finite number of derivatives. We will not discuss them here. S nadard references are \cite{27,28} We note that there are many physically natural reasons that force the existence of resonances. For example, preservation of volume requires that the product of eigenvalues is 1 and the Hamiltonian structure of equations imposes much more severe restrictions. \demo Exercise A linea map $A$ is called Hamiltonian (or symplectic) if $A^{-1} = J A J$ where $J$ is the antisymmetric matrix $J = \left( \pmatrix{ & 0 & 1 \cr &-1 &0} \right)$. Show that if $\lambda$ is an eigenvalue of a real Hamiltonian matrix, so are $ \bar \lambda$ and $1/\lambda$. We also point out that even if the existence of resonances is atypical for a map, it may be unavoidable for families of maps. That is, if a map is resonant, we can make an small perturbation that makes it non-resonant. If we have a family one of whose eigenvalues starts in $-0.8$ and ends up in $-1.2$, all nearby families are going to contain a point that has an eigenvalue $-1$. Notice also that these normal forms are also discontinuous. An small modification in the eigenvalues may make a particular term become non-resonant. (This is quite similar to the Jordan blocks.) Since the structure of the terms depends on exact comparisons calculations in which roundoff is introduced cannot be used to compute normal forms. One has to use exact arithmetic. If one is using inexact arithmetic, the only possible course is not to eliminate the terms for which the product of eigenvalues is close to another eigenvalue. We refer to \cite{24} for programs for the exact calculation of normal forms using symbolic manipulators. We warn that the inductive method we have sketched before is not very efficient and that one has to use a more efficent method. The thory of normal forms for Hamiltonian maps is quite interesting. The last exercise shows that there are neccesarily lots of resonances. Also, it shows that for Hamiltonian maps, there can be no linearly stable equilibrium points. The best that one can hope is to have linearly neutral equilibrium points. For these points, one has to study the stability using higher order terms so that the problem has a considerable interest. It is conjectured that most fixed points for maps are eventualy unstable when the number of variables is four or greater -- this is related to the famous Arnol'd diffusion -- but neverthelss, one can study the normal forms in great detail and obtain estimates for practical times of stability. Since the fixed points are the places where one tends to leave satelites, this is of considerable practical interest \cite{29}. \section{5. Invariant manifolds} By Hartman-Grobman theorem, in a neighborhood of a hyperbolic point, we expect that there is a continuous image of a subspace that converges to the origin, under iteration and another one of points $\Delta$ that converge to the origin under inverse iteration. It is surprising, but true that these sets are manifolds as smooth as the map. The precise statement is: \proclaim Theorem. Let $f:\real^n \to \real^n$. Assume $f(0) =0$, $f\in C^r$, $r\in\real \cup \{\omega,\infty\}$, $$\eqalign{&\real^n = E^s\oplus E^u\ , \quad Df(0) E^s = E^s\ Df (0)E^u = E^u \cr &\|Df(0)^N |_{E^s}\| \le \lambda <1\ , \|Df(0)^{-N} |_{E^u}\| \le \lambda <1\ , \cr}$$ Then, in a neighborhood of zero we can find a $C^r$ function $W: B_\delta \subset E^s\to E^u$ whose graph is invariant under $f$ $f(0)=0$, $Df(0)=0$ (analogously for $E^u$). Moreover, if $x$ is such that $|f^n(x)|\le\delta$ then $x\in \hbox{\rm graph}(W)$). In \cite{30} one can find a very readable proof reproduced in \cite{31}. We will denote these points as $W^\delta$, the ``local stable manifold''. Using this local stable manifold, one can characterize the set of all points that converge to the origin as the union of the inverse images of the local stable manifold. These stable and unstable manifolds tend to wander around large regions of the phase space of the system so that they organize a lot of the behavior. Many calculations of the unstable manifold --- the stable manifold is the unstable manifold for the inverse --- are based on the following result, usually called the $\lambda$-lemma. (See \cite{6} for a proof of an slightly weaker result. \proclaim Theorem. Let $D$ be a $C^r$ disk of the same dimension of $E^u$, sufficiently small, and transversal to $E^s$. Then, $f^n (D) \cap B_\delta$ converges to $W^{u,\delta}$ in the $C^r$ sense. The reason for this convergence is that the direction along $E^u$ stretching and the direction along $E^s$ is contracting. The non-linear effects are negligible in sufficiently small neighborhoods. A practical algorithm in the case that the manifolds are one-dimensional is to pick an small segment, fill it with points, iterate the map, discard the points that go out of the neighborhood and interpolate with the remaining ones --- e.g., using splines --- to make up for those lost. The algorithm to compute the local manifold stops when there is no change in these points. Once this local manifold has been computed, we can interpolate more densely --- refine if desired by iterating again --- and iterate and plot without discarding the points that move out of the neighborhood. This seems to be an extraordinarily effective and reliable algorithm. It, of course fails in situations in which the stable direction is very close to the unstable one (this happens often in Hamiltonian systems near to integrable where the angles are often exponentially small in the perturbation parameter). The only drawback is that one has to select the length of the segment, the number of points. It seems difficult to devise rules that select these parameters without human intervention. Even if the mathematical justification of the $\lambda$-lemma works in any number of dimensions, this algorithm is hard to adapt to two or bigger dimensions. The main source of trouble is that if the two unstable directions have very different eigenvalues, a large number of iterations tends to cluster the points along the one with largest modulus. If one knew how many iterations were necessary, one could counteract this problem by iterating an ellipse which is much larger in the weak unstable direction than in the strong unstable one. Another complication is that interpolation of two dimensional functions is notoriously more difficult than for one-dimensional functions. Of course, the problem becomes much harder as the dimension increases. Even if it is quite possible to compute two dimensional invariant manifolds for some concrete systems, it seems that there is no generally available implementation that is robust enough to be left in the hands of general users. It would most welcome. Another completely different algorithm for 1-dimensional manifolds is trying to find out the power series expansion for the function $\varphi : \real\to \real^n$ that satisfies $$\eqalign{&\varphi (\lambda t) = f\varphi (t)\cr &\varphi (t) = \sum \varphi_n t^n\cr}$$ We see that $\varphi_0=0$, $Df(0)\varphi_1 = \lambda\varphi_1$ so that $\varphi_1$ is an eigenvalue and that to order $n$, the expansion becomes $\lambda^n \varphi_n = Df(0) \varphi_n +R_n$ where $R_n$ is an expression involving the previously computed $\varphi$. Since $\lambda^n$ is not an eigenvalue this equation can be solved for all $\varphi_n$. It turns out that one can estimate the recursion in the series under the assumption that $f$ is analytic. This was done in \cite{32} which is, presumably the first proof of the stable manifold theorem. The algorithm is easy to implement in the case the $f$ is a polynomials. A very nice modern discussion of the implementation and discussion of computational errors can be found in \cite{33}. I have found it more computationally efficient and more stable to use an accelerated convergence method that I will describe now. If we have $\varphi (\lambda t) = f(\varphi (t)) +R$ where $R=O(t^N)$ we try to find $\Delta (t) = \sum_N^{2N}\Delta_N t^n$ in such a way that $(\varphi +\Delta) (\lambda t) = f((\varphi +\Delta)t) +\tilde R$ with $\tilde R= O(t^{2N})$. The equations for $\Delta$ are $$\Delta (\lambda t) = Df \bigl( \varphi (t)\bigr) \Delta (t) - R^{[\le N]}$$ where $R^{[\le 2N]}$ is $\sum_N^{2N} R_n t^n$. If $Df(\varphi (t)) = \sum_{n=0}^N S_n t^n + O(t^{N+1})$ the equation for $\Delta$ amounts to $$\Delta_n \lambda^n = Df (0) \Delta_n + \sum_{k=1}^n S_n \Delta_{n-k}+R_n$$ We see that this can be solved much faster than the non-linear recursions that involve recomputation of the remainders. More importantly, it is much more stable numerically. Similar tricks can be used to accelerate the calculation of normal forms. Notice that the only thing that is used in this algorithm is that $\lambda$ is an eigenvalue but that $\lambda^n$ is not. In those cases, one can construct invariant manifolds which are tangent to the eigenspaces just the same. The algorithm can also be extended easily to the case that $\lambda$ is a diagonal matrix. With some more complication it can also be extended to other cases. \section{6. Bifurcations} Many physical systems involve external parameters and it is a natural question to try to study the behavior for all possible values of the parameter. The implicit function theorem guarantees that if the derivative of the return map at a periodic orbit does not have an eigenvalue 1, the periodic orbit will also exist for sufficiently close values of the parameter and, moreover can be written as an smooth function of these parameters. Hartman-Grobman theorem says that, unless some eigenvalues cross the unit circle the topological class of the periodic point will not change. >From the numerical point of view, this following of non-degenerate periodic orbits is relatively easy to implement. If we know periodic orbits for a value of the parameter, we can take it as initial guess for a Newton method that tries to find the periodic orbit for an slightly bigger value (how to choose it?). Since the solution $x(s)$ of $f(x,s)=x$ satisfies $${d\over ds}x(x)= - \left( {\partial f\over\partial x} (x,s)-\Id \right)^{-1} {\partial f\over\partial s} (x,s)$$ we can try to use a high quality O.D.E.\ solver (with adaptive step!) to obtain very good approximate solutions that can then be refined by an equation solver. These methods --- usually called ``homotopy'' or ``continuation'' are very effective and frequently one finds solutions by connecting the physically relevant problem to another problem that can be solved exactly. Of course, finding suitable systems that approximate the required one is an art. Notice also that the considerations above apply to fixed point equations that involve parameters. Since we reduced many equations for invariant objects to fixed point equations, these methods are widely applicable. We will, nevertheless, discuss only periodic orbits. Severe complications arise when some eigenvalue of the periodic orbit cross the unit circle and, hence the local picture changes quite drastically --- even the topological picture changes. There are three main simple ways in which this can happen: \item{A)} one eigenvalue becomes 1 \item{B)} one eigenvalue becomes $-1$ \item{C)} a pair of complex eigenvalues crosses the unit circle. Even if there are other possibilities such as several eigenvalues coming onto the unit circle at the same value of the parameter, we will ignore them for the moment for the sake of simplicity. Notice also that, for a ``typical'' family, one does not expect that these things happen since arbitrarily close families can break the degeneracy. Of course this is a rather poor justification for restricting our attention to these cases. In Physics, one often has systems with symmetries or geometric properties that force these situations or we may be particularly interesting in studying the system in the situation where this happens --- analogous to putting a thermodynamic system at the tricritical point! The true argument for restricting ourselves to these cases is just lack of space. In these lectures we can only underline the main results. We refer to standard treatises such as \cite{34} for full proofs and precise formulations --- including the explicit form of the ``non-degeneracy'' conditions. A good treatment of the elementary bifurcations suitable for practitiones is \cite{35}. The behavior in A) can be understood by looking at what happens for the family $f_\lambda (x) = X^2 +\lambda$ at $\lambda = 1/4$. We see that the fixed points are $X= {1\over2} \pm \sqrt{(1/4-\lambda) (1/4+\lambda)}$. We see that at $\lambda =1/4$ we obtain a double root and that, for $\lambda >1/4$ there is no root an for $\lambda$ smaller there are two roots at a distance commensurate with $(\lambda -1/4)^{1/2}$. Notice that the fixed point start to change very quickly with the parameter $({d\over d\lambda} x_\lambda\approx {1\over\sqrt{\lambda -1/4}})$ but with a well defined law. It is possible --- but non-trivial --- to show that qualitatively similar things happen for ``typical systems'' in higher dimensions. The proof roughly goes as follows: First, we show that, by a suitable change of variables --- that depends on the parameter $\lambda$ --- we can reduce our system to: $$f\pmatrix{ s\cr u\cr x\cr} = \pmatrix{A_\lambda (s)\cr B_\lambda (u)\cr C_\lambda + (1+\lambda) x+ D_\lambda x^2\cr} + o (|x| + |u| + |s|)^3$$ where $s,u$ correspond to the stable and unstable directions of $Df(0)$ and $A_\lambda (s)$, $B_\lambda (u)$ are contractive and expansive mappings respectively. This change of variables can be computed exactly as in the normal form case, but we just leave the terms that we cannot eliminate. Second, we study the dynamics of the main part of the normal form above. Third, we show that the features of existence of periodic orbits etc. are not modified by the high order terms. If we believe that the perturbation argument can be carried out, we see that the conclusions we wanted to argue are true provided that $C_0\ne0$, $D_\lambda \ne0$. This is what we meant by a ``typical system''. The fact that the stable or unstable components do not play any role in the bifurcation is quite a general phenomenon. It can be understood at the level of normal forms as we have sketched here (this is the so-called Lyapunov-Schmidt reduction) or one can invoke a high powered invariant manifold theorem to show that all the interesting behavior happens on a two dimensional set. (See \cite{31} for a detailed discussion of both methods.) \demo Exercise. In population dynamics it is very natural to assume that zero population in one generation implies zero population for future generations so that one is lead to consider models of the form $f_\lambda (x) = xg_\lambda (x)$. Show that these models violate the genericity assumptions and study the bifurcations that one gets. (An often used example is $f_\lambda (x) = \lambda x(1-x)$. For this example the critical value is $\lambda =1$.) It is very tempting to say that B) reduces to A) just by considering $f^2$. Unfortunately this does not work as can be seen by noticing, by the implicit function theorem we know that there is a family of fixed points for $f$ --- hence for $f^2$ --- on both sides of the critical point. Nevertheless A) predicts that there are no fixed points on one side. The reason why this does not work is that when $f_{\lambda_0}^1 (0)=-1$ $f_{\lambda_0}^2$ violates the non-degeneracy conditions. \demo Exercise. Show that if $f_{\lambda_0} (0)=0$, $f'_{\lambda_0}(0)=-1$ then $f_{\lambda_0}^2{}'' (0)=0$. Nevertheless, it is possible to show that, given some non-degeneracy conditions in $f_\lambda$ we can obtain normal forms up to one order more and the results produced by the main part are not affected by the high order remainder terms. \demo Exercise. Compute explicitly the non-degeneracy conditions. Suffice it to say that the main result is that the fixed point changes from stable to unstable but that also a periodic orbit of period 2 appears on one side of the bifurcation. This periodic orbit has an amplitude proportional to the square root of the difference of the parameter to bifurcation and it has opposite stability to that of the fixed point on the branch. Except for that, all the possibilities can occur. This can be seen very simply by studying the family $$f_\lambda (x) = (-1\pm \lambda) x + Ax^2 + Bx^3$$ for different values of $A,B$. \demo Exercise. Study --- algebraically --- or with the help of an small computer the periodic orbits of the example above and get a classification of all the pictures. The most spectacular bifurcation happens in case C). This is the so-called Hopf bifurcation. It is important to note that the non-degeneracy conditions include that the eigenvalues are not roots of unity up to order 4. In that case, one can show that taking polar coordinates, one can reduce the normal form $$f_\lambda (r,\theta) = \bigl( \varphi_\lambda (r),\theta+\omega (r)\bigr) + o(r^5)$$ where $\varphi_\lambda (r) = (1+\lambda) r+A_\lambda r^2 + B_\lambda r^3$ and $\omega$ has a Taylor expansion that starts with the argument of the eigenvalues that cross the unit circle. Again, under genericity assumptions we can see that the equation for $r$ has a fixed point that grows with the parameter as $\sqrt r$. Notice that this corresponds to an invariant circle. Again, it is possible to show that this is not modified by the high order terms in a substantial way. Notice that this is quite a remarkable result since it produces an invariant circle for a map, quite a non-trivial statement. The theory of bifurcation of differential equations can be reduced to the theory we have just developed by taking time-one maps. (Notice that then, we do not get an analogue of case B).) The proofs also can be done directly. In particular, the treatment of case C) is much simpler since a periodic orbit of a differential equation is a much less surprising statement than an invariant circle for a map. Even if only elementary bifurcations are involved, the complete bifurcation diagrams may be quite complicated, periodic orbits may appear and disappear through any combination of those bifurcations. Actually, some bifurcations tend to appear together. For example, after a differential equation has experienced a Hopf bifurcation, the return map of the resulting periodic orbit may experience a Hopf bifurcation giving rise to a solution with two frequencies (this is part of the so-called Ruelle-Takens proposal for the origin of turbulence). It is also known that period doubling bifurcations tend to appear together and to lead to a cascade with remarkable properties. (See the collection \cite{35,7} for reviews on these properties.) The idea that there are several cascades of bifurcations that can lead from ordered states to complicated motions has been discussed in \cite{7,36} where these phenomena have been given the name of ``scenarios for the transition to chaos''. One can understand the existence of scenarios with the help of the following picture. We may think of the set of bifurcating maps as a submanifold --- with some corners corresponding to degenerate cases --- of the set of maps. In some regions in the set of maps, these manifolds accumulate towards the boundary of a region --- somewhat ill defined --- of chaotic systems. Of course a good numerical package for studying solutions as functions of the parameters should include some provisions to try to defeat the complications introduced by the bifurcations and the encounter with critical points. This is somewat delicate because in the neighborhood of the bifircations, the new periodic orbits move very fast with respect to the parameter. So, one has to take special precautions to extrapolate for new values in a way which is consistent with the laws of the bifurcation and not using a linear extratpolation method. A popular method is the arc-length parameterization method \cite{37}. In this method we introduce a fictitious parameter (call it $s$) and we try to write both the real parameter and the solution as a function of $s$ while at the same time we impose some normalization condition. If we were interested in solving $f(x,\lambda)=x$ we would introduce the system $$\eqalign{ & f(x,\lambda) = x\cr &{\partial f\over\partial x} (x,\lambda) \dot x + {\partial f\over\partial \lambda} \dot\lambda = 1\cr}$$ (where $\dot{}$ denotes the derivative with respect to $s$). If we had a solution $x^*\lambda^*$ we would declare it a good solution for the parameter $s=0$ and we could try to find solutions $x,\lambda$ for small $s$ by solving: $$\eqalign{ &f(x,\lambda) -x=0\cr &{\partial f\over\partial x} (x^*,\lambda^*) (x-x^*) + {\partial f\over\partial\lambda} (x^*,\lambda^*) (\lambda-\lambda^*)=s\cr}$$ If we choose $s$ sufficiently small, this will produce a new pair of solutions. There are several publicaly available implementations of these continuation methods. One of them is PITCON whose theory is discussed in detail in \cite{38}. A really useful program is AUTO which not only follows the solutions but also tries to follow the different branches when one that appear as bifurcations. These two programs are relatively easy to use. One just needs to supply the function on is trying to solve and its derivatives as well as some starting guesses and some global parameters. Even if they fail sometimes, they usually provide enough information that one can try to design better ways. The study of bifurcations typically requires to do quite extensive manipulations. Many of them can be done using symbolic manipulators such as Reduce, Macsyma (we have often used the version of DOE MAXIMA, which is quite affordable and comes the the source code), Maple or Mathematica. These manipulators are often invaluable in substituting expansions into each other. Unfortunately, they are quite bade at trying to discuss the existence of real solutions as functions of the parameter. Even in the case that the equations can be solved just using quadratic equations, the existence or not of solutions involves that certain expressions are positive or not. Of course, since the expressions involve common symbols, the assumptions are not independent. It does not seem that there is any symbolic manipulator capable of dealing with such things. MAPLE, MAXIMA and MACSYMA have ways of assigning attributes to the objects --- such as positiveness --- that can be used to control some simplifications. For example, in many calculations it is quite important to disable simplifications such as sqrt~$[x^2]=x$ --- which are wrong if $x$ is negative! --- or $(x-1)A= (x^2-1)\Rightarrow A= (x-1)$ (it misses $x=1$). In MAXIMA this can be done by controlling the application of ratsimp. A very important algorithm to know about when one has to consider systems of algebraic equations with several unknowns is the Grobner-basis algorithm. This algorithm is, roughly, an analogue of gaussian elimination for systems of non-linear polynomial equations. One can operator systematically on the equations in such a way that they appear in a upper-triangular way. That is, we obtain a system of equations of the form $$F_i (x_1,\ldots, x_i) =0$$ The discussion is very similar to gaussian elimination. It could well happen that $F_1\equiv 7$, hence, there is no solution or that $F_1\equiv 0$, and in that case we have more than one solution. In the ideal case, we can get to solve $F_1(x_1)=0$ and obtain a finite set of values, which when substituted in $F_2 (x_1x_2)$ will produce a finite set of values for $x_2$, etc. The symbolic manipulators mentioned above contain implementations of the Grobner basis algorithm. Nevertheless, their speeds are quite different. For a problem that MAXIMA did in 10 minutes, Mathematica took more than 20 hours and it did not finish! In any case, this algorithm seems quite impractical for more than say 6 equations and 6 variables. For a fuller discussion of the Grobner basis algorithm we refer to \cite{39}. One should note that in this area, even small modifications of the algorithm may make eenourmous differneces in speed for concrete problems. If one knows which normal forms one wants to compute there are several FORTRAN or C programs available. A nice one is discussed in \cite{40}. Unfortunately, the inexact arithmetic makes it impossible to detect degenerate normal forms since this identification depends deciding whether exact equalities take place. \section{6. References.} \item{1.} A. Dou, {\it Ecuaciones Diferenciales Ordinarias} (Dossat, Madrid, 1968). \item{2.} S. Wiggins, {\it Global Bifurcations and chaos} (Springer, New York, 1988). \item{3.} F. C. Moon, {\it Chaotic and Fractal dynamics} (Adison Wesley, New York, 1992) \item{4.} C. Sparrow, {\it The Lorenz equations} (Springer, New York, 1982). \item{5.} V. I. Arnol'd, A. Avez, {\it Ergodic Problem of Classical Mechanics} (Benjamin, New York, 1968). \item{6.} J. Palis, W. de Melo, {\it Geometric Theory of Dynamical Systems} (Springer, New York, 1982). \item{7.} J.--P. Eckmann, D. Ruelle, Ergodic theory of chaos and Strange attractors. {\it Rev. Mod. Phys.}, {\bf 57}, 617--656 (1985). \item{8.} D. Ruelle, {\it Chaotic evolution and strange attractors} (Cambridge Univ. Press, Cambridge, 1989). \item{9.} P. Hartman, {\it Ordinary Differential Equations} (Birkhauser, Boston, 1983). \item{10.} E. Hairer, G. Wanner, E. N\O rset, {\it Solving Ordinary Differential Equations I} (Springer Verlag, New York, 1987). \item{11.} E. Hairer, G. Wanner, {\it Solving Ordinary Differential Equations II} (Springer Verlag, New York, 1991). \item{12.} D. Ruelle, Ergodic theory of differentiable dynamical systems, {\it Pub. Mat. I.H.E.S. }, {\bf 50}, 27--58, (1979). \item{13.} L.S. Young, Dimension, entropy and Lyapunov exponents in differentaible dynamical systems, {\it Physica}, {\bf 124A}, 639-645 (1983). \item{14.}R. A. Johnson, K. Palmer, G. Sell, Ergodic properties of linear dynamical systems. {\it SIAM Jour. Math. Anal.} {\bf 18 } 1--33, (1987). \item{15.} H. Poincar\'e, {\it Les m\'ethodes nouvelles de la m\'ecanique c\'eleste}, (Gauthier Villars, Paris 1891-1897). \item{16.} {R. Ma\~n\'e, An ergodic closing lemma, {\it Ann. of Math.},{\bf 116},{503-541}(1982). \item{17.} A. Katok, Lyapounov exponents, entropy and periodic orbits for diffeomorphisms, {\it Pub. Mat. I.H.E.S.},{\bf 51}, {137--173} (1980). \item{18.} C. Falcolini, R. de la Llave, {\it A rigorous partial justification of Greene's criterion} {\it Jour. Stat. Phys.},{\bf 67},{609--643}, (1992). \item{19.} R. DeVogelere, On the structure of symmetric periodic solutions of conservative systems, with applications, { \it Contributions to the theory of Nonlinear oscillations IV} {S. Lefschetz ed.} (Princeton Univ. Press, Princeton, 1958). \item{20.} H-T. Kook, J.D. Meiss, Periodic orbits for reversible symplectic mappings, {\it Physica}{\bf 35D}, 65--86 (1988). \item{21.} J. Stoer, R. Burlisch, {\it Introduction to numerical analysis} (Springer Verlag, New York, 1980). \item{22.} R. G. Helleman, Variational Solutions of Non-integrable Problems, {\it Topics in Noonlinear Dynamics} (S. Jorna, ed) (A.~I.~P. New York, 1978). \item{23.} J. D. Meiss, Symplectic maps, variational principles and transport. {\it Rev. Mod. Phys.}, {\bf 64}, 795--847, (1992). \item{24.} R. H. Rand, D. Armbruster, {\it Perturbation methods, bifurcation theory and computer algebra}, (Springer--Verlag, New York, 1987). \item{25.} Z. Nitecki, {\it Differentiable Dynamics}, (M.I.T. Press, Cambridge, 1971). \item{26.} P. Hartman, On local homeomorphisms of Euclidean spaces, {\it Bol. Soc. Mat. Mex.}, {\bf 5}, 220--241, (1960). \item{27.} E. Nelson, {\it Topics in Dynamics: I Flows}, (Princeton Univ. Press, Princeton, 1968). \item{28.} S. Sternberg, {\it Celestial Mechanics}, (Benjamin, New York, 1969). \item{29.}A. Giorgilli, A. Delshams, E. Fontich, L. Galgani, C. Sim\'o, Effective stability for a hamiltonian system near an elliptic equilibrium point with an application to the restricted three body problem {\it Jour. Diff. Eq.},{\bf 77}, 167--198, (1989). \item{30.} O. E. Lanford III, Bifurcation of periodic solutions into invariant tori: the work of Ruelle and Takens, {\it Lect. Notes in Math.},{\bf 322}, (1973) \item{31.} J. Marsden, M. McCracken, {\it The Hopf bifurcation and certain of its applications.} ( Springer, New York, 1976). \item{32.} H. Poincar\'e, Sur une classe nouvelle de transcendantes uniformes, {\it Jour. de Math.},{\bf 6},313--365, (1890). \item{33.} W. Fraceschini, L. Russo, Stable and unstable manifolds of the Henon mapping, {\it Jour. Stat. Phys.} {\bf 25}, (1981). \item{34.} S.--N. Chow, J. Hale, {\it Methods of bifurcation theory}, (Springer Verlag, New York, 1982). \item{35.} P. Cvitanovic, {\it Universality in Chaos}, (Adam Hilger, Bristol, 1984). \item{36.} J.--P. Eckmann, Roads to turbulence in dissipative dyanamical systems, {\it Rev. Mod. Phys.}, {\bf 53}, {643--654}, (1981). \item{37.} H. Keller, {\it Numerical methods in bifurcation problems}, (Springer Verlag, New York, 1987). \item{38.} W. C. Rheinboldt, {\it Numerical analysis of parametrized nonlinear equations}, (Wiley, New York, 1986). \item{39.} A. Akritas {\it Elements of computer algebra with applications}, (Wiley, New York, 1989). \item{40.}B.D. Hassard, N.D. Kazarinoff, Y.--H. Wan, {\it Theory and applications of Hopf bifurcation} (Cambridge University Press, Cambridge, 1981). \end