%% 10. Systems of equations % % \setcounter{page}{114} % %% % % Up to this point, our discussion has concerned scalar ODE\kern .5pt s.\ \ In this % chapter we begin to talk about systems of equations, which % involve more than one dependent variable, or % equivalently, a dependent variable that is a vector rather than % a scalar. % %% % % All our vector problems will be of first order. A general first-order % ODE can be written % $${\bf y}'(t) = {\bf f}(t,{\bf y}) \eqno (10.1)$$ % for a linear or nonlinear function $\bf f$, where for each % $t$, ${\bf y}(t)$ % represents a vector of dimension $n$, % $${\bf y}(t) = \pmatrix{y_1^{}(t)\cr \vdots \cr y_n^{}(t) }, \quad % {\bf f}(t,{\bf y}) = % \pmatrix{f_1^{}(t,{\bf y})\cr \vdots \cr f_n^{}(t,{\bf y}) }.$$ % The dependent variables are thus $y_1^{},\dots,y_n^{}$, and % we could write the system out explicitly like this: % $$y_1'(t) = f_1^{}(t, y_1^{}(t),\dots, y_n^{}(t)) ,$$ % $$\vdots$$ % $$y_n'(t) = f_n^{}(t, y_1^{}(t),\dots, y_n^{}(t)) .$$ % For systems of small size we will usually find it % convenient to use variable names without subscripts, such as $u,v,w$. % The reason we can restrict the discussion to first order is that any % higher-order ODE can be written as a system of first-order ODE\kern .5pt s by % introducing additional variables, as outlined in % Exercise 1.5.\ \ In {\sf FLASHI} nomenclature, any % {\sf f\kern 1pt S$\,\cdots$} can be turned into % an {\sf Fs$\,\cdots.$}\footnote{This doesn't % mean that in practice, higher-order systems are necessarily reduced % to first-order form. For example, the Solar % System Dynamics Group at the Jet Propulsion % Laboratory in California solves systems of hundreds % of ODE\kern .5pt s numerically to track the moon and % planets and other bodies: when the media tell you there's going to % be an eclipse, % the computations probably came from JPL.\ \ The methods used % by the SSD Group are based on a second-order formulation, without reduction % to first order.} % %% % % To fully specify an ODE problem from (10.1), one normally % needs $n$ additional conditions. In an initial value problem, % we prescribe the value of $\bf y$ at a fixed time such as $t=0$, % $${\bf y}(0) = {\bf a} . \eqno (10.2)$$ % This amounts to specifying each of the $n$ components, % $$y_1^{}(0) = a_1^{}, ~\dots, \quad y_n^{}(0) = a_n^{} .$$ % In a boundary value problem, different components of $\bf y$ or % combinations thereof may be % specified at different points, for a total normally of $n$ % conditions all together. Exactly what specifications may lead % to a well-posed BVP for a system of equations % is not always a straightforward matter. % For IVPs, however, existence and uniqueness are guaranteed so % long as $\bf f$ is continuous with respect to $t$ and $\bf y$ and Lipschitz % continuous with respect to $\bf y$, as we shall prove in the next % chapter (Theorem~11.2). % %% % % For example, consider the two-variable autonomous system % $$u' = -v - u(u^2+v^2), \quad v' = u - v(u^2+v^2). \eqno (10.3)$$ % Here are the trajectories of $u$ and $v$ for $t\in[\kern .3pt 0,30\kern .3pt]$ % with $u(0) = 1$, $v(0) = 0$.\footnote{For % plots like this showing multiple components, we revert to % Matlab's built-in color tables rather than the light/dark green and blue % colors usually used in this book for linear/nonlinear IVPs and BVPs.} % ODEformats, N = chebop(0,30); N.lbc = [1; 0]; N.op = @(t,u,v) [diff(u)+v+u*(u^2+v^2); diff(v)-u+v*(u^2+v^2)]; [u,v] = N\0; plot([u v]), title('Fig.~10.1:~~Equation (10.3)',FS,11) %% % \vskip 1.02em %% % % \noindent % We see two interlocking oscillations, decaying rapidly in amplitude at % first but then much more slowly. % %% % % Since (10.3) is autonomous, the whole future % of the system is determined by the values of $u$ and $v$ at % any given time. Thus it makes sense % to plot this solution in the {\bf phase plane}, % with $u$ on the horizontal axis and $v$ on the vertical axis. % This is a straightforward generalization of % the phase plane of the last chapter, where the two % axes represented $y$ and $y'$. % Here, however, trajectories need not point % rightward in the upper half-plane. % arrowplot(u,v,CO,ivpnl,LW,.5,MS,4) title('Fig.~10.2.~~Solution of (10.3) in the phase plane',FS,11) axis equal, ylim([-.4 .6]), hold on, plot(0,0,'.k',MS,8) xlabel('$u$',FS,10,IN,LT), ylabel('$v$',FS,10,IN,LT) %% % \vskip 1.02em %% % % \noindent % The dot at the center marks the fixed point of the % system, corresponding to the constant solution % $u(t) = v(t) = 0$. In general, a {\bf fixed point} of % an autonomous ODE ${\bf y}' = {\bf f}({\bf y})$ % is any point ${\bf y}_0^{}$ such that % ${\bf f}({\bf y}_0^{}) = {\bf 0}$, implying that % the constant ${\bf y}(t) = {\bf y}_0^{}$ is a solution. % The dynamics of this problem is angularly symmetric, as % can be confirmed by converting $\bf f$ to polar coordinates. % %% % As in the last chapter, we may think of the trajectory as % following the arrows of a direction field: quiver(N,[-.7 1.1 -.4 .6],'k','scale',2,'xpts',12,'ypts',8) axis equal, ylim([-.4 .6]), hold off title('Fig.~10.3.~~Same with quiver arrows added',FS,11) %% % % \noindent % One way to understand the motion is to % interpret the nonlinear terms of (10.3) as linear terms % with varying coefficients. For a positive constant % $\kappa$, the system % $$u' = -v - \kappa u, \quad v' = u - \kappa v$$ % describes a simple oscillation with damping. Locally, near % any particular point $(u(t), v(t))$, (10.3) will behave % in just this way with $\kappa = u(t)^2+v(t)^2$. This % explains the behavior of the two images above: we start % with a damping coefficient $\kappa = 1$, but at the % end of the trajectory, $\kappa$ is about 60 times smaller: % kappa = u(end)^2 + v(end)^2 %% % % We emphasized in the last chapter that phase plane or more % generally phase space analysis (for $n\ge 3$) % is appropriate only for a system that is autonomous, % so that the state $\bf y$ at a fixed time % determines fully the future solution ${\bf y}(t)$.\footnote{A % nonautonomous system of size $n$ can, however, be written as an autonomous % one of size $n+1$ by introducing a new dependent variable $\tau$ % equivalent to $t$, governed by the equation $\tau' = 1$ and % the initial value $\tau(0) = 0$. See Figure 9.9.} A fundamental property % of phase space analysis is that % {\em distinct trajectories can never meet} % (assuming $\bf f$ is continuous with respect to $t$ and $\bf y$ and % Lipschitz continuous with respect to $\bf y$). % This follows from the uniqueness of solutions, for % if an autonomous ODE has two solutions ${\bf y}_1^{}(t)$ and ${\bf y}_2^{}(t)$ % such that ${\bf y}_1^{}(t_1^{}) = {\bf y}_2^{}(t_2^{})$ % for some $t_1^{}$ and $t_2^{}$, then if solutions are unique, % we must have ${\bf y}(t_1^{}+t) = {\bf y}_2^{}(t_2^{}+t)$ for % all $t$, positive and negative. So the trajectories corresponding % to $\bf y_1^{}$ and $\bf y_2^{}$ must be the same. % %% % % For another example, let us consider % the two-variable system known as the {\bf predator--prey} or % {\bf Lotka--Volterra} equations. % We suppose an environment in which there is a prey species, say rabbits, % whose population as a function of time is $u(t)$, and a predator % species, say foxes, whose population as a function of time is % $v(t)$.\footnote{Volterra's original application of 1926 % involved fish populations in the Adriatic % Sea. Danby writes in {\em Computer Modeling: From Sports % to Spaceflight$\dots$ From Order to Chaos:} Presentations of % this model include many pairs: birds and worms, wolves and deer, % sharks and tuna, to name a few. I have found that my students prefer % foxes and rabbits, so I shall refer to these. After all, who % identifies with a worm?''} % In the absence of foxes, the rabbit population grows % exponentially thanks to abundant vegetation, but it is decreased by % fatal encounters with foxes at a rate proportional to % the product $uv$, % $$u'= u - uv. \eqno (10.4)$$ % The fox population, on the other hand, decays exponentially in the % absence of rabbits to eat but expands when they are available, % $$v'= -\textstyle{1\over 5}v + uv. \eqno (10.5)$$ % (The constant $1/5$ is included to make the dynamics nontrivial. % In a serious application % one would have dimensional constants in front of each term, based % on measurements of field data.) % Combining the two equations with initial conditions, % we obtain the Lotka--Volterra IVP % $$u' = u-uv, ~ v'= -\textstyle{1\over 5}v + uv, \quad t\ge 0, ~ u(0) = u_0^{}, % ~v(0) = v_0^{}. \eqno (10.6)$$ % %% % % \noindent % A quiver plot looks like this: % \noindent % N = chebop(0,60); N.lbc = [1;1]; N.op = @(t,u,v) [diff(u)-u+u*v; diff(v)+.2*v-u*v]; quiver(N,[0 1.5 0 3],'k','xpts',12,'ypts',8) xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT) title('Fig.~10.4.~~Quiver plot for Lotka--Volterra equations (10.6)',FS,11) %% % % This problem has periodic solutions, % where the rabbit population oscillates between large and (very) small and % the fox population oscillates between large and (not so) small, % with a phase lag since a large rabbit population leads to a growing % fox population. Here is the periodic solution % that evolves from the initial populations $u(0)=v(0)=1$. % [u,v] = N\0; plot([u v]) title(['Fig.~10.5.~~Lotka--Volterra ' ... '(= predator--prey) equations'],FS,11) text(40,1.8,'foxes',IN,LT,FS,9) text(36.5,0.7,'rabbits',IN,LT,FS,9) hold on, plot([0 0;60 60],[.2 1; .2 1],'--k',LW,.6), hold off %% % \vskip 1.02em %% % % \noindent % The period is about 17.5. % [~,m] = max(v,'local'); period = diff(m) %% % % To understand the Lotka--Volterra dynamics, again it is helpful to % interpret the nonlinear terms of (10.4)--(10.5) as linear terms % with varying coefficients. Equation (10.4) can be written % $u' = (1-v) u$, showing that locally, % $u$ will decrease exponentially for $v>1$ % (the fox curve lies above the upper dashed line) % and increase exponentially for $v<1$ (it lies below % the upper dashed line). % Thus the rabbit population has an exponential % growth tendency that shuts off as $v$ rises above~$1$. % Similarly, equation (10.5) can % be written $v' = (u-1/5) v$, showing that the fox population has % an exponential decay tendency but grows % whenever there are enough rabbits, $u>1/5$ (the rabbit curve lies % above the lower dashed line). % %% % % Here is the orbit just computed in the % phase plane, winding around counterclockwise several times. % arrowplot(u,v,CO,ivpnl), hold on quiver(N,[0 1.5 0 3],'k'), axis([0 1.5 0 3]) title('Fig.~10.6.~~Lotka--Volterra cycle in the phase plane',FS,11) xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT) plot([.2 0;.2 1.5],[0 1; 3 1],'--k',LW,.6), hold off %% % \vskip 1.02em %% % % \noindent % Let us superimpose several orbits in a single plot, corresponding to % the solutions arising % from various initial populations of rabbits. % for s = .4:.2:1.4 N.lbc = [s; 1]; [u,v] = N\0; arrowplot(u,v,LW,.8,CO,ivpnl,MS,4,YS,.5), hold on end title('Fig.~10.7.~~Superposition of six orbits',FS,11) axis([0 1.5 0 3]) xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT) plot([.2 0;.2 1.5],[0 1; 3 1],'--k',LW,.6), hold off %% % \vskip 1.02em %% % % \noindent % The figures show that in the Lotka--Volterra model, the % effects of the initial conditions last forever: % some fox-rabbit oscillations % have great amplitude while others move little. At the % center is obviously a fixed point, which we may confirm % by setting $u'=v'=0$ in (10.6). This gives the equations % $$u = uv, \quad \textstyle{1\over 5} v = uv,$$ % implying that these equations in fact % have two fixed points, one at $(1/5,1)$ and the other at $(0,0)$. % Let us mark them in the plot. % hold on, plot([0 .2],[0 1],'.k',MS,8), hold off title('Fig.~10.8.~~Same with fixed points added',FS,11) %% % \vskip 1.02em %% % % The Lotka--Volterra equations are perhaps the most famous % system of two ODE\kern .5pt s, appearing in every ODE textbook. % Turning now to three dimensions, it is only % fitting that we should begin with the most famous system of % three ODE\kern .5pt s, the {\bf Lorenz equations.} % In 3D, a new possibility opens up for reasons we shall % explain in a moment: the dynamics of an ODE % may be {\em chaotic}. % The Lorenz equations are the archetypical chaotic system, % consisting of three coupled first-order equations, % $$u' = 10(v-u), \quad v' = u(28-w)-v, \quad w' = uv -(8/3)w. % \eqno (10.7)$$ % The coefficients $10$, $28$, and $8/3$ are arbitrary, and it is % not necessary to use exactly these values. However, so much % study of the Lorenz equations has been based on % these parameter choices that they have become standard. % %% % % Equations (10.7) may look complicated at first % glance, but in fact, the nonlinearities here are rather % simple: the equation for $u'$ is linear, and the % equations for $v'$ and $w'$ are quadratic, their only % nonlinearities consisting of the term $-uw$ in the equation for $v'$ % and the term $uv$ in the equation for $w'$. % To make an IVP, we impose initial conditions such as % $$u(0)= v(0) = -15, ~~ w(0) = 20. \eqno (10.8)$$ % Here is the solution of (10.7)--(10.8) up to time % $t=30$, showing the three components together. % 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 v w]) title('Fig.~10.9.~~Lorenz equations (10.7)--(10.8)',FS,11) %% % \vskip 1.02em %% % % \noindent % To focus on just one component of the solution, here is a plot of $v(t)$. % Note the characteristic property of chaos: apparent % randomness in a deterministic system. We shall discuss % chaos more fully in Chapter~13. % plot(v,CO,ivpnl,LW,.7) title('Fig.~10.10.~~Component $v(t)$',FS,11) %% % \vskip 1.02em %% % The Lorenz equations are autonomous, so they can be % represented in a phase space, which will have three dimensions. % Here is a phase space plot of the trajectory just computed. % The plot includes dots showing the fixed points, which are the % solutions of the equations we obtain by setting % $u'=v'=w'=0$ in (10.7), % $$v = u, \quad 28u - uw =v, \quad uv = 8w/3.$$ % There are three fixed points, one at % $(u,v,w) = (0,0,0)$ and the others at % $$(u,v,w) = (\pm 6\sqrt 2, \kern 1pt \pm 6\sqrt2, \kern 1pt 27).$$ %% % \vskip -2em plot3(u,v,w,CO,ivpnl,LW,.4) xlabel('$u$',FS,10,IN,LT), ylabel('$v$',FS,10,IN,LT) zlabel('$w$',FS,10,IN,LT) title('Fig.~10.11.~~3D phase space for the Lorenz eqs.',FS,11) c = 6*sqrt(2); hold on, plot3([0 c -c],[0 c -c],[0 27 27],'.k',MS,6) view(-35,30), axis([-17 19 -25 25 0 45]), hold off set(gcf,'position',[380 220 540 300]) %% % \vskip 1.02em %% % % Ideas related to phase spaces and their trajectories are % at the heart of the rich subject of {\em dynamical % systems,} and a great deal is known about the % mathematics and application of these notions. % The traditional method of this analysis is to begin % by finding fixed points and then examine the behavior near % these points by approximating the nonlinear problem locally % by a linear problem. As we shall see in Chapter 14, % linear analysis makes sense at any point of a phase % space, not just a fixed point; locally, any smooth nonlinear % ODE might as well be linear. What makes the fixed points % particularly interesting is that trajectories interact in % nontrivial ways at such points, and to see the % overall dynamics in a phase space, it is often fruitful to % start from the fixed points and see how they may be connected % by various trajectories. See Chapter 15. % %% % % One of the fundamental things one can learn from this % kind of geometrical reasoning is that trajectories can get tangled up in % 3D, but not in 2D.\ \ In 2D, since trajectories cannot cross, tangling % is not possible. Therefore, chaos can appear in autonomous % systems with $n\ge 3$ but never with $n=2$.\footnote{A % precise formulation of this idea is provided by the % {\em Poincar\'e--Bendixson theorem.}} Another point % is that since a linear system of ODE\kern .5pt s behaves in essentially % the same way at every point of phase space, its trajectories look essentially % the same everywhere. So they can't get tangled up either, and % chaos can never appear in a linear problem, no matter % how many variables there are. Since we have not % given a mathematical definition of chaos, these observations % must be regarded as heuristic, but they can be made precise. % %% % % Let us finish with an example of a three-variable problem % that is not chaotic, analyzed in pp.~202--204 of % Bender and Orszag, {\em Advanced Mathematical Methods for % Scientists and Engineers,} % $$u' = vw, ~~ v' = -2uw, ~~ w' = uv. \eqno (10.9)$$ % It can be shown that along any trajectory, the quantity % $u(t)^2 + v(t)^2 + w(t)^2$ is constant, and thus trajectories % lie on spheres. Here is a plot showing some of them. % close all, N = chebop(0,10); N.op = @(t,u,v,w) [diff(u)-v*w; diff(v) + 2*u*w; diff(w)-u*v]; for theta = -1.5:.2:1.5 N.lbc = [cos(theta); sin(theta); 0]; [u,v,w] = N\0; plot3(u,v,w,CO,purple,LW,.5), hold on plot3(-u,v,w,LW,.5,CO,blue), plot3(w,v,u,CO,green,LW,.5) plot3(w,v,-u,'k',LW,.5) end plot3([0 0 0 0 1 -1],[0 0 1 -1 0 0],[1 -1 0 0 0 0],'.k',MS,8) xlabel('$u$',FS,5,IN,LT), ylabel('$v$',FS,5,IN,LT) zlabel('$w$',FS,5,IN,LT), set(gca,FS,6), axis equal, xlim([-2 2]) title('Fig.~10.12.~~Phase space for (10.9)',FS,11), hold off set(gcf,'position',[380 220 540 300]) %% % \vskip 1.02em %% % % \noindent % As explained by Bender and Orszag, this example is related to % the mechanical phenomenon known as {\em Eulerian wobble} and can % be used to explain why if you toss a book % in the air, you can make it spin stably around % the shortest axis or the longest, but not the middle one: it will % somersault. This is also called the {\em tennis racket theorem} % or the {\em Dzhanibekov effect.} % %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: SIR model for epidemics}\\[-3pt] % \hrulefill % \end{center} % %% % % A major advance in epidemiology dates to a 1927 study by W. O. Kermack and % A. G. McKendrick % of a 1906 outbreak of bubonic plague in Bombay, now Mumbai % (A contribution to the mathematical % theory of epidemics,'' {\em Proceedings of the Royal Society of London A}.) % Kermack and McKendrick proposed % what is now called the {\em SIR model\/} for the spread of % disease. Although the SIR model % is now considered too simplistic for quantitative % tracking of real epidemics, it exhibits several important properties % and remains the foundation for more comprehensive models. % %% % % This is a {\em compartment model,} in which one imagines % members of a population making the transition continuously from a box % labeled S for % {\em susceptible} to another box~I for {\em infected} and eventually % to a box R for {\em recovered} or {\em removed} (possibly by death). % In nondimensional form the system is % $$S' = -\beta S I, \quad I' = (\kern .5pt % \beta S-\gamma)I, \quad R' = \gamma I, % \eqno (10.10)$$ % where $\beta,\gamma>0$ are parameters % and $S(t)$, $I(t)$, and $R(t)$ denote the proportions of the population % in the three compartments. % These equations model the dynamics of the % transitions ${\rm S}\to {\rm I}$ and ${\rm I}\to {\rm R}$. % First, $S$ decreases and $I$ increases at the rate % $\beta S I$ (so $S$ is monotonically nonincreasing). % This product corresponds to infected people % coming into contact with susceptible ones with % infection rate~$\beta$ per contact (high for measles, low for AIDS). % Second, $R$ increases and~$I$ decreases % at the rate $\gamma I$ (so $R$ is monotonically % nondecreasing). This represents the approximation that % infected people come to the end of their infections, whether % through death or recovery, in a random fashion at a steady rate (again high % for measles, low for AIDS). % %% % % Note that (10.10) is nonlinear because of the product % $SI$. Note also how easily we can spot factors that % are omitted from the model. For example, (10.10) says nothing % about the course of the infection; an infected patient is % simply assumed to recover or die at random with a % probability that does not vary from day to day. Also, the % model assumes the populations are perfectly mixed, % with no spatial or social aspect to the epidemic. In actuality, % spatial and socioeconomic factors can be very important, % as was made famous by John Snow's study of the cholera outbreak % linked to a certain public water pump in a crowded quarter % of London in 1854. % %% % The equations (10.10) have two obvious steady states. % If $S=1$ and $I=R=0$, nobody is infected and the disease % never gets going. If $S=I=0$ and $R=1$, the whole % population is recovered or removed and % the epidemic has passed into history. % The nontrivial dynamics emerges when we begin with a nonzero % infected fraction $I_0^{}$, % $$S(0)=1-I_0^{}, \quad I(0)=I_0^{}, \quad R(0)=0.$$ % Here we show the course of the epidemic with $\beta = 2$, % $\gamma = 1$, and $I_0^{} = 0.0001$. beta = 2; gamma = 1; I0 = 0.0001; N = chebop([0 30]); close all N.op = @(t,S,I,R) [ diff(S) + beta*S*I diff(I) - beta*S*I + gamma*I diff(R) - gamma*I ]; N.lbc = [1-I0; I0; 0]; [S,I,R] = N\0; plot([S,I,R]), hold on, plot(S+I+R,'--k'), hold off xlabel('time $t$',IN,LT,FS,9), ylabel('proportion of population',FS,9) title('Fig.~10.13.~~An SIR epidemic (10.10)',FS,11) text(19.3,.07,'infected, $I(t)$',FS,9,IN,LT) text(19.3,.27,'susceptible, $S(t)$',FS,9,IN,LT) text(19.3,.86,'recovered or removed, $R(t)$',FS,9,IN,LT) text(19.3,1.06,'total $S(t)+I(t)+R(t)$',FS,9,IN,LT) %% % \vskip 1.02em %% % % The model may be simple, but there are % some remarkable things to see in these curves. % First of all, note how a disease % may be active in a small fraction of a population for % a long time, exponentially growing but at a rate so % low that it may go unnoticed.\footnote{What do you % suppose is happening in a carton of milk in your % refrigerator when its best before'' date is still a week away?} % %% % Why does it grow at all? As mathematicians, we may say % that the fixed point $(S,I,R) = (1,0,0)$ is % unstable. Looking at (10.10), we see that the equation % giving exponential growth is the second one, % $I' = (\beta S-\gamma) I$. If $\beta S > \gamma$, the % coefficient is positive and we can expect $I$ to grow exponentially % at the rate $\beta S - \gamma$. For our choice of parameters, this is indeed % the case at $t=0$, with $\beta S - \gamma = 1$. % As epidemiologists, we can interpret this condition as follows: % each newly infected person, on average, infects more than 1 additional % person before ceasing to be infective. %% % Eventually, on the other hand, the epidemic shuts off. As soon % as $S$ falls to the point where $\beta S\le \gamma$, $I(t)$ begins % to diminish and the epidemic starts fading out. % In this example this occurs at the point % $t_{\rm c}$ where $S(t_{\rm c}^{}) = 0.5$. tc = roots(S-0.5) %% % % \noindent % Note that although the infected population begins % to diminish after this point, the epidemic has by % no means yet run its course. In the end, most of % the population has had the disease: % Rinfty = R(end) %% % % \noindent % But not all. This is a fascinating and important aspect % of epidemiology. For a population to be safe from a certain % pathogen, it is not necessary that everybody be immune. It suffices % for a high enough fraction to be immune so that further infections % will not spread. This is the celebrated effect % of {\em herd immunity.} And here we realize that the two obvious'' % fixed points of this system % noted earlier are not the only ones. In fact, % any values $S_\infty^{}$ and % $R_\infty^{} = 1-S_\infty^{}$ in $[\kern .3pt 0,1]$ % may be combined % with $I_\infty^{} = 0$ to make a fixed point. % %% % As $\beta$ increases, the peak of $I(t)$ comes both sooner and % higher: for beta = [1.3 1.6 2 2.8] N.op = @(t,S,I,R) [ diff(S) + beta*S*I diff(I) - beta*S*I + gamma*I diff(R) - gamma*I ]; [S,I,R] = N\0; plot(I), hold on, [val,pos] = max(I); text(pos,val+.02,['$\beta =$ ' num2str(beta)],FS,9,IN,LT,HA,CT) end xlabel('time $t$',FS,9,IN,LT), ylabel('infected proportion',FS,9) ylim([0 .35]) title('Fig.~10.14.~~SIR epidemic for different infection rates',FS,11) %% % \vskip 1.02em %% % % Epidemiological equations are not just used to track biological diseases. % Similar compartment models have been used to study the spread of % computer malware as well as other items that travel across % networks --- like a rumor, or a meme. See our favorite reference % at the end of Chapter~15. % %% % % \smallskip % {\sc History.} % Systems of differential equations were studied by % d'Alembert and Lagrange in the mid-18th century, % but vector notation was not yet % available in those days. Vector formulations % became customary with the work of Peano and others near % the end of the 19th century. % %% % % \smallskip % {\sc Our favorite reference.} J. D. Murray's wide-ranging book % {\em Mathematical Biology\/} % did much to create the field of mathematical % biology as we know it today. In its later extension to two volumes % (Springer, 2001), the book focuses in volume 1 on systems of % ODE\kern .5pt s modeling biological systems, including insect % populations, fishery management, tumor cell growth, % divorce rates, biological oscillators, % propagation of neural signals, and much more. The biology is % varied, and the dynamical effects in the ODE\kern .5pt s are equally varied, combining % to give a deeply satisfying illustration of the power of % mathematics in understanding our world. % \smallskip % %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 10.} Any system of ODE\kern .5pt s can be % reduced to a first-order system of the form % ${\bf y}'(t) = {\bf f}(t,{\bf y})$. % An $n$-variable autonomous system % ${\bf y}'(t) = {\bf f}({\bf y})$ has an $n$-dimensional % phase space. If $n\ge 3$ and the equations are % nonlinear, there is the possibility of chaos. % \vspace{2pt}}} % \end{displaymath} % %% % \smallskip\small\parskip=1pt\parindent=0pt % {\em \underline{Exercise $10.1$}. Four bugs on a rectangle.} % In Chapter 3, using the complex arithmetic trick, we % considered a scalar ODE in which a lion chases an % antelope, always moving directly towards it. % Now consider a problem of four bugs starting at four % corners of a rectangle, each chasing the one counterclockwise % from it at speed $1$. % {\em (a)} Suppose the bugs begin at the corners of % the square $(\pm 1,\pm 1)$. % Draw a figure showing the evolution to time $t=1.99995$, and % report the locations at that time of the four bugs. % (This part of the problem could be reduced % to a scalar ODE --- indeed it can easily be % solved analytically --- but the next part is necessarily a system.) % {\em (\kern .7pt b)} Now suppose the bugs begin at the four corners % of the rectangle $(\pm 1.01, \pm 0.99)$. Answer the % same questions as in part {\em (a)}. To learn about % the extraordinary paths followed by the bugs as $t\to 2$, see % Chapman, Lottes, and Trefethen, Four bugs on a rectangle,'' % {\em Proceedings of the Royal Society,} 2011. % \par % {\em \underline{Exercise $10.2$}. Oregonator.} % The Oregonator, or Field--Noyes model, is an approximation to the % dynamics of the Belousov--Zhabotinskii chemical reaction, which exhibits % periodic oscillations of its reagents. (You can see a video of the reaction at % {\tt https://youtu.be/wxSa9BMPwow?t=3m9s}.) In nondimensional % form the equations are % $$\varepsilon u' = qv - uv + u(1-u),\quad % \delta v' = -qv - uv + w, \quad w' = u-w,$$ % where $q$, $\varepsilon$, and $\delta$ are parameters. % For this exercise use the values $q=8\times 10^{-4}$, $\varepsilon=0.05$, % and $\delta=0.02$. (These are not physically realistic for the BZ % reaction. For realistic values, the system is considered stiff'' and % requires a change in the numerical solver Chebfun uses, as illustrated % in equation 72 of Appendix~B.)\ \ Solve % the system with Chebfun for $u(0)=w(0)=0$, $v(0)=1$, and plot % the three components as functions of time. (A particularly interesting way % to see the dynamics is {\tt comet3(u,v,w)}.) % \par % {\em Exercise $10.3$. Nonautonomous.} % Suppose that a first-order, two-variable ODE system has unique % solutions for any initial data and one of its solutions is % $u(t)= \sin(t),$ $v(t) = \sin(2t)$. % Sketch this orbit in the phase plane and explain % why this ODE cannot be autonomous. % \par % {\em \underline{Exercise $10.4$}. Matrix eigenvalues via ODE isospectral flow.} % Suppose ${\bf A} = [\, a~b\,;~b~c\, ]$ is a $2\times 2$ real symmetric matrix. % The eigenvalues of ${\bf A}$ can be determined from % $\lim_{t\to\infty}^{}a(t)$ and $\lim_{t\to\infty}^{}c(t)$ % in the ODE $a' = 2\kern .3pt b^2$, $b' = b(c-a)$, $c' = -2\kern .3pt b^2$. % {\em (a)} Try this for the matrix $[\, 3~2\,;~2~1\, ]$, plotting % deviations from their limiting values % as a function of $t\in [\kern .3pt 0,2\kern .2pt ]$. % How close are $a(t)$ and $c(t)$ at $t = 1$ and $2$ to their % final values? % {\em (\kern .7pt b)} Do the same for the % matrix $[1~2\,;~2~3]$. What is qualitatively different in this case? % {\em (c)} What are the fixed points of this ODE?\ \ (If % you are curious to learn more about isospectral % flows in linear algebra, including the case of $n\times n$ matrices % and connections with standard algorithms for computing matrix % eigenvalues, % see Deift, Nanda, and Tomei, Ordinary differential equations and % the symmetric eigenvalue problem,'' {\em SIAM Journal on Numerical % Analysis,} 1983.) % \par % {\em Exercise $10.5$. Energy conservation for Eulerian wobble.} % Show that the quantity $u^2+v^2+w^2$ is independent of time for % solutions of the system (10.9). % \par % {\em \underline{Exercise $10.6$}. Four kinds of fixed points} % (from Hairer, N\o rsett, and Wanner). % {\em (a)} Draw a quiver plot of the system $u' = (u-v)(1-u-v)/3$, % $v' = u(2-v)$. Determine the four fixed points analytically and % include them in the plot. % {\em (\kern .7pt b)} Calculate the eigenvalues of the Jacobian at these points % and show that one is a sink, one is a source, one is a saddle, % and one is a spiral. % \par % {\em \underline{Exercise $10.7$}. Phase plane for SIR model.} % {\em (a)} Suppose a solution of (10.10) satisfies % $S(0)+I(0)+R(0)=1$. Show that % $S(t)+I(t)+R(t) = 1$ for all $t$. % Using this conservation property, we can eliminate $S$ or $I$ or % $R$ from the system and thus obtain a two-variable system that can % be analyzed in the phase plane. % {\em (\kern .7pt b)} Elimination of $R$ is particularly trivial. % Write down the resulting two-variable system. Solve it % for the same parameters % as in Figure~10.13 and make an arrowplot of the trajectory in % the $S$-$I$ plane. Where are the fixed points in this plane? %