\chapter{Chaotic Systems and Fractals}
\label{chp:chaos}
\begin{flushright}
\parbox{5.5in}{An apparent paradox is that chaos is deterministic, generated by fixed rules which do not themselves involve any elements of change. We even speak of deterministic chaos. In principle, the future is completely determined by the past; but in practice small uncertainties, much like minute errors of measurement which enter into calculations, are amplified, with the effect that even though the behavior is predictable in the short term, it is unpredictable over the long run.~\cite[page~11]{peitgen92}}
\end{flushright}

\bigskip

\section{Chaotic Systems Are Dynamical Systems With Wild Behavior}
The objective of this chapter is to show the reader ``the wild things that simple nonlinear equations can do'' \cite[page~459]{may76}. These simple nonlinear models with very complicated dynamics\index{simple nonlinear models!with very complicated dynamics} are examples of {\em chaotic systems\/},\index{systems!chaotic}\index{chaotic systems} which are deterministic nonlinear dynamical systems that exhibit a random-like behavior.\index{deterministic nonlinear dynamical systems!with a random-like behavior} A discrete-time model of population dynamics discussed in the next section, Newton's method for finding roots of certain polynomials, or a jetting balloon after its release are examples of dynamical chaotic systems. Trajectories of chaotic dynamical systems are sensitive to initial conditions in the sense that starting from slightly different initial conditions the trajectories diverge exponentially. Another property of a chaotic dynamical system is {\em loss of information\/} about initial conditions. We follow Moon~\cite[page~4]{moon87} to explain this property. Suppose that we are given a state-plane portrait of a chaotic dynamical system in the $(x,y)$-plane. Assume that we specify the $x$-coordinate with accuracy $\Delta x$ and the $y$-coordinate with accuracy $\Delta y$. Then, we could form a grid in the $(x,y)$-plane, shown in Figure~\ref{fig:inf_loss}, where each area would be of size $\Delta x\,\Delta y$. Thus, if we select initial conditions to the stated accuracy, then we know that these initial conditions are located somewhere in the shaded box in the $(x,y)$-plane corresponding to the initial time $t_0$. For a chaotic system this uncertainty grows in time to, say, $N$ boxes at a later time $t_1$, as illustrated in Figure~\ref{fig:inf_loss}. The uncertainty growth,\index{uncertainty growth} as mentioned above, is exponential\index{chaotic systems!uncertainty growth exponential} and for continuous-time systems can be described by
\begin{equation}
N(t)\approx e^{\mu (t-t_0)}N(t_0),
\label{eq:infolossinchs}
\end{equation}
where the constant $\mu$ is a measure of sensitivity with respect to initial conditions.\index{sensitivity!with respect to initial conditions} Lyapunov exponents, that we discuss in a subsequent section, can be used to measure the initial condition sensitivity.\index{Lyapunov exponents!measure the initial condition sensitivity}
\begin{figure}
\centering
\includegraphics[scale=0.5]{lossinfo.eps}
%\vspace{3in}
%\hskip1.25in\special{eps:c:/mybook/lossinfo.eps y=3in}
%\centerline{\psfig{figure=lossinfo.eps height=2.5in}}
%\centerbmp{4in}{3in}{c:/mybook/lossinfo.bmp}
%\hskip.7in\special{bmp:c:/mybook/lossinfo.bmp x=4.57in y=3.5in}
\caption{An illustration of the growth of uncertainty in a chaotic dynamical system.}
\label{fig:inf_loss}
\end{figure}
Chaotic behavior can arise in low-order nonlinear dynamical systems. In the case of continuous-time systems, chaotic behavior can only occur for systems of dimension three or higher. However, in discrete-time systems, chaos can occur even in one-dimensional systems, as we show in the following section.

\section{Chaotic Behavior of the Logistic Equation}

\subsection{The Logistic Equation---An Example From Ecology}
The so called logistic equation, or logistic map,\index{logistic equation}\index{logistic map|see{logistic equation}} whose behavior we analyze in this section, arises, for example, in the analysis of the growth of a population.
\begin{defn}
A map\index{map} is a function whose domain (input) space and range (output) space are the same. Let $\b x(0)$ be a point in $\mb{\Bb R}^n$ and $\b f$ be a map on $\mb{\Bb R}^n$. The orbit,\index{orbit|see{trajectory}} or trajectory,\index{trajectory} of $\b x(0)$ under $\b f$ is the set of points
\[
\left\{ \b x(0),\b f(\b x(0)),\b f^2(\b x(0)),\ldots,\b f^k(\b x(0)),\ldots \right\},
\]
where $\b f^k$ is the $k$-th iterate of $\b f$. The starting point, $\b x(0)$, of the orbit is called the initial value of the orbit.\index{initial value!of the orbit}
\end{defn}\markdefn

We first analyze a situation where there are no predators and there is an unlimited supply of food. For the sake of argument, we consider a population of rabbits. Assume that on the average each rabbit gives birth to $\sigma$ new rabbits in one unit of time, and that no rabbit dies. Let $N(k)$ be the number of rabbits at time $t=k$, then the change in the population over a time period, between $k$ and $k+1$ is
\begin{equation}
N(k+1)-N(k)=\sigma N(k).
\label{eq:over_period}
\end{equation}
Thus,
\begin{equation}
N(k+1)=(1+\sigma)N(k).
\label{eq:over_per}
\end{equation}
Equation~(\ref{eq:over_per}) is a very simple linear model of population growth.\index{linear model!of population growth} Given the number $N(0)$ of rabbits at time $k=0$, the number of rabbits at time $k$ is
\begin{equation}
N(k)=(1+\sigma)^k N(0).
\label{eq:grows_exp}
\end{equation}
One can see from~(\ref{eq:grows_exp}) that the linear model given by~(\ref{eq:over_per}) predicts the exponential grow of population.\index{linear model!predicts the exponential grow of population} We now try to make the above model more realistic. We assume, as before, that the number of births is proportional to the size of the population. Let $b$ be the birth rate, then the number of births at time $k$ is $b N(k)$. We also assume that the number of deaths is proportional to the population size. Let $d$ be the death rate. Then, the number of deaths at time $k$ is $d N(k)$. The change in population in a time period between times $k$ and $k+1$ is
\begin{equation}
N(k+1)-N(k)=bN(k)-dN(k).
\label{eq:over_period2}
\end{equation}
We rewrite equation~(\ref{eq:over_period2}) as
\begin{equation}
N(k+1)=(1+b-d)N(k).
\label{eq:over_per2}
\end{equation}
Let $\sigma=b-d$. Then, equation~(\ref{eq:over_per2}) has the same form as equation~(\ref{eq:over_per}). In both cases, the population grows exponentially. We refer to the constant $\sigma$ as the {\em unrestricted growth rate\/}.\index{unrestricted growth rate} Over short periods of time, the above equation may be sufficient for purposes of an analysis of a population. Over long periods of time, the growth rates may change as the size of a population changes.

We alter the above linear models by introducing the concept of saturation.\index{saturation} We first assume that the environment of the population can only support a certain number of species. We denote this number by $K$, and call it the {\em carrying capacity\/}.\index{carrying capacity} Thus, if $N(k)>K$, then there will not be enough food or space available, and more species will die than will be born. In such a situation the growth rate should be negative. If the population at time $k$ is $N(k) < K$, then the growth rate should be positive. Finally, if $N(k)$ is much smaller than the carrying capacity, then the growth rate should be close to the unrestricted growth rate $\sigma$. A possible function that satisfies the above requirements has the form
\begin{equation}
\sigma \left( 1-\frac{N(k)}{K}\right).
\end{equation}
The change in the size of population then can be described as
\begin{equation}
N(k+1)-N(k)=\sigma \left( 1-\frac{N(k)}{K}\right)N(k),
\end{equation}
or equivalently
\begin{equation}
N(k+1)=(1+\sigma)N(k)-\frac{\sigma}{K}N^2(k).
\label{eq:1d_logistic}
\end{equation}
Equation~(\ref{eq:1d_logistic}) is called the {\em logistic equation\/}.\index{logistic equation} The term $\frac{\sigma}{K}N^2(k)$ is referred to as the {\em damping term\/},\index{damping term!in logistic equation} because it dampens the growth of the population. To simplify~(\ref{eq:1d_logistic}), we scale the variable $N(k)$ by introducing the new variable
\[
x(k)=\frac{\sigma}{\rho\,K}N(k),
\]
where
\[
\rho=1+\sigma.
\]
Then, equation~(\ref{eq:1d_logistic}) becomes
\begin{equation}
\fbox{$\displaystyle x(k+1)=\rho x(k)(1-x(k))$}
\label{eq:1d_log_work}
\end{equation}
We will show that if $0<\rho \le 4$ and $0\le x(0) \le 1$, then $0\le x(k)\le 1$ for $k=1,2,\ldots$. Let
\[
f(x)=\rho x(1-x).
\]
Then,
\[
\frac{d f(x)}{d x}=\rho-2\rho x=0
\]
for $x=1/2$ and $f_{max}(x)=\rho/4$. Hence, for $k=1,2,\ldots$,
\[
0\le x(k)\le 1\quad\mb{if}\quad 0<\rho\le 4\quad\mb{and}\quad 0\le x(0)\le 1.
\]

\subsection{Stability Analysis of the Logistic Map}

We begin our investigation of the logistic equation~(\ref{eq:1d_log_work}) by seeking its equilibrium states (fixed points).\index{equilibrium states|see{fixed points}}\index{fixed points} These can be found by setting $x(k+1)$ equal to $x(k)$ in~(\ref{eq:1d_log_work}) to obtain
\[
x(k)\left(\rho x(k)+\left(1-\rho\right)\right)=0.
\]
It follows that the equilibrium solutions are
\[
x^{*}=0\quad\mb{and}\quad x^{*}=\frac{\rho-1}{\rho}.
\]
We study attracting properties of the fixed points\index{fixed points!attracting properties of} as a function of the parameter $\rho$. Roughly speaking, a fixed point $x^{*}=f(x^{*})$ of the system $x(k+1)=f(x(k))$ is attracting if, for a starting point close to $x^{*}$, the next iterate is still closer. Formally, a fixed point $x^{*}$ is attracting if\index{attracting fixed point}\index{fixed point!attracting} there is an $\varepsilon >0$ such that when $\left|x(0)-x^{*}\right|<\varepsilon$, then $\lim_{k \to \infty}x(k)=x^{*}$.

