%% 15. Stable and unstable fixed points

%%
% <latex>
% \setcounter{page}{185}
% \parindent=.6cm
% \parskip=-1em
% </latex>

%%
% The last chapter dealt with linearization of ODE\kern .5pt s.  One of the
% most important applications of linearization is the analysis
% of fixed points.

%%
% <latex>
% Consider an autonomous equation
% $$ {\bf y}'(t) = {\bf f}({\bf y}(t)). \eqno (15.1) $$
% As defined in Chapters~9 and~10, a {\bf fixed point} of
% (15.1) is a vector ${\bf y}_*^{}$ such that ${\bf f}({\bf y}_*^{}) = 
% {\bf 0}$.   We saw in those chapters that a fruitful way of
% understanding the behavior of an autonomous ODE is
% to begin by examining its fixed points in the phase plane
% (if $n=2$) or more generally in phase space.
% If ${\bf y}_*^{}$ is a fixed point, then the term ${\bf f}_*^{}$
% in Theorem 14.4 vanishes.  Equation (14.17) becomes
% $$ \delta {\bf y}'(t) = {\bf J}_*^{} \delta {\bf y}(t) +
% O(\|\delta {\bf y}(t)\|^2), \eqno (15.2) $$
% with $\delta{\bf y}(t) = {\bf y}(t) - {\bf y}_*^{}$ as before.
% Near a fixed point, therefore, an ODE behaves approximately like the
% equation ${\bf y}' = {\bf A}{\bf y}$ with which
% we began the last chapter --- not just
% affine, but linear, with ${\bf A}$ being the Jacobian matrix ${\bf J}_*^{}$. 
% </latex> 

%%
% <latex>
% To illustrate the structure of some linearizations at
% fixed points in the phase plane,
% here are two plots showing solutions of
% (15.2) (without the ``O'' term)
% corresponding to the diagonal matrices
% $$ {\bf J}_*^{} = \pmatrix{-1 & \hphantom{-}0 \cr \hphantom{-}0 & -1} ,~
% \pmatrix{-2 & \hphantom{-}0 \cr \hphantom{-}0 & -1} . 
% \eqno (15.3) $$
% Each plot shows trajectories emanating from 
% 16 equally spaced initial points on the unit circle, which
% is drawn in black.
% </latex>
ODEformats, th = (pi/8)*(1:16) + .0001; u0 = cos(th); v0 = sin(th);
L = chebop(0,1.4);
op = @(J) @(t,u,v) [diff(u)-J(1,1)*u-J(1,2)*v; ...
                    diff(v)-J(2,1)*u-J(2,2)*v];
subplot(1,2,1), plot(0,0,'.k',MS,6), hold on
c = chebfun('exp(1i*pi*x)'); plot(c,'k',LW,.7)
J = [-1 0; 0 -1]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-1.25 1.25 -1.1 1.4]), set(gca,XT,-1:1,YT,-1:1)
pp = get(gca,'pos'); pp(1) = pp(1)+.015; set(gca,'pos',pp)
text(0,1.15,'stable node',FS,9,HA,CT)
axis square, hold off
subplot(1,2,2), plot(0,0,'.k',MS,6), hold on, plot(c,'k',LW,.7)
J = [-2 0; 0 -1]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-1.25 1.25 -1.1 1.4]), set(gca,XT,-1:1,YT,-1:1)
pp = get(gca,'pos'); pp(1) = pp(1)-.015; set(gca,'pos',pp)
text(0,1.15,'stable node',FS,9,HA,CT)
title(['Fig.~15.1.~~Fixed points for the matrices ' ...
       '(15.3)\kern -.40in'],FS,11,HA,RT)
