%% In case the "geom" style is not available, one can modify this file as
%% indicated in the comments, and then print the figures separately.
%% Erase next two lines if not using "geom" style.
%% Replace by the line following them (erase the leading "%%").
\documentstyle[12pt,geom]{article}
\proofingfalse
%%\documentstyle[12pt]{article}

\newcommand{\st}{\cal W^s}
\newcommand{\unst}{\cal W^u}
\newcommand{\stab}{\cal W^s(p)}
\newcommand{\unstab}{\cal W^u(p)}
\newcommand{\lstab}{\cal W_{loc}^s(p)}
\newcommand{\lunstab}{\cal W_{loc}^u(p)}
\newtheorem[{}{\em}]{theorem}{Theorem}


\begin{document}
\title{Visualization of Stable and Unstable Manifolds of Diffeomorphisms}
\author{Eduardo Tabacman\\The Geometry Center\\University of Minnesota}
\date{November 1993}

\maketitle
\section{Introduction}

Recently there has been interest in the visualization of objects in three
dimensional space, using computer graphics. Some programs that allow the user
to do this in a simple way have appear. (An example of which is GEOMVIEW, the
program used for all the figures in this paper). They require the user to provide
the program with some points on the surface to represent, usually in the form 
of a mesh, or grid of points, and the program takes care of the rest. This is 
consistent with the view geometers take of surfaces. One thinks of a surface as 
the image of a {\em parametrization}, that is, a diffeomorphism $\gamma:
\Bbb R^2\to \Bbb R^d$, or even as the 
map itself. A grid on the surface is then the image of a grid (like $\Lambda=\{
(j,k) : 0\le j,k\le N\}$, for some $N\in \Bbb N$)  in $\Bbb R^2$ by 
this map. 

Knowing an explicit parametrization is then very convenient. However, in some cases one wants to visualize a surface that is defined implicitly, one for which
no explicit parametrization is known. Such is the case of the manifolds 
considered here. The purpose of this paper is to show that there is a parametrization 
for these stable and unstable manifolds, and that it can be approximated by 
means of the diffeomorphism, in a way that is appropiate for using with 
software as the one mentioned above. We pay particular attention to
manifolds of dimension larger than one, mainly because we are interested
in dealing with problems that do not exist for one dimensional
manifolds.

Several people have studied the problems that arise when computing these
kind of manifolds. Most of them look at one dimensional manifolds
(see~\cite{Y:1} and~\cite{S:1}, for example). The paper~\cite{G:1}
considers higher dimensional manifolds. Although they correspond to
flows, and not to maps, some of the remarks made there apply to the
cases we consider. In particular, the problems encountered in computing
the manifolds are similar.

\section{The setup}

When studying the iteration of an invertible map $f:\Bbb R^d\to\Bbb R^d$, one is often interested in the long time behaviour of orbits. The simplest of all orbits in this regard is that of a {\em fixed point}, that is a point
$p$ such that $f(p)=p$, which consists of just that point.

(Next in complexity is the {\em periodic orbit}, i.e. the orbit of a point $p$ such 
that, for some $n\in \Bbb N$, $f^n(p)=p$. This orbit consists of finitely many
points: it is $\{p,f(p),\dots,f^{n-1}(p)\}$.
But since a periodic point of period $n$ can be regarded as a fixed point of the map $f^n$, we will only consider fixed points.)

The behaviour of other orbits can be much more complex, but there are points whose orbits are related to fixed points in a simple way: under iteration by $f$ or $f^{-1}$,
they approach the fixed point.

More precisely, if $p$ is a fixed point of $f$, let
\begin{eqnarray*}
\stab&=&\{x\in\Bbb R^d : f^n(x)\to p \text{ as }n\to\infty\}\\
\unstab&=&\{x\in\Bbb R^d : f^{-n}(x)\to p \text{ as }n\to\infty\}
\end{eqnarray*}

A simple consequence of these definitions is the invariance of these two sets, 
$f(\stab)=\stab\text{, and }f(\unstab)=\unstab$.

$\unstab$ is called the {\em unstable manifold\/} of $p$, while $\stab$ is the 
{\em stable\/} manifold.
Altough these two sets are not always embedded manifolds  in $\Bbb R^d$ (see figure 1), there are local versions that are:  Given a neighborhood
 $\cal U$ of $p$, we can define 