A fixed point $x^{*}$ is repelling if for any $\varepsilon >0$ no matter how small, if $0<\left|x(0)-x^{*}\right|<\varepsilon$, then $\left|x(k)-x^{*}\right|>\varepsilon$ for some $k$.
\begin{thm}
\label{thm:atractrepel}
Suppose $x^{*}=f(x^{*})$ is a fixed point of $x(k+1)=f(x(k))$, where $f$ is continuously differentiable.\index{continuously differentiable} Then, $x^{*}$ is attracting if
\[
\left|f^{'}(x^{*})\right|<1,
\]
and is repelling if\index{fixed point!repelling}
\[
\left|f^{'}(x^{*})\right|>1.
\]
The case when $\left|f^{'}(x^{*})\right|=1$ is inconclusive.
\end{thm}
{\bf Proof}\quad We first show that if $x^{*}=f(x^{*})$ is such that $\left|f^{'}(x^{*})\right|<1$, then $x^{*}$ is attracting. Since by assumption $\left|f^{'}(x^{*})\right|<1$ and $f^{'}$ is continuous, there exists an interval $(x^{*}-\varepsilon,x^{*}+\varepsilon)$ about $x^{*}$ such that for some $L>0$,
\[
\left|f^{'}(x^{*})\right|<L<1\quad\mb{and}\quad\left|f^{'}(x)\right|<L\quad \mb{for every}\quad x\in (x^{*}-\varepsilon,x^{*}+\varepsilon).
\]
Let $x(k) \in (x^{*}-\varepsilon,x^{*}+\varepsilon)$. We will show that $x(k+1)$ is closer to $x^{*}$ than $x(k)$. We have
\[
\left|x(k+1)-x^{*}\right|=\left|f(x(k))-f(x^{*})\right|
\]
because $x^{*}=f(x^{*})$ is a fixed point. Using the mean value theorem\index{mean value theorem} yields
\[
\left|x(k+1)-x^{*}\right|=\left|f(x(k))-f(x^{*})\right|=\left|f^{'}(\tilde{x})\right|\left|x(k)-x^{*} \right|
\]
for some $\tilde{x} \in (x(k),x^{*})$ or $\tilde{x} \in (x^{*},x(k))$ if $x^{*} < x(k)$. Since $x(k) \in (x^{*}-\varepsilon,x^{*}+\varepsilon)$, then $\tilde{x} \in (x^{*}-\varepsilon,x^{*}+\varepsilon)$ also. Thus, $\left|f^{'}(\tilde{x})\right|<L<1$, and therefore
\[
\left|x(k+1)-x^{*}\right|<L\left|x(k)-x^{*}\right|.
\]
From the above, we conclude that if $x(0)\in (x^{*}-\varepsilon,x^{*}+\varepsilon)$, then
\[
\left|x(k)-x^{*}\right|<L^k\left|x(0)-x^{*}\right|.
\]
Because $L<1$, we have $\lim_{k \to \infty}\left|x(k)-x^{*}\right|=0$, that is,
\[
\lim_{k\to \infty}x(k)=x^{*},
\]
which means that $x^{*}$ is an attracting fixed point.

Similarly, if $x^{*}=f(x^{*})$ is such that $\left|f^{'}(x^{*})\right|>1$, then it is a repelling fixed point.\index{repelling fixed point}\index{fixed point!repelling} Indeed, if $\left|f^{'}(x^{*})\right|>1$, then because $f^{'}$ is continuous, there is an interval $(x^{*}-\varepsilon,x^{*}+\varepsilon)$ about $x^{*}$ and a constant $L>1$ such that $\left|f^{'}(x)\right|>L$ for all $x\in (x^{*}-\varepsilon,x^{*}+\varepsilon)$. We will now show that if $x(k) \in (x^{*}-\varepsilon,x^{*}+\varepsilon)$, then $x(k+1)$ will be further away from $x^{*}$ than $x(k)$. Using the mean value theorem,\index{mean value theorem} we obtain
\begin{eqnarray}
\left| x(k+1)-x^{*}\right|&=& \left|f(x(k))-f(x^{*}) \right|\nonumber\\
&=& \left|f^{'}(\tilde{x}) \right|\left|x(k)-x^{*} \right|\nonumber\\
&>& L\left|x(k)-x^{*} \right|.
\label{eq:rateofcovod}
\end{eqnarray}
If $x(k+1)$ is no longer in $(x^{*}-\varepsilon,x^{*}+\varepsilon)$, we are done. Otherwise, we repeat the process to get
\begin{eqnarray*}
\left| x(k+2)-x^{*}\right|&=& \left|f(x(k+1))-f(x^{*}) \right|\\
&>& L\left|x(k+1)-x^{*} \right|\\
&>& L^2\left|x(k)-x^{*} \right|.
\end{eqnarray*}
If $x(k+2)$ is no longer in the interval $(x^{*}-\varepsilon,x^{*}+\varepsilon)$, we are done. If not, we repeat the process to obtain
\[
\left|x(k+m)-x^{*} \right|>L^m\left| x(k)-x^{*}\right|.
\]
Since $(x^{*}-\varepsilon,x^{*}+\varepsilon)$ is a finite interval and $L^m \to \infty$ as $m \to \infty$, at some time, $x(k+m)$ must leave the interval, which means that $x^{*}$ is a repelling fixed point. The proof of the theorem is complete.\\
\hspace*{5.5in}$\Box$

Equation~(\ref{eq:rateofcovod}) tells us that the trajectory $\{x(k)\}$ diverges exponentially from $x^{*}$ as $k\to \infty$. Thus, a trajectory with an initial condition close to $x^{*}$ for the case when $\left|f^{'}(x^{*})\right|=L>1$ will separate from $x^{*}$ at a multiplicative rate of approximately $L>1$ per iteration until it moves significantly away from $x^{*}$. In other words, the distance between $f(x(k))$ starting from $x(0)$ close to $x^{*}$ will grow exponentially, that is, this distance will be magnified approximately by $L$ for each iteration. To quantify this multiplicative rate of separation\index{multiplicative rate!of separation}\index{separation!multiplicative rate of} of points close to $x^{*}$, the Lyapunov exponent\index{Lyapunov exponent} term is introduced that we discuss in Section~\ref{sec:lyapunexpon}.

We now consider the logistic equation~(\ref{eq:1d_log_work}). Recall that $f=\rho x(1-x)$ denotes the right hand side of the logistic equation. The derivative of $f$ with respect to $x$ is
\[
f^{'}(x)=\rho(1-2x).
\]
Therefore, $f^{'}(0)=\rho$. Since $\left|f^{'}(0)\right|<1$ when $\rho <1$, by Theorem~\ref{thm:atractrepel}, the fixed point $x^{*}=0$ is attracting for $\rho <1$ and repelling for $\rho >1$. The fixed point $x^{*}=\frac{\rho-1}{\rho}$ is outside of the range of $x$ for $\rho<1$. Thus, if $0<\rho<1$, we have just one attractive fixed point. For $\rho>1$, the fixed point $x^{*}=0$ is repelling and $x^{*}=\frac{\rho-1}{\rho}$ comes in within the range of the variable $x$. The derivative of $f$ evaluated at $x^{*}=\frac{\rho-1}{\rho}$ is
\[
f^{'}\left( \frac{\rho-1}{\rho}\right)=2-\rho.
\]
By Theorem~\ref{thm:atractrepel}, the fixed point $x^{*}=\frac{\rho-1}{\rho}$ is attracting if $\left| 2-\rho\right|<1$, that is, if
\[
1<\rho<3.
\]
Thus, for $\rho=1$ there is an exchange of attracting properties\index{exchange of attracting properties} from one fixed point to the other. This change in the qualitative behavior at a certain parameter value is called a {\em bifurcation}.\index{bifurcation!change in the qualitative behavior} In Figure~\ref{fig:bifursm}, the fixed points $x^{*}=0$ and $x^{*}=\frac{\rho-1}{\rho}$ are shown and the ranges of $\rho$ for which they are attracting.
\begin{figure}
\centering
\includegraphics[scale=0.5]{bifursm.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/bifursm.bmp}
%\hskip.7in\special{bmp:c:/mybook/bifursm.bmp x=4.67in y=3.5in}
\caption{Exchange of attracting properties of the fixed points for $0<\rho\le 3.44$.}
\label{fig:bifursm}
\end{figure}
For $\rho>3$, there are no attracting fixed points, and the solutions of the logistic equation become very complex. We discuss the logistic equation dynamics for $\rho >3$ after displaying its solutions for $\rho=0.7$, $\rho=1.4$, $\rho=2.8$, and $\rho=3.6$ in Figure~\ref{fig:smorho}. The initial condition in all four cases was $x(0)=0.9$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{smorho.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/smorho.bmp}
%\hskip.7in\special{bmp:c:/mybook/smorho.bmp x=4.54in y=4in}
\caption{Solutions of the logistic equation: (a) $\rho=0.7$, (b) $\rho=1.4$, (c) $\rho=2.8$, and (d) $\rho=3.6$.}
\label{fig:smorho}
\end{figure}
Another way of displaying a solution of a first-order difference equation is a {\em cobweb}.\index{cobweb!a way of displaying a solution} To construct a cobweb for the given first-order discrete system, we first draw the function $y=f(x)$ representing the right-hand side of the difference equation. Then, we draw the straight line $y=x$. The fixed points correspond to the points of intersection of these two graphs. The solution sequence, represented by the cobweb, staring at $x(0)$ is obtained as follows. Starting from $x(0)$, draw upward vertical line to intersect the graph of $f(x)$. This corresponds to the calculation of $x(1)$. The value of $x(1)$ is transferred to the straight line $y=x$. The $x$-coordinate of the intersection is $x(1)$ which is then used to generate $x(2)$ using the above procedure. The process is repeated over and over again. In Figure~\ref{fig:cobwsro}, the cobwebs corresponding to the solutions displayed in Figure~\ref{fig:smorho} are presented.
\begin{figure}
\centering
\includegraphics[scale=0.5]{cobwsro.eps}
%\vspace{4.5in}
%\centerbmp{5.1in}{4.5in}{c:/mybook/cobwsro.bmp}
%\hskip.2in\special{bmp:c:/mybook/cobwsro.bmp x=5.1in y=4.5in}
\caption{Cobwebs for the logistic equation: (a) $\rho=0.7$, (b) $\rho=1.4$, (c) $\rho=2.8$, and (d) $\rho=3.6$.}
\label{fig:cobwsro}
\end{figure}

As we mentioned before, for $\rho>3$ there are no attracting fixed points. As $\rho$ increases, the solutions of the logistic equation exhibit increasing complexity. For $\rho$ just slightly larger than 3, the solution sequence settles down into a steady oscillation of period~2 as shown in Figure~\ref{fig:twoperc}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{twoperc.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/twoperc.bmp}
%\hskip.5in\special{bmp:c:/mybook/twoperc.bmp x=4.54in y=4in}
\caption{A two-cycle solution and its cobweb for $\rho=3.3$ and $x(0)=0.1$.}
\label{fig:twoperc}
\end{figure}
Then, as $\rho$ is further increased, periodic solutions of period $4, 8, 16,\ldots$ appear. In Figure~\ref{fig:forperc}, a solution of period~4 is shown corresponding to $\rho=3.5$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{forperc.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/forperc.bmp}
%\hskip.6in\special{bmp:c:/mybook/forperc.bmp x=4.54in y=4in}
\caption{A period-four solution and its cobweb for $\rho=3.5$ and $x(0)=0.1$.}
\label{fig:forperc}
\end{figure}
For $\rho=3.9$, the solution is shown in Figure~\ref{fig:chaos2}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{chaos2.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/chaos2.bmp}
%\hskip.6in\special{bmp:c:/mybook/chaos2.bmp x=4.54in y=4in}
\caption{A periodic solution and its cobweb for $\rho=3.9$ and $x(0)=0.1$.}
\label{fig:chaos2}
\end{figure}

We next analyze the solutions for $\rho>3$. To simplify the analysis, we first make the substitution
\begin{equation}
x(k)=z(k)+\frac{\rho-1}{\rho},
\label{eq:movetoo}
\end{equation}
which is equivalent to shifting the fixed point\index{shifting the fixed point}\index{fixed point!shifting} $x^{*}=\frac{\rho-1}{\rho}$ to the origin in the $(\rho,z)$-coordinates. Note that
\begin{equation}
x(k+1)=z(k+1)+\frac{\rho-1}{\rho}.
\label{eq:movetoof}
\end{equation}
Substituting~(\ref{eq:movetoo}) and~(\ref{eq:movetoof}) into the logistic equation~(\ref{eq:1d_log_work}) yields
\begin{equation}
z(k+1)+\frac{\rho-1}{\rho}=\rho\left( z(k)+\frac{\rho-1}{\rho}\right)\left(1- z(k)-\frac{\rho-1}{\rho}\right).
\label{eq:nowinz}
\end{equation}
After performing some manipulations, we represent~(\ref{eq:nowinz}) as
\begin{equation}
\fbox{$\displaystyle z(k+1)=(2-\rho)z(k)-\rho z^2(k)$}
\label{eq:inzsimple}
\end{equation}
The fixed points of~(\ref{eq:inzsimple}) are $z^{*}=0$ corresponding to $x^{*}=\frac{\rho - 1}{\rho}$, and $z^{*}=\frac{1-\rho}{\rho}$ corresponding to $x^{*}=0$. We investigate the attracting properties of the fixed point $z^{*}=0$ for $\rho \ge 0$. We see from Figure~\ref{fig:twoperc} that for $\rho$ just greater than 3, the solution of the logistic equation goes to an attracting period-two cycle.\index{attracting!period-two cycle} Two numbers $p$ and $q$ form a period-2 cycle\index{period-2 cycle!attracting}\index{attracting!period-2 cycle} for a first-order discrete systems if, when $x(k)=p$, then $x(k+1)=q$, $x(k+2)=p$, and so on. A period-2 cycle is attracting\index{attracting!2-cycle} if there are intervals $(p_1,p_2)$ and $(q_1,q_2)$ about $p$ and $q$, respectively, such that
\begin{itemize}
\item[(i)] if $x(0)\in (p_1,p_2)$, then $\lim_{k \to \infty}x(2k)=p$ and $\lim_{k \to \infty}x(2k+1)=q$;
\item[(ii)] if $x(0)\in (q_1,q_2)$, then $\lim_{k \to \infty}x(2k)=q$ and $\lim_{k \to \infty}x(2k+1)=p$.
\end{itemize}
The fixed points $p$ and $q$ of the period-2 cycle of equation~(\ref{eq:inzsimple}) can be found by solving the equation
\begin{equation}
z=g(g(z)),
\label{eq:solforcyctwo}
\end{equation}
where $g=(2-\rho)z-\rho z^2$. We can use MATLAB's function \verb!compose!\index{MATLAB's function!compose} to form~(\ref{eq:solforcyctwo}), and then use the function \verb!solve! to solve~(\ref{eq:solforcyctwo}). We obtain four roots. Two of them are the fixed points that we analyzed above. The remaining two roots are the fixed points of the period-2 cycle. They are:
\begin{equation}
p=-\frac{1}{2\rho}\left( \rho-3+\sqrt{\rho^2-2\rho - 3}\right),
\label{eq:pptofcyc}
\end{equation}
and
\begin{equation}
q=-\frac{1}{2\rho}\left( \rho-3-\sqrt{\rho ^2-2\rho - 3}\right).
\label{eq:qptofcyc}
\end{equation}
Once we have found the 2-cycle, we next determine the range of $\rho$ for which it is attracting. A 2-cycle $p$ and $q$ of the system $z(k+1)=g(z(k))$ is attracting if the fixed points $p$ and $q$ of the map $y(k+1)=h(y(k))$ are attracting, where $h(y)=g(g(y))$. By Theorem~\ref{thm:atractrepel}, the fixed points $p=h(p)$ and $q=h(q)$ are attracting if $\left|h^{'}(p)\right|<1$ and $\left|h^{'}(q)\right|<1$. Using the chain rule, we obtain
\begin{equation}
h^{'}(y)=\frac{d}{dy}g(g(y))=g^{'}(g(y))g^{'}(y).
\label{eq:chruletoh}
\end{equation}
Since $g(p)=q$ and $g(q)=p$, we have
\begin{equation}
h^{'}(p)=g^{'}(g(p))g^{'}(p)=g^{'}(q)g^{'}(p).
\label{eq:chrulrp}
\end{equation}
Similarly,
\begin{equation}
h^{'}(q)=g^{'}(p)g^{'}(q).
\label{eq:chrulrq}
\end{equation}
Therefore, $\left|h^{'}(p)\right|<1$ and $\left|h^{'}(q)\right|<1$ when
\begin{equation}
\left|g^{'}(p)g^{'}(q)\right|<1.
\label{eq:gpgqlone}
\end{equation}
Using Theorem~\ref{thm:atractrepel} and the arguments similar to the ones above, we conclude that the 2-cycle represented by $p$ and $q$ is repelling if
\begin{equation}
\left|g^{'}(p)g^{'}(q)\right|>1.
\label{eq:gprepel}
\end{equation}
We next substitute~(\ref{eq:pptofcyc}) and~(\ref{eq:qptofcyc}) into $g^{'}(p)g^{'}(q)$, and evaluate the resulting expression using MATLAB's function \verb!subs!,\index{MATLAB's function!subs} to obtain
\begin{equation}
g^{'}(p)g^{'}(q)=\left(2-\rho-2\rho p \right)\left(2-\rho-2\rho q \right)=4-\rho^2+2\rho.
\label{eq:gpgqprime}
\end{equation}
Thus, the 2-cycle represented by the pair $(p,q)$ is attracting if
\begin{equation}
-1<4-\rho^2+2\rho<1,
\label{eq:monepone}
\end{equation}
where $\rho>3$. Solving~(\ref{eq:monepone}) yields
\begin{equation}
3<\rho<1+\sqrt{6}.
\label{eq:trhoos}
\end{equation}
When $\rho>1+\sqrt{6}$, an attracting 4-cycle appears.\index{attracting!4-cycle} As can be seen from the bifurcation diagram\index{bifurcation diagram} shown in Figure~\ref{fig:lebifurd}, that as the value of $\rho$ increases the period doubling process takes place.
\begin{figure}
\centering
\includegraphics[scale=0.5]{lebifurd.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/lebifurd.bmp}
%\hskip.8in\special{bmp:c:/mybook/lebifurd.bmp x=4.54in y=4in}
\caption{Bifurcation diagram of the logistic equation.}
\label{fig:lebifurd}
\end{figure}
Each of the branches of the bifurcation diagram produces a period-2 cycle of its own.\index{branches of bifurcation diagram!produce period-2 cycles}\index{branches of bifurcation diagram branches!produce period-2 cycles} The process continues as the parameter $\rho$ approaches the value $\rho=3.5699$, until no periodic cycles are created. For $\rho=3.5699$ the behavior is not cyclic anymore. The behavior of the solutions of the logistic equation for $\rho \ge 3.5699$ is complex indeed, as can be seen, for example, by inspecting the bifurcation diagram shown in Figure~\ref{fig:lebifurd}. Various windows of ranges of $\rho>3.5699$ correspond to different periodicities of solutions as well as to aperiodic solutions.\index{aperiodic solutions} These aperiodic solutions are called {\em chaotic solutions\/} or {\em chaotic motions\/}.\index{chaotic solutions}\index{chaotic motions|see{chaotic solutions}} Here, we quote a descriptive definition of chaotic motion due to Thompson and Stewart~\cite[page~197]{thompson86}. ``... one way to define chaos is in the negative---recurrent behavior that is not an equilibrium, a cycle or even a quasi-periodic motion---there is more to be said about chaos. Chaotic motion has some aspect that is provably as random\index{chaotic motion!aspect of randomness} as a coin toss. The randomness arises from sensitive dependence\index{randomness!from sensitive dependence on initial conditions} on imperfectly known initial conditions, resulting for example in broadband noise in the power spectra of observable time series. This seems remarkable because the dynamical system needs no stochastic input to achieve this''. A precise definition of chaotic motion is given in Section~\ref{sec:lyapunexpon}.

\subsection{Period Doubling to Chaos}

We describe one of the ways in which chaotic motion can arise---by repeated flip bifurcations,\index{repeated flip bifurcations!lead to chaos} that is, through period doubling.\index{chaotic motion!by repeated flip bifurcations} Such a behavior arises for $3\le \rho \le 3.5699$. We show that there is a sequence of numbers
\begin{equation}
c_1<c_2<\cdots <c_n <\cdots
\label{eq:perdoubseq}
\end{equation}
such that when
\[
c_n<\rho<c_{n+1},
\]
then the logistic equation has an attracting $2^n$-cycle. Thus, for example, if $c_8<\rho < c_9$, then the solution will settle into an attracting $2^8$-cycle. We already have $c_1=3$ and $c_2=1+\sqrt{6}=3.4495$. We are interested in finding the remaining elements of the sequence~(\ref{eq:perdoubseq}). The sequence~(\ref{eq:perdoubseq}) is convergent, where
\[
\lim_{n \to \infty} c_n=c_{\infty}=3.5699.
\]
To find the elements of the sequence~(\ref{eq:perdoubseq}), we could, in principle, form composite maps\index{composite maps} $f^{2^n}(x)$ and use stability results from the previous subsection. We have observed that if the parameter $\rho$ is slightly bigger than $1+\sqrt{6}$, then the 2-cycle is repelling and the solution goes to an attracting 4-cycle. A cobweb of an attracting 4-cycle\index{attracting!4-cycle}\index{cobweb!of attracting 4-cycle} is shown in Figure~\ref{fig:cobwfcy}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{cobwfcy.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/cobwfcy.bmp}
%\hskip1in\special{bmp:c:/mybook/cobwfcy.bmp x=3.97in y=3.5in}
\caption{Cobweb of an attracting 4-cycle for $\rho=3.5$.}
\label{fig:cobwfcy}
\end{figure}
We could find the four points corresponding to an attracting 4-cycle by solving the algebraic equation
\begin{equation}
x=f(f(f(f(x)))),
\label{eq:fourcommap}
\end{equation}
where $f(x)=\rho x(1-x)$. A graph of the function $y=f^4(x)$ is shown in Figure~\ref{fig:fcommap}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{fcommap.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/fcommap.bmp}
%\hskip1in\special{bmp:c:/mybook/fcommap.bmp x=3.97in y=3.5in}
\caption{A graphical solution of $x=f^4(x)$, where $f(x)=3.5x(1-x)$. Marked points, corresponding to slopes less than one, form an attracting 4-cycle.}
\label{fig:fcommap}
\end{figure}
The solution to~(\ref{eq:fourcommap}) gives eight points. Four of these points, corresponding to $y=f^4(x)$ having negative slopes when intersecting the line $y=x$, form an attracting 4-cycle. Using results from the previous subsection, we could find the range of $\rho$ resulting in attracting 4-cycle solutions. We could then proceed to find the range of $\rho$ corresponding to solutions going to attracting 8-cycles, then 16-cycles, and so forth. However, this approach is only feasible for $n=1$ and $n=2$. We have to take a different approach to compute the elements of the sequence~(\ref{eq:perdoubseq}). We use an iteration technique consisting of four steps:
\begin{itemize}
\item[(1)] Finding the second iterate, $z(k+1)=g(g(z(k))$, of the logistic map~(\ref{eq:inzsimple}).
\item[(2)] Shifting one of the points of the 2-cycle of the second iterate to the origin.
\item[(3)] Neglecting cubic and quartic terms in the shifted map.
\item[(4)] Scaling the map obtained in step~(3).
\end{itemize}
Performing iterations consisting of the above four steps will allow us to find the elements of the sequence~(\ref{eq:perdoubseq}).

The second iterate of~(\ref{eq:inzsimple}) is
\begin{eqnarray}
z(k+2) &=& (2-\rho) z(k+1)-\rho z^2(k+1)\nonumber\\
&=& (2-\rho)\left((2-\rho)z(k)-\rho z^2(k) \right)-\rho \left((2-\rho)z(k)-\rho z^2(k) \right)^2\nonumber\\
&=&(2-\rho)^2 z(k)-\left((2-\rho)\rho+\rho(2-\rho)^2 \right)z^2(k)\nonumber\\
&{}&\quad + 2\rho ^2(2-\rho) z^3(k)-\rho ^3 z^4(k).
\label{eq:suarez}
\end{eqnarray}
The points $p$ and $q$, given by~(\ref{eq:pptofcyc}) and~(\ref{eq:qptofcyc}), are fixed points of the map~(\ref{eq:suarez}). We select one of the points, $p$ or $q$, to work with. Suppose that we choose the point $p$. Consider the transformation
\begin{equation}
v(k)=z(k)-p.
\label{eq:ztovtrans}
\end{equation}
Note that
\begin{equation}
v(k+2)=z(k+2)-p.
\label{eq:ztovtrans2}
\end{equation}
Substituting~(\ref{eq:ztovtrans}) and~(\ref{eq:ztovtrans2}) into~(\ref{eq:suarez}), performing manipulations, and neglecting terms of order three and four,\index{neglecting terms!of order three and four} we obtain
\begin{equation}
v(k+2)=\left(4+2\rho -\rho^2 \right)v(k)-\rho \left(\rho^2-2\rho-3 - 3\sqrt{\rho^2-2\rho -3} \right)v^2(k).
\label{eq:simpvmap}
\end{equation}
One can easily verify the above result using MATLAB's function \verb!subs!\index{MATLAB's function!subs} applied to~(\ref{eq:suarez})--(\ref{eq:ztovtrans2}). We will be studying fixed points of~(\ref{eq:simpvmap}). We replace $k+2$, the argument of $v$ on the left-hand side of~(\ref{eq:simpvmap}), with $k+1$. Let
\begin{equation}
C_1(p)=\rho^2-2\rho-3 - 3\sqrt{\rho^2-2\rho -3}.
\label{eq:c1terminv}
\end{equation}
Then,~(\ref{eq:simpvmap}) can be represented as
\begin{equation}
v(k+1)=\left(4+2\rho -\rho^2 \right)v(k)-\rho C_1(p)v^2(k).
\label{eq:almostthef}
\end{equation}
Note that we if we used $q$ in~(\ref{eq:ztovtrans}) instead of $p$, we would get
\begin{equation}
v(k+1)=\left(4+2\rho -\rho^2 \right)v(k)-\rho C_1(q)v^2(k),
\label{eq:qinstaedp}
\end{equation}
where now $C_1(q)=\rho^2-2\rho-3 + 3\sqrt{\rho^2-2\rho -3}$. In the next step, we transform the $v$-map given by~(\ref{eq:almostthef}) into the same format as the $z$-map\index{z-map} given by~(\ref{eq:inzsimple}). This allows us to apply stability properties of the $z$-map to the $v$-map.\index{v-map} We introduce an iterated parameter $\rho_1$ such that
\begin{equation}
2-\rho_1=4+2\rho_0 - \rho_0^2,
\label{eq:rhooneando}
\end{equation}
and represent~(\ref{eq:almostthef}) as
\begin{equation}
v(k+1)=\left(2-\rho_1\right)v(k)-\rho_0 C_1(p)v^2(k).
\label{eq:afterrho1inv}
\end{equation}
Finally, we scale~(\ref{eq:afterrho1inv}) using the transformation
\begin{equation}
s(k)=\frac{\rho_0 C_1(p)}{\rho_1}v(k).
\label{eq:scalevtos}
\end{equation}
Substituting~(\ref{eq:scalevtos}) into~(\ref{eq:afterrho1inv}) gives
\begin{equation}
\frac{\rho_1}{\rho_0 C_1(p)}s(k+1)=\frac{\rho_1(2-\rho_1)}{\rho_0 C_1(p)}s(k)-\frac{\rho_1^2\rho_0C_1(p)}{\rho_0^2 C_1^2(p)}s^2(k).
\label{eq:r1roco}
\end{equation}
Simplifying~(\ref{eq:r1roco}) yields
\begin{equation}
\fbox{$\displaystyle s(k+1)=\left(2-\rho_1\right)s(k)-\rho_1s^2(k)$}
\label{eq:alrsm}
\end{equation}
The above $s$-map\index{s-map} has the same format as the $z$-map given by~(\ref{eq:inzsimple}) and analyzed in the previous subsection. It follows from these results that the $s$-map given by~(\ref{eq:alrsm}) has an attracting 2-cycle\index{attracting!2-cycle} for $3<\rho_1<1+\sqrt{6}$, which corresponds to an attracting 4-cycle\index{attracting!4-cycle} of the $z$-map\index{z-map} given by~(\ref{eq:inzsimple}). Using~(\ref{eq:rhooneando}), we conclude that $3<\rho_1<1+\sqrt{6}$ corresponds to
\[
1+\sqrt{6}<\rho_0<1+\sqrt{4+\sqrt{6}}=3.5396.
\]
Indeed, representing~(\ref{eq:rhooneando}) as
\[
\rho_0^2-2\rho_0-(2+\rho_1)=0,
\]
and using the quadratic formula\index{quadratic formula} gives
\[
\rho_0=1+\sqrt{3+\rho_1},
\]
where we took the solution that gives $\rho_0>0$. The above concludes the first iteration. The result of this iteration is the element $c_3=1+\sqrt{4+\sqrt{6}}=3.5396$ of the sequence~(\ref{eq:perdoubseq}) and the scaling parameter $C_1(p)$ whose properties will be investigated in the following subsection.

We then proceed with the next iteration. In particular, executing the steps of the iteration method on the $s$-map yields the new map that we call the $t$-map.\index{t-map} This map has the form
\begin{equation}
t(k+1)=\left(4+2\rho_1-\rho_1^2\right)t(k)-\rho_1C_2 t^2(k).
\label{eq:newtmap}
\end{equation}
Introducing the iterated parameter,\index{iterated parameter} $\rho_2$, such that
\begin{equation}
\rho_2=\rho_1^2-2\rho_1-2,
\label{eq:recrelofu}
\end{equation}
and scaling the $t$-map using the transformation
\[
u(k)=\frac{\rho_1 C_2}{\rho_2}t(k)
\]
yields the $u$-map\index{u-map} of the form
\begin{equation}
u(k+1)=\left(2-\rho_2\right)u(k)-\rho_2 u^2(k).
\label{eq:newumap}
\end{equation}
The above map has an attracting 2-cycle for $3<\rho_2<1+\sqrt{6}$ which corresponds to an attracting 4-cycle of the $s$-map, which in turn corresponds to an attracting 8-cycle of the $z$-map. We calculate the range of the parameter $\rho$ corresponding to attracting 8-cycles of the $z$-map by working the recursive relation~(\ref{eq:recrelofu}) to get
\[
c_4=1+\sqrt{3+\rho_1}=1+\sqrt{3+1+\sqrt{3+\rho_2}}=1+\sqrt{4+\sqrt{4+\sqrt{6}}}=3.5573.
\]
In general, after $i+1$ iterations we obtain
\begin{equation}
\rho_{i+1}=\rho_i^2-2\rho_i-2,
\label{eq:afterkit}
\end{equation}
and work backwards to find the term $c_{i+3}$ of the sequence~(\ref{eq:perdoubseq}). We can also compute the elements of the sequence~(\ref{eq:perdoubseq}) using the recursive relation
\begin{equation}
\rho_{i}=\rho_{i+1}^2-2\rho_{i+1}-2,\quad i=1,2,\ldots.
\label{eq:forwrec}
\end{equation}
Solving the above, using the quadratic formula,\index{quadratic formula} yields
\[
\rho_{i+1}=1+\sqrt{3+\rho_i}.
\]
This means that the elements of the sequence~(\ref{eq:perdoubseq}) satisfy the difference equation\index{difference equation}
\begin{equation}
c_{i+1}=1+\sqrt{3+c_i},\quad c(1)=3.
\label{eq:formforci}
\end{equation}
The equation~(\ref{eq:formforci}) has two equilibrium points that are obtained by solving the algebraic equation
\[
c=1+\sqrt{3+c}.
\]
The equilibrium point\index{equilibrium point}
\[
c_{eq}=\frac{3+\sqrt{17}}{2}
\]
is attracting, and we have
\[
\lim_{i\to \infty}c_i=c_{\infty}=c_{eq}=3.5616.
\]
An analysis without truncation when computing second iterates of the logistic equation yields $c_{\infty}=3.5699$. Such an analysis can be found in the paper by Feigenbaum~\cite{feigenbaum83}. After the parameter $\rho$ reaches the limiting value $\rho=3.5699$, no stable solution patterns of the logistic equation emerge. We obtain a chaotic behavior.

\subsection{The Feigenbaum Numbers}

The sequence~(\ref{eq:perdoubseq}) and the scaling factors
\begin{equation}
\frac{\rho_i C_{i+1}}{\rho_{i+1}}
\label{eq:saclfc}
\end{equation}
have certain properties that are universal for all quadratic maps.\index{quadratic maps!universal properties of}\index{universal properties!of quadratic maps} Here, we perform an approximate analysis of these properties following Kahn~\cite[Chapter~16]{kahn90}.

We begin by rewriting equation~(\ref{eq:forwrec}) as
\begin{equation}
\left(\rho_{i+1}-1\right)^2=3+\rho_i.
\label{eq:feigdelo}
\end{equation}
Hence,
\begin{equation}
\rho_{i+1}-1=\sqrt{3+\rho_i},
\label{eq:rhoipo}
\end{equation}
and
\begin{equation}
\rho_{i}-1=\sqrt{3+\rho_{i-1}}.
\label{eq:rhhoi}
\end{equation}
Subtracting~(\ref{eq:rhhoi}) from~(\ref{eq:rhoipo}) gives
\begin{equation}
\rho_{i+1}-\rho_i=\sqrt{3+\rho_i}-\sqrt{3+\rho_{i-1}}.
\label{eq:subonefa}
\end{equation}
Summing~(\ref{eq:rhhoi}) and~(\ref{eq:rhoipo}) yields
\begin{equation}
\rho_{i+1}+\rho_i-2=\sqrt{3+\rho_i}+\sqrt{3+\rho_{i-1}}.
\label{eq:sumotoa}
\end{equation}
Multiplying~(\ref{eq:subonefa}) and~(\ref{eq:sumotoa}), we get
\begin{equation}
\left( \rho_{i+1}-\rho_i\right)\left(\rho_{i+1}+\rho_i - 2\right) = \rho_i - \rho_{i-1}.
\label{eq:mulfbys}
\end{equation}
We rewrite~(\ref{eq:mulfbys}) as
\begin{equation}
\rho_{i+1}+\rho_i - 2=\frac{\rho_i - \rho_{i-1}}{\rho_{i+1}-\rho_i}.
\label{eq:afterrw}
\end{equation}
Let
\begin{equation}
\fbox{$\displaystyle \delta_i=\frac{\rho_i - \rho_{i-1}}{\rho_{i+1}-\rho_i}$}
\label{eq:getfnu}
\end{equation}
The limiting value of $\delta_i$ as $i$ tends to infinity is called the {\em Feigenbaum number\/}~$\delta$.\index{Feigenbaum number $\delta$} To find its value, we use~(\ref{eq:getfnu}) to represent equation~(\ref{eq:afterrw}) as
\begin{equation}
\delta_i=\rho_{i+1}+\rho_i - 2.
\label{eq:newfwd}
\end{equation}
We obtain an approximate value for the Feigenbaum number $\delta$ because of the truncation procedure that we used to obtain~(\ref{eq:saclfc}). We have
\begin{equation}
\delta_{approx}=\lim_{i\to \infty}\left(\rho_{i+1}+\rho_i - 2\right) = 2c_{\infty}-2=5.1231
\label{eq:fndapp}
\end{equation}
The exact value of the Feigenbaum number $\delta$ is $\delta_{exact}=4.669202\ldots$.

The other parameter that characterizes our iteration method is the scaling factor given by~(\ref{eq:saclfc}). The limiting value of~(\ref{eq:saclfc}) as $i$ tends to infinity is called the {\em Feigenbaum number\/}~$\alpha$. Its approximate value is
\begin{eqnarray}
\alpha_{approx}&=&\lim_{i\to \infty}\frac{\rho_i C_{i+1}(p)}{\rho_{i+1}}\nonumber\\
&=&\frac{c_{\infty} C_{\infty}(p)}{c_{\infty}}\nonumber\\
&=& c_{\infty}^2-2c_{\infty}-3 - 3\sqrt{c_{\infty}^2-2c_{\infty}-3}\nonumber\\
&=&-2.24
\label{eq:apfal}
\end{eqnarray}
The exact value of the Feigenbaum number $\alpha$\index{Feigenbaum number $\alpha$} is $\alpha_{exact}=-2.50\ldots$.

\section{Fractals}

This section is devoted to modeling of irregular shapes. These non-Euclidean objects are called {\em fractals\/}.\index{fractal!non-Euclidean object} The term {\em fractal\/}\index{fractal!term coined by Mandelbrot} was coined by Mandelbrot~\cite{mandelbrot83} and it comes from the Latin adjective {\em fractus\/}, meaning ``irregular''.\index{fractus!irregular} The corresponding Latin verb is {\em frangere\/}, which means ``to break to create irregular fragments.''\index{frangere!to break to create irregular fragments} An example of a simple fractal\index{fractal!von Koch snowflake} is the so called the von~Koch snowflake\index{von Koch snowflake} curve that is shown in Figure~\ref{fig:vonkoch}. In this Figure, an iterative procedure for constructing the curve is also shown.
\begin{figure}
\centering
\includegraphics[scale=0.5]{vonkoch.eps}
%\vspace{4in}
%\centerbmp{5.33in}{4in}{c:/mybook/vonkoch.bmp}
%\hskip.2in\special{bmp:c:/mybook/vonkoch.bmp x=5.33in y=4in}
\caption{The first four stages in constructing the von~Koch snowflake curve.}
\label{fig:vonkoch}
\end{figure}
We start with dividing a line segment into thirds and then replacing the middle segment with two equal segments that form two sides of an equilateral triangle. The outcome of this stage is a shape consisting of four segments. In the next stage, each of the four segments is divided into thirds and the middle segments are replaced with two equal segments that form two sides of equilateral triangle. The process is repeated over and over again. The resulting curve is characterized by {\em self-similarity\/} or {\em scaling\/}.\index{self-similarity}\index{scaling}

One of the characteristics of a fractal is its dimension.\index{fractal!dimension of} A line segment is a one-dimensional Euclidean object\index{Euclidean object} and we can divide this line segment into $N$ identical parts. In such a case, the individual parts can be viewed as scaled down, by the factor $r=1/N$, versions of the original object. Similarly a square, a two-dimensional Euclidean object, can be divided into $N$ self-similar parts. These parts can be viewed as scaled down copies by the factor $r=1/\sqrt{N}$ of the original square. In general, a $D$-dimensional object can be divided into $N$  copies of itself. Each copy can be viewed as a scaled down version by the factor
\begin{equation}
r=\frac{1}{N^{1/D}}
\label{eq:dimfactor}
\end{equation}
of the parent object. Equation~(\ref{eq:dimfactor}) can be represented as
\begin{equation}
r^DN=1.
\label{eq:aftermulr}
\end{equation}
We now generalize the above arguments to the objects that are not necessarily Euclidean. Taking logarithms of both sides of~(\ref{eq:aftermulr}) yields the expression for the {\em fractal dimension\/}\index{fractal dimension} of the form
\begin{equation}
\fbox{$\displaystyle D=\frac{\log N}{\log\left( 1/r\right)}$}
\label{eq:fracdimm}
\end{equation}
We use~(\ref{eq:fracdimm}) to compute the dimension of the von~Koch curve\index{dimension!of the von Koch curve}\index{von Koch curve!dimension of}. Any segment of the curve is composed of four sub-segments, that is, $N=4$. Each sub-segment is a scaled down copy of its parent by the factor $r=1/3$. Hence, the dimension of the von~Koch curve is
\begin{equation}
D=\frac{\log N}{\log\left( 1/r\right)}=\frac{\ln 4}{\ln 3}=1.26
\label{eq:fdimvko}
\end{equation}
We can see from the above that the fractal dimension, unlike the dimension of a common Euclidean pattern, need not be an integer.

Another example of a fractal is the {\em Cantor set\/}\index{Cantor set}, sometimes called {\em Cantor's ternary set\/}.\index{Cantor's ternary set|see{Cantor set}} The Cantor set consists of the points of the unit interval $[0,1]$ that remain after ``middle third'' intervals have been successively removed. This set is illustrated in Figure~\ref{fig:cantors}. It is an example of so called {\em disconnected fractal\/}\index{disconnected fractal} or {\em fractal dust\/}.\index{fractal dust}
\begin{figure}
\centering
\includegraphics[scale=0.5]{cantors.eps}
%\vspace{3in}
%\centerbmp{3.4in}{3in}{c:/mybook/cantors.bmp}
%\hskip1.25in\special{bmp:c:/mybook/cantors.bmp x=3.4in y=3in}
\caption{The first four iterations in constructing Cantor's ternary set.}
\label{fig:cantors}
\end{figure}
For Cantor's set, we have $N=2$ and $1/r=3$. Hence, the fractal dimension of the Cantor set\index{fractal dimension!of the Cantor set}\index{Cantor set!fractal dimension of} is
\begin{equation}
D=\frac{\log N}{\log\left( 1/r\right)}=\frac{\ln 2}{\ln 3}=0.631
\label{eq:cansdim}
\end{equation}




\subsection{The Mandelbrot Set}

The Mandelbrot set,\index{Mandelbrot set} denoted ${\cal M}$, is the set of all complex $c\in \mb{\Bb C}$ for which the sequence
\begin{equation}
0, c, c^2+c, (c^2+c)^2+c, \ldots
\label{eq:quacoseq}
\end{equation}
does not tend to $\infty$ as $n$ tends to $\infty$. Thus, the Mandelbrot set is generated by the map in the complex plane,
\begin{equation}
z \mapsto z^2 +c,\quad c\in\mb{\Bb C},
\label{eq:qudmapm}
\end{equation}
or equivalently,
\begin{equation}
z(k+1)=z^2(k)+c,\quad c\in\mb{\Bb C}.
\label{eq:eqquacom}
\end{equation}
Let
\[
z(k)=x(k)+jy(k)\quad\mb{and}\quad c=a+jb.
\]
Then,
\begin{eqnarray*}
z(k+1)&=& x(k+1)+jy(k+1)\\
&=& x^2(k)-y^2(k)+a+j\left( 2x(k)y(k)+b \right)\\
&=& x^2(k)+2jx(k)y(k)+(jy(k))^2+a+jb\\
&=& z^2(k)+c.
%\label{eq:eqbetzab}
\end{eqnarray*}
Thus, we can express the first-order quadratic complex dynamical system~(\ref{eq:eqquacom}) as a second-order real dynamical system of the form
\begin{equation}
\left.\begin{array}{lll}
x(k+1)&=&x^2(k)-y^2(k)+a\\
y(k+1)&=&2x(k)y(k)+b
\end{array}\right\}
\label{eq:reldyesys}
\end{equation}
The representations~(\ref{eq:eqquacom}) and (\ref{eq:reldyesys}) are indeed equivalent.

In terms of the dynamical system~(\ref{eq:reldyesys}), the Mandelbrot set is the set of parameter values, $a$ and $b$, for which the solutions of~(\ref{eq:reldyesys}) starting at $x(0)=y(0)=0$ remain bounded. The Mandelbrot set is shown in Figure~\ref{fig:mandset}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{mandset.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/mandset.bmp}
%\hskip1in\special{bmp:c:/mybook/mandset.bmp x=3.97in y=3.5in}
\caption{The Mandelbrot set for $x(0)=y(0)=0$.}
\label{fig:mandset}
\end{figure}

\subsection{Julia Sets}

A filled Julia set\index{filled Julia set} corresponding to the given values of the parameters $a$ and $b$ of the dynamical system~(\ref{eq:reldyesys}) is the set of initial conditions $(x(0),y(0))$ resulting in bounded system trajectories. For each value of $c=a+jb$ there is a different Julia set. In order to construct a Julia set, pick first a point $(a,b)$, say in the Mandelbrot set. Then, pick an initial condition $(x(0), y(0))$ and compute the resulting trajectory of~(\ref{eq:reldyesys}) starting from the chosen initial condition $(x(0), y(0))$. Plot only the initial conditions $(x(0), y(0))$ that give rise to bounded trajectories. The set of plotted points is the filled Julia set corresponding to the particular parameter values $(a, b)$. The boundary of the filled Julia set is the Julia set.\index{Julia set!the boundary of the filled Julia set} Two examples of filled Julia sets are shown in Figures~\ref{fig:juliar} and~\ref{fig:juliab}. The filled Julia set shown in Figure~\ref{fig:juliar} is called by Keen~(\cite[page~61]{devaney89}), ``Douady's Rabbit''.\index{Julia set!Douady's Rabbit} It is observed by Peitgen in~\cite[page~199]{peitgen88} that in a sense the Mandelbrot set can be considered as a one-page-dictionary of all Julia sets. This means that the Mandelbrot set stores infinitely many different Julia sets.\index{Mandelbrot set!one-page-dictionary of all Julia sets}
\begin{figure}
\centering
\includegraphics[scale=0.5]{juliar.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/juliar.bmp}
%\hskip.7in\special{bmp:c:/mybook/juliar.bmp x=4.67in y=3.5in}
\caption{Filled Julia set for $c=-0.12256117+j0.74486177$.}
\label{fig:juliar}
\end{figure}
\begin{figure}
\centering
\includegraphics[scale=0.5]{juliab.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/juliab.bmp}
%\hskip.7in\special{bmp:c:/mybook/juliab.bmp x=4.67in y=3.5in}
\caption{Filled Julia set for $c=-0.75$.}
\label{fig:juliab}
\end{figure}


\subsection{Iterated Function Systems}

An iterated function system~(IFS)\index{iterated function system (IFS)}\index{IFS|see{iterated function system}} can be used to model natural objects such as clouds, smoke, leaves, or ferns. An IFS uses affine transformations that we discuss next. A two-dimensional affine transformation\index{affine transformation} $\b w:\mb{\Bb R}^2 \to \mb{\Bb R}^2$ from $\mb{\Bb R}^2$ into itself is defined by
\begin{equation}
\b w (\b x)=\b A\b x+\b b,
\label{eq:aftradeifs}
\end{equation}
where $\b A\in\mb{\Bb R}^{2\times 2}$ and $\b b\in \mb{\Bb R}^2$. The affine transformation~(\ref{eq:aftradeifs}) satisfies
\begin{equation}
\left\|\b w(\b x)-\b w(\b y) \right\|\le L\|\b x-\b y\|,
\label{eq:licoofafftr}
\end{equation}
where the Lipschitz constant\index{Lipschitz constant} $L$ is
\[
L=\|\b A\|.
\]
The affine transformation $\b w$ is contractive\index{affine transformation!contractive}\index{contractive!affine transformation} if $L<1$. It is a symmetry if\index{affine transformation!symmetry}\index{symmetry!affine transformation}
\[
\left\|\b w(\b x)-\b w(\b y) \right\|=\|\b x-\b y\|.
\]
The affine transformation $\b w$ is expansive if $L>1$.\index{affine transformation!expansive}\index{expansive!affine transformation}

We now define an IFS. A two-dimensional IFS consists of a set of $N$ affine transformations denoted
\[
\{\b w_1,\b w_2,\ldots,\b w_N  \},
\]
each transformation mapping $\mb{\Bb R}^2$ into itself, together with a set of probabilities
\[
\{p_1,p_2,\ldots,p_N \},
\]
where each $p_i\ge 0$ and
\[
p_1+p_2+\cdots +p_N=1.
\]
An IFS is the set
\[
\{\b w_n,p_n:n=1,2,\ldots,N\}.
\]
Let $L_n$, $n=1,2,\ldots,N$, denote the Lipschitz constant for $\b w_n$. Then, we say that the IFS obeys the average contractivity condition if\index{average contractivity!condition of IFS}
\begin{equation}
L_1^{p_1}L_2^{p_2} \cdots L_N^{p_N}<1.
\label{eq:aveconcoifs}
\end{equation}
Interesting mathematical properties of the IFS that satisfy the condition~(\ref{eq:aveconcoifs}) are reported by Barnsley~\cite{barnsley88}.

An IFS operates as follows. A set of $N$ affine transformations is chosen and an initial point is picked. Then, the affine transformations are iteratively applied over and over again in a random order with probabilities $p_n$, for example, proportional to the determinants of $\b A_n$. The result is a sequence of points from $\mb{\Bb R}^2$,
\[
\b x(0), \b x(1),\ldots
\]
As an example, suppose that $N=4$, where
\begin{eqnarray*}
\b A_1&=&\left[\begin{array}{cc}
0.02 &0\\
 0 &0.21
\end{array}\right]\quad \b b_1=\left[\begin{array}{c}
-0.15\\
 -0.8
\end{array}\right]\\
\b A_2&=&\left[\begin{array}{cc}
0.8496 &0.0255\\
 -0.0255 &0.8496
\end{array}\right]\quad \b b_2=\left[\begin{array}{c}
0\\
0.15
\end{array}\right]\\
\b A_3&=&\left[\begin{array}{cc}
-0.1554& 0.235\\
0.1958 &0.1865
\end{array}\right]\quad \b b_3=\left[\begin{array}{c}
0.05\\
-0.55
\end{array}\right]\\
\b A_4&=&\left[\begin{array}{cc}
0.1554 &-0.235\\
0.1958& 0.1865
\end{array}\right]\quad \b b_4=\left[\begin{array}{c}
-0.35\\
-0.4
\end{array}\right].
\end{eqnarray*}
The above data comes from Bunde and Havlin~\cite{bunde94}. The elements of the set of probabilities
\[
\{p_1,p_2,p_3,p_4\}
\]
were constructed as
\[
p_n=\frac{\left|\det(\b A_n)\right|}{\sum_{i=1}^N \left|\det(\b A_i)\right|}\quad n=1,2,3,4.
\]
The initial point was chosen to be $\b x(0)=\left[\begin{array}{cc}
-0.1 & -0.2
\end{array}\right]^{\top}$. Then, the roulette wheel method\index{roulette wheel method} was used to select the affine transformation to operate on $\b x(0)$. The roulette wheel method operates as follows. Each affine transformation is assigned a slice of a wheel, where the slice is proportional to the probability $p_n$ of the particular transformation $\b w_n$. A random number from the interval $[0,1]$ is generated, which corresponds to spinning the wheel. The transformation under the wheel's marker is selected. Suppose, $\b w_k$ was selected. Then $\b w_k$ operates on $\b x(0)$  and the result is
\[
\b x(1)=\b A_k\b x(0)+\b b_k .
\]
The roulette wheel is spun again to select the affine transformation to act upon $\b x(1)$ in order to obtain $\b x(2)$, etc.

The result of $30,000$ iterations is shown in Figure~\ref{fig:fern}. The object is called the Barnsley fern\index{Barnsley fern}\index{IFS!Barnsley fern} because it was M.~Barnsley~\cite{barnsley88} who introduced the method.
\begin{figure}
\centering
\includegraphics[scale=0.5]{fern.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/fern.bmp}
%\hskip.7in\special{bmp:c:/mybook/fern.bmp x=4.54in y=4in}
\caption{The Barnsley fern.}
\label{fig:fern}
\end{figure}

\section{Lyapunov Exponents}
\label{sec:lyapunexpon}
Trajectories of a chaotic\index{chaotic trajectories!sensitive dependence on initial conditions} dynamical system are characterized by a sensitive dependence upon initial conditions, which means that the distance between the adjacent trajectories grows exponentially with time as illustrated in Figure~\ref{fig:lya_exp}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{lya_exp.eps}
%\vspace{4in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/lya_exp.bmp}
%\hskip.7in\special{bmp:c:/mybook/lya_exp.bmp x=4.67in y=3.5in}
\caption{The change in distance between adjacent chaotic trajectories that is used to compute the largest Lyapunov exponent.}
\label{fig:lya_exp}
\end{figure}
For a discrete dynamical system, if $x^{*}$ is a fixed point of a one-dimensional map $f$ and if $\left\arrowvert f^{\prime}\left( x^*\right)\right\arrowvert=a>1$, then by Theorem~\ref{thm:atractrepel}, trajectories starting near $x^{*}$ will separate from $x^{*}$ at a multiplicative rate\index{multiplicative rate!of separation} of approximately $a$ per iteration until they move significantly far away from $x^{*}$. Let $x_1$ be a point that constitutes the $k$-cycle. Then, for a $k$-cycle,\index{k-cycle} the trajectory starting from a point close to $x_1$, will separate from $x_1$ at a rate approximately equal to the product of the derivatives of the $k$ points of the $k$-cycle. Let this product be denoted $A$ and suppose that $A>1$. Then, it takes $k$ iterations to separate a trajectory starting close to $x_1$ by a distance $A$ from $x_1$. We can say that in the case of a $k$-cycle the average rate of separation is $A^{1/k}$ per iteration. We now are ready to introduce the term {\em Lyapunov number\/}\index{Lyapunov number} to describe the average multiplicative rate of separation\index{average multiplicative rate!of separation} of points close to $x_1$. For a $k$-cycle, the Lyapunov number $L(x_1)$ is
\begin{eqnarray}
L(x_1)&=& \left( \left\arrowvert \left( f^k\right)^{\prime}(x_1) \right\arrowvert\right)^{1/k}\nonumber\\
&=& \left( \left\arrowvert f^{\prime}\left( x_1\right) \right\arrowvert \left\arrowvert f^{\prime}\left( x_2\right) \right\arrowvert \cdots \left\arrowvert f^{\prime}\left( x_k\right) \right\arrowvert\right)^{1/k}
\label{eq:onedlyapnu}
\end{eqnarray}
The concept of Lyapunov number is now extended to non-periodic trajectories.\index{non-periodic trajectories}
\begin{defn}
Let $f$ be a differentiable map on the real line $\mb{\Bb R}$ and let $\{x(0), x(1), \ldots \}$ be a trajectory starting from $x(0)$. Then, the Lyapunov number $L(x(0))$ of the trajectory starting from $x(0)$ is
\[
L(x(0))=\lim_{N\to \infty}\left( \left\arrowvert f^{\prime}(x(0))\right\arrowvert \cdots  \left\arrowvert f^{\prime}(x(N))\right\arrowvert\right)^{1/N}
\]
if the limit exists.
\end{defn}\markdefn

The Lyapunov exponent, $\mu(x(0))$, is then defined as the natural logarithm of the Lyapunov number, that is,
\begin{equation}
\mu(x(0)) = \ln L(x(0))= \lim_{N\to \infty}\frac{\ln\left\arrowvert f^{\prime}(x(0))\right\arrowvert +\cdots +\ln \left\arrowvert f^{\prime}(x(N))\right\arrowvert}{N}
\label{eq:lysexpodi}
\end{equation}
if the limit exists.

We now look at the Lyapunov exponent\index{Lyapunov exponent} from another perspective. Let $d_0$ denote the distance between the two initial states of a continuous dynamical system. Then, for a chaotic motion, at a later time,
\begin{equation}
d(t)=e^{\mu t}d_0.
\label{eq:divochatr}
\end{equation}
The measure of the trajectories' divergence\index{trajectories' divergence} is obtained by averaging the exponential growth\index{averaging!the exponential growth} at many points along a trajectory. This is because the chaotic trajectories' divergence can only be measured locally because if the trajectories are bounded, $d(t)$ in~(\ref{eq:divochatr}) cannot go to infinity. To define the measure of the divergence of chaotic trajectories, first a reference trajectory is obtained. Next a point on an adjacent trajectory\index{adjacent trajectory} is selected and $d(t)/d_0$ is measured. This is illustrated in Figure~\ref{fig:lya_exp}. Then, a new $d_0(t)$ is selected on new adjacent trajectory. Using the collected data, the Lyapunov exponent is computed as
\begin{equation}
\mu(x(0))=\lim_{N\to \infty} \frac{1}{t_N-t_0}\sum_{k=1}^N \ln \frac{d(t_k)}{d_0(t_{k-1})}.
\label{larlyaex}
\end{equation}
There are $n$ Lyapunov exponents for an $n$-dimensional nonlinear dynamical system, just as there are $n$ eigenvalues that govern behavior of an $n$-dimensional linear system. To define these Lyapunov exponents for an $n$-dimensional nonlinear dynamical system model, first choose an $n$-dimensional sphere centered at a point on a reference trajectory. Let the radius of this sphere be $\delta_0$. We can view this sphere as an $n$-dimensional ellipsoid with $n$ semi-axes of equal length $\delta_i(0)=\delta_0$. Next, at a later time $t$ construct an $n$-dimensional ellipsoid with the property that all the trajectories emanating from the previously chosen sphere, pass through this ellipsoid. Let the $n$ semi-axes of the ellipsoid\index{ellipsoid!semi-axes of}\index{semi-axes!of ellipsoid} be denoted $\delta_i(t)$. Then, we can compute the Lyapunov exponents as
\begin{equation}
\fbox{$ \displaystyle \mu_i(\b x(0))=\lim_{t \to \infty}\left(\lim_{\delta_i(0)\to 0}\left\{ \frac{1}{t}\ln \frac{\delta_i(t)}{\delta_i(0)}\right\} \right), \quad i=1,2,\ldots,n  $}
\label{eq:orideolyaex}
\end{equation}
if the limits exist. It is common to order the Lyapunov exponents so that
\[
\mu_1\ge \mu_2\ge \ldots\ge \mu_n.
\]
Thus, the Lyapunov exponents are numbers that are used to measure the exponential attraction or separation\index{exponential attraction!measure of}\index{exponential separation!measure of} in time of two neighboring trajectories of a dynamical system.
\begin{example}\rm
This example was adapted from Moon~\cite[page~193]{moon87}. Suppose that we are given a one-dimensional map
\begin{equation}
x(k+1)=f(x(k)).
\label{eq:odimatile}
\end{equation}
We demonstrate again that in the regions where $f$ is differentiable, the Lyapunov exponent\index{Lyapunov exponent} is determined by $|f^{\prime}|$. Indeed, let $x(k)$ and $x(k)+\varepsilon$ be two initial conditions. Then, we have
\[
d_0(k-1)=\varepsilon\quad\mb{and}\quad d(k) = \left|f(x(k)+\varepsilon)-f(x(k))\right|\approx \left|\left. f^{\prime}\right\arrowvert_{x=x(k)}\right|\varepsilon.
\]
Using~(\ref{larlyaex}), we obtain
\begin{eqnarray*}
\mu(x(0))&=&\lim_{N\to \infty} \frac{1}{t_N-t_0}\sum_{k=1}^N \ln \frac{d(t_k)}{d_0(t_{k-1})}\\
&=&\lim_{N \to \infty}\frac{1}{N}\sum_{k=1}^N \ln\left|f^{\prime}\left(x(k)\right) \right|,
\end{eqnarray*}
which is the same as~(\ref{eq:lysexpodi})
\end{example}\markex

We use the notion of Lyapunov exponent to formally define chaotic motion.\index{chaotic motion!defintion of}\index{chaotic motion!Lyapunov exponent} The other term that we use in this definition is an asymptotically periodic trajectory that we define next. This definition is adapted from Alligood et al.~\cite[page~108]{alligood97}.
\begin{defn}
\label{def:asymper}
A trajectory $\{\b x(0), \ldots, \b x(N),\ldots \}$ of a smooth map $\b f$ is called {\em asymptotically periodic\/} if it converges to a periodic trajectory as $N\to \infty$; that is, there exists a periodic trajectory $\{ \b y_1, \b y_2, \ldots, \b y_k, \b y_1, \b y_2, \ldots\}$ such that
\[
\lim_{N\to \infty} \left\Arrowvert \b x(N)-\b y(N)\right\Arrowvert =0.
\]
\end{defn}\markdefn

Using the definition of Lyapunov exponent and the above definition, we can define a chaotic trajectory. This definition is adapted from Alligood et al.~\cite[page~110]{alligood97}.
\begin{defn}
Let $f$ be a map of the real line $\mb{\Bb R}$, and let $\{ x(0), x(1),\ldots \}$ be a bounded trajectory of the map $f$. This trajectory is chaotic if
\begin{enumerate}
\item the trajectory is not asymptotically periodic,
\item the Lyapunov exponent $\mu(x(0))$ is greater then zero.
\end{enumerate}
\end{defn}\markdefn

For a map on $\mb{\Bb R}^n$, each trajectory has $n$ Lyapunov numbers. The Lyapunov numbers, in this case, measure the rates of separation\index{Lyapunov numbers!measure the rates of separation} from the current point of the trajectory along $n$ orthogonal directions. The first direction is the one where the separation between the nearby points is the greatest. The second direction is the one in which the separation is the greatest from all directions perpendicular to the first direction. The third direction will be the one where the separation is the greatest among directions perpendicular to the first two directions, and so on. We illustrate the above in two dimensions in Figure~\ref{fig:twodilyanum}, where a unit disk is centered at $\b x(0)$. After $N$ iterates of the two-dimensional map $\b f$, the disk evolves into an ellipse. The per-iteration changes of the ellipse axes are the Lyapunov numbers of $\b x(0)$. Thus, the Lyapunov numbers measure the contraction or expansion near the trajectory of $\b x(0)$ during the first $N$ iterations. The natural logarithm of each Lyapunov number is a Lyapunov exponent.\index{Lyapunov exponent!natural logarithm of Lyapunov number}
\begin{figure}
\centering
\includegraphics[scale=0.5]{tdillye.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/tdillye.bmp}
%\hskip1.25in\special{bmp:c:/mybook/tdillye.bmp x=4in y=3in}
\caption{The unit disk is mapped into an ellipse by a two-dimensional map.}
\label{fig:twodilyanum}
\end{figure}

Let $D\b f(\b x(0))$ be the first derivative of $\b f$ evaluated at $\b x(0)$. Let $J_N=D\b f^N(\b x(0))$ be the first derivative of the $N$-th iterate of $\b f$. We denote by $B$ the unit sphere centered at $\b x(0)$. With the above notation, we are ready to give the definition of Lyapunov exponents for a map on $\mb{\Bb R}^n$. This definition is adapted from Alligood et al.~\cite[page~195]{alligood97}.
\begin{defn}
Let $\b f$ be a smooth map on $\mb{\Bb R}^n$ and let $r_k^N$ be the length of the $k$-th longest orthogonal axis of the ellipsoid $J_NB$ for a trajectory with the initial state $\b x(0)$. The $k$-th Lyapunov number is
\[
L_k(\b x(0))=\lim_{N\to \infty}\left( r_k^N\right)^{1/N}
\]
if the limit exists. The $k$-th Lyapunov exponent of $\b x(0)$ is
\[
\mu_k(\b x(0))=\ln L_k(\b x(0)).
\]
\end{defn}\markdefn

Using the above definition, we can extend the definition of chaotic trajectories for maps\index{chaotic trajectories!for maps} on $\mb{\Bb R}^n$. This definition is adapted from Alligood et al.~\cite[page~196]{alligood97}
\begin{defn}
Let $\b f$ be a map on $\mb{\Bb R}^n$, $n\ge 1$, and let $\{\b x(0),\b x(1), \b x(2),\ldots \}$ be a bounded trajectory of the map $\b f$. The trajectory is chaotic if
\begin{enumerate}
\item it is not asymptotically periodic,
\item no Lyapunov exponent is equal to zero,
\item $\mu_1(\b x(0))>0$.
\end{enumerate}
\end{defn}\markdefn

For continuous systems, the above definition of a chaotic trajectory has to be modified---see Alligood at al.~\cite[pages~385--386]{alligood97}.

To numerically compute Lyapunov exponents, a reference trajectory of a given nonlinear system is constructed and the ellipsoid semi-axes lengths $r_k$-s are found using linearized models. We now analyze this process in some detail. For this, suppose that we are given an $n$-th order continuous dynamical system model
\begin{equation}
\dot{\b x}=\b f(\b x),
\label{eq:tiinnosymole}
\end{equation}
where $\b x(t) \in \mb{\Bb R}^n$ is the system's state at time $t$, and $\b f:\mb{\Bb R}^n \to \mb{\Bb R}^n$ is called the vector field.\index{vector field} The solution to~(\ref{eq:tiinnosymole}) subject to initial condition $\b x(t_0)=\b x_0$  is called a trajectory and is denoted $\b{\phi}_t(\b x_0)$. For every $t$, the mapping
\[
\b{\phi}_t:\mb{\Bb R}^n \to \mb{\Bb R}^n,
\]
maps the state space $\mb{\Bb R}^n$ into itself. The mapping $\b{\phi}_t$ is called the {\em flow\/} of the system~(\ref{eq:tiinnosymole}).\index{flow!of the system} For example, if~(\ref{eq:tiinnosymole}) is linear, that is, $\b f(\b x)=\b A\b x$, where $\b A\in \mb{\Bb R}^{n\times n}$, then
\[
\b{\phi}_t(\b x_0)=e^{\b A(t-t_0)}\b x_0.
\]
%In the case of a time-varying $n$-th order model
%\begin{equation}
%\dot{\b x}=\b f(t,\b x),
%\label{eq:tvarndinm}
%\end{equation}
%the vector field $\b f$ depends on time. The solution %to~(\ref{eq:tvarndinm}) subject to the initial condition $\b x(t_0)=\b x_0$ %is denoted $\b{\phi}_t(t_0,\b x_0)$.

Consider now an equilibrium state $\b x_{e}$ of~(\ref{eq:tiinnosymole}). The local behavior of~(\ref{eq:tiinnosymole}) in a neighborhood of $\b x_{e}$ is determined by the linearized model of~(\ref{eq:tiinnosymole}) about the equilibrium state $\b x_{e}$. The linearized model\index{linearized model} is
\begin{equation}
\frac{d}{dt}\Delta \b x=D\b f(\b x_{e})\Delta \b x,
\label{eq:limooftiinle}
\end{equation}
where $\Delta \b x=\b x-\b x_{e}$. We can also linearize the system~(\ref{eq:tiinnosymole}) about its trajectory $\b{\phi}_t(\b x_0)$. The linearized model about the trajectory $\b{\phi}_t(\b x_0)$ is
\begin{equation}
\frac{d}{dt}\Delta \b x=D\b f(\b{\phi}_t(\b x_0))\Delta \b x,\quad \Delta \b x(t_0)=\Delta \b x_0.
\label{eq:liabatra}
\end{equation}
The solution to~(\ref{eq:liabatra}) is
\[
\Delta \b x(t)=\b{\Phi}_t(\b x_0)\Delta \b x_0.
\]
We call the $n\times n$ matrix $\b{\Phi}_t(\b x_0)$ the {\em state transition matrix\/}\index{state transition matrix} of the system~(\ref{eq:liabatra})---see Section~\ref{sec:sol_LTV_uncon} for more information on the state transition matrix properties. The transition matrix satisfies the equation
\[
\frac{d}{dt}{\b{\Phi}}=D\b f(\b{\phi}_t(\b x_0))\b{\Phi}, \quad \b{\Phi}_0=\b I_n.
\]
Let
\[
\left\{m_i(t)\right\}_{i=1}^n
\]
be the eigenvalues of the transition matrix $\b{\Phi}_t(\b x_0)$. Then, the Lyapunov exponents of $\b x_0$ are
\begin{equation}
\fbox{$\displaystyle \mu_i\left(\b x_0\right)=\lim_{t \to \infty}\frac{1}{t}\ln\left| m_i(t)\right|,\quad i=1,2,\ldots,n $}
\label{eq:deolyaex}
\end{equation}
if the limits exist. The above is an alternative to~(\ref{eq:orideolyaex}) representation of the Lyapunov exponents.
\begin{example}\rm
The Lyapunov exponents of an equilibrium state $\b x_{e}$ of the system~(\ref{eq:tiinnosymole}), with $t_0=0$, can be found by first linearizing~(\ref{eq:tiinnosymole}) about the equilibrium state to obtain a linearized model of the form~(\ref{eq:limooftiinle}). The transition matrix of~(\ref{eq:limooftiinle}) is
\[
\b{\Phi}_t(\b x_{e})=e^{D\b f(\b x_{e})t}.
\]
The $i$-th eigenvalue of the transition matrix $\b{\Phi}_t(\b x_{e})$ is $m_i(t)=e^{\lambda_i t}$, where $\lambda_i$ is the $i$-th eigenvalue of the matrix $D\b f(\b x_{e})$. Therefore, the Lyapunov exponents of $\b x_{e}$ are
\begin{eqnarray}
\mu_i\left(\b x_e\right) &=& \lim_{t \to \infty}\frac{1}{t}\ln\left| m_i(t)\right|\nonumber\\
&=& \lim_{t \to \infty}\frac{1}{t}\ln\left| e^{\Re \lambda_i t}e^{j\Im \lambda_i t}\right|\nonumber\\
&=& \lim_{t \to \infty}\frac{1}{t}\ln\left(\left| e^{\Re \lambda_i t}\right|\left|\cos\left( \Im \lambda_i t\right) + j \sin \left(\Im\lambda_i t\right) \right|\right)\nonumber\\
&=& \lim_{t\to \infty }\frac{1}{t}t\,\Re \lambda_i\nonumber\\
&=& \Re \lambda_i.
\label{eq:lyaforextinc}
\end{eqnarray}
Thus, the Lyapunov exponents in this case are the real parts of the eigenvalues\index{Lyapunov exponents!real parts of the eigenvalues of linearized system} of the linearized system about the equilibrium state. Note that if $\b x_0$ is in the basin of attraction of $\b x_{e}$, that is, if
\[
\lim_{t \to \infty}\b{\phi}_t(\b x_0)=\b x_{e},
\]
then the Lyapunov exponents of $\b x_0$ are the same as the Lyapunov exponents of $\b x_{e}$.
\end{example}\markex
\begin{example}\rm
Suppose that we are given a discrete-time dynamical system model,
\begin{equation}
\b x(k+1)=\b f(\b x(k)),\quad \b x(0)=\b x_0.
\label{eq:distnomole}
\end{equation}
Let $\b x_{e}$ be the equilibrium state of the model~(\ref{eq:distnomole}), that is,
\[
\b x_{e}=\b f\left( \b x_{e}\right).
\]
The linearized model is
\[
\Delta \b x(k+1)=D\b f(\b x_{e})\Delta x(k).
\]
The state transition matrix of the linearized model is
\begin{equation}
\b{\Phi}_k(\b x_{e})=\left(D\b f\left(\b x_{e} \right) \right)^k
\label{eq:sttrfodis}
\end{equation}
If $\lambda_i$ is the $i$-th eigenvalue of the matrix $D\b f\left(\b x_{e} \right)$, then $\lambda_i^k$ is the $i$-th eigenvalue of the transition matrix~(\ref{eq:sttrfodis}). The Lyapunov exponents of $\b x_{e}$, in this case, are
\begin{eqnarray}
\mu_i(\b x_e) &=& \lim_{k \to \infty}\frac{1}{k}\ln\left| m_i(k)\right|\nonumber\\
&=& \lim_{k \to \infty}\frac{1}{k}\ln\left| \lambda_i\right|^k\nonumber\\
&=&\ln \left|\lambda_i \right|.
\label{eq:lyaefoditsy}
\end{eqnarray}
\end{example}\markex

In summary, Lyapunov exponents can be used to characterize a chaotic motion.\index{Lyapunov exponents!characterize a chaotic motion} Specifically, a necessary condition for a chaotic motion to occur in a dynamical system with bounded trajectories is that the largest Lyapunov exponent be positive. An efficient numerical method for evaluating the largest Lyapunov exponent can be found, for example, in the paper by Parker and Chua~\cite{parker87}, see also their book~\cite{parker89} on practical numerical algorithms for chaotic systems.


In Figure~\ref{fig:le_loma}, we show a plot of the Lyapunov exponent versus the parameter $\rho$ for the logistic equation~(\ref{eq:1d_log_work}).
\begin{figure}
\centering
\includegraphics[scale=0.5]{le_loma.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/le_loma.bmp}
%\hskip1in\special{bmp:c:/mybook/le_loma.bmp x=4.54in y=4in}
\caption{Lyapunov exponent versus parameter $\rho$ for the logistic equation~(\ref{eq:1d_log_work}).}
\label{fig:le_loma}
\end{figure}

\section{Discretization Chaos}

In this section, we consider examples of chaotic systems that are discrete-time analogs of continuous-time systems, where the latter systems do not exhibit any chaotic behavior. We begin by considering a method of obtaining discrete models from continuous models.\index{discretization chaos}\index{discrete models!from continuous models} Suppose that we are given a scalar differential equation
\begin{equation}
\dot{x}=ax,\quad x(t_0)=x_0.
\label{eq:simfodeco}
\end{equation}
The solution to the above differential equation is
\begin{equation}
x(t)=e^{a(t-t_0)}x(t_0).
\label{eq:sotofodeq}
\end{equation}
Let $\Delta t=t-t_0$. Then, we can express~(\ref{eq:sotofodeq}) as
\begin{equation}
x(1)=e^{a \Delta t}x(0),
\label{eq:soforwdelt}
\end{equation}
where $x(1)=x(t_0+\Delta t)$ and $x(0)=x(t_0)$. More generally,
\begin{equation}
x(k+1)=e^{a \Delta t}x(k),
\label{eq:dismoofsied}
\end{equation}
where $x(k+1)=x((k+1)\Delta t)$ and $x(k)=x(k\Delta t)$. The discrete model~(\ref{eq:dismoofsied}) of the continuous differential equation~(\ref{eq:simfodeco}) is exact in the sense that the values of~(\ref{eq:sotofodeq}) and~(\ref{eq:dismoofsied}) are the same at the sampling times. We generalize the above method to a class of nonlinear differential equations. Consider a scalar differential equation
\begin{equation}
\dot{x}(t)=x(t) h(x(t)).
\label{eq:nonldieq}
\end{equation}
Dividing both sides of the above by $x$ and then integrating yields
\begin{equation}
\int\frac{dx}{x}=\int h(x(t)) dt.
\label{eq:dibxaint}
\end{equation}
Suppose now that we hold $h(x(t))$ constant on the interval $\left[\begin{array}{ll}
k\Delta t, &(k+1)\Delta t
\end{array}\right)$. Then, evaluating~(\ref{eq:dibxaint}) over the above interval gives
\[
\ln\left( \frac{x((k+1)\Delta t)}{x(k\Delta t)}\right)=h(x(k\Delta t))\Delta t.
\]
Using the the notation introduced after~(\ref{eq:dismoofsied}), we represent the above as
\begin{equation}
\fbox{$ \displaystyle x(k+1)=e^{h(x(k)) \Delta t}x(k) $}
\label{eq:evafmkas}
\end{equation}
\begin{example}\rm
Consider the system of differential equations
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1 &=& \gamma x_1\left( 1-x_1-0.5 x_2\right)\\
\dot{x}_2 &=& 2\gamma x_2\left( 1-0.5x_1-x_2\right)
\end{array}\right\}
\label{eq:granaathcos}
\end{equation}
where $\gamma=2.5$. The above system models a pair of competing species.\index{competing species!population dynamics} For a discussion of modeling of competing species population dynamics, the reader is referred to Subsection~\ref{subsec:popugromo}. The system~(\ref{eq:granaathcos}) was analyzed by Grantham and Athalye~\cite{grantham90}. A trajectory of the system~(\ref{eq:granaathcos}) for $x_1(0)=x_2(0)=0.5$ is shown in Figure~\ref{fig:co_specp}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{co_specp.eps}
%\vspace{4in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/co_specp.bmp}
%\hskip1in\special{bmp:c:/mybook/co_specp.bmp x=4.54in y=4in}
\caption{A trajectory of the system~(\ref{eq:granaathcos}).}
\label{fig:co_specp}
\end{figure}
We discretize~(\ref{eq:granaathcos}) using the formula~(\ref{eq:evafmkas}). For $\Delta t=1$, we obtain the following discrete-time analog\index{discrete-time analog}
\begin{equation}
\left.\begin{array}{ccl}
x_1(k+1) &=& x_1(k)\exp \left(\gamma\left( 1-x_1(k)-0.5 x_2(k)\right)\right)\\
x_2(k+1) &=& x_2(k)\exp \left(2\gamma\left( 1-0.5x_1(k)-x_2(k)\right)\right)
\end{array}\right\}
\label{eq:disanofcosex}
\end{equation}
The trajectory of the discrete model~(\ref{eq:disanofcosex}) starting from the same initial condition as its continuous-time analog is shown in Figure~\ref{fig:dico_spp}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{dico_spp.eps}
%\vspace{4in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/dico_spp.bmp}
%\hskip1in\special{bmp:c:/mybook/dico_spp.bmp x=4.54in y=4in}
\caption{A trajectory of the system~(\ref{eq:granaathcos}).}
\label{fig:dico_spp}
\end{figure}
The continuous-time model trajectory asymptotically converges to the equilibrium state $\left[\begin{array}{cc}
2/3 & 2/3
\end{array}\right]^{\top}$, while the discrete-time trajectory wanders around chaotically.
\label{ex:twokiocosp}
\end{example}\markex

\section{Controlling Chaotic Systems}

A trajectory of a chaotic system behaves in a random-like fashion. The set on which the chaotic motion is taking place attracts trajectories starting from a set of initial conditions. Roughly speaking, this attracting set is called a {\em chaotic attractor\/}.\index{chaotic attractor} It is not possible to predict when a trajectory will pass through a given point in the chaotic attractor. However, any given point in the chaotic attractor will be visited by the chaotic system's trajectory in a finite time. This property of a chaotic system is referred to as {\em ergodicity\/}.\index{ergodicity!property} Using the ergodicity property one can propose a control algorithm whose task it is to bring the chaotic system's trajectory to a prespecified state in the system's state space. The idea behind the chaotic control algorithm,\index{chaotic control algorithm} that we discuss next, is to set up ``fishing nets'' in the chaotic attractor in a neighborhood of an operating point to ``catch'' the system trajectory. Once the trajectory is ``caught'', a control law is used to bring the trajectory to the desired operating point. The controller is constructed using linear quadratic regulator~(LQR)\index{linear quadratic regulator (LQR)} technique. The chaotic algorithm\index{chaotic algorithm|see{chaotic control algorithm}} is illustrated with two examples; the H\'{e}non map and the bouncing ball system.
\subsection {Ingredients of Chaotic Control Algorithm}
Chaotic motion is characterized by extreme sensitivity to initial conditions. This might suggest that it would be impossible, for example, to stabilize a chaotic system with respect to a desired operating point, because any perturbation resulting from applying a control action would lead to a chaotic motion's exponential grow in time. However, it turns out that precisely this property of a chaotic motion can be used to convert a chaotic attractor to a periodic motion. A chaotic attractor,\index{chaotic attractor} as observed by Ott, Grebogi, and Yorke~\cite{ott90} typically has imbedded within it an infinite number of unstable periodic orbits. These periodic orbits are dense in the attractor in the sense that they pass through any neighborhood, however small, of any point in the attractor. Ott et al.~\cite{ott90} noticed that these properties could be advantageous in achieving control objectives. They proposed a method, now called the OGY method, which involves converting, using small time-dependent perturbations of a system parameter, a chaotic attractor into a periodic orbit. The idea behind this method is to exploit the ergodicity property of a chaotic system. Specifically, all of the unstable periodic orbits of the chaotic attractor will eventually be visited by the system trajectory. Once the trajectory gets to a neighborhood of the desired unstable periodic orbit, a small perturbation is applied to stabilize the orbit of interest.

The chaotic control algorithm, that we discuss next, can be viewed as a variant of the OGY method.\index{Ott, Grebogi, and Yorke (OGY) method}\index{OGY method|see{Ott, Grebogi, and Yorke (OGY) method}} This algorithm was presented by Vincent~\cite{vincent97} and is characterized by two ingredients. They are: a chaotic attractor and a controllable target. These components of the chaotic algorithm are illustrated in Figure~\ref{fig:cha_algo}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{cha_algo.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/cha_algo.bmp}
%\hskip1.25in\special{bmp:c:/mybook/cha_algo.bmp x=4in y=3in}
\caption{A pictorial representation of the chaotic algorithm.}
\label{fig:cha_algo}
\end{figure}
We assume that either the system under consideration is chaotic or that it can be made chaotic using an open-loop control. Next, we select a fixed point of interest and assume that a feedback law can be constructed such that the closed-loop system is non-chaotic and the chosen fixed point is asymptotically stable. The region of attraction\index{region of attraction} of the fixed point should have a non-empty intersection with the chaotic attractor. A controllable target\index{controllable target!chaotic control algorithm} can be any subset of the region of attraction of the fixed point that has a non-empty intersection with the chaotic attractor.\index{chaotic attractor} Note that the fixed point does not have to be located in the chaotic attractor. If we start the system from any initial state in the chaotic attractor, the resulting chaotic trajectory will at some time arrive at the intersection region of the chaotic attractor and the controllable target. The chaotic algorithm keeps track when the trajectory reaches this intersection region. When this happens, the open-loop control, that was used to create chaotic motion, is turned off and the feedback control is turned on to drive the system towards the fixed point. The region of attraction intersecting the chaotic attractor plays the role of a ``fishing net'' for the chaotic trajectory.


\subsection{Chaotic Control Algorithm}

We consider two classes of single-input dynamical systems. The first class is modeled by nonlinear difference equation,
\begin{equation}
\b x(k+1)=\b f(\b x(k),u(k)),
\label{eq:chacondysy}
\end{equation}
and the other class is modeled by nonlinear differential equation,
\begin{equation}
\dot{\b x}(t)=\b f(\b x(t),u(t)),
\label{eq:chaconcondysy}
\end{equation}
The control $u$, in both cases, will be bounded. The control $u$ will be either a specified function of time or a specified function of the system state. In the first case, we have an open-loop control, in the latter case we have a closed-loop control. We assume that for a given system there exists an open-loop control\index{open-loop control!chaotic control algorithm} such that the system has a chaotic attractor. In addition, we assume that for a specified constant control $u=u_{e}$ there is a corresponding equilibrium state located in a vicinity of the chaotic attractor. The equilibrium state satisfies the equation
\[
\b x_{e}=\b f(\b x_{e},u_{e})
\]
for the discrete-time system, and it satisfies the relation
\[
\b 0=\b f(\b x_{e},u_{e})
\]
for the continuous-time system.

We now briefly describe the way we obtain a controllable target for the system of interest. We first linearize the system about the selected equilibrium state. The linearized discrete-time system\index{linearized!discrete-time system} has the form,
\[
\Delta \b x(k+1)=\b A\Delta \b x(k)+\b b\Delta u(k),
\]
where $\Delta \b x(k)=\b x(k)-\b x_{e}\quad\mb{and}\quad \Delta u(k)=u(k)-u_{e}$. The linearized continuous-time system has the form\index{linearized!discrete-time system}
\[
\frac{d}{dt}{\Delta \b x}(t)=\b A\Delta \b x(t)+\b b\Delta u(t),
\]
where $\Delta \b x(t)=\b x(t)-\b x_{e}\quad\mb{and}\quad \Delta u(t)=u(t)-u_{e}$. Then, the LQR method,\index{LQR method} described in Chapter~\ref{chpt:lqr}, is used to construct a state variable feedback control law, $\Delta u(\b \Delta \b x)=-\b k\Delta \b x$, that guarantees the origin of the linearized system,
\begin{equation}
\Delta \b x(k+1)=(\b A-\b b\b k)\Delta \b x(k),
\label{eq:optiliclosy}
\end{equation}
or
\begin{equation}
\frac{d}{dt}\Delta \b x(t)=(\b A-\b b\b k)\Delta \b x(t),
\label{eq:opticontclsy}
\end{equation}
to be asymptotically stable. This implies that for the corresponding nonlinear system driven by the state feedback controller, the equilibrium state about which the system was linearized, is asymptotically stable. We use a Lyapunov function obtained from the linearized system to estimate the region of asymptotic stability~(RAS) of the closed-loop nonlinear system. The Lyapunov function\index{Lyapunov function!chaotic control algorithm}\index{chaotic control algorithm!Lyapunov function} has the form, $V(\b x)=\b x^{\top}\b P_l\b x$, where $\b P_l$ is obtained by solving the corresponding Lyapunov equation for the optimal closed-loop linear system given by~(\ref{eq:optiliclosy}) or by~(\ref{eq:opticontclsy}). This Lyapunov function is then used to determine a controllable target for the nonlinear system driven by a bounded controller\index{bounded controller} of the form,
\begin{equation}
u=\mb{sat}(\b k\b x)+u_{e},
\label{eq:lqrfonnlsat}
\end{equation}
where
\[
\mb{sat}(v)=\left\{\begin{array}{lcl}
u_{\max}-u_{e}&\mb{if}&v>u_{\max}-u_{e}\\
v&\mb{if}& u_{\min}-u_{e}\le v\le u_{\max}-u_{e}\\
u_{\min}-u_{e}&\mb{if}&v<u_{\min}-u_{e},
\end{array}\right.
\]
and the bounds $u_{\min}$ and $u_{\max}$ are specified by a designer. Using the above control law, we then seek the largest level curve $V(\b x)=V_{\max}$ such that for the points $\{\b x:V(\b x)\le V_{\max} \}$, for the discrete-time system,
\begin{eqnarray*}
\Delta V(k) &=& V(\b x(k+1))-V(\b x(k))\\
&=& (\b x(k+1)-\b x_{e})^{\top}\b P_l(\b x(k+1)-\b x_{e})-(\b x(k)-\b x_e)^{\top}\b P_l(\b x(k)-\b x_e)\\
&=&\left( \b f(\b x(k),u(k))-\b x_{e}\right)^{\top}\b P_l\left( \b f(\b x(k),u(k))-\b x_{e}\right)-\left( \b x(k)-\b x_{e}\right)^{\top}\b P_l\left(\b x(k)-\b x_{e} \right)\\
&<& 0,
\end{eqnarray*}
or for the continuous-time system,
\begin{eqnarray*}
\dot{V}(\b x) &=& 2\left( \b x-\b x_{e}\right)^{\top} \b P_l\dot{\b x}\\
&=& 2\left( \b x-\b x_{e}\right)^{\top} \b P_l\b f(\b x,u)\\
&<& 0.
\end{eqnarray*}
We are now ready to describe the chaotic algorithm.\index{chaotic algorithm} When the given nonlinear system trajectory is in the region $\{\b x:V(\b x)>V_{\max} \}$, an open-loop control is applied to keep the system in the chaotic attractor until the system reaches a point in the set $\{\b x:V(\b x)\le V_{\max} \}$, where the open-loop control is turned off and the feedback saturating control law is applied to force the system trajectory to asymptotically approach the equilibrium state. More formally, we state the chaotic algorithm in the form of the following algorithm:\\
\hspace*{1in}IF $V(\b x)>V_{\max}$\\
\hspace*{1in}THEN $u=$ open-loop control to keep the system in the chaotic attractor\\
\hspace*{1in}ELSE $u=\mb{sat}(\b k \b x)+u_{e}$\\
\hspace*{1in}END\\
We illustrate the above algorithm when applied to the H\'{e}non map.\index{H\'{e}non map!chaotic algorithm}
\begin{example}\label{ex:henon_cont}\rm
The dynamical system known as the H\'{e}non map is described by the difference equations,
\begin{equation}
\left[\begin{array}{c}
x_1(k+1)\\
x_2(k+1)
\end{array}\right]=\left[\begin{array}{l}
1+x_2(k)-1.4 x_1^2(k)\\
0.3x_1(k)
\end{array}\right].
\label{eq:henonattr}
\end{equation}
It was named after the French astronomer M.~H\'{e}non who discovered chaotic behavior of the above system in 1974. The orbit of $\b x(0) = \left[\begin{array}{cc}
1 & 0
\end{array}\right]^{\top}$ of the H\'{e}non map is shown in Figure~\ref{fig:henon_ma}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{henon_ma.eps}
%\vspace{3.5in}
%\centerline{\psfig{figure=c:/mybook/henon_ma.eps,height=3in}}
%\centerbmp{4in}{3in}{c:/mybook/henon_ma.bmp}
%\hskip1in\special{bmp:c:/ee675/henon_ma.bmp x=4.67in y=3.5in}
\caption{An orbit of the H\'{e}non map.}
\label{fig:henon_ma}
\end{figure}
We modify the H\'{e}non map by adding the control input\index{chaotic algorithm!H\'{e}non map} into the first equation so that the map becomes,
\begin{equation}
\left[\begin{array}{c}
x_1(k+1)\\
x_2(k+1)
\end{array}\right]=\left[\begin{array}{l}
-1.4 x_1^2(k)+x_2(k)+1+u\\
0.3x_1(k)
\end{array}\right].
\label{eq:henonattrcon}
\end{equation}
There are two fixed states corresponding to $u=u_{e}=0$. They are obtained by solving two algebraic equations,
\[
\left\{\begin{array}{ccc}
x_{1, eq} &=& -1.4 x^2_{1, e}+x_{2, e}+1\\
x_{2, eq} &=& 0.3 x_{1, e}
\end{array}\right.
\]
Solving the above gives
\[
\left[\begin{array}{cc}
0.6314 & 0.1894
\end{array}\right]^{\top}\quad\mb{and}\quad \left[\begin{array}{cc}
-1.1314 & -0.3394
\end{array}\right]^{\top}.
\]
We select to work with the first of the above fixed states. We linearize the system~(\ref{eq:henonattrcon}) about the first fixed state to obtain
\[
\b A=\left.\left[\begin{array}{cc}
-2.8 x_1 & 1\\
0.3 & 0
\end{array}\right]\right\arrowvert_{x_1=0.6314}=\left[\begin{array}{cc}
-1.7678 & 1\\
0.3 & 0
\end{array}\right],\quad \b b=\left[\begin{array}{c}
1\\
0
\end{array}\right].
\]
We then use the discrete LQR design\index{discrete!LQR design} method to find the state feedback. The weight matrices of the associated quadratic performance index are chosen to be $\b Q=\b I_2$ and $r=1$. Solving the discrete algebraic Riccati equation\index{discrete!algebraic Riccati equation}\index{algebraic Riccati equation!discrete} yields
\[
\b P=\left[\begin{array}{cc}
3.9351 & -1.5008\\
-1.5008 & 1.7974
\end{array}\right]\quad\mb{and}\quad \b k=\left[\begin{array}{cc}
-1.5008 & 0.7974
\end{array}\right].
\]
Hence,
\[
\b A_c=\b A-\b b\b k=\left[\begin{array}{cc}
-0.2670 & 0.2026\\
0.3000 & 0
\end{array}\right].
\]
Solving the discrete Lyapunov equation,\index{discrete!Lyapunov equation}
\[
\b A_c^{\top}\b P_l\b A_c-\b P_l=-\b I_2
\]
gives
\[
\b P_l=\left[\begin{array}{cc}
1.1902 & -0.0686\\
-0.0686 & 1.0489
\end{array}\right].
\]
The next step is to calculate $V_{\max}$ corresponding to the largest level curve such that the Lyapunov first difference $\Delta V <0$ on this level curve for the nonlinear model. We obtain level curves by trial and error as follows. We first evaluate
\begin{eqnarray}
V &=& \left( \b x-\b x_{e}\right)^{\top}\b P_l\left( \b x-\b x_{e}\right)\nonumber\\
&=& \left[\begin{array}{cc}
x_1 -0.6314 & x_2-0.1894
\end{array}\right]\left[\begin{array}{cc}
1.1902 & -0.0686\\
-0.0686 & 1.0489
\end{array}\right]\left[\begin{array}{c}
x_1 -0.6314\\
x_2-0.1894
\end{array}\right]\nonumber\\
&=& 1.1902(x_1-0.6314)^2-0.1372(x_1-0.6314)(x_2-0.1894)\nonumber\\
&{}&\quad +1.0489(x_2-0.1894)^2.
\label{eq:expfovzer}
\end{eqnarray}
We next choose an initial state $\b x(0)$ such that $V(\b x(0))=V_0$, where $V_0$ is a constant selected by the designer. Let $x_2(0)=0.1894$, which is the value of the second component of the equilibrium state of interest. Substituting this value into~(\ref{eq:expfovzer}) allows us to compute the first component of $\b x(0)$ such that $V(\b x(0))=V_0$. We have
\[
\b x(0)=\left[\begin{array}{cc}
\sqrt{0.8402 V_0}+0.6314 & 0.1894
\end{array}\right]^{\top}.
\]
Then, a level curve\index{level curve} is constructed by solving the differential equation,
\[
\nabla ^{\top} V \dot{\b x}=0,
\]
subject to the initial condition $\b x(0)$. To satisfy the above equation, we can select
\begin{equation}
\left\{\begin{array}{ccc}
\dot{x}_1 &=& \frac{\partial V}{\partial x_2}\\
\dot{x}_2 &=& -\frac{\partial V}{\partial x_1}
\end{array}\right.
\label{eq:solflcurv}
\end{equation}
After constructing a level curve, we evaluate $\Delta V$ at these points and check if $\Delta V < 0$ at every point on the level curve. For the control law
\begin{equation}
u=0.3\,\mb{sat}(\b k\b x/0.3),
\label{eq:satcoforhenma}
\end{equation}
applied to the nonlinear system~(\ref{eq:henonattrcon}), we found $V_{\max}=0.0109$. Plots of $x_1$ and $x_2$ versus time for the closed-loop system~(\ref{eq:henonattrcon}), (\ref{eq:satcoforhenma}) are shown in Figures~\ref{fig:henon_x1} and~\ref{fig:henon_x2}, respectively, for the initial condition $\b x(0)=\left[\begin{array}{cc}
-0.5 & 0
\end{array}\right]^{\top}$.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=henon_x1.eps,height=3in}}
%\centering\psfig{file=c:/ee675/henon_x1.ps,height=3.5in}
%\hskip1in\special{ps:c:/ee675/henon_x1.ps y=3.5in}
\caption{A plot of $x_1$ versus time for the controlled H\'{e}non map of Example~\ref{ex:henon_cont}.}
\label{fig:henon_x1}
\end{figure}
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=henon_x2.eps,height=3in}}
%\centering\psfig{file=c:/ee675/henon_x2.ps,height=3.5in}
%\hskip1in\special{ps:c:/ee675/henon_x2.ps y=3.5in}
\caption{A plot of $x_2$ versus time for the controlled H\'{e}non map of Example~\ref{ex:henon_cont}.}
\label{fig:henon_x2}
\end{figure}
The orbit of $\b x(0) = \left[\begin{array}{cc}
-0.5 & 0
\end{array}\right]^{\top}$ of the controlled H\'{e}non map\index{controlled H\'{e}non map}\index{H\'{e}non map!controlled} is shown in Figure~\ref{fig:henx1x2}, where the controllable target is also depicted as the set of points inside the ellipse $V(\b x-\b x_e)= 0.0109$.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=henx1x2.eps,height=3in}}
%\centering\psfig{file=c:/ee675/henx1x2.ps,height=3.5in}
%\hskip1in\special{ps:c:/ee675/henx1x2.ps y=3.5in}
\caption{An orbit of the controlled H\'{e}non map of Example~\ref{ex:henon_cont}, where the initial condition is denoted by ``$\circ$'' and labelled ``IC''. The subsequent iterations prior to hitting the controllable target are denoted by ``$\star$''. The iterations in the controllable target are denoted by ``$\cdot$''.}
\label{fig:henx1x2}
\end{figure}
\label{ex:henonex}
\end{example}\markex

\begin{example}\rm {\bf Controlling the bouncing ball on a vibrating plate}\\
The bouncing ball on a vibrating plate\index{bouncing ball!on a vibrating plate} is an another example of a chaotic system.\index{chaotic system!bouncing ball on a vibrating plate} In this example, we use the chaotic algorithm to design a control law for the vibrating plate, as shown in Figure~\ref{fig:bch_ball}, so that the ball, starting at rest on the plate can be bounced up to and then maintained bouncing at a prescribed height. The amplitude of the vibrating plate will be kept fixed. The control input is the frequency of the plate vibrations.\index{control input!frequency of the plate vibrations}
\begin{figure}
\centering
\includegraphics[scale=0.5]{bch_ball.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/bch_ball.bmp}
%\hskip1.25in\special{bmp:c:/mybook/bch_ball.bmp x=4in y=3in}
\caption{Bouncing ball on a vibrating plate.}
\label{fig:bch_ball}
\end{figure}
To proceed, we need a mathematical model of the ball-plate system, which we derive next.
\begin{center}
{\bf Modeling a Bouncing Ball on a Vibrating Plate}
\end{center}
Let the displacement of the plate in the vertical direction be
\[
y(t)=A\sin(\omega t),
\]
where $A$ is the amplitude of the plate displacement, $\omega$ is the angular frequency of vibrations and $t$ is time. Then, the velocity of the plate is
\begin{equation}
W(t)=\frac{d y(t)}{dt}=A\omega\cos(\omega t).
\label{eq:veloftheplat}
\end{equation}
We denote by $U$ the velocity of the ball just before impact and let $V$ denote the velocity of the ball just after impact. In our modeling process, we use the principle of conservation of momentum,\index{principle!of conservation of momentum}\index{conservation of momentum!principle of} which states that momentum of the system before and after collisions is constant, that is,
\[
\sum(\mb{masses} \times \mb{velocities})_{before}=\sum (\mb{masses} \times \mb{velocities})_{after}.
\]
Applying the conservation of momentum principle to our system consisting of the bouncing ball and the vibrating plate yields
\begin{equation}
m_{ball}U+m_{plate}W_{before}=m_{ball}V+m_{plate}W_{after}.
\label{eq:resofpconmo}
\end{equation}
To proceed, we define the so called {\em coefficient of restitution}\index{coefficient of restitution}
\begin{equation}
e=-\frac{V-W_{after}}{U-W_{before}}.
\label{eq:decoeore}
\end{equation}
Performing simple manipulations on~(\ref{eq:decoeore}), we obtain
\begin{equation}
W_{after}=Ue-W_{before}e+V.
\label{eq:eqreofrpom}
\end{equation}
Substituting the above into~(\ref{eq:resofpconmo}) and rearranging gives
\begin{equation}
m_{ball}U+m_{plate}W_{before}=\left( m_{ball}+m_{plate}\right)V+e m_{plate}U-e m_{plate}W_{before}.
\label{eq:mbaumpwb}
\end{equation}
Dividing both sides of the above by $(m_{ball}+m_{plate})$ and rearranging, we get
\begin{equation}
V=\frac{m_{ball}-e m_{plate}}{m_{ball}+m_{plate}}U+\frac{m_{plate}+em_{plate}}{m_{ball}+m_{plate}}W_{before}.
\label{eq:afdivbsbmbpmp}
\end{equation}
Let
\[
M=\frac{m_{ball}}{m_{plate}}.
\]
Substituting the above into~(\ref{eq:afdivbsbmbpmp}) gives
\begin{equation}
V=\frac{M-e}{M+1}U+\frac{e+1}{M+1}W_{before}.
\label{eq:uscamdeab}
\end{equation}
The ball is acted upon by a constant acceleration $g$ and therefore $V(t)=gt$. Hence, from one bounce to the next,
\begin{equation}
t(k+1)=t(k)+\frac{2V(k)}{g}.
\label{eq:frontooth}
\end{equation}
We assume that the amplitude, $A$, of the plate vibrations\index{amplitude!of the plate vibrations} is negligible compared with the height of the bounce. Referring to Figure~\ref{fig:fromkto}, we obtain
\begin{equation}
U(k+1)=-V(k).
\label{eq:ulatisvno}
\end{equation}
\begin{figure}
\centering
\includegraphics[scale=0.5]{fromkto.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/fromkto.bmp}
%\hskip1.25in\special{bmp:c:/mybook/fromkto.bmp x=4in y=3in}
\caption{Velocities of the bouncing ball from one bounce to the next.}
\label{fig:fromkto}
\end{figure}
Let
\[
a_1=\frac{e+1}{M+1}\quad\mb{and}\quad a_2=\frac{M-e}{M+1}.
\]
Taking into account the above notation and~(\ref{eq:uscamdeab}), we obtain a relation between velocities of the bouncing ball\index{bouncing ball} and the vibrating plate\index{vibrating plate} at the bounce at a time $t(k+1)$,
\begin{equation}
V(k+1)=a_2 U(k+1)+a_1 W_{before}(k+1).
\label{eq:uvandwbefore}
\end{equation}
We use~(\ref{eq:ulatisvno}) to represent (\ref{eq:uvandwbefore}) as
\begin{equation}
V(k+1)=-a_2V(k)+a_1W(k+1).
\label{eq:usaoaat}
\end{equation}
Furthermore, by~(\ref{eq:veloftheplat}) and~(\ref{eq:frontooth}),
\[
W(k+1)=A\omega \cos(\omega t(k+1))=A\omega \cos\left(\omega \left( t(k)+\frac{2V(k)}{g}\right) \right).
\]
Hence, we can represent~(\ref{eq:usaoaat}) as
\begin{equation}
V(k+1)=-a_2V(k)+a_1A\omega \cos\left(\omega \left( t(k)+\frac{2V(k)}{g}\right) \right).
\label{eq:vkploww}
\end{equation}
Let
\begin{equation}
\phi(k)=\omega t(k)\quad\mb{and}\quad \psi(k)=2\omega V(k)/g,
\label{eq:defphiapsi}
\end{equation}
where $\phi(k)$ is the current phase angle of the plate and $\psi(k)$ is the phase angle change between the current bounce of the ball and the next.\index{phase angle change!between current bounce of the ball and the next} Then, we can write
\begin{equation}
\fbox{$\displaystyle\begin{array}{ccl}
\phi(k+1) &=& \phi(k)+\psi(k)\\
\psi(k+1) &=& -a_2 \psi(k)+\hat{a}_1\omega^2\cos(\phi(k)+\psi(k)),
\end{array} $}
\label{eq:bamausino}
\end{equation}
where $\hat{a}_1=2Aa_1/g$. The above dynamical system is known as the {\em ball map}.\index{ball map} In our simulations, we use the same values of the ball map parameters as those used by Vincent~\cite{vincent97}, which are given in Table~\ref{tab:bounce_ball}.
\begin{table}
\caption{Simulation parameters of the ball map.}
\begin{center}
\begin{tabular}{cc}\hline
$M$ & 1/26\\
$e$ & 0.8\\
$A$ & 0.013 m\\ \hline
\end{tabular}
\end{center}
\label{tab:bounce_ball}
\end{table}
Using these parameters, we compute $a_1= 1.73333$, $\hat{a}_1= 0.004594 \;\mb{sec}^2$, and $a_2= - 0.733333$.

There are many possible periodic equilibrium solutions to the ball map~(\ref{eq:bamausino})\index{periodic equilibrium solutions!to the ball map}\index{ball map!periodic equilibrium solutions to}. For instance, the ball can bounce to a fixed height at every $n$ cycles of the plate, or it can bounce to $m$ different heights at $nm$ cycles of the plate. Here, we will only be concerned with the case when $m=1$. There are bands of frequencies of the vibrating plate that will not support periodic motion for $m=1$. Examples of such frequency bands, as listed by Vincent~\cite{vincent97}, are:
\begin{eqnarray*}
32.3966&<\omega<& 33.0782\\
36.4575&<\omega<& 38.1954\\
40.5278&<\omega<& 42.7038\\
44.4249&<\omega<& 46.7796
\end{eqnarray*}
A frequency of the plate vibration from one of the above ranges will either result in a chaotic motion of the ball or a chaotic transient to a periodic motion with $m>1$. In our example, we use $\omega = 45$~rad/sec to generate chaotic motion.\index{chaotic motion generation}

The ball map~(\ref{eq:bamausino}) may produce negative values for $\psi$ corresponding to negative bounce heights,\index{bounce heights!of the ball map} which is physically impossible. To overcome this problem, whenever the map produces a negative value for $\psi$, this value is multiplied by a random negative number between zero and one in order to simulate the fact that the ball cannot go through the plate. In Figure~\ref{fig:bball_ma}, an orbit of the modified ball map for the initial conditions $\phi(0)=0$ and $\psi(0)=9.1743$~rads and $\omega=45$~rad/sec is shown. Note that the values of $\psi$ are divided by $2\pi$ before plotting.
\begin{figure}
%\centering
%\includegraphics[scale=0.5]{fromkto.eps}
%\vspace{3.5in}
\centerline{\psfig{figure=bball_ma.eps,height=3in}}
%\centerbmp{4in}{3in}{c:/ee675/bball_ma.bmp}
%\hskip1in\special{bmp:c:/ee675/bball_ma.bmp x=4.67in y=3.5in}
\caption{An orbit of the ball map.}
\label{fig:bball_ma}
\end{figure}
The chaotic attractor\index{chaotic attractor!of the ball map} is located between the upper and lower bounds as depicted in the figure. The upper bound is given by
\[
\psi_{upper}=2\hat{a}_1\omega^2+\hat{a}_1\omega^2\cos(\phi(k)+\psi(k)),
\]
while the lower bound is
\[
\psi_{lower}=\hat{a}_1\omega^2 \cos(\phi(k)+\psi(k)),
\]
When $\psi_{lower}$ goes negative it is set equal to zero.

\begin{center}
{\bf Controller Construction for a Bouncing Ball on a Vibrating Plate}
\end{center}

To proceed with the controller design, we first find fixed points of the ball map.\index{fixed points!of the ball map}\index{ball map!fixed points of} They are found by solving the equations,
\begin{equation}
\left\{\begin{array}{ccl}
\bar{\phi}&=&\bar{\phi}+\bar{\psi}\\
\bar{\psi} &=& -a_2 \bar{\psi}+\hat{a}_1 \omega^2 \cos(\bar{\phi}+\bar{\psi}).
\end{array}\right.
\label{eq:eqfofifipts}
\end{equation}
The fixed points are obtained by evaluating the right-hand side of the first above equation modulo $2n\pi$, where $n$ is the number of cycles of the plate between bounces. We obtain,
\begin{equation}
\left\{\begin{array}{ccl}
\bar{\psi}&=&2n \pi\\
\cos(\bar{\phi}) &=& \frac{2n\pi(1+a_2)}{\hat{a}_1 \omega^2}.
\end{array}\right.
\label{eq:soleqfofifipts}
\end{equation}
Linearizing the ball map\index{linearizing!the ball map} about~(\ref{eq:soleqfofifipts}) gives
\begin{equation}
\left[\begin{array}{c}
x_1(k+1)\\
x_2(k+1)
\end{array}\right]=\left[\begin{array}{cc}
1&1\\
a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}
x_1(k)\\
x_2(k)
\end{array}\right],
\label{eq:linbalma}
\end{equation}
where
\begin{eqnarray*}
x_1(k) &=& \phi(k)-\bar{\phi}\\
x_2(k) &=& \psi(k)-\bar{\psi}\\
a_{21} &=& -\hat{a}_1 \omega^2 \sin(\bar{\phi})\\
a_{22} &=& -a_2+a_{21}.
\end{eqnarray*}
We next construct a control law for selecting frequencies of the vibrating plate in such a way that the steady-state constant height of the bouncing ball is reached. Vincent~\cite{vincent97} showed that for $n=1$ the range of frequencies for which there exists an asymptotically stable equilibrium\index{asymptotically stable equilibrium} is $\left[19.0977~\mb{rad/sec}, 28.9505~\mb{rad/sec}\right]$. The corresponding range of the bounce heights is $\left[0.0578~\mb{m}, 0.1327~\mb{m}\right]$. Thus, for example the height, $h$, of 0.1~m is in the above range and from the energy conservation law,\index{energy conservation law}
\[
\frac{1}{2}m_{ball}V^2=m_{ball}gh,
\]
we obtain
\[
V=1.3994~\mb{m/sec}^2.
\]
Using the above,~(\ref{eq:defphiapsi}), and~(\ref{eq:soleqfofifipts}), we get
\[
\bar{\psi}=2\pi=\frac{2\omega V}{g}\quad\mb{and}\quad \cos(\bar{\phi})=\frac{2\pi(1+a_2)}{\hat{a}_1\omega^2},
\]
from which we compute,
\[
\omega=22~\mb{rad/sec}\quad\mb{and}\quad\bar{\phi}=41.1^{\circ}=0.7173~\mb{rad}.
\]
The linearized ball map about the above equilibrium state\index{linearized ball map!about equilibrium state} is
\[
\b x(k+1)=\left[\begin{array}{cc}
1 & 1\\
-1.4617 & -0.7283
\end{array}\right]\b x(k),
\]
which has its poles in the open unit circle.\index{open unit circle} Solving the discrete Lyapunov equation,\index{discrete!Lyapunov equation} $\b A^{\top}\b P\b A-\b P=-\b I_2$, gives
\[
\b P=\left[\begin{array}{cc}
8.6414 & 4.6978\\
4.6978 & 5.9598
\end{array}\right].
\]
Employing the same method as in the previous example, we obtain $V_{\max}=5.006$. The controllable target is
\[
\left\{\left[\begin{array}{c}
\phi\\
\psi
\end{array}\right]:V(\phi-\bar{\phi},\psi-\bar{\psi})=\left[\begin{array}{cc}
\phi -0.7173 &\psi -2\pi
\end{array}\right] \b P\left[\begin{array}{c}
\phi-0.7173\\
\psi-2\pi
\end{array}\right] \le V_{max}\right\}.
\]
The controllable target is depicted in Figure~\ref{fig:bb_cntr}. We use $\omega=45$~rad/sec when the system trajectory is outside of the controllable target,\index{controllable target!of the ball map} and apply $\omega=22$~rad/sec when the system trajectory hits the controllable target, that is, we use the algorithm:\\
\hspace*{1in}IF $V(\phi-\bar{\phi},\psi-\bar{\psi})>V_{\max}$\\
\hspace*{1in}THEN $\omega=45$~rad/sec\\
\hspace*{1in}ELSE $\omega=22$~rad/sec\\
\hspace*{1in}END\\
An orbit of $\phi(0)=0$, $\psi(0)=9.1743$ of the controlled modified ball map\index{controlled modified ball map}\index{modified ball map}\index{ball map!modified} using the above algorithm is shown in Figure~\ref{fig:bb_cntr}.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=bb_cntr.eps,height=3in}}
%\centerbmp{4.67in}{3.5in}{c:/ee675/bb_cntr.bmp}
%\hskip1in\special{bmp:c:/ee675/bb_cntr.bmp x=4.67in y=3.5in}
\caption{An orbit of the controlled modified ball map.}
\label{fig:bb_cntr}
\end{figure}
A plot of $\phi(k)$ versus time is shown in Figure~\ref{fig:bb_maphi}.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=bb_maphi.eps,height=3in}}
%\centerbmp{4.67in}{3.5in}{c:/ee675/bb_maphi.bmp}
%\hskip1in\special{bmp:c:/ee675/bb_maphi.bmp x=4.67in y=3.5in}
\caption{A plot of $\phi(k)$ versus time of the controlled modified ball map.}
\label{fig:bb_maphi}
\end{figure}
A plot of $\psi(k)$ versus time is shown in Figure~\ref{fig:bb_mapsi}.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=bb_mapsi.eps,height=3in}}
%\centerbmp{4.67in}{3.5in}{c:/ee675/bb_mapsi.bmp}
%\hskip1in\special{bmp:c:/ee675/bb_mapsi.bmp x=4.67in y=3.5in}
\caption{A plot of $\psi(k)$ versus time of the controlled modified ball map.}
\label{fig:bb_mapsi}
\end{figure}
\end{example}\markex

\section{An Application to Secure Communication}

The purpose of secure communication\index{secure communication} is to limit the availability of the transmitted information to a third party even if the communication channel\index{communication channel} is public or has been compromised. With an ever increasing reliance on wireless and internet communication\index{wireless communication}\index{internet communication} of sensitive information, including data related to individuals, electronic commerce and other financial transactions, secure communication is of utmost importance.


One approach to secure transfer of signals is to hide them in broad spectrum noise-like signals.\index{broad spectrum!noise-like signals} Chaotic signals are characterized by spread spectrum,\index{chaotic signals!characterized by spread spectrum} where the power of a chaotic signal is spread across a broad range of frequencies.\index{chaotic signal!power spread across a broad range of frequencies} Thus, chaotic signals are particularly suitable for secure communication because the transmitted signal masked by a chaotic signal may appear to a third party as a noise. The concept of synchronized systems,\index{synchronized systems} that we define next, plays important role when using chaotic systems in secure communication.
\begin{defn}
Two dynamical systems,
\[
\dot{\b x}=\b f(\b x)\quad\mb{and}\quad \dot{\tilde{\b x}}=\b g(\tilde{\b x}), \quad \b x, \tilde{\b x}\in \mb{\Bb R}^n,
\]
are said to be synchronized if
\[
\b e(t)=(\b x(t) - \tilde{\b x}(t)) \to \b 0\quad \mb{as}\quad t\to \infty .
\]
\end{defn}\markdefn

Cuomo, Oppenheim, and Strogatz~\cite{cuomo93} used the concept of synchronized chaotic systems (SCS's) to secure communications with the Lorenz circuit implemented in both the transmitter and receiver. Specifically, a chaotic masking signal was added at the transmitter to the message. Then, at the receiver, the masking was regenerated and subtracted from the received signal.

Synchronization is widely used in electronics, for example, in television, radar, and electronic watches~\cite{tang83}. The general idea behind synchronization is the adaptive generation of a reference signal that is used to estimate the frequency and phase of the received signal. The problem of designing a dynamical system whose behavior mimics that of a given chaotic system is referred to as chaotic synchronization.\index{chaotic synchronization!definition of} This is a challenging design problem because chaotic systems are characterized by a sensitive dependence on initial conditions. Thus, one may be tempted to conclude that chaotic synchronization is not feasible because even infinitesimal errors in initial conditions or system parameters will result in the exponential divergence of the systems' trajectories.


Although chaotic systems would seem to defy synchronization, Pecora and Carroll~\cite{pecora90,carroll91} demonstrated that chaotic system synchronization is possible. Morg\"{u}l and Solak~\cite{morgul96} proposed a systematic approach to chaotic system synchronization using observers.\index{synchronization!using observer} Other observer-based approaches to this problem were presented by Grassi and Masscolo~\cite{grassi97}, Azemi and Yaz~\cite{azemi00}, Lin et al.~\cite{lin05}, and more recently by Jiang, Tang, and Chen~\cite{jiang06}. Lian, Chiu, Chaing, and Liu~\cite{lian01} proposed a systematic design methodology using fuzzy observer-based approach.\index{synchronization!using fuzzy observer}\index{fuzzy observer!chaotic synchronization}

In this Section, we show how to combine fuzzy observers, chaos, and LMIs  to develop systems for secure communication. Specifically, we construct a fuzzy observer-based secure communication system using a chaotic system. The design is accomplished with LMIs. Our exposition is based on the paper by Lian et al.~\cite{lian01}.

\subsection{Fuzzy Representation of Chaotic Systems}

In our construction of a secure communication system, we employ a chaotic system represented as a fuzzy system.\index{chaotic system!represented as a fuzzy system} We use
the Lorentz's system described by a system of three differential equations,
\begin{eqnarray*}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array}\right] &=& \left[\begin{array}{c}
a(-x_1+x_2)\\
bx_1-x_2-x_1x_3\\
-cx_3+x_1x_2
\end{array}\right]\\
y&=& x_1,
%\label{eq:lorenz_sysy}
\end{eqnarray*}
where $a=10$, $b=28$, and $c=8/3$; see also Exercise~\ref{eq:lorentz_system}. To design Lorentz's system fuzzy representation,\index{Lorentz's system!fuzzy representation} we use $x_1$ to construct fuzzy rules. We note that $x_1(t) \in [-d, \; d]$, where $d=30$. We use two fuzzy rules:
\begin{eqnarray*}
\mb{Rule i}: &{}& \mb{ IF } x_1(t)\mb{ is } F_i(x_1(t)) \mb{ THEN }\\
&{}& \left\{\begin{array}{lcl}
\dot{\b x}(t)&=&\b A_i\b x(t)+\b b_i\\
\b y(t)&=&\b c_i\b x(t),\quad i=1,2,
\end{array}\right.
\end{eqnarray*}
where
\[
F_1(x_1(t))=\frac{1}{2}\left(1 + \frac{x_1(t)}{d} \right) \quad \mb{and} \quad F_2(x_1(t))=\frac{1}{2}\left(1 - \frac{x_1(t)}{d} \right),
\]
\[
\b A_1=\left[\begin{array}{ccc}
-10 & 10 & 0\\
28 & -1 & -d\\
0 & d & -\frac{8}{3}
\end{array}\right],\quad \b A_2=\left[\begin{array}{ccc}
-10 & 10 & 0\\
28 & -1 & d\\
0 & - d & -\frac{8}{3}
\end{array}\right],
\]
\[
\b c_1 = \b c_2 = \left[\begin{array}{ccc}
1 & 0 & 0
\end{array}\right],\quad \mb{and}\quad \b b_1=\b b_2=\b 0 .
\]
Note that on the interval $[-d, d]$, we have $F_i(x_1(t))\ge 0$, $i=1,2$. In addition,
\[
F_i(x_1)+F_2(x_1)=1.
\]
The fuzzy Lorentz system\index{fuzzy Lorentz system} takes the form
\begin{eqnarray}
\dot{\b x}(t) &=& \left(\frac{F_1}{F_1+F_2}\b A_1 + \frac{F_2}{F_1+F_2}\b A_2\right)\b x(t)\nonumber\\
&=& \left(F_1\b A_1 + F_2\b A_2\right)\b x(t) .
\label{eq:fuzzy_lorentz}
\end{eqnarray}
In Figure~\ref{fig:fuzzy_lorentz}, we show plots of the fuzzy Lorentz system over the time inerval $[0, 100]$ for the initial condition $\b x(0)=\left[\begin{array}{ccc}
-8 & -8 & 24
\end{array}\right]^{\top}$ .
\begin{figure}
\centerline{\psfig{figure=fuzzy_lorentz.eps,height=4.3in}}
\caption{Plots of trajectories of the fuzzy Lorentz system~\eqref{eq:fuzzy_lorentz}.}
\label{fig:fuzzy_lorentz}
\end{figure}

\subsection{Constructing Secure Communication System}

We first construct a transmitter. Our approach is to hide the message to be transmitted in a noise-like signal\index{message to be transmitted!hidden in a noise-like signal} being generated by a fuzzy chaotic system\index{fuzzy chaotic system} of the form
\begin{equation}
\left. \begin{array}{ccl}
\dot{\b x} &=& \sum_{i=1}^r \mu_i\left(\b A_i \b x + \b b_i + \b l_i m \right)\\
y &=& \b c\b x + m,
\end{array} \right\}
\label{eq:lo_transmitter}
\end{equation}
where
\[
\mu_i=\frac{F_i}{\sum_{i=1}^r F_i},
\]
$m$ is the message to be transmitted, and $\b l_i$, $i=1,2,\ldots,r$, are gain vectors selected in such a way that the matrices
\[
\b A_i - \b l_i\b c,\quad i=1,2,\ldots,r,
\]
are asymptotically stable.

The receiver is a fuzzy observer of the transmitter,\index{receiver!as a fuzzy observer of the transmitter}
\begin{equation}
\left. \begin{array}{ccl}
\dot{\tilde{\b x}} &=& \sum_{i=1}^r \mu_i\left(\b A_i \tilde{\b x} + \b b_i + \b l_i(y-\tilde{y}) \right)\\
\tilde{y}&=&\b c\tilde{\b x},
\end{array} \right\}
\label{eq:lo_receiver}
\end{equation}

Note that
\begin{eqnarray*}
y-\tilde{y}&=& \b c\b x + m - \b c\tilde{\b x}\\
&=& \b c(\b x-\tilde{\b x}) + m.
\end{eqnarray*}
Let
\[
\b e=\b x-\tilde{\b x},
\]
be the estimation error of the transmitter's state\index{estimation error!of the transmitter's state} $\b x$ by the receiver's state $\tilde{\b x}$. Then
\begin{equation}
y-\tilde{y}=\b c\b e + m
\label{eq:differ_y_and_tide_y}
\end{equation}
Thus if
\[
\b e(t) \to 0\quad\mb{as}\quad t\to \infty,
\]
then we can recover the message because
\[
y-\tilde{y} \to m \quad\mb{as}\quad t\to \infty .
\]

