%% 14. Linear systems and linearization %% % % \setcounter{page}{172} % \parindent=.6cm % \parskip=-1em % %% % % 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. % %% % % 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. % %% % % 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. % %% % % \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.$ % } % %% % % \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 % %% % % 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. % %% % % \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)$$ % } % %% % % \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 % %% % % 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. % %% % % \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)$$ % } % %% % % \def\qed{\vrule height4pt width3pt depth4pt} % {\em Proof.} Again the formulas are readily verified and % uniqueness comes from Theorem 11.2. \qed % %% % % \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. % %% % % 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$. % %% % % 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. % %% % 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. %% % % \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)$$ % } % %% % % \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 % %% % % 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. % %% % 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) %% % \vskip 1.02em %% % % \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$. % 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 %% % \vskip 1.02em %% % % \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. % %% % % 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_*^{}$. % %% % % \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)$$ % } % %% % \vskip -0.5em %% % 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)$$ %% % \vskip -1.5em 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) %% % \vskip 1.02em %% % % \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$. % 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 %% % \vskip 1.02em %% % % \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! % %% % % 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. % 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) %% % \vskip 1.02em %% % % \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. % 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 %% % \vskip 1.02em %% % % \newpage % \begin{center} % \hrulefill\\[1pt] % {\sc Application: linearized Lorenz trajectories}\\[-3pt] % \hrulefill % \end{center} % %% % % 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. % %% % % 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}. % %% % % 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 ]$. % 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]; %% % \vskip 1.02em %% % % \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. % J(4) %% % \vskip 1.02em %% % % \noindent % Now here is something curious. We look at the eigenvalues of this matrix, % eig(J(4)) %% % \vskip 1.02em %% % % \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? % %% % % 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$. % 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) %% % \vskip 1.02em %% % % \noindent % Most of the time, the spectral abscissa is positive, and its mean % is positive, % mean(spec_absc) %% % % \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. % %% % % 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. % 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) %% % \vskip 1.02em %% % % \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.'' % %% % % \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}$. % %% % % \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. % %% % % \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} % %% % \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)}. %