\begin{eqnarray*}
\lstab&=&\{x\in\Bbb R^d : f^n(x)\in\cal U\text{ for all $n\in \Bbb N$, and }f^n(x)\to p \text{ as }n\to\infty\}\\
\lunstab&=&\{x\in\Bbb R^d : f^{-n}(x)\in\cal U\text{ for all $n\in \Bbb N$, and }f^{-n}(x)\to p \text{ as }n\to\infty\}\text{ ,}
\end{eqnarray*}
and these two are embedded in $\Bbb R^d$. ($\lstab$ and $\lunstab$ depend on the choice
of $\cal U$, but we will omit this dependence in the notation). 

Also, $\unstab=\bigcup_{i=1}^\infty f^i(\lunstab)$
, and $\stab=\bigcup_{i=1}^\infty f^{-i}(\lstab)$. 

%% Unstable manifold for Henon map.
%% Erase next line if not using the "geom" style. 
\fig{h2-.51-1600-3200}{One dimensional manifold for Hen\'on's map}{A piece of a one dimensional unstable manifold. Note that the manifold accumulates onto itself, so it is not embedded in $\Bbb R^2$.}
 

A simple example :

Let $\bold A$ be an invertible $d\times d$ matrix, all of whose eigenvalues are real and distinct, and let 
\begin{eqnarray*}
f:\Bbb R^d&\to&\Bbb R^d\\ 
x&\to& f(x)=\bold A\,x
\end{eqnarray*}
Then $0\in\Bbb R^d$ is a fixed point of $f$.

What are $\unst$ and $\st$ ?

Since all the eigenvalues are distinct, there is a basis of eigenvectors for
$\Bbb R^d$, $\{v_1,\dots,v_d\}$, and we can write any $x\in\Bbb R^d$ as 
$x=c_1v_1+\cdots+c_dv_d$, for some $c_i\in\Bbb R$.

Then $$f^n(x)=\bold A^n\,x=\sum_{i=1}^dc_i\bold A^n\,v_i=\sum_{i=1}^dc_i\lambda_i^nv_i.$$
Since the $v_i$ are linearly independent, to have $f^n(x)\rightarrow 0$ we must
have\ \linebreak[4]$c_i\lambda_i^n\rightarrow 0\quad\forall i=1,\dots,d$, as $n\rightarrow\infty$.
Therefore, either $|\lambda_i|<1 \text{ or } c_i=0$.

In other words, 
$$\st(0)=\text{\sl span}\{v_i: \text{\,corresponding $\lambda_i$ have }|\lambda_i|<1\}.$$
Similarly, 
$$\unst(0)=\text{\sl span}\{v_i: \text{\,corresponding $\lambda_i$ have }|\lambda_i|>1\}.$$
We see then that in this case (as is true for any linear map), $\st$ and $\unst$
are planes, whose dimensions are given by the number of eigenvalues that have absolute value less than $1$ for $\st$ and greater than $1$ for $\unst$.

Some of this structure carries on to the nonlinear case:
\begin{theorem}
Suppose that $f:\Bbb R^d\rightarrow\Bbb R^d$ is a diffeomorphism, with fixed point $p$.

If $Df(p)$ (the Jacobian matrix of $f$ at $p$) has stable and unstable spaces 
$E^u$ and $E^u$ (as in the previous example), then $f$ has $\lstab$ and $\lunstab$ that are manifolds as smooth as $f$, of the same dimensions as $E^s$ and $E^u$, respectively.

Moroever, the affine plane $p+E^s$ at $p$ is tangent to $\lstab$, and  the same is true for $p+E^u$ and $\lunstab$.
\end{theorem}

For a proof, see \cite{H:ODE}. 

So, near $p$, and to a first order, $\lstab$ and $\lunstab$ are given by the corresponding manifolds of the linearization of the map $f$.

Before continuing, note that $\stab$ is the unstable manifold for $f^{-1}$,
so that we can restrict our treatment to one of the two, say $\unstab$.
Also, we will assume that the manifolds have dimension 2, and that the fixed point $p$ is {\em hyperbolic}, that is, no eigenvalue of $Df(p)$ has
modulus 1. Moreover, we will assume that all the eigenvalues are real and 
distinct.
Putting all this together, we assume that $Df(p)$ has two eigenvalues larger
than 1: $1<\lambda_1<\lambda_2$, with corresponding eigenvectors $\{v_1,v_2\}$
, which form a basis for  $E^u$.
None of this is essential, but it makes the exposition more clear.

