%% 14. Linear systems and linearization

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

%%
% <latex>
% Much of mathematics springs
% from two topics:\ calculus and linear algebra.  The subject
% of ODE\kern .5pt s is obviously rooted in calculus, and it has a fundamental
% linear algebra side too.  If an ODE is linear, then it is
% ``all'' linear algebra.  If it is nonlinear, then linear
% algebra still determines its local behavior
% near each point of a trajectory.
% </latex>

%%
% <latex>
% Let us start with the linear case.  If $\bf A$ is an
% $n\times n$ matrix, a first-order linear ODE with
% $n$ variables can be defined by
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t), \eqno (14.1) $$
% a special case of the general system (10.1).  Conversely,
% any first-order linear, autonomous, homogeneous $n$-variable system of ODE\kern .5pt s
% (``{\sf FLAsH}'')
% can be written in this form for some matrix ${\bf A}$.
% If ${\bf y}_0^{}$ is
% an $n$-vector, then one solution of (14.1) is
% $$ {\bf y}(t) = e^{t{\bf A}} {\bf y}_0^{}, \eqno (14.2) $$
% where $e^{t{\bf A}}$ is the {\bf matrix exponential,} defined by
% $$ e^{t{\bf A}} = {\bf I} + t{\bf A} + {1\over 2!} (t{\bf A})^2 + {1\over 3!} (t{\bf A})^3
% + \cdots . \eqno (14.3) $$
% This series always converges, providing a well-defined matrix
% exponential for any ${\bf A}$ and $t$.\footnote{A proof can be based
% on the observation that for any $t$, the powers $(t{\bf A})^k$ grow at most
% geometrically as $k\to\infty$ whereas the factorials $k!$ grow
% faster than geometrically.}  In fact, (14.3) is the
% Taylor series of $e^{t{\bf A}}$ about $t=0$, and from the theory
% of analytic functions it is known that the series can be differentiated
% term by term, yielding
% $$ {d\over dt} \kern .7pt e^{t{\bf A}} = {\bf A} e^{t{\bf A}}. $$
% This formula confirms that (14.2) is a solution of (14.1).
% Equations (14.1)--(14.2) are matrix generalizations of
% equations (2.1)--(2.2) of Chapter~2.
% </latex>

%%
% <latex>
% So linear, autonomous, homogeneous systems of ODE\kern .5pt s
% are all about exponentials of matrices.
% We record this conclusion as a theorem formulated for IVPs.
% </latex>

%%
% <latex> 
% \smallskip
% {\em
% {\bf Theorem 14.1.  Solution of first-order linear autonomous
% homogeneous IVP system (\textsf{\textbf{FLAsHI}}).}
% The problem
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t),
% \quad {\bf y}(0) = {\bf y}_0^{}, \eqno (14.4) $$
% where\/ ${\bf y}_0^{}$ is an $n$-vector and\/ $\bf A$ is an
% $n\times n$ matrix, has the unique solution
% $$ {\bf y}(t) = e^{t{\bf A}} {\bf y}_0^{} \eqno (14.5) $$
% valid for all $t$, $-\infty < t < \infty.$
% }
% </latex> 

%%
% <latex>
% \smallskip
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  We have just outlined an argument
% for existence, and uniqueness follows from Theorem 11.2. A suitable
% Lipschitz constant for (14.4) is $K = \|{\bf A}\|$; as
% explained in Chapter 11, the choice of norm does not matter.  \qed
% \smallskip
% </latex>

%%
% <latex>
% A small extension of (14.1) is to make the ODE inhomogeneous
% (``{\sf FLAsh}''),
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t) + {\bf g}, $$
% where ${\bf g}$ is a fixed $n$-vector.  We might call this
% an {\bf affine} autonomous system of equations, though most of the
% time we will just say ``linear.''
% Following Theorem 2.4 for
% the scalar case, we readily derive the solution as follows.
% </latex>