Our task then is to formulate conditions guaranteeing that $\b e(t) \to 0\quad\mb{as}\quad t\to \infty$. Let us analyze the dynamics of the estimation error,\index{dynamics!of the estimation error}
\begin{eqnarray}
\dot{\b e} &=& \dot{\b x} - \dot{\tilde{\b x}}\nonumber\\
&=& \sum_{i=1}^r \mu_i\b A_i \b x + \sum_{i=1}^r \mu_i\b b_i + \sum_{i=1}^r \mu_i\b l_i m - \sum_{i=1}^r \mu_i\b A_i \tilde{\b x} - \sum_{i=1}^r \mu_i\b b_i - \sum_{i=1}^r \mu_i\b l_i(y-\tilde{y})\nonumber\\
&=& \sum_{i=1}^r \mu_i \b A_i \b e  + \sum_{i=1}^r \mu_i\b l_i m - \sum_{i=1}^r \mu_i\b l_i (\b c\b x + m -\b c\tilde{\b x})\nonumber\\
&=& \sum_{i=1}^r \mu_i \b A_i \b e - \sum_{i=1}^r \mu_i\b l_i\b c (\b x  - \tilde{\b x})\nonumber\\
&=& \sum_{i=1}^r \mu_i (\b A_i - \b l_i \b c) \b e .
\label{eq:lo_sec_tr_err}
\end{eqnarray}
We are guaranteed that $\b e(t) \to 0\quad\mb{as}\quad t\to \infty$ if can select the gains $\b l_i$, $r=1,2,\ldots, r$ such that there exists a common $\b P=\b P^{\top} > 0$ so that
\begin{equation}
(\b A_i - \b l_i \b c)^{\top} \b P + \b P(\b A_i - \b l_i \b c) < 0,\quad r=1,2,\ldots, r .
\label{eq:enough for exist_sec_comm}
\end{equation}
The problem of solving simultaneously the above $r$ Lyapunov inequalities\index{solving simultaneously!Lyapunov inequalities} can be represented as an LMI feasibility problem.\index{LMI feasibility problem} Indeed, performing multiplications in~\eqref{eq:enough for exist_sec_comm} yields
\begin{equation}
\b A_i^{\top}\b P - \b c^{\top}\b l_i^{\top}\b P + \b P\b A_i - \b P\b l_i\b c < 0
\label{eq:after_mul_of_enough}
\end{equation}
Let
\[
\b M_i=\b P\b l_i.
\]
Then we also have that $\b M_i^{\top}=\b l_i^{\top}\b P$. Taking the above into account, we represent~\eqref{eq:enough for exist_sec_comm} in the form of $r$ LMIs,
\begin{equation}
\b A_i^{\top}\b P + \b P\b A_i - (\b c^{\top}\b M_i^{\top} + \b M_i \b c) <0, \quad r=1,2,\ldots, r .
\label{eq:lmi_form_of_enough}
\end{equation}
We summarize the above development in the form of the following proposition that appears in~\cite[p.~214]{lian01}.
\begin{proposition}
If there exist $\b P=\b P^{\top}>0$ and $\b l_i$, $r=1,2,\ldots, r$ such that for $r=1,2,\ldots, r$,
\begin{eqnarray*}
\b A_i^{\top}\b P + \b P\b A_i - (\b c^{\top}\b M_i^{\top} + \b M_i \b c) <0,\quad \b M_i=\b P\b l_i,
\end{eqnarray*}
then
\[
\b l_i = \b P^{-1}\b M_i, \quad r=1,2,\ldots, r
\]
and
the estimation error $\b e(t) \to 0$ as $t\to \infty$.
\end{proposition}