\section{A First Try, and some Problems}
The fact that the affine plane $p+E^u$ (which we can compute using linear algebra) is an approximation to $\lunstab$, together with the relation 
between the local versions of the manifolds to the global manifolds, suggest a straightforward way of computing them:  
Take a small piece of $E^u$ near $p$ (which can be found in terms of the 
eigenvectors of $Df(p)$), so small that it can be regarded as a good approximation to $\lunstab$. Then we can get $\unstab$ by iterating this little
piece under $f$. Similarly, we can get $\stab$ by using $E^s$ and $f^{-1}$.
To use this in a way that is meaningful for a computer, we can put a grid on
$E^u$, and iterate the points on it.

An algorithm for this procedure may look like this :
\begin{itemize}
\bf
\item Choose a small $\epsilon>0$, and $N,n\in\Bbb N$.
\item Compute  $f^n(\epsilon (j v_1+k v_2))$, for $j=1,\dots,N, k=1,\dots,N$.
\end{itemize}
One problem that appears immediately, even when $\unstab$ is one dimensional, is that points tend to separate from each other under iteration, exponentially with $n$ (the number of iterations). The rate of this separation will be appproximately $\lambda$, the eigenvalue corresponding to $E^u$.

If we want to graph a large portion of $\unstab$, we will have to start with a 
grid with many points, so that the initial distance between them is small, and so that the final one ($\approx\lambda^n$ of the original distance) is acceptable.

But if $\unstab$ has dimension higher than one (two, as we are considering here) then there may be other problems. The most common is that, when the eigenvalues are different, some directions will be stretched more than others.

$E^u=\text{\sl span}\{v_1,v_2\}$, and if we start with a nice little square near $p$, after several iterations we will end up with something more resembling of a line in the $v_2$ direction than of any two dimensional being. (See figure 2).

%% Figures of rectangle & badly stretched one
%% Erase next line if not using the "geom" style.
\doublefig{linear5x5}{Effect of a linear map}{A nice grid, {\em before\/} and {\em after\/} four iterations of a linear map. (The scales of the
two pictures are different: see Notes on Figures (Section 5)).}

Four iterations is a rather modest number. For this map, we would need more than 
$10$ iterations to get a decent portion of the unstable manifold. After so many iterations, the final grid would be undistinguishable from a line. (See figure 3).
 
%% Figures of a simple (and bad) approach
%% Erase next line if not using the "geom" style.
\doublefig{2usp5-3-10x10-12}{An unsuccessful try}{Starting with a small uniform 
grid (left), we iterated it $12$ times, to get the almost curve on the right.
(Again, the scales are different in both pictures. See Section 5 for more details).}
 
It is in trying to avoid this situation that we look for a nice parametrization
of $\unstab$, that is a function $\gamma:\Bbb R^2\rightarrow\Bbb R^d$ whose 
image gives $\unstab$.

\section{A different way}

The construction of the parametrization of the manifold is based in the following theorem:
\begin{theorem}[Hartman-Grobman]
Suppose $f:\Bbb R^d\rightarrow\Bbb R^d$ is a diffeomorphism with a hyperbolic fixed point $p$.
Then there is a neighborhood $U$ of $p$ in $\Bbb R^d$ and a homeomorphism (that 
is, a continous map with continous inverse) $h:U\to\Bbb R^d$ with $h(p)=0$ conjugating $f$ to $Df(p)$ : 
\begin{equation}
f=h\circ Df(p)\circ h^{-1}. 
\label{conj}
\end{equation}
\end{theorem}

Again, we refer to \cite{H:ODE} for a proof.

Since $p$ is hyperbolic, 
we can split $\Bbb R^d=E^u\times E^s$, where $E^u$ is the unstable subspace  
for $Df(p)$, and $E^s$ is its stable subspace.
The map $h$ has to be thought as $h:U\to h(U)\subset E^u\times E^s$.
If we now define $\gamma: h(U)\cap E^u\to U$ by $\gamma(v)=h^{-1}(v)$, we have
that $$\gamma^{-1}\circ f\circ\gamma=Df(p).$$ From this, one gets  that 
$$\gamma=f^n\circ\gamma\circ Df^{-n}(p)\qquad\forall n\in\Bbb Z.$$

