%% 13. Chaos % % \setcounter{page}{156} % %% % % Having just considered functions that are truly random, % at least in principle, we now turn to the celebrated phenomenon of % solutions of ODE\kern .5pt s that look random'' % even though they are not. Besides this property, % two other hallmarks of chaos % are {\em sensitive dependence on initial conditions} and % {\em strange attractors.} In this % chapter we explore these notions, connecting % them with the number known as the Lyapunov exponent. % %% % % We begin with the Lorenz equations as given in equation (10.7), % $$u' = 10(v-u), \quad v' = u(28-w)-v, \quad w' = uv -(8/3)w. % \eqno (13.1)$$ % In Figure 10.11 we plotted the trajectory for $t\in [\kern .3pt 0,30\kern .3pt]$ % emanating from the initial data % $$u(0)= v(0) = -15, ~ w(0) = 20. \eqno (13.2)$$ % Here is an image of this trajectory in phase space, seen from % the angle that amounts to a projection onto the $u$--$w$ plane. % ODEformats, N = chebop(0,30); N.op = @(t,u,v,w) [diff(u)-10*(v-u); ... diff(v)-u*(28-w)+v; diff(w)-u*v+(8/3)*w]; N.lbc = [-15; -15; 20]; [u,v,w] = N\0; plot(u,w,CO,ivpnl,LW,0.4), axis([-25 25 0 50]) c = 6*sqrt(2); hold on, plot([0 c -c],[0 27 27],'.k',MS,8), hold off xlabel('$u$',FS,10,IN,LT), ylabel('$w$',FS,10,IN,LT) title(['Fig.~13.1.~~Lorenz trajectory projected '... 'onto $u$--$w$ plane'],FS,11) %% % \vskip 1.02em %% % % \noindent % This image reveals the famous butterfly'' structure of % the {\bf strange attractor} for the Lorenz equations.\footnote{Not % to be confused with the {\em butterfly effect\/} that also % originates with Lorenz (not his famous 1963 paper, but a later one % from 1969). % If a butterfly flaps its wings in the Amazon, thanks % to sensitive dependence of the earth's weather on initial conditions, % might that cause a hurricane in Texas?} % What it means to be a strange attractor is that typical % orbits rapidly approach this % set, which is not a simple curve or surface but % has the form of a fractal. For the Lorenz equations, the % fractal dimension is approximately $2.06$, so % the butterfly is just a little bit thicker than a % 2-dimensional manifold.\footnote{The {\em fractal\/} or % {\em Hausdorff dimension\/} of a set is a precisely % defined mathematical quantity that may take any % nonnegative real value, not just an integer.} % %% % Now let us look at dependence on initial conditions. Here % the function $v(t)$ just computed is plotted in green, as in % Chapter 10. In addition, another function $\tilde v(t)$ is plotted % in red that arises from the very slightly perturbed initial data % $$u(0)= v(0) = -15, ~ w(0) = 20.00001.\eqno (13.3)$$ %% % \vskip -1.5em N.lbc = [-15; -15; 20.00001]; [u2,v2,w2] = N\0; plot(v2,'r'), hold on, plot(v,CO,ivpnl), hold off title(['Fig.~13.2.~~Lorenz equations: ' ... 'sensitive dependence on initial data'],FS,11) %% % \vskip 1.02em %% % % \noindent % The two initial conditions differ by less than % one part in $10^6$, so one would expect $v$ and $\tilde v$ to % be close to one another, and for about 13 time units, so they are. % But throughout that time, the trajectories % are steadily (on average) separating, and % eventually the differences are of size $O(1)$ and % the functions become completely uncoupled. % An initial difference of size $10^{-5}$ has grown % by six orders of magnitude in 15 time units. % %% % Everything in this problem is smooth, differentiable, and % indeed analytic. For % example, the derivative of the value $v(30)$ with respect % to the third initial condition $w(0)$ is a perfectly well-defined % number --- it is just that it is a very big one! Here is % an approximation to % the derivative of $v(5)$ with respect to $w(0)$, (v2(5)-v(5))/(w2(0)-w(0)) %% % % \noindent % Here is the analogous approximation to % the derivative of $v(10)$ with respect to $w(0)$, % (v2(10)-v(10))/(w2(0)-w(0)) %% % % \noindent % Note that these numbers are rapidly growing with $t$. % If we compute the corresponding finite difference at $t=30$, we get % (v2(30)-v(30))/(w2(0)-w(0)) %% % % \noindent % but this greatly underestimates the actual derivative % of $v(30)$ with respect to $w(0)$, % which is on the order of $10^{13}$. We can explain the % underestimation by noting that the % initial perturbation by $0.00001$ is not nearly small enough to % grow by 13 orders of magnitude before reaching size~$O(1)$, % at which point it stops growing. % %% % % To see more of this effect, it is interesting to plot % $\|\tilde{{\bf y}}(t) - {\bf y}(t)\|$ on a log scale, % where ${\bf y}$ and $\tilde {\bf y}$ are the % 3-component vector solutions resulting % from the original and the perturbed initial conditions. % (The norm $\|\cdot\|$ we use is % the square root of the sum of the squares of the components % $u(t), v(t), w(t)$.) % Here is such a plot, and it reveals two distinct regimes. % For $t<15$, the perturbations grow exponentially. % For $t>15$, they reach size $O(1)$ (numbers roughly in the % range $10$--$100$) and stop growing. % tt = linspace(0,30,400); err = sqrt((u2(tt)-u(tt)).^2+(v2(tt)-v(tt)).^2+(w2(tt)-w(tt)).^2); semilogy(tt,err,'k'), ylim([1e-6 1e3]) title(['Fig.~13.3.~~Exponential growth of perturbations ' ... 'to scale $O(1)$'],FS,11) xlabel('$t$',FS,10,IN,LT) ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,10,IN,LT) %% % \vskip 1.02em %% % % \noindent % The growth rate of this initial phase is known as the % {\bf Lyapunov exponent,}~$\lambda$, for this ODE.\ \ We % can approximate it % by a least-squares fit to the data for $t\in[\kern .3pt 0,15]$. % ii = find(tt<=15); c = polyfit(tt(ii),log(err(ii)),1); e0 = c(2); lambda = c(1) %% % % \noindent % The significance of this number is that, approximately % speaking, perturbations % in the Lorenz trajectories typically grow initially at the rate % $$\| \tilde {\bf y}(t) - {\bf y}(t) \| \approx C e^{0.91 t}.$$ % We can plot this fit to the data as a dashed line. % hold on, semilogy([0 18],exp(e0+lambda*[0 18]),'--r',LW,1.5) hold off title('Fig.~13.4.~~Slope of line = Lyapunov exponent $0.91$',... FS,11,IN,LT) xlabel('$t$',FS,10,IN,LT) ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,9,IN,LT) %% % \vskip 1.02em %% % % In this experiment, we started from a finite perturbation, % which grew exponentially over a finite time span. The mathematical % definition of a Lyapunov exponent is based on % infinitesimal perturbations, which may grow forever. % To be precise, $\lambda$ is defined by this formula, involving % a supremum over all initial perturbations $\delta {\bf y}(0)$: % $$\lambda = \lim_{t\to\infty} \limsup_{\|\delta {\bf y}(0) \|\to 0} t^{-1} % \log {\|\delta {\bf y}(t)\|\over \|\delta {\bf y}(0)\|}. \eqno (13.4)$$ % In principle this number might vary from one initial point ${\bf y}(0)$ % to another, in which case one would typically be interested in the % supremum over all ${\bf y}(0)$.\footnote{A further complication is that % what we have called the Lyapunov exponent should more properly be termed % the {\em maximal\/} Lyapunov exponent. For an $n$-dimensional dynamical % system, imagine initial conditions in the form of an infinitesimal % spherical cloud or blob centered at a point ${\bf y}(0)$. As $t\to\infty$, % the blob will elongate and compress into an ellipsoid, with the % lengths of the axes of the ellipsoid growing or decaying at % different exponential rates on average % as $t\to \infty$ in a manner that can be made precise by the tool known % as the singular value decomposition. These growth and decay % rates are the Lyapunov exponents of the system.} % %% % The Lyapunov exponent for the Lorenz equations % is known to be about 0.9056. By chance, our % experiment has matched this number more closely than should really be % expected, statistically speaking (see Exercise 13.1). %% % % The article that launched the study of chaos % was published by Edward Lorenz in 1963 with the % title Deterministic nonperiodic flow,'' and is one of % the most important scientific publications of the 20th % century.\footnote{The term chaos'' came later, coined by Jim % Yorke in 1975.} % The fact that it took until after 1963 for the phenomenon to be % widely recognized can be attributed to the invention of computers. % Chaos is not an unusual behavior at all in nonlinear systems of % ODE\kern .5pt s, but it could not be easily seen in the closed form solutions % that were available in the days of hand calculation. % %% % % The historical side of our next example is particularly % interesting. Among the most important ODE\kern .5pt s of all scientifically % are the equations of the % {\boldmath\bf $n$-body problem,} by % which we mean the equations that govern the % trajectories of $n$ point masses attracting % each other gravitationally and moving according to Newton's laws. % In its simplest form this is a set of $3n$ ODE\kern .5pt s of % second order: for each body, there are three spatial coordinates, % and the differential equations express the fact that the % acceleration of each body is equal to % the sum of the inverse-square attractions % to the others. % For problems confined to a plane, the count reduces to % $2n$ variables, an $x$ and a $y$ coordinate for each body. We are % going to look at such a problem, and for convenience we will employ % our usual trick of encoding $x$ and $y$ as a complex variable % $z=x+iy$, as first introduced % in Chapter~3. This brings us to a second-order system of % ODE\kern .5pt s in $n$ complex unknowns: % $$z_j'' = - \sum_{i\ne j} {z_j^{} - z_i^{}\over | z_j-z_i|^3}, % \quad 1 \le j \le n. \eqno (13.5)$$ % The analytic solution to the 2-body problem goes back to Kepler % and Newton, and it involves stable ellipses (as in Exercise 4.1). % But there is no analytic solution for the 3-body problem, let alone % for the $n$-body problem with $n>3$. % %% % % For a long time, mathematicians and physicists % have asked, is our solar system stable? Might % its regular orbits one day break down? This problem proved % theoretically intractable % but has led to a great deal of important mathematics along the way. % Even today, after extensive computational investigations, % there is some controversy % over whether or not the solar system % is stable or chaotic, but the difficulty of this question is a % result of a special circumstance: % the mass of the sun is much greater than that % of the planets, and acts as a strong regulating influence. % The sun's mass is so dominant that to first approximation, % each planet's orbit is just the solution of a 2-body problem involving % itself and the sun. % By contrast, if you look at the % trajectories of $n$ attracting bodies of equal masses % in the absence of a heavy sun to act as a policeman, chaos is the rule. % Before computers, nobody could look at such trajectories. % %% % The idealized configuration we shall examine % consists of three planets in a plane with no sun, % initially positioned with zero velocity at the % vertices of a 3-4-5 right triangle. % The following commands compute the orbits up to $t=150$. u0 = 0; v0 = 3; w0 = 4i; N = chebop(0,150); N.op = @(t,u,v,w) [ ... diff(u,2) + (u-v)/abs(u-v)^3 + (u-w)/abs(u-w)^3; ... diff(v,2) + (v-u)/abs(v-u)^3 + (v-w)/abs(v-w)^3; ... diff(w,2) + (w-u)/abs(w-u)^3 + (w-v)/abs(w-v)^3]; N.lbc = @(u,v,w) [u-u0; v-v0; w-w0; diff(u); diff(v); diff(w)]; [u,v,w] = N\0; %% % % \noindent % To begin with we plot the orbits just up to $t=20$. % plot(u,'interval',[0 20]), hold on, plot(v,'r','interval',[0 20]) plot(w,CO,ivpnl,'interval',[0 20]), hold off title('Fig.~13.5.~~Orbits of three planets to $t=20$',FS,11) axis equal, ylim([-2.5 5]) %% % \vskip 1.02em %% % % \noindent % In this image the red, green, % and blue curves show paths of the three planets % as~$t$ increases. At first they are well separated, % but soon two close red-blue flybys % are observed. Here are the distance and time of the closest approach: % [closest_distance,closest_time] = min(abs(u{0,20}-v{0,20})) %% % % Plotting the orbits up to $t=50$ gives a picture hard % to interpret in detail. It is clear, however, that % the planets are swinging around each other in an % effectively random fashion. % plot(u,'interval',[0 50]), hold on, plot(v,'r','interval',[0 50]) plot(w,CO,ivpnl,'interval',[0 50]), hold off title('Fig.~13.6.~~Orbits continued to $t=50$',FS,11) axis equal, ylim([-2.5 5]) %% % \vskip 1.02em %% % % One might expect this behavior % to continue forever, but in fact, it doesn't. % Around $t=86$, the three-body system breaks up (self-ionizes''), % with the red planet shooting off to the northwest and % the green and blue ones spiraling off together to the southeast. % This curious effect is % revealed in the plot below of the orbits up to $t=100$, looking like % ribbons on a birthday present. % Thus the 3-body problem for an initial 3-4-5 triangle turns % out to exemplify what is known as {\em transient chaos.} % plot(u,'interval',[0 100]), hold on plot(v,'r','interval',[0 100]), plot(w,CO,ivpnl,'interval',[0 100]) title(['Fig.~13.7.~~Orbits to $t=100$, ' ... 'with self-ionization near $t=86$'],FS,11) axis equal, ylim([-2.5 5]), hold off %% % \vskip 1.02em %% % % \noindent % Here, for example, is the time at which the red planet % passes through coordinate $x=-2$: % t_escape = roots(real(v)+2) %% % Just as we did with the Lorenz equations, let us examine % the effect of a small perturbation in one component of % the initial data. % Changing the initial coordinate $w(0)$ from $4i$ to % $3.9999i$ gives the following orbit up to $t=100$, % which is similar to the former one for a while but then % begins to differ completely. % This time there is no self-ionization. N.lbc = @(u,v,w) [u-u0; v-v0; w-3.9999i; diff(u); diff(v); diff(w)]; [u2,v2,w2] = N\0; plot(u2,'interval',[0 100]), hold on plot(v2,'r','interval',[0 100]) plot(w2,CO,ivpnl,'interval',[0 100]) hold off, axis equal, ylim([-2.5 5]) title('Fig.~13.8.~~Slightly perturbed orbits to $t=100$',FS,11) %% % \vskip 1.02em %% % % \noindent % A plot shows $\|\tilde {\bf y}(t)-{\bf y}(t)\|$ as a function of $t$ % and the estimated Lyapunov exponent $\lambda\approx 0.15$. % tt = linspace(0,150,400); err = sqrt(abs(u2(tt)-u(tt)).^2 + abs(v2(tt)-v(tt)).^2 + ... abs(w2(tt)-w(tt)).^2); semilogy(tt,err,'k'), ylim([1e-5 1e3]), ii = find(tt<=70); title(['Fig.~13.9.~~Exponential growth of ' ... ' perturbations to scale $O(1)$'],FS,11) c = polyfit(tt(ii),log(err(ii)),1); e0 = c(2); lambda = c(1) hold on, semilogy([0 80],exp(e0+lambda*[0 80]),'--r',LW,1.5) hold off, xlabel('$t$',FS,9,IN,LT) ylabel('$\|\tilde{\bf y}(t)-{\bf y}(t)\|$',FS,9,IN,LT) %% % \vskip 1.02em %% % % \noindent % Note the transition from the chaotic exponential % phase to smooth algebraic growth at a rate $O(t)$ % after the bodies have separated. % %% % % Our third example of a chaotic system is the % {\bf R\"ossler equations,} % $$u' = -v-w, ~~ v' = u+\textstyle{1\over 5} v, ~~ % w' = {1\over 5} + w(u-c), \eqno (13.6)$$ % where the constant $1/5$ has been fixed % arbitrarily and $c$ is a parameter. % This system is even simpler than the Lorenz % equations (13.1) in that only one of the three equations is nonlinear. % The R\"ossler equations illustrate % the phenomenon of {\em period doubling} as a route to chaos. % First, here we solve (13.6) for $t\in [\kern .3pt 0,300\kern .3pt]$ % with $c=2$ with initial conditions % $u(0)=2$ and $v(0) = w(0) = 0$. % The image on the left shows the trajectory in the three-dimensional % $u$-$v$-$w$ phase space, which settles down to a regular % oscillation, a limit cycle. The image on the right shows % the same trajectory projected onto the $u$--$v$ plane, just % the part for $t>200$. The initial transient % has died away, and thus all one sees is the 2D projection % of the limit cycle. % N = chebop(0,300); N.lbc = [2; 0; 0]; N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-2)]; [u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,LW,.3,CO,ivpnl,LW,.5) xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10) subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.7) axis([-10 10 -10 10]), axis square text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10) title(['Fig.~13.10.~~R\"ossler system (13.6), $c=2$ ' ... '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT) %% % \vskip 1.02em %% % % \noindent % As the parameter $c$ is increased, a change takes place, a bifurcation % (see Chapter~17). % At around $c=2.8$, the limit cycle undergoes {\bf period doubling}, % with the trajectory unfolding into a double loop. % The image for $c=3.5$ shows this period 2'' solution. % N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-3.5)]; [u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.3) xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10) subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.6) axis([-10 10 -10 10]), axis square text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10) title(['Fig.~13.11.~~R\"ossler system (13.6), $c=3.5$ ' ... '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT) %% % \vskip 1.02em %% % % \noindent % As $c$ increases further, more bifurcations take place. % A second period doubling to a period 4'' % solution occurs around $c=3.7$, and here is % a pair of images for $c=4$. % N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-4)]; [u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.3) xlabel('$u$',FS,10), ylabel('$v$',FS,10), zlabel('$w$',FS,10) subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.5) axis([-10 10 -10 10]), axis square text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10) title(['Fig.~13.12.~~R\"ossler system (13.6), $c=4$ ' ... '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT) %% % \vskip 1.02em %% % % \noindent % The phenomenon of period doubling was made famous in the 1970s by % Mitchell Feigenbaum, who showed that it is % a route to chaos in many problems. As $c$ increases further (the % next transition is around $c=4.1$), % the period doubles again and again, % infinitely often. Each time interval from one doubling to % the next is asymptotically $4.6692\dots$ times shorter than the last; this % is {\em Feigenbaum's constant.} Finally, for $c$ greater than % about $4.2$, % the system is chaotic. The following images for $c=5$ show % the chaotic regime, with the orbits settling down not to a limit % cycle but to a strange attractor. % N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-5)]; [u,v,w] = N\0; subplot(1,2,1), plot3(u,v,w,CO,ivpnl,LW,.2) subplot(1,2,2), plot(u{200,300},v{200,300},CO,ivpnl,LW,.4) axis([-10 10 -10 10]), axis square text(8.7,8,'$t>200$',FS,9,HA,RT), xlabel('$u$',FS,10), ylabel('$v$',FS,10) title(['Fig.~13.13.~~R\"ossler system (13.6), $c=5$ ' ... '\vrule height 16pt width 0pt depth 7pt'],FS,11,HA,RT) %% % \vskip 1.02em %% % % Our final example of chaos is governed by a % nonautonomous equation due to Moon and Holmes, % a chaotic {\em nonlinear forced oscillator}. % The equation is % $$y'' + \textstyle{1\over 4} y' - y + y^3 = 0.4 \cos(t). % \eqno (13.7)$$ % As always, the constants here are somewhat arbitrary and % could be adjusted. % This is just a second order ODE, so in the absence of the % nonautonomous forcing function, the trajectories could not % be chaotic. The forcing function, however, makes this % system equivalent to a first-order autonomous system in three % variables $y$, $y'$, and $t$ (see footnote 3 of % Chapter~10). Thus the dimension is % great enough for chaos to be a possibility, and a computation % confirms its appearance. % t = chebfun('t',[0 300]); N = chebop(0,300); N.op = @(t,y) diff(y,2) + .25*diff(y)-y+y^3; N.lbc = [0;0]; rhs = 0.4*cos(t); y = N\rhs; clf, plot(y,CO,ivpnl), ylim([-3 3]) title('Fig.~13.14.~~Nonlinear forced oscillator (13.7)',FS,11) %% % \vskip 1.02em %% % We have not defined chaos! This is because there is % no universally accepted definition, though the features of sensitive % dependence and strange attractors are important. Note % that exponential divergence of trajectories alone cannot % be enough, since the solutions of $y' = y$ diverge at % the rate $e^{\kern .3pt t}$ but nobody would call this equation chaotic. % One needs the exponential divergence to be combined with % some kind of global boundedness associated with nonlinearity, % and this is where the strange attractors come in. %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: chaos in a food web}\\[-3pt] % \hrulefill % \end{center} % %% % % In Chapter 10 we saw that a simple model of the interactions between % rabbits and foxes led to limit cycles in their populations. If we % introduce a third species into the fray, this can % transform the dynamics completely. This possibility was discovered % by Bob May (later Lord May) % in the 1970s, a story made famous in the % book {\em Chaos:\ Making a New Science} by James Gleick. % %% % % Rabbits need to eat too! Suppose in addition % to the populations of rabbits, $u(t)$, % and foxes, $v(t)$, we consider the population of carrots, $c(t)$. We % suppose that, left on their own, these tend to grow logistically % (see Exercises~3.15 and~3.16). Now % rabbits consume carrots, and foxes consume rabbits, and neither animal % species can be sustained without food. A reasonable model of the system is % $$c'=c(1-c)-f_1^{}(c)u, ~~ u'=f_1^{}(c)u - f_2^{}(u)v- d_1^{}u, ~~ v' % = f_2^{}(u)v - d_2^{} v, \eqno (13.8)$$ % with $f_i^{}(z)=a_i^{} z/(1+b_i^{}z)$ for $i=1,2$. The consumption interactions % are little different from those in the Lotka--Volterra equations. % They include a % saturation effect for the consumer species, accounting for extra % competition as populations grow large. This model was proposed by % Hastings and Powell in Chaos in a three-species food chain,'' % {\em Ecology}, 1991, where the parameter choices were % $a_1^{}=5$, $a_2^{}=0.1$, $b_2^{}=2$, $d_1^{}=0.4$, $d_2^{}=0.01$. The parameter % $b_1^{}$ was allowed to vary. % %% % Increasing $b_1^{}$ increases the effect of competition between rabbits. % Here is what happens to the rabbits with $b_1^{}=2.5$, % a nonchaotic regime. a1 = 5; a2 = 0.1; b2 = 2; d1 = 0.4; d2 = 0.01; f1 = @(z,b1) a1*z./(1+b1*z); f2 = @(z) a2*z./(1+b2*z); N = chebop(0,3000); b1 = 2.5; N.op = @(t,c,u,v) [ diff(c)-(c*(1-c)-f1(c,b1)*u); diff(u) - (f1(c,b1)*u-f2(u)*v-d1*u); diff(v) - (f2(u)*v-d2*v) ]; N.lbc = [.4;1;9]; [c,u,v] = N\0; plot(u{2000,3000},CO,ivpnl) title(['Fig.~13.15.~~Nonchaotic rabbit population, ' ... '$b_1=2.5$ in (13.8)'],FS,11,IN,LT) xlabel('$t$',FS,10,IN,LT), ylabel('rabbits, $u(t)$',FS,9) %% % \vskip 1.02em %% % % \noindent % At first glance this may look somewhat irregular. But in phase space we % see an ordinary limit cycle of doubled period, just % as in Figure 13.11. % plot3(c{2000,3000},u{2000,3000},v{2000,3000},LW,.7,CO,ivpnl), view(63,68) xlabel('carrots, $c(t)$',FS,9), ylabel('rabbits, $u(t)$',FS,9) zlabel('foxes, $v(t)$',FS,9) title('Fig.~13.16.~~In phase space, a limit cycle',FS,11) set(gcf,'position',[380 220 540 300]) %% % \vskip 1.02em %% % % \noindent % If we set $b_1^{}=3.5$, increasing % the competition between rabbits, the picture changes fundamentally. % b1 = 3.5; N.op = @(t,c,u,v) [ diff(c)-(c*(1-c)-f1(c,b1)*u); diff(u) - (f1(c,b1)*u-f2(u)*v-d1*u); diff(v) - (f2(u)*v-d2*v) ]; [c,u,v] = N\0; plot3(c{2000,3000},u{2000,3000},v{2000,3000},LW,.4,CO,ivpnl), view(63,68) xlabel('carrots, $c(t)$',FS,9), ylabel('rabbits, $u(t)$',FS,9) zlabel('foxes, $v(t)$',FS,9) title(['Fig.~13.17.~~Chaotic rabbit population, ' ... '$b=3.5$ in (13.8)'],FS,11,IN,LT) set(gcf,'position',[380 220 540 300]) %% % \vskip 1.02em %% % % \noindent % Looking again at the rabbit population, we see repeated but aperiodic % crashes in the population, separated by a variable number and size of % shorter cycles of boom and bust. % close all, plot(u{2000,3000},CO,ivpnl) title('Fig.~13.18.~~Chaotic rabbit population',FS,11) xlabel('$t$',FS,10,IN,LT), ylabel('rabbits, $u(t)$',FS,9) %% % \vskip 1.02em %% % Obviously real food webs can be far more complex than this. But this % example shows that the change % from two species to three can have a quantum effect, and for reasons % that are mathematical, not biological. %% % % \smallskip % {\sc History.} The roots of chaos are often traced to % Poincar\'e and Hadamard in the 1890s, long before the days of % computers, but the subject did % not become widely studied until the work of Lorenz, May, Yorke, % Feigenbaum % and others in the 1960s and 1970s. Chaos became known % to the public with James Gleick's 1987 book {\em Chaos:\ Making % a New Science} mentioned earlier and the 1993 movie {\em Jurassic Park}. % %% % % \smallskip % {\sc Our favorite reference.} % One of the best mathematical textbooks ever written is % S. H. Strogatz, {\em Nonlinear Dynamics and Chaos: With Applications % to Physics, Biology, Chemistry, and Engineering,} Westview % Press, 2015 (first published in 1994). Strogatz combines marvelous % clarity with a precise presentation of the mathematics % of this fascinating subject. % \smallskip % %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 13.} % Many nonlinear ODE\kern .5pt s with three or more variables are chaotic, though % this fact was not widely recognized % until the last third of the % twentieth century. Chaos is characterized by the % property that perturbations may grow % exponentially with time, yet global orbits remain bounded. % The maximal rate of exponential growth is known as the Lyapunov exponent. % \vspace{2pt}}} % \end{displaymath} % %% % \smallskip\small\parskip=1pt\parindent=0pt % {\em \underline{Exercise $13.1$}. Lucky Lyapunov exponents.} % It was mentioned in the text that the estimate $0.9052$ of the % Lyapunov exponent was fortuitously close to the true value. Confirm % this by computing the name number based on data for % $t\in [\kern .3pt 0,14]$ and $t\in [\kern .3pt 1,15]$ instead of % $t\in [\kern .3pt 0,15]$. % \par % {\em \underline{Exercise $13.2$}. Smaller perturbation of the Lorenz % equations.} % In the text we perturbed the Lorenz equations by changing the value % of $w(0)$ from $20$ to $20+ 10^{-5}$. {\em (a)} % Reproduce Figures 13.2 and 13.3 % with the smaller perturbation $20 + 10^{-9}$ and comment on the results. % {\em (\kern .7pt b)} Also repeat the finite-difference estimates of the derivatives % of $v(5)$, $v(10)$, and $v(30)$ with respect to $w(0)$, and % comment on the results. % \par % {\em \underline{Exercise $13.3$}. Alternative choices of % the Lorenz coefficient\/ $28$.} % In the Lorenz equations, let $r$ denote the parameter that traditionally % takes the value $28$. Starting from our usual initial conditions, make plots of % $u(t)$ against $w(t)$ as in Figure 13.1 % for $t\in [\kern .3pt 0, 100\kern .3pt]$ with % $r = 20, 22, 24$; also make plots in each case of $u(t)$ against % $t$ as in Figure 13.2. Which case seems to be chaotic? % Which one gives the clearest example of transient chaos? % \par % {\em \underline{Exercise $13.4$}. Lorenz equations with a breeze.} % In the Lorenz equations (13.1), let the first equation be % changed to $u' = 10(v-u) - a$, where $a$ is a parameter. % As in the last exercise, make plots of % $u(t)$ against $w(t)$ and of $u(t)$ against $t$ % for $t\in [\kern .3pt 0, 30\kern .3pt]$ with % $a = 20, 25, 30$. Comment on the solutions. % \par % {\em \underline{Exercise $13.5$}. Two electrons and a nucleus.} % Consider the highly idealized problem of two electrons % of mass 1 and charge $-1$ orbiting a nucleus of % charge $+2$ fixed at the origin of the $x$-$y$ plane. Let $z(t)$ and % $w(t)$ be the positions of the electrons represented with the % usual complex variable trick. Consider trajectories for % $t\in [\kern .3pt 0,20]$ starting from initial conditions % $z = i$, $z'=1$, $w=-i$, $w' = a$. % {\em (a)} Write down the ODE governing the evolution of $w(z)$ % and $z(t)$ % assuming an inverse-square electrostatic force law with constant~$1$. % The electrons repel each other while being attracted to the nucleus. % {\em (\kern .7pt b)} Plot the two trajectories in the case $a=1$. The % configuration is symmetric, and the symmetry should be % preserved for all $t$. % {\em (c)} Now make similar plots for $a = 0.5, 0.6, 0.7, 0.8, 0.9,$ % and $0.99$. Discuss the results. % \par % {\em \underline{Exercise $13.6$}. Chaos and cellular automata.} % Discrete processes, in which no real numbers are involved, can also % exhibit chaotic properties. Illustrate this by simulating the % Rule 30'' cellular automaton described in Stephen Wolfram's 2002 book % {\em A New Kind of Science.}\footnote{On p.~27 of the original % edition Wolfram calls this probably % the single most surprising scientific discovery I have ever % made.''} To do this, let $n$ be a % positive even integer and initialize an % $n\times n$ matrix ${\bf A}$ to zero apart from the value 1 in the middle % of the top row, $a_{1,n/2}^{} = 1$. Now sweep from % rows $2$ to $n$, updating each entry according to the following % rule: $a_{ij}^{}$ is changed to the the value $1$ if the three % entries above it, $a_{i-1,j-1:j+1}$, correspond to the binary % representation of one of the numbers $1,2,3,4$ rather than % $0,5,6,7$. (The numbers in columns $1$ and $n$ can be % left unchanged.) Take $n=200$ and use the Matlab {\tt spy} command to % plot the matrix at the end. Note that the structure is % completely deterministic yet has apparently % random properties. For another view of the randomness, plot % $a_{i,n/2}^{}$ as a function of $i$. % \par % {\em \underline{Exercise $13.7$}. Random Fibonacci sequence.} % Here, on the other hand, is an example of a discrete process with % true randomness. % Suppose $x_0^{} = x_1^{} = 1$ and for $k\ge 1,$ % $x_{k+1}^{} = x_{k-1}^{}+x_k^{}$. Then it is well known % that $x_k^{}$ grows asymptotically at a rate determined by % $\phi^k$, where $\phi$ is the golden ratio $(1+\sqrt 5\kern .5pt)/2 % \approx 1.618.$ Consider the % {\em random Fibonacci sequence} defined by % $x_{k+1}^{} = \pm x_{k-1}^{} \pm x_k^{}$, where at each step each sign is % independently $+$ or $-$ with probability $0.5$. % On a semilogy scale, make plots of $|x_k^{}|$ vs.\ $k$ % up to the maximum values % $k = 100$ and $5000$. Numerically estimate the {\em Lyapunov constant} for % this process, that is, the constant $C$ such that with probability~1, % the sequence grows in absolute value at a rate $C^k$ as $k\to\infty$. % (If you are curious to learn more, see D. Viswanath, % {\em Mathematics of Computation\/} 69 (2000), pp.~1131--1155.) %