\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
2007 Conference on Variational and Topological Methods: Theory, Applications,
Numerical Simulations, and Open Problems.
{\em Electronic Journal of Differential Equations},
Conference 18 (2010),  pp. 45--56.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{45}
\title[\hfilneg EJDE-2010/Conf/18/\hfil
Invariant sets for dynamical systems]
{A geometric approach to invariant sets for dynamical systems}

\author[D. Medina, P. Padilla\hfil EJDE/Conf/18 \hfilneg]
{David Medina, Pablo Padilla}  % in alphabetical order

\address{David Medina \newline
Instituto Tecnol\'ogico Superior de Perote\\
Carretera Perote-M\'exico km. 2.5\\
Centro, Perote, Veracruz, C. P. 91270, M\'exico}
\email{medina@math.unam.mx}

\address{Pablo Padilla \newline
Departamento de Matem\'aticas y Mec\'anica, IIMAS-UNAM,
Apartado Postal 20-726, C. P. 01000 M\'exico,  M\'exico}
\email{pablo@mym.iimas.unam.mx}

\thanks{Published July 10, 2010.}
\subjclass[2000]{37L05}
\keywords{Invariant sets; dynamical systems; area functional;
\hfill\break\indent steepest descent method}

\begin{abstract}
 In this article, we present a geometric framework to study invariant
 sets of dynamical systems associated with differential equations.
 This framework is based on properties of invariant sets for
 an area functional. We obtain existence results
 for heteroclinic and periodic orbits. We also implement this
 approach numerically by means of the steepest descent method.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}

\section{Introduction}

In this work, we consider the area generated by a curve under the
action of the flow defined by an autonomous differential equation.
To fix ideas we work in $\mathbb{R}^2$ although the method is
 completely general.

This approach has been considered by Smith \cite{Smith} as well as
by Li and Muldowney \cite{Li-Muldowney0, Li-Muldowney1,
Li-Muldowney2, Li-Muldowney3, Li-Muldowney4}. The first author
establishes bounds on the Hausdorff dimension of $\omega$-limit
sets. The last two authors extend the Bendixson's criterion for
nonexistence of periodic solutions and give stability results for
some periodic orbits based on area functionals.


The main goal of this paper is to obtain invariant sets considering
a functional that describes the area that generates a curve under a
flow and observing that the functional vanishes on invariant curves
and, conversely, whenever its infimum is achieved at a curve of
finite length, then this curve is invariant (see Theorem \ref{T:1}).
This variational principle is used in order to obtain a numerical
implementation that will give us an approximation to this curve by
employing the steepest descent method.
This implementation is done in the cases of heteroclinic orbits and
limit cycles.


\section{The area functional and invariant sets}

Given a dynamical system defined on a set $X$, its evolution is
described by a monoparametric family of operators $\{S(t)\}_{t\geq
0}$, that maps $X$ into itself and possesses the standard
semigroup properties:
\begin{gather*}
S(t+s)=S(t)\cdot S(s),\quad \forall s,t\geq 0, \\
S(0)=I.
\end{gather*}
The basic property of this family is that $S(t)$ is a continuous
operator  from $X$ to $X$.

A set $Y \subset X$ is \emph{positively invariant} for $S(t)$ if
\begin{equation}
S(t)Y \subset Y, \quad \forall t>0\,,
\end{equation}
and \emph{negatively invariant} if
\begin{equation}
S(t)Y \supset Y, \quad \forall t>0.
\end{equation}
When the set is both positively and negatively invariant we say
that it is an invariant set.
Thus, a set $Y \subset X$  is an invariant set for the semigroup
$\{S(t)\}_{t\geq 0}$ if
\begin{equation}
S(t)Y=Y, \quad \forall t\geq 0.
\end{equation}
The simplest examples of invariant sets are  equilibrium
points, heteroclinic orbits and limit cycles.

A \emph{heteroclinic orbit} is a solution in phase space that joins
two different equilibrium points. One of them the $\omega$-limit set
of this orbit and the other its $\alpha$-limit set.

A \emph{stable limit cycle} is a closed trajectory in phase space
such that any  trajectory nearby spirals into it when $t$
increases.

To describe the characterization of invariant sets through
an area functional, let us consider the case where $S(t)$ is
generated by the differential equation
\begin{equation}
\dot{x}=f(x),\label{E:1}
\end{equation}
where \emph{f} is a vector field in $\mathbb{R}^2$.

Thus, for a curve $\gamma(s)$, where $\gamma$ is assumed to be
parameterized by arc length $\gamma:[s_0,s_1] \to \mathbb{R}^2$,
$S(t)\gamma (s)$, for some fixed $t\in [0,T]$ represents the
\emph{displaced curve} under the flow generated by $S$ at time
$t$. Analogously, if we consider $S(t)\gamma(s)$ for all $t\in
[0,T]$ and all $s\in [s_0,s_1]$ we obtain a strip on the phase
plane, as shown in figure  \ref{fig1}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig1} 
\caption{Strip generated by a curve under a flow.}\label{fig1}
\end{center}
\end{figure}