We can use this to extend the domain of definition of $\gamma$ to all of $E^u$, as follows :

If $v\in E^u$, $Df^{-n}(p)(v)\to 0$ as $n\to\infty$. Therefore, there is an 
$n\in \Bbb N$ such that $Df^{-n}(p)(v)\in h(U)\int E^u$. Since this is in the 
domain of $\gamma$, we can define 
$$\gamma(v)=f^n\circ\gamma\circ Df^{-n}(p)(v).$$

It is easy to see that this is a well defined homemorphism, and that $\gamma(E^u)=\stab$. In other words, this $\gamma$ is the parametrization we want. This
last formula gives the key to computing $\gamma$:

Let $i:E^u\to\Bbb R^d$ be $i(v)=p+v$ (that is, translation by $p$). The image of
this map, the plane $p+E^u$, is tangent to $\unstab$, so that near $p$ the two 
sets are close. Therefore, we can regard $\gamma\approx i$ near 0. Then we can 
compute $\gamma(v)$ by taking $n$ so large that $Df^{-n}(p)(v)$ is close to 0,
and using the approximation $$\gamma(v)\approx f^n\circ i \circ Df^{-n}(p)(v).$$
To be more concrete, an algorithm of this procedure for generating a grid on 
$\unstab$ may look something like this :
\begin{itemize}
\bf
\item Choose a small $\epsilon>0$, and $N\in\Bbb N$.
\item Choose $n$ so that $Df^{-n}(p)(N v_1+N v_2)\in B_{\epsilon}(p)$.
\item Compute $f^n\circ i\circ Df^{-n}(p)(j v_1+ k v_2)$, for $j=1,\dots
,N, k=1,\dots,N$.
\end{itemize}
(See figure 4 for an example).
Here $\{v_1,v_2\}$ is a basis for $E^u$, not necesarily of eigenvectors.
The size of $\epsilon$ be should be determined by experimenting. One can try
a value and then repeat the procedure with a smaller $\epsilon$, until the 
results do not change noticeably.


%% Nice pictures of the good way of doing things
%% Erase next line if not using the "geom" style.
\triplefig{2sp5-3-9x9}{The three stages}{From the original grid, get a new one by iteration of $Df(p)$, and then the final result, after iteration under $f$.
For details (as what is $f$, or the different scales of the pictures), see 
Section 5)}

This is really not so different from the idea tried in the previous section.\newline
$Df^{-n}(p)$ is a linear map, so that if $\Lambda$ is a rectilinear grid, 
$Df^{-n}(p)\,.\,\Lambda$ is another rectilinear grid, only much smaller. Therefore, we are iterating a small rectilinear grid in $E^u$ near $p$, several times under $f$.
The difference is that while the grid $\Lambda$ is a nice uniform one, after 
iterating under $Df^{-1}(p)$ we will have distorted it, stretching more in some directions than in others, pretty much in the opposite way as to what $f$ does. This way, when we iterate this small grid under $f$, we approximately undo the distortions, and end up with a grid that looks more uniform, like the original one, $\Lambda$.

For example, we use a grid $\Lambda$ with sides in the direction of the
eigenvectors in $E^u$ (so that $x\in\Lambda\ \Longrightarrow\  x=c_1 v_1+c_2 v_2,\quad $for some $c_1,c_2\in\Bbb Z$), we will have $Df^{-n}(p)\,.\,x=c_1\lambda_1^{-n}v_1+c_2\lambda_2^{-n}v_2$. So, the net effect is to iterate under $f$ a small grid, with sides in the same direction as the original $\Lambda$, but of different lengths, according to the eigenvalues in each direction.

We can use the eigenvectors (or generalized eigenvectors) as coordinate system in $E^u$ (and then $Df(p)$ will be in Jordan form), but this is not necessary. If we can see what $E^u$ is, we don't need to find the eigenvalues or eigenvectors.