axis square, hold off
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% In both figures, all trajectories are converging to the origin.
% This is because both eigenvalues of ${\bf J}_*^{}$, which in this
% case are simply the diagonal entries, are negative.
% If all the eigenvalues of the Jacobian at a fixed point are
% negative, or more generally satisfy
% $\hbox{Re\kern 1pt}\lambda<0$ and thus lie in the open
% left half of the complex plane, then ${\bf y}_*^{}$ is called
% a {\bf sink}.  This implies
% that all orbits starting sufficiently close to ${\bf y}_*^{}$ converge to 
% ${\bf y}_*^{}$ at an exponential rate.
% The image on the left above may be the first one that
% comes to mind when one thinks about sinks, but the image on
% the right is more typical: most trajectories approach a
% sink along special directions corresponding to the eigenvectors
% associated with eigenvalues of least negative real part.
% In this example the eigenvectors are $(1,0)^T$, the direction
% of exponential decay at the rate $e^{-2t}\kern 1pt$, and
% $(0,1)^T$, the direction of slower decay at the rate
% $e^{-t}\kern 1pt$.  As $t$ increases, the component in
% the $(1,0)^T$ direction becomes negligible compared with the component
% in the $(0,1)^T$ direction, so trajectories approach the
% origin along the latter, vertical axis.
% </latex>

%%
% But still we have diagonal matrices, so these
% pictures do not show the general behavior.
% To illustrate some further possibilities, here are figures for
% two nondiagonal matrices,
% $$ {\bf J}_*^{} = \pmatrix{-1 & \hphantom{-}1 \cr -1 & -1} ,~
% \pmatrix{-2 & 2.5 \cr \hphantom{-}0 & -1} , \eqno (15.4) $$
% whose eigenvalues are $\{-1+i, -1-i\}$ and
% $\{-2, -1\}$, respectively.
% Again the eigenvalues are in the left half-plane, so
% $\bf 0$ is again a sink, but the images are quite different,
% showing a combination of rotation mixed with decay.
subplot(1,2,1), plot(0,0,'.k',MS,6), hold on, plot(c,'k',LW,.7)
J = [-1 1; -1 -1]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; 
  [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-1.35 1.35 -1.2 1.5]), set(gca,XT,-1:1,YT,-1:1)
pp = get(gca,'pos'); pp(1) = pp(1)+.015; set(gca,'pos',pp)
text(0,1.2,'stable spiral',FS,9,HA,CT)
axis square, hold off
subplot(1,2,2), plot(0,0,'.k',MS,6), hold on, plot(c,'k',LW,.7)
J = [-1 2.5; 0 -1/2]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; 
  [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-1.35 1.35 -1.2 1.5]), set(gca,XT,-1:1,YT,-1:1)
