%%% 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