%% 13. Chaos
% <latex>
% \setcounter{page}{156}
% </latex>

%%
% <latex>
% Having just considered functions that are truly random,
% at least in principle, we now turn to the celebrated phenomenon of
% solutions of ODE\kern .5pt s that ``look random'' 
% even though they are not.  Besides this property,
% two other hallmarks of chaos
% are {\em sensitive dependence on initial conditions} and
% {\em strange attractors.}  In this
% chapter we explore these notions, connecting
% them with the number known as the Lyapunov exponent.
% </latex>

%% 
% <latex>
% We begin with the Lorenz equations as given in equation (10.7),
% $$ u' = 10(v-u), \quad v' = u(28-w)-v, \quad w' = uv -(8/3)w.
% \eqno (13.1) $$ 
% In Figure 10.11 we plotted the trajectory for $t\in [\kern .3pt 0,30\kern .3pt]$
% emanating from the initial data
% $$ u(0)= v(0) = -15, ~ w(0) = 20. \eqno (13.2) $$
% Here is an image of this trajectory in phase space, seen from
% the angle that amounts to a projection onto the $u$--$w$ plane.
% </latex>
ODEformats, N = chebop(0,30);
N.op = @(t,u,v,w) [diff(u)-10*(v-u); ...
    diff(v)-u*(28-w)+v; diff(w)-u*v+(8/3)*w];