%%
% <latex> 
% \smallskip
% {\em
% {\bf Theorem 14.2.  Solution of first-order linear autonomous
% IVP system with constant inhomogeneity (\textsf{\textbf{FLAshI}}).}
% The problem
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t) + {\bf g},
% \quad {\bf y}(0) = {\bf y}_0^{}, \eqno (14.6) $$
% where\/ ${\bf y}_0^{}$ and\/ $\bf g$ are\/ $n$-vectors and\/ $\bf A$ is an
% $n\times n$ matrix, has the unique solution
% $$ {\bf y}(t) = e^{t{\bf A}} {\bf y}_0
% + \int_0^t e^{(t-s){\bf A}} {\bf g} \kern .5pt ds . \eqno (14.7) $$
% If $A$ is nonsingular this reduces to
% $$ {\bf y}(t) = e^{t{\bf A}} {\bf y}_0
% + {\bf A}^{-1} (e^{t{\bf A}} - {\bf I}) {\bf g}. \eqno (14.8) $$
% }
% </latex> 

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  It can be verified explicitly that (14.7)
% and (14.8) are solutions of (14.6), and 
% uniqueness again follows from Theorem 11.2. \qed
% \smallskip
% </latex>

%%
% <latex>
% Theorem 14.2 assumes a constant inhomogeneity, but
% the variable case is also straightfoward
% (again ``{\sf FLAsh}''),
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t) + {\bf g}(t), $$
% where ${\bf g}(t)$ is a $n$-vector function of\/ $t$.
% Here is a systems generalization of a result midway between
% Theorems 2.3 and 2.4 (since $\bf g$ but not $\bf A$ has
% been made $t$-dependent).  As in those theorems, the
% ODE is understood to hold everywhere except at any
% points of discontinuity.
% </latex>

%%
% <latex> 
% \smallskip
% {\em
% {\bf Theorem 14.3.  Solution of first-order linear autonomous
% inhomogeneous IVP system (\textsf{\textbf{FLAshI}}).}
% The problem
% $$ {\bf y}'(t) = {\bf A} {\bf y}(t) + {\bf g}(t),
% \quad {\bf y}(0) = {\bf y}_0^{}, \eqno (14.9) $$
% where\/ ${\bf y}_0^{}$ is an\/ $n$-vector,
% ${\bf g}(t)$ is an\/ $n$-vector piecewise continuous function
% of\/ $t$, and\/ $\bf A$ is an
% $n\times n$ matrix, has the unique solution
% $$ {\bf y}(t) = e^{t{\bf A}} {\bf y}_0
% + \int_0^t e^{(t-s){\bf A}} {\bf g}(s) \kern .5pt ds . \eqno (14.10) $$
% }
% </latex> 

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  Again the formulas are readily verified and
% uniqueness comes from Theorem 11.2. \qed
% </latex>

%%
% <latex>
% \smallskip
% The discussion above is confined to autonomous problems,
% but nonautonomous linear equations of the
% form ${\bf y}' = {\bf A}(t) \kern .3pt {\bf y} + {\bf g}(t)$
% are important too.  
% In this situation the matrix exponential $e^{t{\bf A}}$
% generalizes to a so-called {\em fundamental matrix} ${\bf Y}(t)$,
% discussed in equations (19.17)--(19.18) of Chapter~19.
% The special case ${\bf g}={\bf 0}$ also appears
% in equation (14.32) below.
% </latex>

%%
% <latex>
% Some ODE\kern .5pt s are linear to start with,
% but the greatest importance of
% linear ODE\kern .5pt s comes from the fact that a nonlinear ODE can be
% {\em linearized.}  By this we mean that locally, near
% a given time $t_0^{}$, the solution of a nonlinear ODE
% $$ {\bf y}'(t) = {\bf f}(t,{\bf y}(t)) \eqno (14.11) $$
% will evolve approximately like that of a linear one.
% In fact, every autonomous ODE looks locally like an affine
% autonomous system of the form (14.6), so long as
% the coefficients are twice differentiable.
% We now work out the details of such approximations,
% which are a matter of multivariate calculus.
% As usual, we keep the formulas simple by assuming $t_0^{} = 0$,
% even though in an application we may linearize about
% values $t_0^{} \ne 0$.
% </latex>