Our fuzzy Lorentz-based secure communication system\index{fuzzy Lorentz-based!secure communication system} consists of the transmitter~\eqref{eq:lo_transmitter} and the receiver~\eqref{eq:lo_receiver}. The message to be transmitted\index{message!to be transmitted} is
\[
m(t)=0.01\,\mb{square}(2\pi t),
\]
where $\mb{square}(\cdot)$ generates the square wave. In Figure~\ref{fig:lo_sec_comm3} we show the transmitted signal, $y = \b c\b x + m$.
\begin{figure}
\centerline{\psfig{figure=lo_sec_comm3.eps,height=3.3in}}
\caption{Transmitted signal, $\tilde{y}$.}
\label{fig:lo_sec_comm3}
\end{figure}

In Figure~\ref{fig:lo_sec_comm2}, we compare the message signal $m$ against its recovered version. As can be seen, after transients decay the recovered signal\index{recovered signal} matches the message signal perfectly.
\begin{figure}
\centerline{\psfig{figure=lo_sec_comm2.eps,height=3.3in}}
\caption{Plots of the message signal, $m$, and its recovered version at the receiver's end obtained from the difference, $y-\tilde{y}$; see~\eqref{eq:differ_y_and_tide_y}.}
\label{fig:lo_sec_comm2}
\end{figure}

Finally, in Figure~\ref{fig:lo_sec_comm1}, we show a state-space trajectory of of the transmitted signal,\index{transmitted signal} that is, the signal generated by the transmitter.
\begin{figure}
\centerline{\psfig{figure=lo_sec_comm1.eps,height=3.3in}}
\caption{State-space trajectory of the transmitter~\eqref{eq:fuzzy_lorentz}.}
\label{fig:lo_sec_comm1}
\end{figure}

We add that the initial condition of the transmitter was $\b x(0)=\left[\begin{array}{ccc}
1& 1 & 1
\end{array}\right]^{\top}$. We set the zero initial condition for the receiver.

\section{Notes}
Sandefur~\cite{sandefur90} is a lucid elementary introduction to chaotic
discrete-time dynamical systems.\index{chaotic!discrete-time dynamical systems} Peitgen, J\"{u}rgens and Saupe~\cite{peitgen92} and Alligood, Sauer and Yorke~\cite{alligood97} are superb introductions to chaos and
dynamical systems. Papers by May~\cite{may76} and Feigenbaum~\cite{feigenbaum83} are
recommended for their originality and clarity. Introductory texts to nonlinear
dynamics and chaos for scientists and engineers are the books by Moon~\cite{moon87}
and Thompson and Stewart~\cite{thompson86}. Numerical algorithms for chaotic systems
can be found in the book by Parker and Chua~\cite{parker89}. Fractals\index{fractals} are covered in
Mandelbrot~\cite{mandelbrot83}, Peitgen and Saupe~\cite{peitgen88}, Devaney and
Keen~\cite{devaney89}, and Bunde and Havlin~\cite{bunde94}.
Barnsley~\cite{barnsley88}, and Peitgen and Saupe~\cite{peitgen88} are recommended
for their coverage of iterated function systems.\index{iterated function systems} A reader who is interested in a
mathematical approach to chaotic behavior of dynamical systems should consult the
book by Arrowsmith and Place~\cite{arrowsmith90}. Peterson~\cite{peterson93} gives a
remarkable account of the development of models of celestial mechanics leading to
the discovery of chaos in the solar system.\index{chaos!in the solar system} The use of fractals in education, art,
music, fashion, chess and medicine is described in an edited volume by
Pickover~\cite{pickover96}. Other more recent papers on controlling chaotic systems,
in addition to the ones already referenced in the previous section,
are:~\cite{yu97}, \cite{zeng97}, \cite{yang98}. Applications of chaos in information
processing are discussed in the papers by Dmitriev, Panas and
Starkov~\cite{dmitriev91} and Andreyev et al.~\cite{andreyev96}. Aihara, Takabe and
Toyoda~\cite{aihara90} and Adachi and Aihara~\cite{adachi97} discuss the presence of
chaos in biological neurons\index{chaos!in biological neurons} and propose mathematical models of chaotic neural
networks. Applications of chaotic neural networks in the complex memory search for
patterns with a sought feature are discussed by Nara, Davis and
Totsuji~\cite{nara93}; see also~\cite{chan_chaos99}.

Edward Lorenz was a meteorologist at the Massachusetts Institute of Technology~(MIT). His discovery of chaotic behavior of his mathematical model of weather is a fascinating one. As reported by Gleick~\cite[Chapter~1]{gleick88} and Alligood et al~\cite[Chapter~9]{alligood97}, to speed up his computer output which was simulating a weather model, Lorenz used just three significant digits to print out the results of simulations although the calculations were performed using several more digits of accuracy. One day in the winter of 1961, Lorenz came across a simulation that he decided to repeat. He typed in an initial condition for the new simulation from the printout using the three digit format. He then performed the simulation expecting the same behavior as before of the repeated simulation. Instead, ``he saw something unexpected, something that planted a seed for a new science.''~\cite[page~16]{gleick88}. The resulting simulation was totally different from the previous one. The difference between the two solutions was the result of the Lorenz's weather model\index{Lorenz's weather model} sensitive dependence on initial conditions. In the computer memory, according to Gleick~\cite{gleick88}, the initial condition for the first run was 0.506127, while to save space, Lorenz entered 0.506 for the second run. This seemingly inconsequential difference caused the two simulations to flow apart until there was no resemblance at all between them. Later, Lorenz greatly simplified his weather model to just three nonlinear equations that are given in Exercise~\ref{ex:lorenzb}. The phase portrait of these famous equations resemble butterfly's wings,\index{butterfly's wings!symbol for early activity in chaotic systems} which became a symbol for the early activity in chaotic systems.