N.lbc = [-15; -15; 20]; [u,v,w] = N\0;
plot(u,w,CO,ivpnl,LW,0.4), axis([-25 25 0 50])
c = 6*sqrt(2);
hold on, plot([0 c -c],[0 27 27],'.k',MS,8), hold off 
xlabel('$u$',FS,10,IN,LT), ylabel('$w$',FS,10,IN,LT)
title(['Fig.~13.1.~~Lorenz trajectory projected '...
       'onto $u$--$w$ plane'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% This image reveals the famous ``butterfly'' structure of 
% the {\bf strange attractor} for the Lorenz equations.\footnote{Not
% to be confused with the {\em butterfly effect\/} that also
% originates with Lorenz (not his famous 1963 paper, but a later one
% from 1969).
% If a butterfly flaps its wings in the Amazon, thanks
% to sensitive dependence of the earth's weather on initial conditions,
% might that cause a hurricane in Texas?}
% What it means to be a strange attractor is that typical
% orbits rapidly approach this
% set, which is not a simple curve or surface but
% has the form of a fractal.  For the Lorenz equations, the
% fractal dimension is approximately $2.06$, so
% the butterfly is just a little bit thicker than a
% 2-dimensional manifold.\footnote{The {\em fractal\/} or
% {\em Hausdorff dimension\/} of a set is a precisely
% defined mathematical quantity that may take any
% nonnegative real value, not just an integer.}
% </latex>

%%
% Now let us look at dependence on initial conditions.  Here
% the function $v(t)$ just computed is plotted in green, as in
% Chapter 10.  In addition, another function $\tilde v(t)$ is plotted
% in red that arises from the very slightly perturbed initial data
% $$ u(0)= v(0) = -15, ~ w(0) = 20.00001.\eqno (13.3) $$
%%
% <latex> \vskip -1.5em </latex>
N.lbc = [-15; -15; 20.00001]; [u2,v2,w2] = N\0;
plot(v2,'r'), hold on, plot(v,CO,ivpnl), hold off
title(['Fig.~13.2.~~Lorenz equations: ' ...
       'sensitive dependence on initial data'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The two initial conditions differ by less than
% one part in $10^6$, so one would expect $v$ and $\tilde v$ to
% be close to one another, and for about 13 time units, so they are.
% But throughout that time, the trajectories
% are steadily (on average) separating, and
% eventually the differences are of size $O(1)$ and
% the functions become completely uncoupled.
% An initial difference of size $10^{-5}$ has grown
% by six orders of magnitude in 15 time units.
% </latex>

%%
% Everything in this problem is smooth, differentiable, and
% indeed analytic.  For
% example, the derivative of the value $v(30)$ with respect
% to the third initial condition $w(0)$ is a perfectly well-defined
% number --- it is just that it is a very big one!  Here is
% an approximation to 
% the derivative of $v(5)$ with respect to $w(0)$,
(v2(5)-v(5))/(w2(0)-w(0))

%%
% <latex>
% \noindent
% Here is the analogous approximation to
% the derivative of $v(10)$ with respect to $w(0)$,
% </latex>
(v2(10)-v(10))/(w2(0)-w(0))

%%
% <latex>
% \noindent
% Note that these numbers are rapidly growing with $t$.
% If we compute the corresponding finite difference at $t=30$, we get
% </latex>
(v2(30)-v(30))/(w2(0)-w(0))

%%
% <latex>
% \noindent
% but this greatly underestimates the actual derivative
% of $v(30)$ with respect to $w(0)$,
% which is on the order of $10^{13}$.  We can explain the
% underestimation by noting that the
% initial perturbation by $0.00001$ is not nearly small enough to
% grow by 13 orders of magnitude before reaching size~$O(1)$,
% at which point it stops growing.
% </latex>

%%
% <latex>
% To see more of this effect, it is interesting to plot
% $\|\tilde{{\bf y}}(t) - {\bf y}(t)\|$ on a log scale,
% where ${\bf y}$ and $\tilde {\bf y}$ are the
% 3-component vector solutions resulting
% from the original and the perturbed initial conditions.
% (The norm $\|\cdot\|$ we use is
% the square root of the sum of the squares of the components
% $u(t), v(t), w(t)$.)
% Here is such a plot, and it reveals two distinct regimes.
% For $t<15$, the perturbations grow exponentially.
% For $t>15$, they reach size $O(1)$ (numbers roughly in the
% range $10$--$100$) and stop growing.  
% </latex>
tt = linspace(0,30,400);
err = sqrt((u2(tt)-u(tt)).^2+(v2(tt)-v(tt)).^2+(w2(tt)-w(tt)).^2);
semilogy(tt,err,'k'), ylim([1e-6 1e3])
title(['Fig.~13.3.~~Exponential growth of perturbations ' ...
       'to scale $O(1)$'],FS,11)
xlabel('$t$',FS,10,IN,LT)
ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The growth rate of this initial phase is known as the
% {\bf Lyapunov exponent,}~$\lambda$, for this ODE.\ \ We
% can approximate it
% by a least-squares fit to the data for $t\in[\kern .3pt 0,15]$.
% </latex>
ii = find(tt<=15);
c = polyfit(tt(ii),log(err(ii)),1); e0 = c(2); lambda = c(1)

%%
% <latex>
% \noindent
% The significance of this number is that, approximately
% speaking, perturbations
% in the Lorenz trajectories typically grow initially at the rate
% $$ \| \tilde {\bf y}(t) - {\bf y}(t) \| \approx C e^{0.91 t}. $$
% We can plot this fit to the data as a dashed line.
% </latex>
hold on, semilogy([0 18],exp(e0+lambda*[0 18]),'--r',LW,1.5)
hold off
title('Fig.~13.4.~~Slope of line = Lyapunov exponent $0.91$',...
    FS,11,IN,LT)
xlabel('$t$',FS,10,IN,LT)
ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,9,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% In this experiment, we started from a finite perturbation,
% which grew exponentially over a finite time span.  The mathematical
% definition of a Lyapunov exponent is based on 
% infinitesimal perturbations, which may grow forever.
% To be precise, $\lambda$ is defined by this formula, involving
% a supremum over all initial perturbations $\delta {\bf y}(0)$:
% $$ \lambda = \lim_{t\to\infty} \limsup_{\|\delta {\bf y}(0) \|\to 0} t^{-1}
% \log {\|\delta {\bf y}(t)\|\over \|\delta {\bf y}(0)\|}. \eqno (13.4) $$
% In principle this number might vary from one initial point ${\bf y}(0)$
% to another, in which case one would typically be interested in the
% supremum over all ${\bf y}(0)$.\footnote{A further complication is that
% what we have called the Lyapunov exponent should more properly be termed
% the {\em maximal\/} Lyapunov exponent.  For an $n$-dimensional dynamical
% system, imagine initial conditions in the form of an infinitesimal
% spherical cloud or blob centered at a point ${\bf y}(0)$.  As $t\to\infty$,
% the blob will elongate and compress into an ellipsoid, with the
% lengths of the axes of the ellipsoid growing or decaying at
% different exponential rates on average
% as $t\to \infty$ in a manner that can be made precise by the tool known
% as the singular value decomposition.  These growth and decay
% rates are the Lyapunov exponents of the system.}
% </latex>

%%
% The Lyapunov exponent for the Lorenz equations
% is known to be about 0.9056.  By chance, our
% experiment has matched this number more closely than should really be
% expected, statistically speaking (see Exercise 13.1).

%%
% <latex>
% The article that launched the study of chaos
% was published by Edward Lorenz in 1963 with the
% title ``Deterministic nonperiodic flow,'' and is one of
% the most important scientific publications of the 20th 
% century.\footnote{The term ``chaos'' came later, coined by Jim
% Yorke in 1975.}
% The fact that it took until after 1963 for the phenomenon to be
% widely recognized can be attributed to the invention of computers.
% Chaos is not an unusual behavior at all in nonlinear systems of
% ODE\kern .5pt s, but it could not be easily seen in the closed form solutions
% that were available in the days of hand calculation.
% </latex>

%%
% <latex>
% The historical side of our next example is particularly
% interesting.  Among the most important ODE\kern .5pt s of all scientifically
% are the equations of the
% {\boldmath\bf $n$-body problem,} by
% which we mean the equations that govern the
% trajectories of $n$ point masses attracting
% each other gravitationally and moving according to Newton's laws.
% In its simplest form this is a set of $3n$ ODE\kern .5pt s of
% second order: for each body, there are three spatial coordinates,
% and the differential equations express the fact that the
% acceleration of each body is equal to
% the sum of the inverse-square attractions
% to the others.
% For problems confined to a plane, the count reduces to 
% $2n$ variables, an $x$ and a $y$ coordinate for each body.  We are
% going to look at such a problem, and for convenience we will employ 
% our usual trick of encoding $x$ and $y$ as a complex variable 
% $z=x+iy$, as first introduced
% in Chapter~3.  This brings us to a second-order system of
% ODE\kern .5pt s in $n$ complex unknowns:
% $$ z_j'' = - \sum_{i\ne j} {z_j^{} - z_i^{}\over | z_j-z_i|^3},
% \quad 1 \le j \le n.  \eqno (13.5) $$
% The analytic solution to the 2-body problem goes back to Kepler
% and Newton, and it involves stable ellipses (as in Exercise 4.1).
% But there is no analytic solution for the 3-body problem, let alone
% for the $n$-body problem with $n>3$.  
% </latex>

%%
% <latex>
% For a long time, mathematicians and physicists
% have asked, is our solar system stable?  Might
% its regular orbits one day break down?  This problem proved 
% theoretically intractable
% but has led to a great deal of important mathematics along the way.
% Even today, after extensive computational investigations,
% there is some controversy
% over whether or not the solar system
% is stable or chaotic, but the difficulty of this question is a
% result of a special circumstance:
% the mass of the sun is much greater than that
% of the planets, and acts as a strong regulating influence.
% The sun's mass is so dominant that to first approximation,
% each planet's orbit is just the solution of a 2-body problem involving
% itself and the sun.
% By contrast, if you look at the
% trajectories of $n$ attracting bodies of equal masses
% in the absence of a heavy sun to act as a policeman, chaos is the rule.
% Before computers, nobody could look at such trajectories.
% </latex>

%%
% The idealized configuration we shall examine
% consists of three planets in a plane with no sun,
% initially positioned with zero velocity at the
% vertices of a 3-4-5 right triangle.
% The following commands compute the orbits up to $t=150$.
u0 = 0; v0 = 3; w0 = 4i; N = chebop(0,150);
N.op = @(t,u,v,w) [ ...
    diff(u,2) + (u-v)/abs(u-v)^3 + (u-w)/abs(u-w)^3; ...
    diff(v,2) + (v-u)/abs(v-u)^3 + (v-w)/abs(v-w)^3; ...
    diff(w,2) + (w-u)/abs(w-u)^3 + (w-v)/abs(w-v)^3];
N.lbc = @(u,v,w) [u-u0; v-v0; w-w0; diff(u); diff(v); diff(w)];
[u,v,w] = N\0;

%%
% <latex>
% \noindent
% To begin with we plot the orbits just up to $t=20$.
% </latex>
plot(u,'interval',[0 20]), hold on, plot(v,'r','interval',[0 20])
plot(w,CO,ivpnl,'interval',[0 20]), hold off
title('Fig.~13.5.~~Orbits of three planets to $t=20$',FS,11)
axis equal, ylim([-2.5 5])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% In this image the red, green,
% and blue curves show paths of the three planets
% as~$t$ increases.  At first they are well separated,
% but soon two close red-blue flybys
% are observed.  Here are the distance and time of the closest approach:
% </latex>
[closest_distance,closest_time] = min(abs(u{0,20}-v{0,20}))

%%
% <latex>
% Plotting the orbits up to $t=50$ gives a picture hard
% to interpret in detail.  It is clear, however, that 
% the planets are swinging around each other in an
% effectively random fashion.  
% </latex>
plot(u,'interval',[0 50]), hold on, plot(v,'r','interval',[0 50])
plot(w,CO,ivpnl,'interval',[0 50]), hold off
title('Fig.~13.6.~~Orbits continued to $t=50$',FS,11)
axis equal, ylim([-2.5 5])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% One might expect this behavior
% to continue forever, but in fact, it doesn't.
% Around $t=86$, the three-body system breaks up (``self-ionizes''),
% with the red planet shooting off to the northwest and
% the green and blue ones spiraling off together to the southeast.
% This curious effect is
% revealed in the plot below of the orbits up to $t=100$, looking like
% ribbons on a birthday present.
% Thus the 3-body problem for an initial 3-4-5 triangle turns
% out to exemplify what is known as {\em transient chaos.}
% </latex>
plot(u,'interval',[0 100]), hold on
plot(v,'r','interval',[0 100]), plot(w,CO,ivpnl,'interval',[0 100])
title(['Fig.~13.7.~~Orbits to $t=100$, ' ...
       'with self-ionization near $t=86$'],FS,11)
axis equal, ylim([-2.5 5]), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here, for example, is the time at which the red planet
% passes through coordinate $x=-2$:
% </latex>
t_escape = roots(real(v)+2)

%%
% Just as we did with the Lorenz equations, let us examine
% the effect of a small perturbation in one component of
% the initial data.  
% Changing the initial coordinate $w(0)$ from $4i$ to
% $3.9999i$ gives the following orbit up to $t=100$,
% which is similar to the former one for a while but then
% begins to differ completely.
% This time there is no self-ionization.
N.lbc = @(u,v,w) [u-u0; v-v0; w-3.9999i; diff(u); diff(v); diff(w)];
[u2,v2,w2] = N\0; plot(u2,'interval',[0 100]), hold on
plot(v2,'r','interval',[0 100])
plot(w2,CO,ivpnl,'interval',[0 100])
hold off, axis equal, ylim([-2.5 5])
title('Fig.~13.8.~~Slightly perturbed orbits to $t=100$',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% A plot shows $\|\tilde {\bf y}(t)-{\bf y}(t)\|$ as a function of $t$
% and the estimated Lyapunov exponent $\lambda\approx 0.15$.
% </latex>
tt = linspace(0,150,400);
err = sqrt(abs(u2(tt)-u(tt)).^2 + abs(v2(tt)-v(tt)).^2 + ...
      abs(w2(tt)-w(tt)).^2);
semilogy(tt,err,'k'), ylim([1e-5 1e3]), ii = find(tt<=70);
title(['Fig.~13.9.~~Exponential growth of ' ...
       ' perturbations to scale $O(1)$'],FS,11)
c = polyfit(tt(ii),log(err(ii)),1); e0 = c(2); lambda = c(1)
hold on, semilogy([0 80],exp(e0+lambda*[0 80]),'--r',LW,1.5)
hold off, xlabel('$t$',FS,9,IN,LT)
ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,9,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Note the transition from the chaotic exponential
% phase to smooth algebraic growth at a rate $O(t)$
% after the bodies have separated.
% </latex>

%%
% <latex>
% Our third example of a chaotic system is the 
% {\bf R\"ossler equations,}
% $$ u' = -v-w, ~~ v' = u+\textstyle{1\over 5} v, ~~
% w' = {1\over 5} + w(u-c), \eqno (13.6) $$
% where the constant $1/5$ has been fixed
% arbitrarily and $c$ is a parameter.
% This system is even simpler than the Lorenz
% equations (13.1) in that only one of the three equations is nonlinear.
% The R\"ossler equations illustrate
% the phenomenon of {\em period doubling} as a route to chaos.
% First, here we solve (13.6) for $t\in [\kern .3pt 0,300\kern .3pt]$
% with $c=2$ with initial conditions
% $u(0)=2$ and $v(0) = w(0) = 0$. 
% The image on the left shows the trajectory in the three-dimensional
% $u$-$v$-$w$ phase space, which settles down to a regular
% oscillation, a limit cycle.  The image on the right shows
% the same trajectory projected onto the $u$--$v$ plane, just
% the part for $t>200$.  The initial transient
% has died away, and thus all one sees is the 2D projection
% of the limit cycle.
% </latex>
N = chebop(0,300); N.lbc = [2; 0; 0];
N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-2)];
[u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,LW,.3,CO,ivpnl,LW,.5)
xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10)
subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.7)
axis([-10 10 -10 10]), axis square
text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10)
title(['Fig.~13.10.~~R\"ossler system (13.6), $c=2$ ' ...
   '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% As the parameter $c$ is increased, a change takes place, a bifurcation
% (see Chapter~17).
% At around $c=2.8$, the limit cycle undergoes {\bf period doubling},
% with the trajectory unfolding into a double loop.
% The image for $c=3.5$ shows this ``period 2'' solution.
% </latex>
N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-3.5)];
[u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.3)
xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10)
subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.6)
axis([-10 10 -10 10]), axis square
text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10)
title(['Fig.~13.11.~~R\"ossler system (13.6), $c=3.5$ ' ...
   '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% As $c$ increases further, more bifurcations take place.
% A second period doubling to a ``period 4''
% solution occurs around $c=3.7$, and here is
% a pair of images for $c=4$.
% </latex>
N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-4)];
[u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.3)
xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10)
subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.5)
axis([-10 10 -10 10]), axis square
text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10)
title(['Fig.~13.12.~~R\"ossler system (13.6), $c=4$ ' ...
   '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The phenomenon of period doubling was made famous in the 1970s by
% Mitchell Feigenbaum, who showed that it is
% a route to chaos in many problems.  As $c$ increases further (the
% next transition is around $c=4.1$),
% the period doubles again and again,
% infinitely often.  Each time interval from one doubling to
% the next is asymptotically $4.6692\dots$ times shorter than the last; this
% is {\em Feigenbaum's constant.}  Finally, for $c$ greater than
% about $4.2$,
% the system is chaotic.  The following images for $c=5$ show
% the chaotic regime, with the orbits settling down not to a limit
% cycle but to a strange attractor.
% </latex>
N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-5)];
[u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.2)
subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.4)
axis([-10 10 -10 10]), axis square
text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10)
title(['Fig.~13.13.~~R\"ossler system (13.6), $c=5$ ' ...
   '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Our final example of chaos is governed by a
% nonautonomous equation due to Moon and Holmes,
% a chaotic {\em nonlinear forced oscillator}.
% The equation is
% $$ y'' + \textstyle{1\over 4} y' - y + y^3 = 0.4 \cos(t).
% \eqno (13.7) $$
% As always, the constants here are somewhat arbitrary and
% could be adjusted.
% This is just a second order ODE, so in the absence of the
% nonautonomous forcing function, the trajectories could not
% be chaotic.  The forcing function, however, makes this
% system equivalent to a first-order autonomous system in three
% variables $y$, $y'$, and $t$ (see footnote 3 of
% Chapter~10).  Thus the dimension is
% great enough for chaos to be a possibility, and a computation
% confirms its appearance.
% </latex>
t = chebfun('t',[0 300]); N = chebop(0,300);
N.op = @(t,y) diff(y,2) + .25*diff(y)-y+y^3;
N.lbc = [0;0]; rhs = 0.4*cos(t);
y = N\rhs; clf, plot(y,CO,ivpnl), ylim([-3 3])
title('Fig.~13.14.~~Nonlinear forced oscillator (13.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% We have not defined chaos!  This is because there is
% no universally accepted definition, though the features of sensitive
% dependence and strange attractors are important.  Note
% that exponential divergence of trajectories alone cannot
% be enough, since the solutions of $y' = y$ diverge at
% the rate $e^{\kern .3pt t}$ but nobody would call this equation chaotic.  
% One needs the exponential divergence to be combined with
% some kind of global boundedness associated with nonlinearity,
% and this is where the strange attractors come in.

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: chaos in a food web}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% In Chapter 10 we saw that a simple model of the interactions between
% rabbits and foxes led to limit cycles in their populations. If we
% introduce a third species into the fray, this can
% transform the dynamics completely.  This possibility was discovered
% by Bob May (later Lord May)
% in the 1970s, a story made famous in the 
% book {\em Chaos:\ Making a New Science} by James Gleick.
% </latex>

%%
% <latex>
% Rabbits need to eat too! Suppose in addition
% to the populations of rabbits, $u(t)$,
% and foxes, $v(t)$, we consider the population of carrots, $c(t)$. We
% suppose that, left on their own, these tend to grow logistically
% (see Exercises~3.15 and~3.16). Now
% rabbits consume carrots, and foxes consume rabbits, and neither animal
% species can be sustained without food. A reasonable model of the system is
% $$ c'=c(1-c)-f_1^{}(c)u, ~~ u'=f_1^{}(c)u - f_2^{}(u)v- d_1^{}u, ~~ v'
% = f_2^{}(u)v - d_2^{} v, \eqno (13.8) $$
% with $f_i^{}(z)=a_i^{} z/(1+b_i^{}z)$ for $i=1,2$. The consumption interactions
% are little different from those in the Lotka--Volterra equations.
% They include a
% saturation effect for the consumer species, accounting for extra
% competition as populations grow large. This model was proposed by
% Hastings and Powell in ``Chaos in a three-species food chain,''
% {\em Ecology}, 1991, where the parameter choices were
% $a_1^{}=5$, $a_2^{}=0.1$, $b_2^{}=2$, $d_1^{}=0.4$, $d_2^{}=0.01$.  The parameter
% $b_1^{}$ was allowed to vary. 
% </latex>

%%
% Increasing $b_1^{}$ increases the effect of competition between rabbits.
% Here is what happens to the rabbits with $b_1^{}=2.5$, 
% a nonchaotic regime.
a1 = 5; a2 = 0.1; b2 = 2; d1 = 0.4; d2 = 0.01;
f1 = @(z,b1) a1*z./(1+b1*z);  f2 = @(z) a2*z./(1+b2*z);
N = chebop(0,3000); b1 = 2.5;
N.op = @(t,c,u,v) [ diff(c)-(c*(1-c)-f1(c,b1)*u); 
                    diff(u) - (f1(c,b1)*u-f2(u)*v-d1*u); 
                    diff(v) - (f2(u)*v-d2*v) ];
N.lbc = [.4;1;9]; [c,u,v] = N\0; plot(u{2000,3000},CO,ivpnl)
title(['Fig.~13.15.~~Nonchaotic rabbit population, ' ...
       '$b_1=2.5$ in (13.8)'],FS,11,IN,LT)
xlabel('$t$',FS,10,IN,LT), ylabel('rabbits, $u(t)$',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% At first glance this may look somewhat irregular. But in phase space we
% see an ordinary limit cycle of doubled period, just
% as in Figure 13.11.
% </latex>
plot3(c{2000,3000},u{2000,3000},v{2000,3000},LW,.7,CO,ivpnl), view(63,68)
xlabel('carrots, $c(t)$',FS,9), ylabel('rabbits, $u(t)$',FS,9)
zlabel('foxes, $v(t)$',FS,9)
title('Fig.~13.16.~~In phase space, a limit cycle',FS,11)
set(gcf,'position',[380 220 540 300])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% If we set $b_1^{}=3.5$, increasing 
% the competition between rabbits, the picture changes fundamentally.
% </latex>
b1 = 3.5;
N.op = @(t,c,u,v) [ diff(c)-(c*(1-c)-f1(c,b1)*u); 
                    diff(u) - (f1(c,b1)*u-f2(u)*v-d1*u); 
                    diff(v) - (f2(u)*v-d2*v) ];
[c,u,v] = N\0; 
plot3(c{2000,3000},u{2000,3000},v{2000,3000},LW,.4,CO,ivpnl), view(63,68)
xlabel('carrots, $c(t)$',FS,9), ylabel('rabbits, $u(t)$',FS,9)
zlabel('foxes, $v(t)$',FS,9)
title(['Fig.~13.17.~~Chaotic rabbit population, ' ...
       '$b=3.5$ in (13.8)'],FS,11,IN,LT)
set(gcf,'position',[380 220 540 300])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Looking again at the rabbit population, we see repeated but aperiodic
% crashes in the population, separated by a variable number and size of
% shorter cycles of boom and bust.
% </latex>
close all, plot(u{2000,3000},CO,ivpnl)
title('Fig.~13.18.~~Chaotic rabbit population',FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('rabbits, $u(t)$',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% Obviously real food webs can be far more complex than this. But this
% example shows that the change
% from two species to three can have a quantum effect, and for reasons
% that are mathematical, not biological. 

%%
% <latex>
% \smallskip
% {\sc History.}  The roots of chaos are often traced to
% Poincar\'e and Hadamard in the 1890s, long before the days of
% computers, but the subject did
% not become widely studied until the work of Lorenz, May, Yorke,
% Feigenbaum
% and others in the 1960s and 1970s.  Chaos became known
% to the public with James Gleick's 1987 book {\em Chaos:\ Making
% a New Science} mentioned earlier and the 1993 movie {\em Jurassic Park}.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.} 
% One of the best mathematical textbooks ever written is
% S. H. Strogatz, {\em Nonlinear Dynamics and Chaos: With Applications
% to Physics, Biology, Chemistry, and Engineering,} Westview
% Press, 2015 (first published in 1994).  Strogatz combines marvelous
% clarity with a precise presentation of the mathematics
% of this fascinating subject.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 13.} 
% Many nonlinear ODE\kern .5pt s with three or more variables are chaotic, though
% this fact was not widely recognized
% until the last third of the
% twentieth century.  Chaos is characterized by the 
% property that perturbations may grow
% exponentially with time, yet global orbits remain bounded.
% The maximal rate of exponential growth is known as the Lyapunov exponent.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em \underline{Exercise $13.1$}.  Lucky Lyapunov exponents.}
% It was mentioned in the text that the estimate $0.9052$ of the
% Lyapunov exponent was fortuitously close to the true value.  Confirm
% this by computing the name number based on data for
% $t\in [\kern .3pt 0,14]$ and $t\in [\kern .3pt 1,15]$ instead of
% $t\in [\kern .3pt 0,15]$.
% \par
% {\em \underline{Exercise $13.2$}.  Smaller perturbation of the Lorenz
% equations.}
% In the text we perturbed the Lorenz equations by changing the value
% of $w(0)$ from $20$ to $20+ 10^{-5}$.  {\em (a)}
% Reproduce Figures 13.2 and 13.3
% with the smaller perturbation $20 + 10^{-9}$ and comment on the results.
% {\em (\kern .7pt b)} Also repeat the finite-difference estimates of the derivatives
% of $v(5)$, $v(10)$, and $v(30)$ with respect to $w(0)$, and
% comment on the results.
% \par
% {\em \underline{Exercise $13.3$}.  Alternative choices of 
% the Lorenz coefficient\/ $28$.}
% In the Lorenz equations, let $r$ denote the parameter that traditionally
% takes the value $28$.  Starting from our usual initial conditions, make plots of
% $u(t)$ against $w(t)$ as in Figure 13.1
% for $t\in [\kern .3pt 0, 100\kern .3pt]$ with
% $r = 20, 22, 24$; also make plots in each case of $u(t)$ against
% $t$ as in Figure 13.2.  Which case seems to be chaotic?
% Which one gives the clearest example of transient chaos?
% \par
% {\em \underline{Exercise $13.4$}.  Lorenz equations with a breeze.}
% In the Lorenz equations (13.1), let the first equation be
% changed to $u' = 10(v-u) - a$, where $a$ is a parameter.
% As in the last exercise, make plots of
% $u(t)$ against $w(t)$ and of $u(t)$ against $t$
% for $t\in [\kern .3pt 0, 30\kern .3pt]$ with
% $a = 20, 25, 30$.  Comment on the solutions.
% \par
% {\em \underline{Exercise $13.5$}.  Two electrons and a nucleus.}
% Consider the highly idealized problem of two electrons
% of mass 1 and charge $-1$ orbiting a nucleus of
% charge $+2$ fixed at the origin of the $x$-$y$ plane.  Let $z(t)$ and
% $w(t)$ be the positions of the electrons represented with the
% usual complex variable trick.  Consider trajectories for
% $t\in [\kern .3pt 0,20]$ starting from initial conditions
% $z = i$, $z'=1$, $w=-i$, $w' = a$.
% {\em (a)} Write down the ODE governing the evolution of $w(z)$
% and $z(t)$
% assuming an inverse-square electrostatic force law with constant~$1$.
% The electrons repel each other while being attracted to the nucleus.
% {\em (\kern .7pt b)} Plot the two trajectories in the case $a=1$.  The
% configuration is symmetric, and the symmetry should be
% preserved for all $t$.
% {\em (c)} Now make similar plots for $a = 0.5, 0.6, 0.7, 0.8, 0.9,$
% and $0.99$.  Discuss the results.
% \par
% {\em \underline{Exercise $13.6$}.  Chaos and cellular automata.}
% Discrete processes, in which no real numbers are involved, can also
% exhibit chaotic properties.  Illustrate this by simulating the
% ``Rule 30'' cellular automaton described in Stephen Wolfram's 2002 book
% {\em A New Kind of Science.}\footnote{On p.~27 of the original
% edition Wolfram calls this ``probably
% the single most surprising scientific discovery I have ever
% made.''}  To do this, let $n$ be a
% positive even integer and initialize an
% $n\times n$ matrix ${\bf A}$ to zero apart from the value 1 in the middle
% of the top row, $a_{1,n/2}^{} = 1$.  Now sweep from
% rows $2$ to $n$, updating each entry according to the following
% rule: $a_{ij}^{}$ is changed to the the value $1$ if the three
% entries above it, $a_{i-1,j-1:j+1}$, correspond to the binary
% representation of one of the numbers $1,2,3,4$ rather than
% $0,5,6,7$.  (The numbers in columns $1$ and $n$ can be
% left unchanged.)  Take $n=200$ and use the Matlab {\tt spy} command to
% plot the matrix at the end.  Note that the structure is
% completely deterministic yet has apparently
% random properties.  For another view of the randomness, plot
% $a_{i,n/2}^{}$ as a function of $i$.
% \par
% {\em \underline{Exercise $13.7$}.  Random Fibonacci sequence.}
% Here, on the other hand, is an example of a discrete process with
% true randomness.
% Suppose $x_0^{} = x_1^{} = 1$ and for $k\ge 1,$
% $x_{k+1}^{} = x_{k-1}^{}+x_k^{}$.  Then it is well known
% that $x_k^{}$ grows asymptotically at a rate determined by
% $\phi^k$, where $\phi$ is the golden ratio $(1+\sqrt 5\kern .5pt)/2
% \approx 1.618.$  Consider the 
% {\em random Fibonacci sequence} defined by 
% $x_{k+1}^{} = \pm x_{k-1}^{} \pm x_k^{}$, where at each step each sign is 
% independently $+$ or $-$ with probability $0.5$.
% On a semilogy scale, make plots of $|x_k^{}|$ vs.\ $k$
% up to the maximum values
% $k = 100$ and $5000$.  Numerically estimate the {\em Lyapunov constant} for
% this process, that is, the constant $C$ such that with probability~1,
% the sequence grows in absolute value at a rate $C^k$ as $k\to\infty$.
% (If you are curious to learn more, see D. Viswanath, 
% {\em Mathematics of Computation\/} 69 (2000), pp.~1131--1155.)
% </latex>