%%
% <latex>
% Suppose we have an autonomous ODE
% $$ {\bf y}'(t) = {\bf f}({\bf y}(t)) \eqno (14.12) $$
% and are interested in the behavior of solutions
% ${\bf y}(t)$ near a particular point~${\bf y}_*^{}$.
% Specifically, we consider the solution of the IVP
% $$ {\bf y}'(t) = {\bf f}({\bf y}(t)), \quad 
% {\bf y}(0) = {\bf y}_0^{} \eqno (14.13) $$
% for some initial value
% ${\bf y}_0^{}$ close to or equal to ${\bf y}_*^{}$.
% We use the abbreviation
% $$ {\bf f}_*^{} = {\bf f}({\bf y}_*^{}) \eqno (14.14) $$
% and suppose that $\bf f$ is twice differentiable at $\bf y={\bf y}_*^{}$.
% The partial derivative $\partial {\bf f} / \partial {\bf y}$
% is an $n\times n$ {\bf Jacobian matrix},
% $$ {\bf J} = \displaystyle{\partial {\bf f} \over \partial {\bf y}} =
% \pmatrix{\displaystyle{\partial f_1^{}\over \partial y_1^{}} & \cdots &
% \displaystyle{\partial f_1^{}\over \partial y_n^{}}
% \cr \vdots & & \vdots \cr
% \noalign{\vskip 5pt}
% \displaystyle{\partial f_n^{}\over \partial y_1^{}} & \cdots &
% \displaystyle{\partial f_n^{}\over \partial y_n^{}} }. $$
% We abbreviate the value of $\bf J$ at ${\bf y} = {\bf y}_*^{}$ by
% $$ {\bf J}_*^{} = {\bf J}({\bf y}_*^{}). \eqno (14.15) $$
% We now calculate from the definition of ${\bf J}_*^{}$
% $$ {\bf f}({\bf y}(t)) = {\bf f}_*^{} +
% {\bf J}_*^{} ({\bf y}(t)-{\bf y}_*^{}) +
% O(\|{\bf y}(t)-{\bf y}_*^{}\|^2). \eqno (14.16) $$
% This equation tells us that up to an error
% $O(\|{\bf y}(t) - {\bf y}_*^{}\|^2)$, any autonomous ODE
% can be regarded as linear, or rather affine, as in (14.6).
% Here $\|\cdot\|$ can denote any norm on the space of
% $n$-vectors, since all norms on a finite-dimensional space 
% are equivalent.
% </latex>