\section{Exercises}
\begin{exercise}\label{eq:lorentz_system}\rm
Simulate the system of differential equations
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array}\right]=\left[\begin{array}{c}
a(-x_1+x_2)\\
bx_1-x_2-x_1x_3\\
-cx_3+x_1x_2
\end{array}\right],\quad \b x(0)=\left[\begin{array}{c}
-8\\
-8\\
24
\end{array}\right],
\label{eq:lorenzode}
\end{equation}
where $a=10$, $b=28$, and $c=8/3$. The set of differential equations~(\ref{eq:lorenzode}) is commonly referred to as the {\em Lorenz equations},\index{Lorenz equations} after the American meteorologist E.~N.~Lorenz who first published these equations in 1963. They model fluid flows in the atmosphere. Plot $x_2$ versus $x_1$, then $x_3$ versus $x_1$, and also $x_3$ versus $x_2$. Use MATLAB's command \verb!plot3! to obtain a 3-D plot of the solution of the Lorenz equations.
\label{ex:lorenzb}
\end{exercise}
\begin{exercise}\rm
Simulate the system of the difference equations
\begin{equation}
\left[\begin{array}{c}
x_1(k+1)\\
x_2(k+1)
\end{array}\right]=\left[\begin{array}{c}
1-x_2(k)-\left| x_1(k)\right|\\
x_1(k)
\end{array}\right],\quad \b x(0)=\left[\begin{array}{c}
-0.1\\
0
\end{array}\right].
\label{eq:gingereq}
\end{equation}
The solution ``fills'' a region that looks like a ``gingerbread man''.\index{gingerbread man!chaotic set} For this reason, Devaney, in~\cite[page~149]{peitgen88}, refers to this example of a chaotic set in the plane as a {\em chaotic gingerbreadman}.\index{chaotic gingerbreadman}
\end{exercise}
\begin{exercise}\rm
Simulate the system of differential equations
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array}\right]=\left[\begin{array}{c}
-x_2-x_3\\
x_1+ax_2\\
b+x_3(x_1-c)
\end{array}\right],
\label{eq:rosslera}
\end{equation}
where $a=0.398$, $b=2$, and $c=4$. The set of differential equations~(\ref{eq:rosslera}) is referred to as the {\em R\"{o}ssler equations},\index{R\"{o}ssler equations} after O.~E.~R\"{o}ssler who reported on their chaotic behavior in 1976. R\"{o}ssler called the chaotic attractor of the above system the ``spiral chaos''.\index{spiral chaos} This type of chaotic attractor is also called a {\em folded band}.\index{folded band}\index{chaotic attractor} Select an initial condition, integrate the equations, then plot $x_2$ versus $x_1$, $x_3$ versus $x_1$, and $x_3$ versus $x_2$. Use MATLAB's command \verb!plot3!\index{MATLAB's command!plot3} to obtain a 3-D plot of the folded band attractor.
\end{exercise}
\begin{exercise}\rm
For the dynamical system
\begin{equation}
x(k+1)=x(k)+r\left( 0.8-\mb{atan}(x(k)) \right),
\label{eq:huiandzakmap}
\end{equation}
plot a bifurcation diagram\index{bifurcation diagram} for the parameter $r$ in the interval $[3, 12]$. The map~(\ref{eq:huiandzakmap}) was analyzed in~\cite{hui94}.
\end{exercise}
\begin{exercise}\rm
For the dynamical system
\begin{equation}
x(k+1)=x(k)+r\left( 0.8-\frac{1-e^{-x(k)}}{1+e^{-x(k)}} \right),
\label{eq:bipologmabd}
\end{equation}
plot a bifurcation diagram\index{bifurcation diagram} for the parameter $r$ in the interval $[10, 55]$.
\end{exercise}
\begin{exercise}\rm
For the dynamical system
\begin{equation}
x(k+1)=x(k)+r\left( 0.8-\frac{1}{1+e^{-x(k)}} \right),
\label{eq:unipolomabf}
\end{equation}
plot a bifurcation diagram for the parameter $r$ in the interval $[10, 45]$.
\end{exercise}
\begin{exercise}\label{exer:chaosegame}\rm
Write a MATLAB program of the following ``chaos game'':\index{chaos game} Select three points that are vertices of an equilateral triangle. Choose at random one of the vertices as a starting point. Then, choose at random another vertex and move halfway to it from the starting vertex. Mark this halfway point with a dot. This halfway point is your new starting point. Choose at random one of the vertices and move halfway to the chosen vertex from the new starting point. Mark the new halfway point with the dot. Repeat the process. The outcome of this chaos game is a fractal called the Sierpi\'{n}ski gasket. Generate a plot of the Sierpi\'{n}ski gasket. Other methods of constructing the Sierpi\'{n}ski gasket\index{Sierpi\'{n}ski gasket} are described in~\cite{peitgen92} and in Exercise~\ref{exer:gasket}. The term ``chaos game'' was coined by M.~F.~Barnsley, the author of~\cite{barnsley88}.
\end{exercise}
\begin{exercise}\label{exer:gasket}\rm
The result of the chaos game described in Exercise~\ref{exer:chaosegame} is the Sierpi\'{n}ski gasket that can also be constructed using the following procedure. Start with an equilateral triangle. Divide it into four equal sub-triangles, quarters, and throw away the middle one, upside-down, triangle. This step leaves three equilateral triangles half the size of the original triangle. In the second step, divide each of the remaining three triangles into four equal quarter-triangles and throw away the middle, upside-down, triangle in each of the three triangles. This step leaves nine triangles one-quarter the size of the original triangle. Repeat the process. The remaining area as the number of iterations tends to infinity is called the Sierpi\'{n}ski gasket.
\begin{itemize}
\item[(i)] Show that the area of the Sierpi\'{n}ski gasket\index{Sierpi\'{n}ski gasket!zero area of} is zero.
\item[(ii)] Determine the fractal dimension of the Sierpi\'{n}ski gasket.\index{Sierpi\'{n}ski gasket!fractal dimension of}
\end{itemize}
\end{exercise}
\begin{exercise}\label{exer:carpet}\rm
In this exercise, we concern ourselves with another fractal, the Sierpi\'{n}ski carpet.\index{Sierpi\'{n}ski carpet} Start with a square. Divide it into nine subsquares, each one-third the size of the original square. Remove the central square. Repeat the procedure. The remaining area as the number of iterations goes to infinity is called the Sierpi\'{n}ski carpet.
\begin{itemize}
\item[(i)] Show that the area of the Sierpi\'{n}ski carpet is zero.\index{Sierpi\'{n}ski carpet!zero area of}
\item[(ii)] Determine the fractal dimension of the Sierpi\'{n}ski carpet.\index{Sierpi\'{n}ski carpet!fractal dimension of}
\end{itemize}
\end{exercise}

\begin{exercise}\label{exer:dbl_pend}\rm
For the double pendulum\index{double pendulum!chaotic behavior} system described in Subsection~\ref{subsection:two_link_robot} set the input torques to zero and simulate the system behavior starting from the initial state $\b x(0)=\left[\begin{array}{cccc}
0 & 0 & \pi & .001
\end{array}\right]^{\top}$.
Verify that one of the equilibrium states corresponding to zero input torques is $\b x_e=\left[\begin{array}{cccc}
0 & 0 & \pi & 0
\end{array}\right]^{\top}$. Linearize the state-space model about this equilibrium point. Then, use the LQR approach to compute the state feedback matrix\index{LQR approach!to compute state feedback matrix} $\b K$ for the weight matrices
\[
\b Q=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 10 & 0\\
0 & 0 & 0 & 10
\end{array}\right]\quad\mb{and}\quad \b R=\left[\begin{array}{cc}
1 & 0\\
0 & 10
\end{array}\right].
\]
Solve the Lyapunov equation\index{Lyapunov equation}
\[
(\b A-\b B\b K )^{\top}\b P_l+\b P_l(\b A-\b B\b K)=-\b Q .
\]
Suppose that the bound on the input torques are: $|T_1|\le 5$ and $|T_2|\le 2$. Use the method of Example~\ref{ex:henonex} to find $V_{\max}$ and apply the chaotic control law to stabilize the closed-loop system\index{chaotic control law!to stabilize the closed-loop system} with respect to $\left[\begin{array}{cccc}
0 & 0 & \pi & 0
\end{array}\right]^{\top}$. Plot $x_1$, $x_2$, and $V=(\b x-\b x_{e})^{\top}\b P_l(\b x-\b x_{e})$ versus time.
\end{exercise}