Assume that  $x_0$ and $x_1$ are fixed points of \eqref{E:1}
(i.e. $f(x_0)=0,\,f(x_1)=0$) and $\gamma(s)$ is a continuous curve
that joins these points, with
$\gamma(s_0)=x_0$ and $\gamma(s_1)=x_1$. Let $\Gamma$ be the set
of these curves.
Thus, ``an element of area'', is given by $|f(\gamma(s))\times
\dot{\gamma}(s)|$, as shown in figure \ref{fig2}.


\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig2} 
\caption{Differential element for the area functional.}\label{fig2}
\end{center}
\end{figure}

Therefore the area, $A_t (\gamma(s))$, is given by:
\begin{equation}
A_t(\gamma(s))=\Big(\int_{s_0}^{s_1}|f(\gamma(s))\times
\dot{\gamma}(s)|\,ds \Big)+o(t).
\end{equation}
It is clear that if the curve is invariant then $A_t(\gamma(s))=0$.
The converse is also true under appropriate assumptions.

Since the functional is bounded from below (clearly $A_t\geq 0$),
then a natural way to study invariant sets is through the
minimization of $A_t$.

Let be  $A^*(\gamma)=\int_{s_0}^{s_1}|\dot{\gamma}\times f|\,ds$,
$l$ fixed and $\Gamma_l=\{\gamma \in \Gamma
|\mbox{length}(\gamma)\leq l\}$. If the infimum of the area
functional is zero in $\Gamma$ and equal to the infimum in $\Gamma_l$
then this infimum will be an invariant curve, as established in the
following theorem.

\begin{theorem}\label{T:1}
If there is an $l$ such that $\inf_{\Gamma}A=\inf_{\Gamma_l}A^*$
and $\inf_{\Gamma}A^*=0$, then $\inf_{\Gamma}A^*$ is attained in
some $\gamma^{*}$ that is invariant under the flow generated by
\eqref{E:1}.
\end{theorem}

\begin{proof}
Without loss of generality, we can work in $\Gamma_l$. ($\Gamma_l$ is
bounded and equicontinuous by hypothesis).
Let $\gamma_n$ be a minimizing sequence in $\Gamma_l$. By Arzela -
Ascoli's theorem there is a subsequence $\gamma_{n_j}$ such that
$\gamma_{n_j}\stackrel{unif} {\longrightarrow} \gamma^{*}$.
By Fatou's lemma
\[
0\leq \int_{s_0}^{s_1}|\dot{\gamma}^{*}\times
f(\gamma^{*})|=\int_{s_0}^{s_1}\liminf|\dot{\gamma}_{n_j}\times
f(\gamma_{n_j})| \leq \lim \int_{s_0}^{s_1}|\dot{\gamma}_{n_j}\times
f(\gamma_{n_j})| \to 0.
\]
Therefore, $|\dot{\gamma}^{*}\times f(\gamma^{*})|=0$, $\forall$ s
$\in [s_0,s_1]$. This implies that $\dot{\gamma}^{*}$ is parallel to
$f(\gamma^{*})$ and so $\gamma^{*}$ is invariant for the flow
generated by \eqref{E:1}.
\end{proof}

\begin{remark} \rm
In general, the equality $\inf_{\Gamma}A=\inf_{\Gamma_l}A^*$  does
not hold, since $\Gamma_l \subset \Gamma$. This fact only implies
that $\inf_{\Gamma}A \leq \inf_{\Gamma_l}A^*$.
\end{remark}