As mentioned before, one can use the same procedure even if there are some repeated eigenvalues, not all of them are real, or if some have modulus 1. The only
thing we need is a basis for $E^u$ (which is the case described is simpler, 
being given by eigenvectors), but if we can see what $E^u$ is, we don't need to
find the eigenvalues or eigenvectors. (See the notes for figures 3 and 4).

If the eigenvectors are not as nice, the effect of $Df(p)$ will not be as easy
to describe, and the image of $\Lambda$ will not be another grid with the sides
in  the same direction as the original one, but the procedure remains the same.

\pagebreak
\section{Notes on the figures}

{\bf Figure 1 :}\ \ The (real) Hen\'on map is 
\begin{eqnarray*}
h:&\Bbb R^2&\to\Bbb R^2\\
&\pmatrix{x\cr y\cr}&\to\pmatrix{x^2-y+a\cr x\cr}.
\end{eqnarray*}

In this picture, $a=0.51\/$, and the fixed point is at 
${\scriptstyle\pmatrix{1.7\cr1.7\cr}}$ (denoted by the cross).
\vspace{1.5pc}

\noindent{\bf Figure 2 :}\ \ The linear map in question is 
\begin{eqnarray*}
L:&\Bbb R^2&\to\Bbb R^2\\
&\pmatrix{x\cr y\cr}&\to\bold A\,.\pmatrix{x\cr y\cr}\ ,
\end{eqnarray*}
where $\bold A=\pmatrix{\frac{2+k}2&\frac{2-k}2\cr\frac{2-k}2&\frac{2+k}2\cr}$.
In the picture, $k=3$.

 The eigenvalues are $k$ and $2$, so it is an expansion. (Actually, the matrix $\bold A$ is $\pmatrix{k&0\cr0&2\cr}$ after a rotation 
by $\frac{\pi}4$).

The two figures are not in the same scale. The second one is much larger than the first one (the map is an expansion, after all). 
\vspace{1.5pc}

\noindent{\bf Figures 3,4 :}\ \ They correspond to the map 
\begin{eqnarray*}
f:&\Bbb R^3&\to\Bbb R^3\\
&\pmatrix{x\cr y\cr z\cr}&\to\pmatrix{\frac{2+k}2x+\frac{2-k}2y+yz\cr
\frac{2-k}2x+\frac{2+k}2y+xz\cr\frac12z-xy+xz\cr},
\end{eqnarray*}
with $k=3$.

Its linear part at $p=\pmatrix{0\cr0\cr0\cr}$ is 
$$
Df(p)=\pmatrix{\frac{2+k}2&\frac{2-k}2&0\cr\frac{2-k}2&\frac{2+k}2&0\cr0&0&\frac12\cr}.
$$
Note that $Df(p)$ leaves the $xy-$plane invariant, and that its restriction to it is the same matrix $\bold A$ of figure 2, so that $E^u=xy-$plane.

In figure 3 we iterated a small (side $=10^{-5}$) uniform grid, $12$ times.

In figure 4 we used the procedure described in the text, asking that the intermediate grid be about the same size as the starting one of figure 3.
(For that we had to take $14$ iterations).

Again, the pictures are not to scale. The initial grid in figure 3 and the intermediate one in figure 4 should be much smaller. Just for reference, the 
axis in the last picture in both figures are $10$ units long, while the grid in the first picture of figure 3 and the second one in figure 4 are about $10^{-5}$ of a unit.
\begin{thebibliography}{9}
\bibitem{G:1} J.Guckenheimer and P.Worfolk, {\em Dynamical Systems: Some
computational problems}, in proceedings of the ``Bifurcations and
Periodic Orbits of vector fields'' conference, D.Schlomiuck (Ed.). NATO
ASI Series, Series C: Math. and Phys. Sciences, 1993, D.Reidel
Publishing Co.
\bibitem{H:ODE} P.Hartman, {\em Ordinary Differential Equations}. Wiley: New York. (1964).
\bibitem{S:1} C.Sim\'o, {\em Modern methods in Celestial Mechanics},
D.Benest and C.Froeschl\'e (Eds.), 285--330. Editions Fronti\'eres,
France, 1990.
\bibitem{Y:1} Z.You, E.J.Kostelich, and J.Yorke, {\em Calculating stable
and unstable manifolds}. Int. J. of Bifurcation and Chaos, Vol 1, No.3, 
1991, 605--623.
\end{thebibliography}
\end{document}