%This is a Revtex file.
\documentstyle[aps,prl,twocolumn]{revtex}
%\documentstyle[preprint,aps]{revtex}
\def\e{\varepsilon}
\def\d{\delta}
\begin{document}
\draft
\title{Obstructions to shadowing \\when a Lyapunov
exponent fluctuates about zero}
\author{ Silvina Dawson}
\address{Center for Nonlinear Studies and Theoretical Division\\Los
Alamos National Laboratory\\
Los Alamos, NM 87545}
\author{Celso Grebogi}
\address{Department of Mathematics \\
University of Maryland\\
College Park, MD 20742}
\author{Tim Sauer}
\address{Department of Mathematical Sciences \\George Mason
University \\Fairfax, VA 22030}
\author{James A. Yorke}
\address{Institute for Physical Science and Technology \\
University of Maryland\\
College Park, MD 20742}
\maketitle
\newpage
\begin{abstract}
We study the existence or nonexistence of the true trajectories of chaotic
dynamical systems that lie close to computer-generated trajectories. The
nonexistence of such shadowing trajectories is caused by finite-time Lyapunov
exponents of the system fluctuating about zero. A dynamical mechanism of the
unshadowability is explained through a theoretical model and identified in
simulations of a typical physical system. The problem of fluctuating Lyapunov
exponents is expected to be common in simulations of higher-dimensional systems.
\end{abstract}
\pacs{05.45.+b, 05.40.+j, 06.50.Dc}
\narrowtext
Physical theory is based on differential equation models. Computer
simulation using the equations is often used to obtain information,
such as long-term statistics, on the system being modeled.
In climate modeling, for example, statistics of temperature and
rainfall might be relevant.
A basic requirement needed to interpret the result of a simulation
is that the behavior of a solution determined by the computation,
which
is afflicted with small errors due to truncation and roundoff, is the
same as the behavior of some true solution of the system under study.
If there is a difference in behavior between computed solutions and
actual solutions, the investigator cannot proceed.
A climate model that continues to repeat winter conditions all year
long because of accumulated
numerical errors will not be useful for computing the mean yearly
temperature.
This problem is especially acute when the system is chaotic. In that
case,
trajectories exhibit sensitive dependence on initial conditions: two
trajectories with initial conditions that are extremely close tend to
diverge exponentially from one another. Because of this effect,
a small truncation or rounding error made at any step during the
computation will tend to be greatly magnified by future evolution
of the system. In view of this, it is natural to ask under what
conditions the computed
trajectory will be close to a true trajectory of the model.
Previous work on this topic \cite{GHY1,GHYS,SY,CP,SS}
has resulted in computational techniques for
``verifying'' computer-generated trajectories for low-dimensional
chaotic systems --
that is, to produce a computer-assisted
proof of the existence of a true trajectory of the system,
called a shadowing trajectory, that
closely tracks the computer-generated pseudotrajectory.
(Even in this case, the statistical properties of the
system may not be decided.)
Despite these positive results, not every pseudotrajectory
can be shadowed. We believe that in systems with
high-dimensional chaos, trajectories with intrinsic noise,
such as computer-generated pseudotrajectories, can be shadowed
only for short times. Consideration of simple examples of
nonlinear maps \cite{SY} illustrates that
there are critical points of trajectories where roundoff error or
other noise, perhaps introduced at a distant part of the trajectory,
can introduce new behavior. At such ``glitches'' all true
trajectories diverge from the computed trajectory.
In this case, when there is
no true trajectory that stays near the computed trajectory, we say
that the computed trajectory is unshadowable.
There is no way known to ensure that a given computer simulation is
representative
of a true trajectory of the system (or even visits the entire phase
space attractor) except when a proof of its shadowability is
available. If the intrinsic noise is being injected by the
environment itself rather than by a truncation error of a computer
simulation, unshadowability raises interesting questions about
the validity of deterministic modeling for the system.
In this Letter we describe a cause of unshadowable pseudotrajectories,
that is likely to occur widely in higher-dimensional
chaotic dynamical systems. We will say that a Lyapunov exponent of a
trajectory ``fluctuates about zero'' if for any positive number $T$,
the time-$T$ Lyapunov exponent spends arbitrarily long stretches of
the trajectory being positive, and arbitrarily long stretches being
negative. The finite-time Lyapunov exponents
quantify the expansion and contraction of phase space along the trajectory
over a time span of $T$. We will show that the existence of oven one
Lyapunov exponent which fluctuates about zero cause computer-generated
trajectories to be unshadowable. Since a positive Lyapunov exponent is
the signature of chaotic dynamics, this is a fact of critical
importance to researchers studying the existence of chaos in computer
models.
The manner in which a Lyapunov exponent fluctuating about zero leads
to unshadowable pseudo-trajectories is illustrated well by a
theoretical model studied by Abraham and Smale \cite{AS} in 1970.
In this example, there is an invariant set containing two fixed
points, one with a single local expanding direction, and one with a
two-dimensional local expanding set. Typical trajectories wandering
through the invariant set spend arbitrarily long times near each of the fixed
points. The second largest Lyapunov exponent of such a trajectory
is positive in trajectory segments near one, and is negative near
trajectory segments near the other, so this exponent fluctuates about zero.
A ball of initial conditions beginning near the first fixed point will
be squeezed into a line segment (with small thickness) under evolution
of the dynamics. A computer-generated trajectory beginning in the
ball, with truncation error $\d$, will be
displaced a distance of $\d$ from the line segment. When the region
around the numerical trajectory develops a second expanding direction
by visiting a neighborhood of the second fixed point,
the numerical trajectory will be pushed away exponentially fast from
the line segment of
true trajectories, resulting in an unshadowable trajectory.
Although this example is non-physical, we have found similar
behavior in models of typical mechanical systems such as the
double rotor \cite{Rom}. For certain parameter settings, the double
rotor has a chaotic attractor whose second-largest Lyapunov exponent
fluctuates about zero. As we discuss below, this effect causes almost every
moderately-long numerical trajectory to be unshadowable.
In order to quantify the phenomenon of unshadowability, we introduce
the ideas of continuous shadowability and brittleness.
A {\it continuously-shadowable} pseudotrajectory is a
computer-generated trajectory that can be
continuously deformed into a true trajectory, in such a way that
the errors at each trajectory point are decreased monotonically to
zero. Although this appears to
be a stronger requirement than for a shadowable pseudotrajectory,
it turns out that the hypotheses of the original
Anosov-Bowen shadowing theorems \cite{A,B,MH}, as well as
computer-assisted shadowing methods mentioned above, imply not only
that pseudo-trajectories are shadowable,
but that they are continuously shadowable. We
argue that continuous deformability to a true trajectory is a minimum
requirement for accepting a computed simulation as meaningful
information.
A fundamental phenomenon connecting continuously-shadowable
pseudo-trajectories to their shadowing (true) trajectories
is the existence of a constant of proportionality between the
error magnitude of the pseudotrajectory and the distance
the pseudotrajectory must move in phase space to be deformed into a
true trajectory. We call this constant of proportionality the {\it
brittleness} of the pseudotrajectory. An obvious necessary condition
for continuous shadowability is that the brittleness times the error
magnitude of the pseudotrajectory is smaller than the extent of the
attractor in phase space.
This leads to a practical algorithm for investigating the shadowability of a
computer-generated trajectory. Using Jacobian information available
from the simulation, it is possible to calculate a first-order
approximation for the brittleness (the {\it test brittleness}),
using small deformations of
the pseudotrajectory. Knowledge of the test brittleness is a useful
diagnostic for continuous shadowability of the
pseudotrajectory, or lack thereof.
Let $f$ denote the map which represents one time step of the dynamics.
For example,
it may represent the time-T map
produced by an ODE-solver with one-step truncation error bounded by
$\d$.
(By one-step error, or noise, we mean the discrepancy between a time-T
step of the
ODE solver and a true time step of the differential equation, starting
from the previous point.) If the
present time is $t_0$ and the present state of the dynamical system is
$p_0$, then the correct state at time $t_0+T$ is $f(p_0)$. The
ODE-solver will produce $p_1$, where $|p_1-f(p_0)| < \d$.
Then $p_0$ and $p_1$ are two points of a
$\d$-{\it pseudotrajectory} of the dynamical system $f$. Further
integration of
the simulation of $f$ results in a $\d$-pseudotrajectory $\{p_0, \ldots,
p_N\}$ of length $N+1$, where $|p_{i+1}-f(p_i)| < \d$ for $i=0, \ldots,
N-1$.
It would be desirable to know the existence of a true trajectory $x_0,
\ldots, x_N$ (that is, $f(x_i)=x_{i+1}$ for $i=1, \ldots, N-1$) that
$\e$-shadows the pseudotrajectory. The true trajectory $\{x_i\}^N_{i=0}$ is
said to be an $\e$-{\it shadowing trajectory} for the pseudotrajectory
$\{p_i\}^N_{i=0}$ if $|x_i-p_i|<\e$ for $i=0,\ldots, N$. Because
of the exponential divergence of trajectories,
if a shadowing trajectory exists, then the initial condition $x_0$ will
differ from $p_0$.
In fact, the first point $x_0$ of the true shadowing trajectory is expected
to be found along the unstable direction emanating from $p_0$.
The shadowing lemma of Anosov \cite{A} and Bowen \cite{B} is a
theoretical result for dynamical systems with hyperbolic structure.
An attractor is called {\it hyperbolic} if the tangent
spaces at each point of the attractor can be decomposed into uniformly
expanding and contracting subspaces, such that the angle between these
subspaces is bounded away from zero. The shadowing lemma states that
for each nonzero distance
$\e$, there exists an error magnitude $\d$ such that each
$\d$-pseudotrajectory can be $\e$-shadowed. Furthermore, under the
hyperbolicity assumption, each such pseudotrajectory can be
continuously shadowed within $\e$. (See, for example, the proof of the
shadowing lemma in \cite{MH}.)
We say that a pseudotrajectory has a {\it glitch} at
point $n$ if $\{p_i\}^n_{i=0}$ can be continuously shadowed, but
$\{p_i\}^{n+1}_{i=0}$ cannot. One type of glitch is caused by a
Lyapunov exponent fluctuating about zero, as described above. A
schematic representation of a second type
of glitch is shown in Fig.~1. Assume that $q$ is a fixed point.
Assume that the pseudotrajectory has no error until
iterate $n$, when error pushes $p_n$ across the stable manifold of
$q$. True
trajectories must follow the unstable manifold of $p_{n-1}$, which
separates exponentially from the stable manifold of $q$, so that no continuous
shadowing trajectory can exist.
When a pseudotrajectory is continuously moved to a true
trajectory by deforming its noise to zero, we will consider the distance
moved by the pseudotrajectory to be the maximum distance any single
trajectory point was moved. We call this the {\it shadowing distance}.
The brittleness is the constant of
proportionality between the shadowing distance of the
pseudotrajectory and
the original magnitude of the one-step error.
The proportionality holds over a large range of noise levels, as long
as the pseudotrajectory itself is not significantly changed.
The brittleness is independent of the error magnitude (for small
one-step errors) but depends on the error directions -- a
different set of one-step errors of the same magnitudes would in
general lead to a different proportionality
constant.
The brittleness should be defined to
be the maximum of this magnification factor found over all
possible error directions.
The brittleness of a pseudotrajectory is a measure of its inability
to be shadowed. If a pseudotrajectory is created with noise level
$10^{-10}$ for a chaotic attractor of unit size, and if its
brittleness is greater than $10^{10}$, then one cannot expect a true trajectory
closely shadowing the pseudotrajectory. For hyperbolic systems,
pseudo-trajectories of infinite length have finite (although possibly
very large) brittleness
\cite{MH}. For nonhyperbolic systems, one typically finds the
brittleness to be an increasing function of the orbit length.
As the length of a
trajectory of a nonhyperbolic chaotic process increases, the
brittleness grows as the trajectory gets increasingly close
to nonhyperbolic regions of the dynamics. The expected length
between glitches is therefore related to the amount of hyperbolicity
possessed by the system.
Although the brittleness cannot be computed exactly without knowing the
true trajectory, the test brittleness is a
first-order approximation to this constant that is explicitly computable.
The test brittleness is determined from the Jacobians (first derivative
matrices) at the points $p_i$. Assuming the error at step $i$ of the
computation is $\d_i$, we want to approximate the correction $c_i$ such
that $f(p_i+c_i)=p_{i+1}+c_{i+1}$ for $i=0,\ldots, N-1$. Write $c_i$
as the sum $c_i=s_i+u_i$ of components in the stable (contracting) and
unstable (expanding)
directions at $p_i$. Set $s_0=u_N=0$, and recursively solve
\begin{equation} \label{app1}
s_{i+1}=S(Df(p_i)s_i+\d_i)
\end{equation}
and
\begin{equation} \label{app2}
u_{i}=U(Df(p_{i+1})^{-1}[u_{i+1}+\d_i]),
\end{equation}
for $s_i$ and $u_i$, $i=0,\ldots,N-1$, where $S$ and $U$ denote
projection onto the stable and unstable directions, respectively.
Then the ratio of the maximum magnitude of correction $c_i$ to the
magnitude of one-step error $\d_i$ is the test brittleness. Since the
computation is linear in $\d_i$, it is possible to choose the vectors $\d_i$ to
have magnitude one. The $\d_i$ are typically chosen to have randomly
varying directions.
The double rotor map \cite{Rom} is a four-dimensional map which
describes the
time evolution of a mechanical system consisting of two connected
massless rods.
The first rod rotates
around a fixed pivot; the second rod pivots around
the opposite end of the first rod. There is a mass at the free end
of the first rod, and equal masses at the ends of the second rod.
A delta-function vertical impulse $f(t)$, of magnitude
$\rho$, is applied to one of the ends at a constant time interval, at
which the system's four phase
variables (the two angular positions and momenta) are recorded.
Interesting dynamical behavior is exhibited by this system for various
settings of the parameter $\rho$.
Application of the test brittleness algorithm to computer simulations
of the double rotor are shown in Fig.~2. The difference between the
shadowable case (parameter $\rho=9$) and the unshadowable case
($\rho=8$)is clear from this figure. In the
vertical axis of Fig.~2(a), we graph the test brittleness of
a length 10000 pseudotrajectory for $\rho=9$, created with one-step error
magnitude $10^{-10}$. The vertical extent of the graph is
about $10^7$, which is the test brittleness for this
trajectory. We should then expect the shadowing distance to be about
$10^{-3}$. Fig.~2(b) shows the same information for a typical
pseudotrajectory of the double rotor with $\rho=8$. The test
brittleness in this case is seen to be greater than $10^{40}$. This
leads us to the prediction that at a minimum, 40 decimal digits of accuracy
will be needed per iteration step in order to shadow a typical length
10000 trajectory of the double rotor with $\rho=8$.
The brittleness of a pseudotrajectory increases with its length.
For an orbit of length $10^5$, our estimate of the test
brittleness is $10^{100}$.
The explanation of the difference in shadowability for the cases
$\rho=8$ and $\rho=9$ lies in the different degrees of
hyperbolicity of the two systems.
A numerical study of the behavior of finite-time Lyapunov exponents
for the two parameter values of the double rotor
is shown in Fig.~3. For $\rho=8$, the finite-time
Lyapunov exponents show fluctuation between one and two positive
exponents. The second largest exponent fluctuates about zero.
For $\rho=9$, there are consistently two positive exponents.
A close examination of the dynamics of the double rotor reveals an
explanation for the fluctuating number of positive finite-time
Lyapunov exponents in the $\rho=8$ case. There are
many periodic orbits embedded in the attractor whose local behavior
varies in
a qualitative way. Some of the periodic points have one expanding
direction and three contracting directions, while others have two
expanding and two contracting \cite{Rom}. As a trajectory (or
pseudotrajectory)
moves densely through the attractor, its number of positive
finite-time Lyapunov exponents varies as it moves among the varying type of
periodic orbits. Although the double rotor with $\rho=8$ is the first
physical example
for which this behavior has been demonstrated, we expect it to be
commonly found in higher-dimensional systems.
\acknowledgments
This research was supported in part by grants from the NSF
(Computational Mathematics and Physics programs), and DOE
(Office of Scientific Computing).
\begin{references}
\bibitem{GHY1}
C. Grebogi, S. Hammel, and J. Yorke,
J. Complexity
{\bf 3}, 136 (1987);
Bulletin A.M.S.
{\bf 19}, 465 (1988).
\bibitem{GHYS}
C. Grebogi, S. Hammel, J. Yorke, and T. Sauer,
Phys. Rev. Lett. {\bf 65}, 1527 (1990).
\bibitem{SY}
T. Sauer and J. Yorke,
Nonlinearity {\bf 4}, 961 (1991).
\bibitem{CP}
S.-N. Chow and K. Palmer,
Dynamics and Diff. Eqns. {\bf 3}, 361 (1991);
J. Complexity {\bf 8}, 398 (1992).
\bibitem{SS}
J.M. Sanz-Serna and S. Larsson,
Appl. Num. Math. {\bf 13}, 181 (1993).
\bibitem{AS}
R. Abraham and S. Smale,
Proc. Symp. Pure Math. (AMS) {\bf 14}, 5 (1970).
\bibitem{Rom}
C. Grebogi, E. Kostelich, E. Ott, J.A. Yorke, Physica D {\bf 25}, 347
(1987);
R. Romeiras, C. Grebogi, E. Ott, W.P. Dayawansa,
Physica D {\bf 58}, 165 (1992).
\bibitem{A}
D. V. Anosov,
Proc. Steklov Inst. Math. {\bf 90}, 1 (1967).
\bibitem{B}
R. Bowen,
J. Diff. Eq. {\bf 18}, 333 (1975).
\bibitem{MH}
K. Meyer and G.R. Hall,
{\it Introduction to Hamiltonian Dynamical Systems and the N-body Problem.}
(Springer-Verlag, New York, 1992).
\end{references}
\begin{figure}
\caption{A near-tangency from a nonhyperbolic system. A small error
in the computation of $f(p_{n-1})$ can push $p_n$ across a stable
manifold, resulting in a glitch.}
\label{fig1}
\end{figure}
\begin{figure}
\caption{First-order approximation of the shadowing distance per unit
one-step error, or brittleness, as a function of trajectory point.
The test brittleness is the vertical range of the graph. Results are
shown for a 10000 point trajectory of the double rotor with parameters
(a) $\rho=9$; (b) $\rho=8$.}
\label{fig5}
\end{figure}
\begin{figure}
\caption{Estimates of the four time-$T$ Lyapunov exponents of the double
rotor, for $T = 100, 200, 300$. A dozen simulations were done for each
$T$. (a) In the $\rho=9$ case, there is consistently only one
positive finite-time Lyapunov exponent. (b) For $\rho=8$,
trajectory segments alternate between one- and two-dimensional
expanding subspaces.}
\label{fig7}
\end{figure}
\end{document}