\begin{remark} \rm
Consider the phase portrait and in particular the orbits  shown in
figure \ref{fig3}. In this case, any curve that joins $x_0$ and
$x_1$ with finite length generates a region of strictly positive area. The infimum is not attained
in the class of curves with finite length, and therefore can not be
an invariant curve. Similar examples can be constructed by
considering the irrational flow on the torus.
\end{remark}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig3} 
\caption{Phase portrait of a nonrectifiable curve wich is not invariant.}
\label{fig3}
\end{center}
\end{figure}

\section{Numerical implementation}

Theorem  \ref{T:1} guarantees that the infimum  of the area
functional is attained at some curve. However, we do not have an
analytic or approximate expression. In what follows, we show that
the variational principle can be used numerically to approximate
this curve.
To obtain this approximation, we use the functional
\begin{equation}
\tilde{A}(\gamma(s))=\int_{s_0}^{s_1}|f(\gamma(s))\times
\dot{\gamma}(s)|^{2}\,ds,\label{E:2}
\end{equation}
because it is convenient from the numerical point of view, but
clearly if $\tilde{A}$ vanishes at a curve, so does $A_t$.

Consider as an application the physical pendulum and corresponding
 separatrices; i. e., the two heteroclinic orbits that join
$(-\pi,0)$ with $(\pi,0)$.
The equation reads
\begin{equation}
\ddot{x}=-\sin x\label{E:3}
\end{equation}
or equivalently,
\begin{equation}
(\dot{x},\dot{y})=(y,-\sin x).\label{E:4}
\end{equation}
Note that \eqref{E:3} and \eqref{E:4} constitute a conservative
system (see \cite{Arnold}), in which the total energy
\[
E=\frac{1}{2}\dot{x}^{2}+U(x),
\]
where
\[
U(x)=\int_0^{x}\sin \eta\, d\eta=1-\cos x,
\]
is conserved. Its phase portrait is shown in figure \ref{fig4}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4} %pendulo.eps}
\caption{Phase portrait of  system \eqref{E:4}.}\label{fig4}
\end{center}
\end{figure}

The separatrix lies on the level curve $E=2$ and is given by
\begin{equation}
\frac{1}{2}y^2-\cos x-1=0.\label{E:5}
\end{equation}
Therefore, the upper part of \eqref{E:5} is parameterized by
\begin{equation}
(s,\sqrt{2\cos s+2}),\; -\pi \leq s \leq \pi.\label{E:6}
\end{equation}

Now, suppose that $(s,v(s))$ is an  approximation to \eqref{E:6},
which  satisfies $v(\pm \pi)=0$. Under this assumption, we calculate
the area functional $\tilde{A}$ given by \eqref{E:2}
generated by the flow \eqref{E:4}.
A straightforward calculation gives that  $\tilde{A}(\gamma(s))$ for
$\gamma(s)=(s,v(s))$ is
\begin{equation}
\tilde{A}(\gamma(s))=\int_{-\pi}^\pi (\sin
s+v(s)v'(s))^2ds.\label{E:30}
\end{equation}
Now, we take
\begin{equation}
v(s)=-(s+\pi)(s-\pi)\sum_{i=0}^n a_is^i,\label{E:7}
\end{equation}
which clearly satisfies $v(\pm \pi)=0$.

Note that in \eqref{E:30}, the $a_i$'s values mentioned in
\eqref{E:7} are taken as \emph{variables}. For some choice of these
values, we obtain a parameterized curve. Thus, we want to find the
choice of these coefficients such that \eqref{E:7} be the
\emph{best approximation} to the invariant curve.


The steepest descent method (or gradient method) is one of the
simplest and the most useful minimization methods for unconstrained
optimization, which is based in successive approximations and it
``follows'' minus gradient of the function, because it represents
the direction in which the function decreases most quickly.
The iteration is:
\begin{equation}
x_{i+1}=x_i-\epsilon \nabla \tilde{A}(x_i),\label{E:8}
\end{equation}
where $x_i$ is the $i$-th iteration , $\epsilon$ is a small value
and $\nabla \tilde{A}(x_i)$ represents the gradient of the area
functional evaluated at the $i$-th
iteration.


To determine the convergence of \eqref{E:8}, we use the
following criterion:
\begin{enumerate}
\item Select a parameter $\delta >0$.
\item Calculate $\mathbf{c}=\nabla f(x_i)$.
If $\|\mathbf{c}\|<\delta$ then
$x_i$ is the sought approximation, else repeat.
\end{enumerate}
For  a more detailed description about this method, see for example
\cite{Bonnans, Peressini}.

Now, we choose $n=2$ in \eqref{E:7}, due to the symmetry of the
separatrix with respect to the y-axis,  and consider the
\emph{initial guess} $(2.03,0,-0.27)$, $\epsilon=0.000009$,
and $\delta=0.0011$.
These values were chosen because they represent a good approximation
 to the exact values for the Taylor expansion of \eqref{E:6}.
Observe that
the corresponding value of the norm in the 2000-th iteration is
$0.001055164525$. This fact shows convergence of this method
according to the criterion previously established.
These results are shown in figures \ref{fig5} to \ref{fig8}. The
invariant curve is below the approximating curves. Apparent changes
is size are due to scaling.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig5} % figure8.eps
\caption{Initial curve versus invariant curve.}\label{fig5}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig6} % figure9.eps
\caption{First iteration versus invariant curve.}\label{fig6}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig7} 
\caption{500-th iteration versus invariant curve.}\label{fig7}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig8} % figure13.eps
\caption{1000-th iteration versus invariant curve.}\label{fig8}
\end{center}
\end{figure}

