%% 9. Second-order equations in the phase plane % % \setcounter{page}{96} % %% % % Consider a second-order autonomous ODE, which may be linear or nonlinear, % $$y''(t) = f(y,y') . \eqno (9.1)$$ % Given a pair of values $y$ and $y'$ at a particular % time $t$, this equation tells us the rate of change of $\kern .3pt y'$ at % $t$, and the rate of change of $y$ itself is by definition~$y'$. % Thus with (9.1) we know the rates of change of both % $y$ and $y'$, and we can % represent this information pictorially by a diagram of the % {\bf phase plane,} in which % the horizontal axis represents $y$ and the vertical axis % represents $y'$. For example, here is a quiver plot'' % of the phase plane % for the simple harmonic oscillator or linear pendulum % equation $y'' = -y$. % ODEformats, L = chebop(@(t,y) diff(y,2) + y); quiver(L,[-2.8 2.8 -1.1 1.1],'k',LW,.3,'scale',.5,'xpts',12,'ypts',8) title(['Fig.~9.1.~~Phase plane for $y''''=-y$' ... ' (simple harmonic motion)'],FS,11,IN,LT) axis equal, ylim([-1.1 1.1]) hold on, plot(0,0,'.k',MS,8), hold off xlabel('$y$',FS,10,IN,LT), ylabel('$y''$',FS,11,IN,LT) %% % \vskip 1.02em %% % % \noindent % At each point $(y,y')$, the arrow shows the direction and % magnitude of the vector $(y',y'')$. % Solutions to the ODE will correspond to trajectories % following the arrows around the plane. % Note that for any equation (9.1), regardless of $f$, the arrows will % always point rightward in the upper half-plane and leftward in % the lower half-plane. % The black dot at $y = y'= 0$ marks a % {\em fixed point} (or {\em steady state}) of this equation, which % means, a point at which $y' = y'' = 0$. % %% % % An image like Figure 9.1, while accurate, is sometimes not % as compelling as one in which all the arrrows are % set to have the same length, giving a plot of % a {\em direction field.} This can be done with the % {\tt quiver} option \verb|'normalize'|. % quiver(L,[-2.8 2.8 -1.1 1.1],'k',LW,.3,'scale',.4,'normalize',1,'xpts',12,'ypts',8) title('Fig.~9.2.~~Same but with arrows all of the same length',FS,11,IN,LT) axis equal, ylim([-1.1 1.1]) hold on, plot(0,0,'.k',MS,8), hold off xlabel('$y$',FS,10,IN,LT), ylabel('$y''$',FS,10,IN,LT) %% % \vskip 1.02em %% % % There is nothing more to the dynamics of an autonomous % ODE than its phase plane. If an % initial point $(y,y')$ is specified, then all of the future % trajectory is determined simply by the vector % field.\footnote{To be precise, this requires the solutions % to be unique, which will be assured if $f$ satisfies a % condition of Lipschitz continuity. See Chapter 11.} % This is the power of phase plane analysis: it reduces dynamics % to geometry. % %% % % For example, let us now specify a\/ $\kern .5pt t\kern .5pt$ domain and an initial condition % for the linear pendulum, % $$y''= -y, \quad t\in [\kern .5pt 0, 1.8\pi],~ % y(0) = 0, ~ y'(0)=1, \eqno (9.2)$$ % or in Chebfun, % L.domain = [0,1.8*pi]; L.lbc = [0;1]; %% % % \noindent % Here we superimpose the trajectory corresponding to the % solution of (9.2) on the vector field just displayed. % The trajectory begins at the top of the unit % circle with initial velocity~$1$ % and initial acceleration~$0$, so the curve is % oriented horizontally to the right. As $y$ increases, $y''$ becomes % negative, and the curve begins to bend around, describing % a clockwise circle. % Since the time interval runs $90\%$ of the way to $2\pi$, % the circle is $90\%$ complete. % y = L\0; hold on, arrowplot(y,diff(y),CO,ivp,YS,3.5) title(['Fig.~9.3.~~A trajectory added to the plot, ' ... 'the solution of (9.2)'],FS,11), hold off %% % \vskip 1.02em %% % Let us dispense with the quiver arrows and collect % trajectories with different initial values on a single % plot. The % next picture takes $b=0.5,0.8,\dots, 2$, producing six circles % of corresponding radii. for b = .5:.3:2 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),CO,ivp,MS,4,YS,3), hold on end axis equal, ylim([-2.2 2.2]) plot(0,0,'.k',MS,8) title('Fig.~9.4.~~Six trajectories for (9.2)',FS,11) %% % \vskip 1.02em %% % The phase plane can be informative about boundary-value % problems, too. For example, suppose we change (9.2) to % $$y''= -y, \quad t\in [\kern .5pt 0, 3], % ~ y(0)=1, ~ y(4) = -1.5. \eqno (9.3)$$ % Here is the solution, shown as a seventh trajectory added % to the previous plot (without an arrowhead since this % is a BVP). L = chebop(0,4); L.op = @(t,y) diff(y,2) + y; L.lbc = 1; L.rbc = -1.5; y = L\0; plot(y,diff(y),CO,bvp) hold off, title('Fig.~9.5.~~Solution of BVP (9.3)',FS,11) %% % \vskip 1.02em %% % % Equations (9.2) and (9.3) describe a linear pendulum without % damping. Following Chapters 4 and 8, let us now add % some damping in the form of a small multiple of $y'$, % $$y''= -y-\varepsilon \kern .3pt y', \quad t\in [\kern .5pt 0, 1.8\pi], % ~ y(0)=0, ~ y'(0) = b, \eqno (9.4)$$ % with $\varepsilon=0.2$. % If we follow the % same six trajectories as before, we see that they now lose amplitude % as they evolve, moving towards the fixed point % $y = y'= 0$ as $t\to\infty$. % This is a {\em stable} fixed point since all nearby trajectories % stay close to it; we shall consider such definitions % systematically in Chapter~15.\ \ The % plot also shows the solution to the BVP variant of % (9.4) with the same boundary conditions $y(0)=1$, % $y(4) = -1.5$ as in (9.3). % Note that the solution starts from a greater velocity $y'$ % than before, necessary to compensate for % the damping while still matching the boundary conditions. % L = chebop(0,1.8*pi); L.op = @(t,y) diff(y,2) + 0.2*diff(y) + y; for b = .5:.3:2 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),CO,ivp,MS,3,YS,3), hold on end axis equal, ylim([-2.2 2.2]) L = chebop(0,4); L.op = @(t,y) diff(y,2) + 0.2*diff(y) + y; L.lbc = 1; L.rbc = -1.5; y = L\0; plot(y,diff(y),CO,bvp) plot(0,0,'.k',MS,8), hold off title('Fig.~9.6.~~Damped linear pendulum (9.4)',FS,11) %% % \vskip 1.02em %% % % \noindent % We stopped the flow at $t=1.8\pi$, but of course this is not necessary. % Here is a single trajectory carried to $9\pi$. % This represents a linear pendulum swinging back and forth $4~1/2$ times, % losing amplitude as it swings. % L = chebop(0,9*pi); L.op = @(t,y) diff(y,2)+0.2*diff(y)+y; L.lbc = [0;2]; y = L\0; plot(y,diff(y),CO,ivp), axis equal, ylim([-2.2 2.2]) hold on, plot(0,0,'.k',MS,8), hold off title(['Fig.~9.7.~~A trajectory for (9.4) ' ... 'carried to longer time'],FS,11) %% % \vskip 1.02em %% % % We mentioned that in % the phase plane, the position of a trajectory at any time % $t$ determines its future, which consists of following the unique path % through that point. This uniqueness property % depends on the equation being autonomous. If % we solve a nonautonomous problem, such as % $$y''= -y-0.2\kern .3pt y'-2\cos(2\kern .3pt t), \quad t\in [\kern .5pt 0, 1.8\pi], % ~ y(0)=0, ~ y'(0) = b, \eqno (9.5)$$ % we can plot the solution in the $y$-$y'$ plane, like this (now % with different colors to help distinguish the curves), % L = chebop(0,1.8*pi); L.op = @(t,y) diff(y,2)+0.2*diff(y)+y+2*cos(2*t); for b = .5:.3:2 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),LW,.5,MS,3), hold on end axis equal, ylim([-3 2.4]), plot(4.5,0,'Xr',MS,30,LW,2) title(['Fig.~9.8.~~Inhomogeneous problem --- not truly '... 'a phase plane!'],FS,11), hold off %% % \vskip 1.02em %% % % \noindent % However, the point of planar analysis has % been lost since the future of a trajectory is not % determined by its current position in the plane, and the % curves cross each other. The plot of trajectories % has become merely a plot of trajectories, no longer % an encapsulation of the dynamics of the system. % Such plots can be interesting and informative; see for example % Exercise 9.2 and the Duffing oscillator in Appendix~B.\ \ They % usually don't belong to phase plane analysis, however. % An alternative in such cases is to include the time variable in % the plot on an additional axis. Here, for example, we repeat the % calculation on the longer interval $[\kern .3pt0 ,40\kern .3pt]$, % for a single trajectory, and plot the result in $t$-$y$-$y'$ space. % L.domain = [0,40]; y = L\0; t = chebfun('t',[0 40]); plot3(t,y,diff(y),CO,ivp,LW,.5), view(-48,27) axis([0 40 -2.5 2.5 -3 2]) axis square, xlabel('$t$',FS,10), ylabel('$y$',FS,10), zlabel('$y''$',FS,10) title('Fig.~9.9.~~Trajectory in $t$-$y$-$y''$ space',FS,11) set(gcf,'position',[380 220 540 300]) %% % \vskip 1.02em %% % % Returning to (9.4), in the language of Chapter 4, the damping % with $\varepsilon = 0.2$ is {\em subcritical\/}: from (4.11) we see that % the critical damping parameter for this equation is $\varepsilon=2$. % Here are images of subcritical, critical, % and supercritical trajectories with % $\varepsilon = 1,2,4$, now for $t\in [\kern .3pt 0, 1.5\pi]$. % Note that it is the % critically damped trio of trajectories that finishes closest % to the origin, confirming that $\varepsilon=2$ gives % the most effective damping. % close all, L.domain = [0,1.5*pi]; subplot(1,3,1), L.op = @(t,y) diff(y,2)+diff(y)+y; for b = 1:-.2:.6 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),MS,3,YS,.8), hold on end axis equal, ylim([-0.4 1.25]) pp = get(gca,'pos'); pp(1) = pp(1)+.03; set(gca,'pos',pp) text(.21,1.09,'subcritical, $\varepsilon=1$',FS,9,HA,CT) hold on, plot(0,0,'.k',MS,8), hold off subplot(1,3,2), L.op = @(t,y) diff(y,2)+2*diff(y)+y; for b = 1:-.2:.6 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),MS,3,YS,.8), hold on end axis equal, ylim([-0.4 1.25]) pp = get(gca,'pos'); pp(1) = pp(1)-.0003; set(gca,'pos',pp) set(gca,YT,0:.5:1,YTL,{' ',' ',' '}) text(.21,1.09,'critical, $\varepsilon=2$',FS,9,HA,CT) hold on, plot(0,0,'.k',MS,8), hold off title(['Fig.~9.10.~~Subcritical, critical, and ' ... 'supercritical damping'],FS,11,HA,CT) subplot(1,3,3), L.op = @(t,y) diff(y,2)+4*diff(y)+y; for b = 1:-.2:.6 L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),MS,3,YS,.8), hold on end axis equal, ylim([-0.4 1.25]) set(gca,YT,0:.5:1,YTL,{' ',' ',' '}) pp = get(gca,'pos'); pp(1) = pp(1)-.03; set(gca,'pos',pp) text(.12,1.09,'supercritical, $\varepsilon=4$',FS,9,HA,CT) hold on, plot(0,0,'.k',MS,8), hold off %% % \vskip 1.02em %% % % Phase plane analysis becomes particularly interesting % for nonlinear ODE\kern .5pt s. For example, the % first nonlinear equation of this book was the % van der Pol equation. With slightly % different coefficients from (1.2), let us write this as % $$y'' + y - \mu (1-y^2)y' = 0, \quad t\in[\kern .5pt 0, 10\kern .5pt], % ~ y(0) = a, y'(0) = 0. \eqno (9.6)$$ % Here are phase plane plots for $\mu = 0.125$ and $\mu = 1.5$. % With the weak damping parameter $\mu=0.0125$, the system is % not far from the linear pendulum, with trajectories winding % slowly in or out to an asymptotic curve known % as a {\bf limit cycle}. With the stronger damping parameter % $\mu = 1.5$, trajectories converge to the limit cycle much faster, % and its shape is far from circular, just as the % van der Pol orbit plotted in Chapter 1 is far from a sine wave. % N = chebop(0,15); clf for j = 1:2 subplot(1,2,j), mu = 0.125; if j==2, mu = 1.5; end N.op = @(t,y) diff(y,2)-mu*(1-y^2)*diff(y)+y; for a = [1 3] N.lbc = [a;0]; y = N\0; arrowplot(y,diff(y),LW,.6,MS,4), hold on end axis equal, ylim([-4.2 4.2]), plot(0,0,'.k',MS,8) text(4.5,3.1,['$\mu = ' num2str(mu) '$'],FS,9,HA,RT) end subplot(1,2,1) title('\kern 2.4 in Fig.~9.11.~~Van der Pol equation',FS,11) %% % \vskip 1.02em %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: nonlinear pendulum}\\[-3pt] % \hrulefill % \end{center} % %% % % The linear pendulum equation $y'' = -y$ describes oscillation in % the context of Hooke's Law, where the restoring force is % proportional to the displacement --- simple harmonic motion. % When the amplitudes of motion are small, % this is the right model for a spring or a % pendulum and for many other vibrating systems, but when % the amplitudes get bigger, the physics always becomes nonlinear. % Different problems have different nonlinearities, but % there is no doubt as to the archetypical problem of this % kind: it is the {\bf nonlinear pendulum}, % corresponding to an idealized point mass moving in a % circle at the end of a rigid weightless bar. % The equation is % $$y''= -\sin(y), \eqno (9.7)$$ % where $y(t)$ represents the angle from the vertical in radians % at time $t$ and constants are set to 1. % For an entire book on the subject, % see Baker and Blackburn, {\em The Pendulum,} % Oxford University Press, 2005. % %% % % This is a perfect example for phase plane analysis. % First we draw a quiver plot together with stable and % unstable fixed points (black and red, respectively). % There are stable fixed points at $(y,y') = (2\pi j,0)$ % for each integer~$j$ and unstable fixed % points at the in-between locations $(2\pi(\kern .7pt j+1/2),0)$. % N = chebop(0,1.8*pi); N.op = @(t,y) diff(y,2)+sin(y); clf quiver(N,[-3 23 -6 6],'k',LW,.3,'scale',.5,'xpts',12,'ypts',8) axis equal, xlim([-3 23]), hold on plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10) plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10), hold off title('Fig.~9.12.~~Nonlinear pendulum (9.7)',FS,11) %% % \vskip 1.02em %% % % To see more, let us plot some trajectories. % We start from $y=0$ and % $y'(0) = b = 1,1.2,\dots, 4$, calculating trajectories % over the interval $t\in [\kern .5pt 0, 1.8\pi]$.\footnote{To us this % looks like the hair of Botticelli's Venus; or is it % the Starbucks logo?} % plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10), hold on plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10) for b = 1:.2:4 N.lbc = [0;b]; y = N\0; arrowplot(y,diff(y),LW,.6,CO,ivp,YS,3,MS,3) end axis equal, xlim([-3 23]), hold off title('Fig.~9.13.~~Nonlinear pendulum (9.7)',FS,11) %% % \vskip 1.02em %% % % \noindent % Physically, these curves can be interpreted as follows. % If $y'(0)=b<2$, the pendulum does not have enough % energy to swing over the top, and the trajectory goes % around and around forever on a single loop in the % phase plane, corresponding to the pendulum swinging % back and forth, with the angle $y(t)$ remaining % bounded. The loops are not circles, but they % approach circles for small amplitude, where the % distinction between $y$ as in (9.2) and $\sin(y)$ % as in (9.7) fades away. % If $y'(0)=b>2$, on the other hand, % the pendulum has enough energy to % swing over, and $y(t)$ keeps increasing monotonically. % In the absence of damping, a pendulum that swings over % once swings over infinitely many times as $t\to\infty$. % The phase plane is $2\pi$-periodic. % %% % Notice the trajectory starting at $y'(0) = 2$. This one appears % to stop at the $y$ axis, and that is exactly what it does. % It has just enough energy to fly up toward the vertical % configuration, with $y(t)$ approaching $\pi$ as $t\to\infty$, but % it never quite hits the top for any finite value of $t$. %% % Of course a real pendulum is sure to have some losses. % Here is the same image corresponding to an equation with subcritical % damping, % $$y''= -\sin(y) - \mu y', \quad t\in [\kern .5pt 0, 1.8\pi], % ~ y(0)=0, ~ y'(0) = b, \eqno (9.8)$$ % with $\mu = 0.1$. % A curve is also added to the plot corresponding to the % BVP defined by $y(0)= 0$, $y(1.95\pi) = 20$. % We see that to reach $y=20$ by the end of % the time interval, a greater initial speed is needed. N.op = @(t,y) diff(y,2)+sin(y)+.1*diff(y); for b = 1:.2:4 N.lbc = [0;b]; y = N\0; arrowplot(y,diff(y),LW,.6,CO,ivp,YS,2,MS,3), hold on end plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10) plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10) axis equal, xlim([-3 23]) title('Fig.~9.14.~~Damped nonlinear pendulum (9.8)',FS,11) N.lbc = 0; N.rbc = 20; y = N\0; plot(y,diff(y),CO,bvpnl), hold off %% % \vskip 1.02em %% % % \noindent % Now, a trajectory % with sufficient initial energy swings over, but it % keeps losing energy as it goes, so it may or may % not swing over a second time. We can see more % if we increase the time interval to $[\kern .5pt 0,8\pi]$. % The trajectory with initial velocity $y'(0) = 4$ swings % over four times before eventually winding down to rest. % As before, the function $f(y,y')$ defining the phase plane % is $2\pi$-periodic, though none of the individual trajectories % are periodic. % N = chebop(0,8*pi); N.op = @(t,y) diff(y,2)+sin(y)+.1*diff(y); plot(pi*[0 2 4 6 8],[0 0 0 0 0],'.k',MS,10), hold on plot(pi*[1 3 5 7 9],[0 0 0 0 0],'.r',MS,10) for b = 1:.2:4 N.lbc = [0;b]; y = N\0; arrowplot(y,diff(y),LW,.6,CO,ivp,MS,3) end axis equal, xlim([-3 29]), hold off title('Fig.~9.15.~~Longer time interval for (9.8)',FS,11) %% % \vskip 1.02em %% % % \medskip % {\sc History.} Four hundred years ago, Galileo understood % the essentials of the linear % pendulum. In his {\em Dialogues Concerning Two New % Sciences\/} the character Salviati % says, As to the times of vibration of bodies suspended by % threads of different lengths, they bear to each other the % same proportion as the square roots of the lengths of % the threads; or one might say the lengths are to each other % as the squares of the times; so that if one wishes to make % the vibration-time of one pendulum twice that of another, % he must make its suspension four times as long.'' % %% % % \smallskip % {\sc Our favorite reference.} When it comes to phase plane % analysis, one of the oldest books is also one of the % nicest: H. T. Davis, {\em Introduction to Nonlinear % Differential and Integral Equations,} Dover, 1962 (first published % in 1960). For example, we enjoy Chapter 10 on The phase plane % and its phenomena,'' Chapter 11 on Nonlinear % mechanics,'' and Chapter 12 on Some particular equations.'' % \smallskip % %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 9.} The phase plane for a % second-order autonomous homogeneous ODE\/ $y'' = f(y,y')$ is % the plane with coordinates\/ $y$ and $y'$. The position of % a solution at a particular time\/ $t$ determines its future trajectory % as the unique curve passing through this point. Plotting % phase plane trajectories gives a quick way to interpret behavior of % ODE\kern .5pt s, including nonlinear ones. % \vspace{2pt}}} % \end{displaymath} % %% % % \small\parindent=0pt\parskip=1pt % {\em Exercise $9.1$. $y'' = y$.} % Draw phase plane trajectories or % quiver arrows for the equation $y'' = y$ (either % by hand or with the computer). % Where in this image does the solution of Fig.~5.3 lie? % (It comes very close to following the {\em stable manifold} of % the fixed point for 30 time units and then the % {\em unstable manifold} for a further 30 times units; see % Chapter~15.) % \par % {\em \underline{Exercise $9.2$}. Moving slowly in % the phase plane.} % Consider the ODE $y'' + y = y^2$. % {\em (a)} What are the fixed points? % {\em (\kern .7pt b)} Draw a quiver plot, choosing axes to make % the plot as informative as possible. Based on this % information, draw a sketch by hand of the key points % of the phase plane dynamics. % Describe qualitatively what orbits will remain bounded, % and how they will behave; also what orbits will diverge % to $\infty$, and how they will behave. % {\em (c)} Find two distinct solutions that satisfy % $y(0) = y(10) = 2$. Plot them both as functions of $t$ and % in the phase plane. % What are the values of $y'(0)$ for these two solutions? % \par % {\em \underline{Exercise $9.3$}. Region-filling orbits.} % (Adapted from Davis, {\em Introduction to Nonlinear % Differential and Integral Equations,} \S 10.4.) % {\em (a)} Use the method of undetermined coefficients to find the % analytical solution of the IVP $y''+2\kern .3pt y = % -2\cos(2\kern .3pt t)$, $y(0)=1$, $y'(0)=\sqrt{2}$. % {\em (\kern .7pt b)} This solution has a curious property in the phase plane: the % curve eventually gets arbitrarily close to every point in a % near-elliptical region. Plot the solution over $t\in[\kern .3pt 0,100\kern .3pt]$ to see the % effect. (This is fun to watch using {\tt comet} as well.) % {\em (c)} The solution just considered is nonperiodic. Periodic % solutions are obtained from the same initial conditions, on the % other hand, if $\cos(2\kern .3pt t)$ is replaced % by $\cos(k\kern .3pt t)$ for certain % values of $k>0$. Find the smallest three such values analytically, % and produce the corresponding plot for $t\in[\kern .3pt 0,100\kern .3pt]$. % \par % {\em \underline{Exercise $9.4$}. Period of the nonlinear pendulum.} % The code % \begin{verbatim} % y = @(s) chebop(@(y) diff(y,2)+sin(y),[0 15-5*log(pi-s)],[s;0],[])\0; % T = @(s) 2*min(diff(roots(y(s)))); % \end{verbatim} % produces two anonymous functions $y$ and $T$ applicable for % values $s\in (0,\pi)$. % Explain what these functions compute and how they do it. % Make a plot of $T(s)$ for $s\in[\kern .3pt 0.1,3.14]$ % and explain its principal features. % Explain in a general way (not necessarily with % mathematical details) the role of the quantity $15-5\log(\pi-s)$. %