%%
% The formulas will be simpler if we work with the difference
% variables
% $$ \delta{\bf y}(t) = {\bf y}(t) - {\bf y}_*^{}, \quad
% \delta{\bf y}_0^{} = {\bf y}_0^{} - {\bf y}_*^{}. $$
% With these definitions, (14.8) and (14.16) 
% give us the following conclusions.

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 14.4.  Linearization of an autonomous
% system of ODE\kern .5pt s (\textsf{\textbf{FlAsHI}}).}\ \ Let
% ${\bf y}_*^{}$ be fixed and assume $\bf f$ is twice differentiable
% at\/ ${\bf y} = {\bf y}_*^{}$ with Jacobian ${\bf J}_*^{}$.
% The solution of\/ $(14.13)$ satisfies
% $$ \delta {\bf y}'(t) = {\bf f}_*^{} + {\bf J}_*^{} \delta {\bf y}(t) +
% O(\|\delta {\bf y}(t)\|^2) \eqno (14.17) $$
% and if\/ ${\bf J}_*^{}$ is nonsingular,
% $$ \delta {\bf y}(t) = e^{t{\bf J}_*^{}} \delta {\bf y}_0
% + ({\bf J}_*^{})^{-1}(e^{t{\bf J}_*^{}} - {\bf I}) \kern .7pt {\bf f}_*^{} +
% O(t\kern .7pt\Delta(t)^2)\eqno (14.18) $$
% if\/ $\Delta(t)$ is an upper bound on $\|\delta{\bf y}(s)\|$ for\/
% $0\le s \le t$.  If\/ ${\bf y}_0^{} = {\bf y}_*^{}$,
% then $\delta{\bf y} = {\bf 0}$, and this estimate becomes
% $$ \delta {\bf y}(t) = ({\bf J}_*^{})^{-1}(e^{t{\bf J}_*^{}}
% - {\bf I})\kern .7pt {\bf f}_*^{} + O(t^3).\eqno (14.19) $$
% }
% </latex>

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Sketch of proof.}
% Equation (14.17) is a restatement of
% (14.16). Equation (14.18) follows by
% applying Theorem 14.2 to (14.17) with ${\bf A}={\bf J}_*^{}$ and
% ${\bf h} = {\bf f}_*^{}$.
% Equation (14.19) follows
% from (14.18) together with the estimate $\Delta(t) =O(t)$
% in the special case $\delta{\bf y}(0) = {\bf 0}$.
% (In this outline we are assuming that the $O(\kern 1pt \cdot\kern 1pt )$ terms
% behave in the obvious matter, with a small perturbation of an ODE
% resulting in a small perturbation of its solution.  A careful
% argument could justify such steps by the use of a lemma known
% as {\em Gronwall's inequality.})
% \qed
% \smallskip
% </latex>

%%
% <latex>
% We now explore three examples of Theorem 14.4 in action, first
% for a scalar first-order problem,
% then for a scalar second-order problem
% analyzed by reduction to a first-order system, then for a 
% two-variable first-order system.
% </latex>

%%
% For the scalar first-order example, we recall the IVP (3.9),
% $$ y'(t) = y(t)^2, \quad y(0) = 1, \eqno (14.20) $$
% whose solution $y(t) = (1-t)^{-1}$ blows up to $\infty$
% at $t=1$.  Here is a solution plotted as in Chapter 3.
ODEformats
N = chebop(0,1); N.op = @(t,y) diff(y) - y^2; N.lbc = 1;
N.maxnorm = 25; y = N\0; plot(y,CO,ivpnl), axis([0 1.1 0 30])
hold on, plot([1 1],[0 30],'--k')
title('Fig.~14.1.~~Solution to blowup equation (14.20)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The Jacobian for this equation is the scalar $J(y) = 2\kern .3pt y$,
% which takes the value $J_*^{} = 2\kern .3pt y_0^{}$ at $y_*^{} = y_0^{}$,
% so Theorem 14.4 gives the following linear approximation for
% $t\approx 0$, $y\approx y_0^{}$:
% $$ y'(t) \approx (y_0^{})^2 + 2\kern .3pt y_0^{}(y(t)-y_0^{}), \quad y(0) = y_0^{}, $$
% or near a possibly nonzero time $t_0^{}$,
% $$ y'(t) \approx (y_0^{})^2 + 2\kern .3pt y_0^{}(y(t)-y_0^{}), \quad y(t_0^{}) = y_0^{}. \eqno (14.21) $$
% Let us add dotted arcs to the plot corresponding
% to solutions of this linear equation emanating from
% $t_0^{} = 0.5$ and $0.8$.
% </latex>
for t0 = [0.5 0.8]
  y0 = y(t0); f0 = y0^2; J0 = 2*y(t0); L = chebop(t0,t0+.22);
  L.op = @(t,u) diff(u) - f0 - J0*(u-y0); L.lbc = y0;
  u = L\0; plot(u,':r',LW,2.5), plot(t0,y0,'.r')
end
title('Fig.~14.2.~~Same with two linearizations added',FS,11)
hold off
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% Each linearized solution is shown over an interval of length $0.22$.
% At the initial point, as asserted by (14.19),
% the linearization matches both the slope and the curvature.
% The approximation starting at $t = 0.5$ tracks the true solution 
% closely over the interval shown,
% whereas the approximation starting at $t = 0.8$ loses
% accuracy more quickly.
% Note that as the solution to a linear problem, it will exist for
% all $t$, unlike the underlying nonlinear problem with its blowup.
% </latex>

%%
% <latex>
% Next we look at an example involving a second-order
% scalar ODE.\ \ This is such an important situation that it is worth
% spelling out the special form Theorem~14.4 takes in this case. 
% Suppose we have a second-order scalar autonomous problem
% $$ u''(t) = F(u(t),u'(t)), \quad u(0) = u_0^{},
% ~~ u'(0) = v_0^{}. \eqno (14.22) $$
% Setting $v(t) = u'(t)$, we can write this as a first-order system
% in the variable
% ${\bf y}(t) = (u(t),v(t))^T$,
% $$ u'(t) = v(t), ~~ v'(t) = F(u(t),v(t)),
% \quad u(0) = u_0^{}, ~~ v(0) = v_0^{}. $$
% The Jacobian ${\bf J}_*^{}$ is
% $$ {\bf J}_*^{} = \pmatrix{0 & 1 \cr\noalign{\vskip 4pt}
% \partial F/\partial u & \partial F/\partial v} , $$
% and (14.17) becomes $u'(t) = v(t)$ together with
% $$ v'(t) = F(u_*^{},v_*^{}) +
% {\partial F\over \partial u} (u(t)-u_*^{}) + 
% {\partial F\over \partial v} (v(t)-v_*^{}). $$
% Translating this back to the original second-order
% scalar context gives this variant of Theorem 14.4,
% with $\delta u(t) = u(t) - u_*^{}$ and
% $\delta v(t) = u'(t) - v_*^{}$.
% </latex>

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 14.5.  Linearization of an autonomous second-order
% scalar ODE (\textsf{\textbf{flASHI}}).}\ \ Let
% the function $F$ defining the scalar second-order IVP\/ $(14.22)$
% be twice differentiable with respect to its first and second
% arguments at $u=u_*^{}$, $u' = v_*^{}$, and define
% $$ F_*^{} = F(u_*^{},v_*^{}), \quad
% a_*^{} = {\partial F\over \partial u} (u_*^{},v_*^{}), \quad
% b_*^{} = {\partial F\over \partial u'} (u_*^{},v_*^{}). $$
% Then $u(t)$ satisfies
% $$ (\delta u(t))'' = F_*^{} + a_*^{} \delta u(t) +
% b_*^{} \delta v(t)
% + O((\delta u(t))^2) + O((\delta v(t))^2).  \eqno (14.23) $$
% If $u_0^{} = u_*^{}$ and $v_0^{} = v_*^{}$, this becomes
% $$ (\delta u(t))'' = F_*^{} + a_*^{} \delta u(t) +
% b_*^{} \delta v(t) + O(t^2) .  \eqno (14.24) $$
% }
% </latex>

%%
% <latex> \vskip -0.5em </latex>

%%
% For an example, we turn to our friend the van der Pol equation, as in (8.8),
% $$ u''=5(1-u^2)u'-u, \quad t\in [\kern .3pt 0,15], \quad u(0)=1, ~u'(0)=0.
% \eqno (14.25) $$
%%
% <latex> \vskip -1.5em </latex>
N = chebop(0,15);
N.op = @(t,u) diff(u,2) - 5*(1-u^2)*diff(u) + u;
N.lbc = [1;0]; u = N\0; plot(u,CO,ivpnl), hold on
title('Fig.~14.3.~~Van der Pol equation (14.25)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% By Theorem 14.5, the linearized approximation near a point $t_0^{}$
% with $u(t_0^{})= u_*^{}$ and $u'(t_0^{})=v_*^{}$ is
% $$ u''(t) \approx 5(1-(u_0^{})^2)v_0^{} - u_0^{}
% (10u_0^{}v_0^{}+1)(u(t)-u_0^{}) + 5(1-u_0^2) (u'(t) - v_0^{}) . $$
% We use this result to add three curves to the plot corresponding
% to times $t_0^{} = 3$, $6.9$, and $11.5$. 
% </latex>
for t0 = [3 6.9 11.5]
  u0 = u(t0); v0 = deriv(u,t0);
  F0 = 5*(1-u0^2)*v0-u0; a0 = -(10*u0*v0+1); b0 = 5*(1-u0^2);
  L = chebop(t0,t0+1.5);
  L.op = @(t,w) diff(w,2) - F0 - a0*(w-u0) - b0*(diff(w)-v0);
  L.lbc = [u0; v0];
  w = L\0; plot(w,':r',LW,2.5), plot(t0,u0,'.r')
end
title('Fig.~14.4.~~Same with three linearizations added',FS,11)
hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% According to (14.24), the linearization of this
% problem captures $u''$ with accuracy $O(t^2)$.  It
% follows that $u$ itself is captured
% with accuracy $O(t^4)$: in the figure,
% the dotted lines match not just $y(0)$,
% $y'(0)$, and $y''(0)$, but also $y'''(0)$ (see Exercise~14.3).
% Yes, the oscillatory middle approximation is correct!
% </latex>

%%
% <latex>
% Now we consider an example to illustrate Theorem 14.5 with
% a generic autonomous system of two equations, that is, one
% not obtained from a second-order scalar problem.  
% Here are the Lotka--Volterra equations (10.6):
% $$ u' = u-uv, ~~ v'= -\textstyle{1\over 5}v + uv, \quad t\ge 0, ~ u(0) = u_0^{},
% ~v(0) = v_0^{}. \eqno (14.26) $$
% In Chapter~10 we plotted a limit cycle solution.
% </latex>
N = chebop(0,20); N.lbc = [1;1];
N.op = @(t,u,v) [diff(u)-u+u*v; diff(v)+.2*v-u*v];
[u,v] = N\0; arrowplot(u,v,CO,ivpnl,MS,4,YS,.5), hold on, axis([0 1.3 0 2.7])
plot([0 .2],[0 1],'.k',MS,8)
title('Fig.~14.5.~~Lotka--Volterra cycle in the phase plane',FS,11)
xlabel('rabbits, $u$',FS,9,IN,LT), ylabel('foxes, $v$',FS,9,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The Jacobian for (14.25) is
% $$ {\bf J}_0^{} = \pmatrix{1-v & -u \cr\noalign{\vskip 4pt}
% v & u - 1/5}. $$
% This enables us to add dotted linearized trajectories as before.
% </latex>
for t0 = [0 2 14]
  u0 = u(t0); v0 = v(t0); f01 = u0-u0*v0; f02 = -v0/5+u0*v0;
  L = chebop(t0,t0+1); L.lbc = [u0; v0];
  L.op = @(t,U,V) [diff(U)-f01-(1-v0)*(U-u0)+u0*(V-v0); ...
                   diff(V)-f02-v0*(U-u0)-(u0-1/5)*(V-v0)];
  [U,V] = L\0; plot(U,V,':r',LW,2.5), plot(u0,v0,'.r')
end
title('Fig.~14.6.~~Same with three linearizations added',FS,11)
hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \newpage
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: linearized Lorenz trajectories}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%% 
% <latex>
% Chaos is the quintessential nonlinear effect, yet at each
% point of time, a chaotic ODE, like any other ODE, behaves linearly.
% Let us explore this effect for the Lorenz equations (13.1),
% $$ u' = 10(v-u), ~~ v'= u(28-w) - v,
% ~~ w' = uv - (8/3) w. \eqno (14.27)  $$
% Suppose at some time $t_0^{}$ a trajectory ${\bf y}(t)$ of
% (14.27) takes
% the values ${\bf y}(t_0^{}) = {\bf y}_*^{} = (u_*^{}$, $v_*^{}$, $w_*^{})^T$.  
% If we define as usual $\delta{\bf y} = {\bf y} - {\bf y}_*^{},$
% then from Theorem 14.4 we have
% $$ \delta {\bf y}' \approx {\bf f}_*^{} + {\bf J}_*^{} \delta {\bf y}, 
% \eqno (14.28) $$
% where we can calculate the Jacobian as
% $$ {\bf J}(u,v,w) = \pmatrix{-10 & 10 & 0 \cr
% 28-w & -1 & -u\cr v & u & -8/3 }. \eqno (14.29) $$
% This gives us a linear model of the local behavior near
% the given time and solution values.
% </latex>

%%
% <latex>
% Let us focus on what the model tells us about perturbations.  If
% ${\bf y}(t)$ is the given trajectory, let ${\bf y}(t) + \Delta {\bf y}(t)$ be
% another nearby trajectory.  Then by subtracting (14.28) applied to one
% solution from (14.28) applied to the other, we obtain
% $$ \Delta {\bf y}' \approx {\bf J}_*^{} \Delta {\bf y}. \eqno (14.30) $$
% This approximation is valid near the given time $t_0^{}$, but
% if we let ${\bf J}$ vary with $t$, so that it always corresponds
% to the local Jacobian along the trajectory, we get
% $$ \Delta {\bf y}' \approx {\bf J}(t) \Delta {\bf y}. \eqno (14.31) $$
% If $\Delta {\bf y}$ is infinitesimal, this becomes an equality,
% $$ \Delta {\bf y}' = {\bf J}(t) \Delta {\bf y} \quad
% \hbox{($\Delta{\bf y}$ infinitesimal)}. \eqno (14.32) $$
% This kind of linearization is applicable to all nonlinear problems.
% {\em Infinitesimal perturbations of a nonlinear ODE evolve
% according to a linear, homogeneous, nonautonomous equation.}
% In other words, infinitesimal perturbations are a problem of type
% {\sf FLasHI}.
% </latex>

%%
% <latex>
% We can look at the Jacobians for the Lorenz equations numerically.   
% Starting from our usual initial values,
% we run the system up to time $t=10$ and
% make a function {\tt J} that evaluates ${\bf J}(t)$ for any
% $t\in [\kern .3pt 0,10\kern .3pt ]$.
% </latex>
N = chebop(0,10);
N.op = @(t,u,v,w) [diff(u)+10*u-10*v; ...
    diff(v)-(28-w)*u+v ; diff(w)-u*v+(8/3)*w];
N.lbc = [-15; -15; 20]; [u,v,w] = N\0;
J = @(t) [-10 10 0; 28-w(t) -1 -u(t); v(t) u(t) -8/3];
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here for example is ${\bf J}(t)$ at
% $t = 4$, a number big enough that initial transients
% have died away and the solution is close to the strange attractor.
% </latex>
J(4)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Now here is something curious.  We look at the eigenvalues of this matrix,
% </latex>
eig(J(4))
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% They all have negative
% real part!  This is a surprise because the Lorenz equations are
% chaotic, meaning that infinitesimal perturbations grow exponentially.
% In fact, we saw in Chapter 13 that on average they grow
% approximately at the rate $\exp(0.91 t)$. 
% Yet eigenvalues with negative real part correspond to exponential
% decay, not growth.  What's going on?
% </latex>

%%
% <latex>
% The explanation is that even though a chaotic system is characterized by
% exponential growth of perturbations on average, there can be points
% along its trajectories where the perturbations are shrinking.  
% This is just what happens with the Lorenz equations.
% To see the effect, we can plot the {\em spectral abscissa} of ${\bf J}(t)$, the
% maximum real part of its eigenvalues, as a function of $t$.
% </latex>
spec_absc = chebfun(@(t) max(real(eig(J(t)))),[4,10],'splitting','on');
plot(spec_absc,'k'), ylim([-6 12])
title(['Fig.~14.7.~~Spectral abscissa of ' ...
       'Lorenz eqs.\ Jacobian'],FS,11)
xlabel('$t$',FS,10), ylabel('spectral abscissa',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Most of the time, the spectral abscissa is positive, and its mean
% is positive,
% </latex>
mean(spec_absc)

%%
% <latex>
% \noindent 
% Nevertheless, it regularly dips below zero.
% If the nonlinear problem were frozen with fixed linear dynamics at such a
% time, then perturbations would eventually decay rather than grow.
% </latex>

%%
% <latex>
% As a fine point of linear algebra, let us say a little more of
% what it means for a matrix $\bf J$ to
% have all its eigenvalues in the complex left half-plane.  This implies
% that solutions to the constant-coefficient problem ${\bf y}' = {\bf J}
% {\bf y}$ 
% must decay as $t\to\infty$.  It does not imply that they must
% decay even at the start, for small $t$.  A different condition of
% linear algebra is needed to ensure this stronger property, namely
% that the {\em numerical abscissa} of the matrix is negative.
% The numerical abscissa of $\bf J$ is the maximum eigenvalue of
% $({\bf J}+{\bf J}^T)/2$, and here we add another curve
% to the plot showing how this quantity varies for our Lorenz trajectory.
% </latex>
numer_absc = chebfun(@(t) max(eig((J(t)+J(t)')/2)),[4,10],'splitting','on');
hold on, plot(numer_absc,'r'), ylim([-6 12]), hold off
title(['Fig.~14.8.~~Numerical abscissa of ' ...
       'Lorenz eqs.\ Jacobian'],FS,11)
xlabel('$t$',FS,10), ylabel('numerical abscissa',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Evidently even the numerical abscissa dips below zero
% intermittently, implying that at certain times
% the behavior of the Lorenz equations is locally dissipative by any definition, 
% though chaos reigns in the large.  Analogous effects can be found
% in atmospheric dynamics and other scientific problems and
% are sometimes quantified by means of a ``local Lyapunov exponent.''
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}
% Perhaps the most extreme case of linearization in the history of
% science was reported in {\em Physical Review Letters} in February
% 2016 in the paper ``Observation of gravitational waves from a
% binary black hole merger,'' which is sure to win a Nobel Prize
% before long for some of its 1011 authors.\footnote{It is one of
% the most compellingly written papers we have seen --- exciting
% reading for anyone interested in science.  The first three authors,
% incidentally, are Abbott, Abbott, and Abbott, and the last three
% are Zucker, Zuraw, and Zweizig.} Einstein's ten coupled
% partial differential equations of general relativity, from 1915, are
% forbiddingly nonlinear.  In 1916, however, Einstein
% proposed that when linearized to small
% amplitudes, the equations predict the possibility of propagation
% of {\em gravitational waves.} Meanwhile, 1.3 billion years earlier,
% two black holes had merged together in an event so cataclysmic
% that for a fraction of a second it radiated energy with greater power than
% all the stars in the observable universe combined.
% For 1.3 billion years the gravitational waves flew outward at the speed of
% light, losing amplitude in the usual inverse-radius fashion,
% and they reached the earth on September 14, 2015, when the
% minute oscillations were detected by the LIGO project simultaneously
% in Louisiana and
% Washington State, USA.\ \ One might fancifully say that the
% merger of the black holes was the most nonlinear event ever
% observed by mankind, and the signal by which it was actually observed
% was the most linear.  The amplitude of that signal
% was almost unimaginably small, as is
% reflected in the astonishing label of the $y$-axis in Figure~1
% of the 2016 paper, informing the reader that the figure shows
% relative deflections in units of $10^{-21}$.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  
% A classic textbook with a great emphasis on the linear
% algebra side of ODE theory was Hirsch and Smale, 
% {\em Differential Equations, Dynamical Systems, and Linear Algebra,}
% published in 1974 by Academic Press.  Later a third author was
% added, Devaney, and the title was changed to
% {\em Differential Equations, Dynamical Systems, and an
% Introduction to Chaos} (third edition Elsevier, 2013).  The books
% are quite different, with less about the linear algebra side in
% the later edition; and both are very interesting.  
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 14.} 
% An autonomous ODE ${\bf y}'(t) = {\bf f}({\bf y}(t))$ with a sufficiently
% smooth function $\bf f$ can be approximated near a particular
% time~$t_0^{}$ and value ${\bf y}_*^{}$ by the affine equation
% $\delta {\bf y}'(t) \approx {\bf f}_*^{} + {\bf J}_*^{} \delta {\bf y}(t),$
% where~${\bf f}_*^{}$ and ${\bf J}_*^{}$ are the function $\bf f$ and
% Jacobian matrix $\partial {\bf f}/ \partial {\bf y}$ frozen at
% $t=t_0^{}$ and ${\bf y} = {\bf y}_*^{}$ and
% $\delta {\bf y}(t) = {\bf y}(t) - {\bf y}_*^{}$.
% Globally over a range
% of values of\/~$t$, infinitesimal perturbations to a nonlinear trajectory
% evolve according to the linear homogeneous nonautonomous ODE 
% $\Delta {\bf y}' = {\bf J}(t) \Delta {\bf y},$ where
% ${\bf J}(t)$ is the Jacobian matrix at time $t$ for the given nonlinear
% trajectory.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em Exercise $14.1$.  Exponential of a matrix and linear systems of BVPs.}
% {\em (a)} Find formulas for the powers ${\bf A}^k$ of the
% $2\times 2$ matrix ${\bf A} = \hbox{\tt [0 1; -1 0]}$ (Matlab notation).
% {\em (\kern .7pt b)} Using the results of (a) and familiar formulas
% for Taylor series, find a formula for $\exp(t{\bf A})$.
% {\em (c)} What is the smallest $t>0$ for which 
% $\exp(t{\bf A})$ is skew-diagonal, that is, having zeros in the
% positions on the diagonal?  Call this number~$T$.
% What is the matrix $\exp(T\kern -1pt {\bf A})$?
% {\em (d)} Show that the BVP ${\bf y}' = {\bf A}{\bf y}$ with boundary conditions
% $y_1(0) = a$, $y_2(T)=b$ is ill-posed, having either no solutions
% or infinitely many solutions depending on the values of $a$ and $b$.
% {\em (e)} Show on the other hand that if~${\bf A}$ is any diagonal matrix,
% this BVP is well-posed, with a unique solution for
% each $a$ and $b$.
% \par
% {\em Exercise $14.2$.  Affine autonomous linear systems.}
% Verify that (14.7) is a solution of (14.6).
% \par
% {\em \underline{Exercise $14.3$}.  Confirming fourth-order accuracy.}
% According to the text, the dashed red lines in Figure 14.4 match
% the green curve they approximate to accuracy $O(t^4)$, or more
% precisely $O((t-t_0^{})^4)$ for an approximation near $t=t_0^{}$.
% Confirm this numerically for the oscillatory middle curve 
% by means of an appropriate log-log plot.
% \par
% {\em \underline{Exercise $14.4$}.  Exponential of a nonnormal matrix.}
% Let $\bf A$ be the $6\times 6$ bidiagonal
% matrix with $-0.5, -0.6, \dots, -1$ on the
% main diagonal, $2,2,\dots,2$ on the first superdiagonal, and $0$ everywhere
% else. 
% {\em (a)} Calculate the spectral abscissa and numerical abscissa
% of~$\bf A$.  
% {\em (\kern .7pt b)} Make a semilogy plot of
% $\|e^{t{\bf A}}\|$ as a function of $t\in [\kern .3pt 0, 40\kern .3pt]$,
% where $\|\cdot\|$ is the 2-norm. (In Matlab the appropriate command
% is {\tt expm}.)  Superimpose appropriate dashed lines to
% indicate how this curve matches the results of {\em (a)}.
% </latex>