On other hand, if we consider (clearly does not satisfy $v(\pm
\pi)=0$)
\begin{equation}
v(s)=\sum_{k=0}^n a_ks^k \label{E:9}
\end{equation}
instead of \eqref{E:7}, we choose $n=11$ (a truncated Taylor
series) and the initial guess:
\begin{equation}\label{E:10}
(2.03,0,-0.27,0,6\times 10^{-3},0,-5\times
10^{-5},0,1.95 \times 10^{-7},0,-6 \times 10^{-10})
\end{equation}
We fix $\epsilon=3\times 10^{-10}$. A graphical representation of
the initial curve is shown in figure \ref{fig9}. Then we obtain the
curve after 10 iterations and compare it with the
corresponding invariant curve (see figure \ref{fig10}).

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig9} % figure3.eps
\caption{Initial curve versus invariant curve. }\label{fig9}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig10} % figure2.eps}
\caption{Graph of 10-th iteration versus invariant
curve.}\label{fig10}
\end{center}
\end{figure}

When considering \eqref{E:9}-\eqref{E:10}, the value of the
gradient of the area functional after 10 iterations is
$2356.113011$. This fact shows that it is not a convenient
choice.

\section{Application to limit cycles}

In this section we shall see how to extend the previous ideas to the
study of limit cycles.
Since these sets are closed curves, we use truncated Fourier series
instead of polynomials in order to approximate them.
Let us consider:
\begin{gather}
\dot{x_1}=-x_2+x_1(x_1^2+x_2^2-1)\label{E:12}\\
\dot{x_2}=x_1+x_2(x_1^2+x_2^2-1).\label{E:13}
\end{gather}
By using polar coordinates, we see that $x_1^2+x_2^2=1$
is a stable limit cycle (see figure \ref{fig11}).

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig11} % estable.eps
\caption{Phase portrait of \eqref{E:12}-\eqref{E:13}.}\label{fig11}
\end{center}
\end{figure}

To characterize this limit cycle by minimization of the
functional $\tilde{A}$ given by \eqref{E:2}, we use the usual
parametrization
\begin{equation}
\gamma(s)=(\cos s,\sin s),\; 0\leq s \leq 2\pi. \label{E:14}
\end{equation}
Assume that
\begin{gather}
u(s)=\sum_{k=0}^2 a_k\cos ks +\sum_{k=1}^2b_k\sin ks,\label{E:15}\\
v(s)=\sum_{k=0}^2 a_k'\cos ks+ \sum_{k=1}^2b_k'\sin ks, \label{E:16}
\end{gather}
is a finite approximation to \eqref{E:14}.
Therefore,
\begin{equation}
\tilde{A}(\gamma(s))=\int_0^{2\pi}
((-v+u(1-u^2-v^2))v'-(u+v(1-u^2-v^2))u')^2\,ds.
\end{equation}
We choose $\epsilon=\mbox{0.005}$ and $\delta=0.0006$ in
order to apply steepest descent to this functional.
The corresponding initial guess is
\begin{equation}
(0,1.3,0.1,0.1,0.1,0.1,0.1,0.01,1.3,0.1).
\label{E:17}
\end{equation}
Since $|\nabla \tilde{A}(x_{110})|<\mbox{0.0006}$, we obtain
convergence as shown in figures \ref{fig12}-\ref{fig15}.