pp = get(gca,'pos'); pp(1) = pp(1)-.015; set(gca,'pos',pp)
text(0,1.2,'stable node',FS,9,HA,CT)
axis square, hold off
title(['Fig.~15.2.~~Fixed points for the matrices ' ...
       '(15.4)\kern -.40in'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Note that in the left image, there is no special direction along
% which trajectories eventually straighten out.
% Such a direction would correspond to an eigenvector,
% but for this matrix, the eigenvectors are complex.
% The eigenvalues
% are $-1+i$ and $-1-i$, both in the left half-plane, which
% explains the decay toward the origin.  In the right image,
% there are real eigenvectors again but
% they are far from orthogonal.  The effect of this is
% that although eventually all the
% trajectories decay to the origin, some of them grow
% for a while before decaying.  This phenomenon is known
% as {\em transient growth.}
% </latex>

%%
% <latex>
% All these plots correspond to the same simple case of
% a $2$-variable problem with a sink.  This is only
% the beginning of the many configurations that can arise
% in linearized analysis of fixed points. If all the
% eigenvalues are in the open right half-plane, that is, with
% $\hbox{Re\kern 1pt}\lambda>0$, then ${\bf y}_*^{}$ is a {\bf source}
% and the arrows are reversed.  If some eigenvalues are in
% the left half-plane and the others are in the right half-plane,
% then ${\bf y}_*^{}$ is a {\bf saddle point}. 
% Here are examples of saddle points corresponding to
% the matrices
% $$ {\bf J}_*^{} = \pmatrix{-1 & 0 \cr \phantom{-}0 & 1} ,~
% \pmatrix{-1 & 2 \cr \phantom{-}0 & 1} .
% \eqno (15.5) $$
% The eigenvalues of both matrices are $\{-1,1\}$.
% </latex>
subplot(1,2,1), plot(0,0,'.k',MS,6), hold on, plot(c,'k',LW,.7)
L = chebop(0,0.9); J = [-1 0; 0 1]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; 
  [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-3 3 -2.7 3.3]), set(gca,XT,-3:3,YT,-3:3)
pp = get(gca,'pos'); pp(1) = pp(1)+.015; set(gca,'pos',pp)
text(0,2.8,'saddle point',FS,9,HA,CT)
axis square, hold off
subplot(1,2,2), plot(0,0,'.k',MS,6), hold on, plot(c,'k',LW,.7)
J = [-1 2; 0 1]; L.op = op(J);
for k = 1:16
  L.lbc = [u0(k); v0(k)]; 
  [u,v] = L\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis([-3 3 -2.7 3.3]), set(gca,XT,-3:3,YT,-3:3)
pp = get(gca,'pos'); pp(1) = pp(1)-.015; set(gca,'pos',pp)
text(0,2.8,'saddle point',FS,9,HA,CT)
axis square, hold off
title(['Fig.~15.3.~~Fixed points for the matrices ' ...
       '(15.5)\kern -.40in'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% There is a general terminology and theory of behavior near
% fixed points that goes beyond these linearized approximations.
% We say that a fixed point of a system of the form (15.1)
% is {\bf Lyapunov stable} if for any neighborhood
% $V$ of~${\bf y}_*^{}$, there is a neighborhood
% $U\subseteq V$ such that every trajectory that starts in~$U$
% remains in $V$ for all $t$.
% This condition does not require decay, just boundedness.
% The fixed point is {\bf asymptotically stable}
% if in addition, $U$ can be chosen such that every trajectory
% that starts in $U$ converges to ${\bf y}_*^{}$ as
% $t\to\infty$.\footnote{Despite this careful terminology,
% most of the time we will be more casual and just say
% a fixed point is stable if all nearby trajectories converge
% to it.}
% If~${\bf y}_*$ is not Lyapunov stable, it is {\bf Lyapunov unstable}, implying
% that some (not necessarily all) trajectories starting near
% ${\bf y}_*$ diverge away as $t$ increases.
% The following theorem, which we give without proof, summarizes
% some of the relationships between these general notions
% and the eigenvalues of~${\bf J}_*^{}$.
% </latex>

%%
% <latex>
% \medskip
% {\em
% {\bf Theorem 15.1.  Stability and eigenvalues of the
% Jacobian (\textsf{\textbf{FlAsHI}}).}\ \ Let\/ ${\bf y}_*^{}$ be a fixed point of an
% autonomous ODE\/ $(15.1)$ where\/ ${\bf f}$ is twice differentiable
% at\/~${\bf y}_*^{}$, and let ${\bf J}_*^{}$ be the associated
% Jacobian matrix. 
% If all the eigenvalues $\lambda$ of\/ ${\bf J}_*^{}$ satisfy
% $\hbox{\rm Re\kern 1pt} \lambda < 0$, then\/ ${\bf y}_*^{}$ is 
% asymptotically stable,
% and if at least one of them 
% satisfies $\hbox{\rm Re\kern 1pt} \lambda > 0$,
% then\/ ${\bf y}_*^{}$ is Lyapunov unstable.
% }
% \medskip
% </latex>

%%
% Note that the theorem leaves open the situation in which
% all eigenvalues satisfy $\hbox{Re\kern 1pt} \lambda\le 0$ but not
% all satisfy $\hbox{Re\kern 1pt} \lambda<0$.  In this case
% ${\bf y}_*^{}$ is unstable if there
% is a defective multiple eigenvalue with $\hbox{Re\kern 1pt}\lambda=0$ (i.e., 
% associated with a Jordan block of size $\ge 2$).  If all
% eigenvalues with $\hbox{Re\kern 1pt}\lambda=0$ are nondefective,
% then linear analysis is not enough
% to determine stability or asymptotic stability; it depends
% on the higher-order nonlinear behavior of $\kern 1pt \bf f$.

%%
% <latex>
% Let us now consider three examples, systems of ODE\kern .5pt s
% we looked at in Chapter~10.  The first of these was
% the Lotka--Volterra equations (10.4)--(10.5),
% $$ u' = u - uv, \quad v' = -\textstyle{1\over 5}v + uv, $$
% whose behavior was plotted in Figs.~10.4--10.8 and again
% in Figs.~14.5--14.6.
% To calculate the fixed points we set $u'=0$, implying $u=0$ or $v=1$,
% and $v'=0$, implying $v=0$ or $u=1/5$.  So the fixed points
% are $(u,v) = (0,0)$ and $(1/5,1)$.  The Jacobian of the system at
% $(u,v)$ is
% $$ {\bf J} (u,v) = \pmatrix{1-v & -u \cr v & u -1/5}, $$
% implying
% $$ {\bf J} (0,0) = \pmatrix{1 & 0 \cr 0 & -1/5}, \quad
% {\bf J} (1/5,1) = \pmatrix{0 & -1/5 \cr 1 & 0}. $$
% The first matrix has eigenvalues $1$ and $-1/5$, so
% $(0,0)$ is a saddle point.  The image on the right
% eigenvalues $\pm i/\sqrt 5$, both imaginary,
% so $(1/5,1)$ is a center, a point of neutral stability.
% This explains the rotation of trajectories that is the conspicuous
% feature of Figs.~10.6--10.8.
% </latex>

%%
% <latex>
% The next example from Chapter 10 is the Lorenz equations (10.7),
% $$ u' = 10(v-u), \quad v' = u(28-w)-v, \quad w' = uv-(8/3)w, $$
% plotted in Figs.~10.9--10.11 and considered further in Chapters~13
% and~14.
% To find fixed points we calculate that $u'=0$ implies
% $v=u$, $v'=0$ then implies $u=v=0$ or $w=27$, and
% $w'=0$ implies $u=v=w=0$ or $u=v = \pm\sqrt{8\cdot 27/3} = \pm 6\sqrt 2$.
% So there are three fixed points, and they are
% $(u,v,w) = (0,0,0)$ and $(\pm 6\sqrt 2,\kern .8pt \pm 6\sqrt 2,\kern .8pt  27)$.
% The Jacobian of the system at $(u,v,w)$ is
% $$ {\bf J} (u,v,w) = \pmatrix{-10 & 10 & 0 \cr 28-w & -1 & -u \cr  v & u & -8/3}, $$
% implying
% $$ {\bf J} (0,0,0) = \pmatrix{-10 & 10 & 0 \cr 28 & -1 & 0 \cr  0 & 0 & -8/3} $$
% and
% $$ {\bf J} (\pm 6\sqrt 2,\kern .8pt \pm 6\sqrt 2,\kern .8pt 27) =
% \pmatrix{-10 & 10 & 0 \cr 1 & -1 & \mp 6\sqrt 2 \cr \pm 6\sqrt 2 &\
% \pm 6\sqrt 2 & -8/3}. $$
% The eigenvalues of $J(0,0,0)$ are about $-22.8$, $-2.7$, and $11.8$, so
% this is an unstable saddle point.  The other fixed points are
% the more interesting ones, with eigenvalues about $-13.9$ and
% $0.9 \pm 10.2i$.  Thus these fixed points are unstable too, mildly so, but
% because of the large imaginary term they have a big rotational component,
% as we know well from the Lorenz trajectories of Figs.~10.11 and~13.1.
% </latex>

%%
% <latex>
% The other main example of Chapter 10 is the system
% of equations (10.10) of an SIR model
% from epidemiology, which with the parameter
% choices $\beta =2$ and $\gamma=1$ takes the form
% $$ S' = -2 SI, \quad I' = (2S-1)I, \quad R' = I. $$
% Behaviors were explored in Figs.~10.13--10.14. 
% For a fixed point analysis we first set $R'=0$, implying $I=0$.  This then
% implies $R'=S'=0$, so we see that every choice of $R$ and $S$ gives a 
% fixed point, so long as $I=0$.  The Jacobian matrix is
% $$ {\bf J} (S,I,R) = \pmatrix{-2I & -2S & 0 \cr 2I & 2S-1 & 0 \cr  0 & 1 & 0}, $$
% which at a fixed point becomes
% $$ {\bf J} (S,I,R) = \pmatrix{0 & -2S & 0 \cr 0 & 2S-1 & 0 \cr  0 & 1 & 0}, $$
% Since this matrix is block lower-triangular, with an upper-left $2\times 2$
% block that is itself upper-triangular, we see that the eigenvalues
% are $0$, $0$, and $2S-1$.  The latter number has
% immediate significance: if $S>1/2$, the matrix has a positive eigenvalue
% and the system is unstable, ready to begin an epidemic as
% soon as any patient gets infected.
% </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: transition to turbulence in a pipe}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% Of all the fixed points in the mathematical sciences, perhaps
% none has received more attention, or caused
% more confusion, than laminar fluid flow in a
% pipe.\footnote{There are competitors, though the equations
% are not so clear-cut, in climate science.  Whenever you hear the
% phrase ``tipping point,'' you can be sure that there is a question of
% stability of a fixed point at hand.  One tipping point of great
% concern involves the melting of the ice in the Arctic: as ice
% melts, the earth reflects less light out to space, and the
% ice melts faster.  See D. Paillard, ``The
% timing of Pleistocene glaciations from a simple
% multiple-state climate model,'' {\em Nature} 391 (1998), pp.~378--381.}
% (Laminar means smooth and steady.)
% This discussion is adapted from
% Trefethen, Trefethen, Reddy, and Driscoll,
% ``Hydrodynamic stability without eigenvalues,'' {\em Science,} 1993.
% </latex>

%%
% <latex>
% The problem was made famous by Osborne Reynolds in 1883.
% Imagine a long circular pipe with a fluid such as water flowing
% through it.  The flow is governed by the set of time-dependent PDE\kern .5pt s 
% known as the Navier--Stokes equations, determining the evolution of the
% velocity field ${\bf v}({\bf x},t)$, and an analytical solution corresponding
% to laminar flow can be written down easily: $\bf v$ is a time-independent field
% pointing along the pipe, corresponding to a steady
% flow with a velocity that is maximal at the centerline and
% decreases quadratically to zero at the wall.  This solution
% is a fixed point of the equations, and it is mathematically stable.
% We have not discussed PDE\kern .5pt s and their stability, but just
% as for ODE\kern .5pt s, the idea is that
% any sufficiently small perturbation of the laminar flow velocity
% field in a pipe must eventually die away.
% So a mathematician would expect that laminar flow of water through a
% pipe should be possible at any speed.
% </latex>

%%
% <latex>
% The paradox is that in practice, this is not what is
% observed.  If the flow in a pipe 
% is fast enough, it is invariably not laminar but turbulent --- complicated,
% apparently chaotic, highly time-dependent.
% Clearly the mathematical solution, laminar flow, has something wrong
% with it in the laboratory.
% What is going on?  How can a flow that is stable
% mathematically be unstable in practice?
% </latex>

%%
% <latex>
% The explanation is that although 
% sufficiently small perturbations of the laminar velocity flow field
% must eventually decay, the threshold that defines
% ``sufficiently small'' is too small
% to be counted upon in practice.
% The mathematics of these high-speed flow problems
% is such that the minimal amplitude of perturbations that do 
% {\em not\/} eventually decay is tiny.
% The slightest imperfection in the pipe or the
% smoothness of the inflow, or the slightest vibration of the laboratory,
% may be enough to kick the system into instability.
% Geometrically, we say that the {\em basin of attraction
% of the laminar state is very narrow.}  This makes
% the laminar state often effectively unobservable in practice.
% </latex>

%%
% <latex>
% A simple ODE model explains how an extremely narrow
% basin of attraction can come about in a set of equations
% that seems far from extreme.
% Let $R>0$ be a parameter, a
% caricature of the Reynolds number, the nondimensional
% centerline velocity of the laminar flow.  Let 
% $u$ and $v$ be two dependent variables, caricatures
% of the field of velocity perturbations of the laminar flow,
% satisfying the equations
% $$ u' = -R^{-1}u + v - v\sqrt{u^2+v^2}, \quad
% v' = -2R^{-1}v + u\sqrt{u^2+v^2}. \eqno (15.6) $$
% Rewriting (15.6) in matrix form reveals the structure more plainly:
% $$ \pmatrix{u\cr v}' =
% \pmatrix{-R^{-1} & 1 \cr 0 & -2R^{-1}} \pmatrix{u\cr v} 
% + \sqrt{u^2+v^2} \pmatrix{0 & -1 \cr 1 & \hphantom{-}0} \pmatrix{u\cr v}.
% \eqno (15.7) $$
% The term on the left is linear, and the term on the
% right is quadratic (if $u$ and $v$ are doubled, it multiplies
% by 4).
% As always, it is the linear term
% that governs behavior for sufficiently small $u$ and $v$.  Since
% the matrix is triangular, its eigenvalues are the diagonal entries,
% $-R^{-1}$ and $-2R^{-1}$, and since these numbers are
% negative, $(0,0)$ is a sink.
% </latex>

%%
% <latex>
% On the left below is an image of the linear part of the
% problem for $R=10$ showing this sink, very much like
% the second panel of Figure 15.2.  Something physically important
% is revealed in this image: a great deal of transient growth
% before the eventual decay.  (In the fluid mechanics problems
% this effect is sometimes called ``lift-up,'' in which vorticity aligned
% with a shear flow excites a growth in local velocity anomalies.)
% On the right is an image of the nonlinear problem for the same
% value $R=10$.   What happens here is that the quadratic term
% catches hold of the linear amplification and moves it onto
% a new track entirely.  Though the initial conditions in this
% nonlinear experiment
% start at a distance of only $0.02$ from the fixed point,
% most of the trajectories spiral up to size $O(1)$ rather than
% decaying to the center.  For this model, this is the caricature
% of transition to turbulence. 
% </latex>
th = (pi/8)*(1:16) + .0001; u0 = cos(th); v0 = sin(th);
clf, subplot(1,2,1), plot(0,0,'.k',MS,6), hold on, plot(.02*c,'k')
N = chebop(0,18); R = 10;
N.op = @(t,u,v) [diff(u) + u/R - v; diff(v) + 2*v/R];
for k = 1:16
  N.lbc = [.02*u0(k); .02*v0(k)];
  [u,v] = N\0; arrowplot(u,v,CO,ivp,MS,3)
end
axis(.07*[-1 1 -1 1]), set(gca,XT,.05*(-1:1),YT,.05*(-1:1))
text(-.054,-.044,'linearly stable',FS,8), axis square
text(-.054,-.0548,'fixed point...',FS,8)
pp = get(gca,'pos'); pp(1) = pp(1)+.015; set(gca,'pos',pp)
subplot(1,2,2), plot(0,0,'.k',MS,6), hold on, plot(.02*c,'k')
N = chebop(0,30); 
N.op = @(t,u,v) [diff(u) + u/R - v + v*sqrt(u^2+v^2)
                 diff(v) + 2*v/R - u*sqrt(u^2+v^2)];
for k = 1:16
  N.lbc = [.02*u0(k); .02*v0(k)];
  [u,v] = N\0; big = (norm([u(end) v(end)])>.02);
  if big, arrowplot(u,v,CO,'r',MS,3)
  else arrowplot(u,v,CO,ivpnl,MS,3), end
end
axis(.08*[-1 1 -1 1]), set(gca,XT,.05*(-1:1),YT,.05*(-1:1))
axis square, hold off
h1 = text(-.076,-.0405,'...but nonlinearly unstable to',FS,8);
h2 = text(-.066,-.052,'small finite perturbations.',FS,8);
h3 = text(-.066,-.0635,'Red trajectories grow to $O(1)$.',...
     FS,8,IN,LT);
pp = get(gca,'pos'); pp(1) = pp(1)-.015; set(gca,'pos',pp)
title(['Fig.~15.4.~~Stable node with narrow ' ...
       'basin of attraction\kern -.6in'],FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here we zoom out the second plot, showing where the red
% trajectories are heading.
% </latex>
axis(1.5*[-1 1 -1 1]), axis square
delete(h1), delete(h2), delete(h3)
text(0.04,-.3,'...with the',FS,8)
text(0.04,-.5,'nonlinear orbits',FS,8)
text(0.04,-.7,'shown in full.',FS,8,IN,LT)
set(gca,XT,-1:1,YT,-1:1)
title(['Fig.~15.5.~~Same dynamics shown on a ' ...
       'larger scale\kern -.5in'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Let us draw a plot to visualize this behavior
% in another way.  For the same equations (15.7) with
% $R=10$, Figure 15.6 shows $(u^2+v^2)^{1/2}$ as a function
% of\/ $t$ for six trajectories emanating from
% initial points $(u,v) = (0, v_0^{})$ with $v_0^{}
% = 0,001, 0.0025, 0.005, 0.01, 0.02, 0.04$.  The three
% initial conditions of lowest amplitude lead to trajectories
% that eventually decay to zero, but the other three, shown in
% red, increase to $O(1)$.  
% </latex>
clf, N = chebop(0,60);
N.op = @(t,u,v) [diff(u) + u/R - v + v*sqrt(u^2+v^2) 
                 diff(v) + 2*v/R - u*sqrt(u^2+v^2)];
for v0 = [.001 .0025 .005 .01 .02 .04]
  N.lbc = [0; v0]; [u,v] = N\0;
  big = (norm([u(end) v(end)])>.2);
  if big, semilogy(sqrt(u.^2+v.^2),'r')
  else semilogy(sqrt(u.^2+v.^2),CO,ivpnl), end, hold on
end
title(['Fig.~15.6.~~ODE model (15.7) ' ...
       'of transition to turbulence'],FS,11)
xlabel('t',FS,10,IN,LT), ylabel('energy',FS,9), ylim([1e-4 10])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Unlike turbulence,
% the model (15.7) is non-chaotic, and indeed it could not 
% possibly be chaotic since it is an autonomous first-order system with
% just two variables.  Similar models with three variables
% instead of two, however,
% combine the narrow basin of attraction of the
% laminar state with chaotic long-time trajectories.
% See Baggett, Driscoll, and Trefethen,
% ``A mostly linear model of transition to turbulence,''
% {\em Physics of Fluids,} 1995.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}  The general theory of stability of
% ODE\kern .5pt s was developed by Aleksandr Mikhailovich
% Lyapunov, a student of Chebyshev,
% who was one of the outstanding Russian mathematicians of the era
% before the 1917 revolution.  Lyapunov's interest in questions of stability
% began with fluid and solid mechanics and was set forth in
% his great work of 1892, {\em The General Problem of
% the Stability of Motion.}  His academic descendants include
% Smirnov, Sobolev, Kantorovich, Ladyzhenskaya, and other
% major figures of twentieth century Russian mathematics.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.} 
% The idea of a tipping point was made famous by 
% Malcolm Gladwell's bestseller
% {\em The Tipping Point: How Little Things Can Make a
% Big Difference,} which first appeared in 2000.
% Gladwell emphasizes tipping points related
% to the mathematics of epidemiology --- which, as he vividly shows,
% applies to many more areas than just epidemiology.
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 15.} The starting point of analysis of an
% autonomous system of ODE\kern .5pt s ${\bf y}' = {\bf f}({\bf y})$
% is determination and classification of its fixed points
% ${\bf y}_*^{}$, defined by the
% condition ${\bf f}({\bf y}_*^{}) = {\bf 0}$.
% For ${\bf y}\approx {\bf y}_*^{}$, the difference $\delta {\bf y}(t) =
% {\bf y}(t) - {\bf y}_*^{}$ evolves approximately according to the
% linear equation $(\delta {\bf y})' = {\bf J}_*^{} (\delta {\bf y})$,
% where ${\bf J}_*^{} = (\partial{\bf f}/\partial{\bf y})({\bf y}_*^{})$ is
% the Jacobian matrix at ${\bf y}_*^{}$.  A fixed point is stable if all
% the eigenvalues of\/ ${\bf J}_*^{}$ lie in the left half of the
% complex plane and unstable if at least one eigenvalue is in the
% right half-plane.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em Exercise $15.1$.  Unstable or stable?}
% Classify the fixed point $y = y' = y''=0$ for the ODE\kern .5pt s
% {\em (a)} $y''+5y' + 4y = 0$,
% {\em (\kern .7pt b)} $y'''+6y''+11y' + 6y = 0$.
% \par
% {\em Exercise $15.2$.  A system with four fixed points.}
% Consider the system
% $u' = {1\over 3}(u^{}-v^{})(1-u^{}-v^{})$, $v' = u^{}(2-v^{})$. 
% {\em (a)} Give a formula for
% the Jacobian matrix {\bf J} as a function of $u$ and $v$.
% {\em (\kern .7pt b)} Determine the fixed points and evaluate the
% Jacobian at these points.
% {\em (c)} Find the eigenvalues analytically, and classify the nature
% of each fixed point.
% \par
% {\em Exercise $15.3$.  Eulerian wobble.}
% Determine the fixed points of (10.9) and analyze their stability.
% \par
% {\em Exercise $15.4$.  Fixed points are unattainable.}
% Let ${\bf y}'(t) = {\bf f}({\bf y})$ be an autonomous ODE satisfying
% the continuity hypotheses of Theorems 11.2 or Theorem 11.3,
% and suppose ${\bf y}_*^{}$ is a fixed point.  Let ${\bf y}(t)$
% be the solution for $t>0$ with
% initial value ${\bf y}(0)$.  Prove that if
% ${\bf y}(0) \ne {\bf y}_*^{}$,
% then ${\bf y}(t) \ne {\bf y}_*^{}$ for all $t$.  
% \par
% {\em \underline{Exercise $15.5$}.  A cyclic system of three ODE\kern .5pt s}
% (adapted from Guckenheimer and Holmes,
% Structurally stable heteroclinic cycles, {\em Mathematical Proceedings of
% the Cambridge Philosophical Society,} 1988).
% Consider the system of ODE\kern .5pt s $u' = u(1-u^2-bv\kern .3pt ^2-c\kern .3pt w^2)$,
% $v' = v(1-v^2-b\kern .3pt w^2-c\kern .3pt u^2)$,
% $w' = w(1-w^2-b\kern .3pt u^2-c\kern .3pt v^2)$,
% where $b$ and $c$ are parameters.
% {\em (a)} Plot the solution $u(t)$ for $t\in [\kern .3pt 0,800\kern .3pt]$
% with $b=0.55$, $c=1.5$ and
% initial conditions $u(0)=0.5$, $v(0)=w(0)=0.49$.
% Make similar plots of $v(t)$ and $w(t)$, and also of the
% whole trajectory in $u$-$v$-$w$ space, and comment on these shapes.
% {\em (\kern .7pt b)} What are the four fixed points of this system that the plots just
% drawn come close to?
% For large $t$, the trajectory moves approximately in a cycle from
% one fixed point, to another, to a third, and then back again.  (It is
% approximating a {\em heteroclinic cycle}.) 
% Which fixed point is the trajectory near at $t=800\kern .5pt$?
% {\em (c)} Find the eigenvalues of the appropriate matrix at one of
% these fixed points.  What does this tell us about the structure of
% this fixed point?  How does that fit with the observed
% trajectory?
% \par
% {\em Exercise $15.6$.  Oregonator.}
% Exercise 10.2 presented the nonlinear equations known as the Oregonator.
% Here, continue with the parameters as presented in that exercise.
% {\em (a)} The origin $u=v=w=0$ is a fixed point. Determine the relevant
% eigenvalues and classify the linearized behaviour of the system there.
% {\em (\kern .7pt b)} There is another fixed point with $u=w>0$. Find it, determine
% its eigenvalues, and classify it. 
% </latex>