\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig12} % stable1.eps
\caption{Initial guess versus stable limit cycle of
\eqref{E:12}-\eqref{E:13}.}\label{fig12}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig13} % stable2.eps
\caption{First iteration versus stable limit cycle.}\label{fig13}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig14} % stable3.eps
\caption{10-th iteration versus stable limit cycle.}\label{fig14}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig15} % stable5.eps
\caption{111-th iteration versus stable limit cycle.}\label{fig15}
\end{center}
\end{figure}

On other hand, if we consider
\begin{gather}
\dot{x_1}=-x_2+x_1(x_1^2+x_2^2-1)^2\label{E:18}\\
\dot{x_2}=x_1+x_2(x_1^2+x_2^2-1)^2\label{E:19}
\end{gather}
instead of \eqref{E:12}-\eqref{E:13}, then $x_1^2+x_2^2=1$ is an
unstable limit cycle (see figure \ref{fig16}).

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig16} %inestable.eps
\caption{Phase portrait of \eqref{E:18}-\eqref{E:19}.}\label{fig16}
\end{center}
\end{figure}

As before, we use \eqref{E:15}-\eqref{E:16},
$\epsilon=0.005,\, \delta=0.0055$ and the initial
guess:
\[
(0,1.2,0.1,0.1,0.1,0.1,0.1,0.01,1.1,0.1)\,.
\]

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig17} % unstable0.eps
\caption{Initial guess versus unstable limit cycle of
\eqref{E:18}-\eqref{E:19}.}\label{fig17}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig18} % unstable1.eps
\caption{First iteration versus unstable limit cycle.}\label{fig18}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig19} %unstable10.eps
\caption{10-th iteration versus unstable limit cycle.}\label{fig19}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.3\textwidth]{fig20} % unstable400.eps
\caption{110-th iteration versus unstable limit cycle.}\label{fig20}
\end{center}
\end{figure}

Since $|\nabla \tilde{A}(x_{400})|<0.0055$, we obtain
convergence as shown in figures \ref{fig17}--\ref{fig20}.

\begin{thebibliography}{00}

\bibitem{Arnold} Arnold V. I., \emph{Mathematical Methods of Classical
Mechanics}, Graduate Text in Mathematics, Springer - Verlag, 1989.

\bibitem{Bonnans} Bonnans J. F., et. al. \emph{Numerical
Optimization: Theoretical and Practical Aspects}, Universitext,
Springer-Verlag, 2006.

\bibitem{Li-Muldowney0} Li Yi M. and Muldowney James S., \emph{On R.
A. Smith's autonomus convergence theorem}, Rocky Mountain Journal of
Mathematics, volume 25, number 1, 365--381, 1995.

\bibitem{Li-Muldowney1} Li Yi M. and Muldowney James S., \emph{On
Bendixson's Criterion}, Journal of Differential Equations 106,
27--39, 1993.

\bibitem{Li-Muldowney2} Li Yi M. and Muldowney James S.,
\emph{Evolution of surface functionals and differential equations},
Ordinary and Delay Differential Equations, Longman Scientific and
Technical, 144--148, 1992.

\bibitem{Li-Muldowney3} Li Yi M. and Muldowney James S., \emph{Lower
bounds for the Haussdorff dimension of attractors}, Journal of
Dynamics and Differential Equations 7, 455--467, 1995.

\bibitem{Li-Muldowney4} Li Yi M. and Muldowney James S.,
\emph{Dynamics of differential equations on invariant manifolds},
Journal of Differential Equations 168, 295--320, 2000.

\bibitem{Peressini} Peressini A. L., et. al., \emph{The Mathematics
of Nonlinear Programming}, Undergraduate Text in Mathematics,
Springer-Verlag, 1988.

\bibitem{Smith} Smith R. A., \emph{Some
Applications of Hausdorff dimension inequalities for ordinary
differential equations}, Proy. Roy. Soc. Edinburg Sect. A 104,
235--259, 1986.

\end{thebibliography}

\end{document}