%% 4. Second-order equations and damped oscillations

%%
% <latex>
% \setcounter{page}{35}
% </latex>

%%
% <latex>
% The most famous of all ODE\kern .5pt s is of order two --- {\em Newton's second
% law of motion,} force equals mass times acceleration:
% $$ F = ma. $$
% What makes this a second-order equation is that acceleration is the
% second derivative of position with respect to time.  Following the
% standard notation of this book, we can rename $a$ as $y''$ and
% solve for it to get
% $$ y'' = {F(t,y)\over m} . \eqno (4.1) $$
% This particular form of Newton's law describes a mass $m$ moving
% along a line with position $y(t)$ at time $t$ and
% subject to a force $F$ that may depend on $t$ and/or $y(t)$.
% </latex>

%%
% <latex>
% Staying with 17th century England for a moment, a famous
% choice of $F$ is {\em Hooke's Law},
% $$ F(t,y) = - b\kern .3pt y , \eqno (4.2) $$
% where $b$ is a positive constant.  This applies for
% example to a spring stretched a distance $y$ beyond its
% rest position.  Combining (4.1) and (4.2) gives the ODE
% $$ y'' = -{b\kern .3pt y\over m} . $$
% We shall write this in a standard form
% that will also be useful in later chapters by defining
% $$ \omega = \sqrt{b\over m} . $$
% The equation becomes
% $$ y'' +\omega^2 y = 0 . \eqno (4.3) $$
% This is the most basic of all second-order ODE\kern .5pt s, the second-order
% analogue of (1.1), and it governs
% {\em simple harmonic motion.}
% As a second-order equation, it has a two-dimensional
% space of solutions spanned by the functions
% $$ \sin (\omega\kern .3pt t) , ~~ \cos (\omega\kern .3pt t) . $$
% Note that these functions are {\em periodic,} with
% period $T$ given by the formula
% $$ T = {2\pi\over \omega }.  \eqno (4.4) $$
% The number $\omega $ is the {\em frequency}.\footnote{One must
% be careful to distinguish this {\em angular frequency,}
% representing radians per unit time,
% from the quantity $1/T$, also called frequency,
% that represents cycles per unit time.
% The musical note A above middle C corresponds to 440 Hz, which
% is an angular frequency of $\approx 2765$ radians per second.}
% The general solution of (4.3) can be written\footnote{When
% the independent variable is
% $x$ rather than $t$, $\omega$ is replaced by $k$, known as
% the {\em wave number}.  In this case $\sin(\omega\kern .3pt t)$ and $\cos(\omega\kern .3pt t)$
% become $\sin(kx)$ and $\cos(kx)$.}
% $$ A \sin (\omega\kern .3pt t) + B\cos (\omega\kern .3pt t) .  \eqno (4.5) $$
% Other representations of this space of solutions are
% also useful, such as
% $$ A \sin(\omega\kern .3pt t + \phi), \eqno (4.6) $$
% where $\phi$ is an phase shift, or the
% complex exponentials expression\footnote{Mathematicians generally
% prefer the complex exponentials representation for
% linear problems, where sines and cosines can be obtained
% by appropriate superpositions.  For nonlinear problems,
% superposition is not available, and it is generally necessary
% to stay with real variables.}
% $$ A e^{i\kern .3pt\omega\kern .3pt t} + B e^{-i\kern .3pt\omega\kern .3pt t} , \eqno (4.7) $$
% where $A$ and $B$ may also
% be complex.  This form is particularly useful when one wishes to consider
% the generalization to
% situations in which $\omega $ is complex, corresponding to a
% constant $\omega^2$ in (4.2) that is negative or complex.
% From (4.7) it can be seen that
% in this case solutions to (4.1) contain components that grow and
% decay exponentially with $t$.
% </latex>

%%
% <latex>
% Here for example is a solution of (4.3) with $\omega =1$ and
% initial conditions $y(0)=1$, $y'(0) = 0$.  We know this solution
% exactly: it is $y(t) = \cos(t)$.  The number displayed in
% the figure title confirms that the period is $2\pi$.
% </latex>
ODEformats, L = chebop(0,60); L.op = @(t,y) diff(y,2) + y;
L.lbc = [1;0]; y = L\0; plot(y,CO,ivp), axis([0 60 -2 2])
[~,maxima] = max(y,'local'); T = maxima(3) - maxima(2);
title(sprintf( ...
  'Fig.~4.1.~~Simple harmonic motion (4.3): $T = %6.4f$',T),FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% In this book we are always free to experiment.
% For example, what if we make the problem nonlinear
% by replacing $y$ by $y^5\kern .7pt$?  The equation becomes
% $$ y'' + y^5 = 0. \eqno (4.8) $$
% A figure shows that the solution
% is again a periodic oscillation of amplitude 1,
% but the period is $34\%$ longer:
L.op = @(t,y) diff(y,2) + y^5;
y = L\0; plot(y,CO,ivpnl), axis([0 60 -2 2])
[~,maxima] = max(y,'local'); T = maxima(3) - maxima(2);
title(sprintf( ... 
  'Fig.~4.2.~~Nonlinear spring law (4.8): $T = %6.4f$',T),FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The reason the period has increased is easy to see:
% $|y^5|$ is smaller than~$|y|$ for $|y|<1$, so this
% nonlinear spring has a weaker restoring force than the linear one and
% oscillates more slowly.  The
% shape of the oscillation has also changed, coming
% closer to a sawtoothed alternation of straight segments.
% We are closer to a ``bang-bang''
% situation, with most of the acceleration, hence
% most of the curvature, appearing for $y\approx \pm 1$.
% </latex>

%%
% <latex>
% As another experiment, let us return to the original
% linear problem (4.3) but add a small term involving the first derivative,
% giving an equation associated with
% {\em damped oscillation,}\footnote{In the van der Pol
% equation (1.2), we saw a similar
% oscillator except with nonlinear damping: positive damping for
% $|y|>1$, negative for $|y|<1$.}
% $$ y'' + 0.1 y' + \omega^2 y = 0. \eqno (4.9) $$
% Physically, this corresponds to an additional applied force
% that is proportional to velocity and in the opposite
% direction, slowing motions down.
% Here is the result, again with $\omega =1$ and $y(0) = 1$,
% $y'(0) = 0$.
% </latex>
L.op = @(t,y) diff(y,2) + 0.1*diff(y) + y;
y = L\0; plot(y,CO,ivp), axis([0 60 -1.5 1.5])
[~,maxima] = max(y,'local'); T = maxima(3) - maxima(2);
title(sprintf( ...
  'Fig.~4.3.~~Damped harmonic motion (4.9): $T = %6.4f$',T),FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Note from the number displayed in the figure
% title that the ``period'' is greater than $2\pi$, but
% only about $1\%$ greater; we shall explain this in a moment.
% Strictly speaking, this solution is not periodic, but
% the maxima are still separated by fixed intervals.  We
% shall see in (4.13) that $y$ is equal to a periodic function times a
% decaying exponential.
% </latex>

%%
% <latex>
% Few subjects are more vital to practical engineering
% than the damping of oscillations.\footnote{Standard texts
% on the analysis of vibration include Meirovitch, {\em Elements of
% Vibration Analysis,} and S. S. Rao, {\em Mechanical
% Vibrations.}}  In cars and aircraft, vibrations
% with frequencies on the order of hundreds or thousands of Hertz 
% must be damped to keep sound levels tolerable.  In bridges and
% buildings, vibrations with frequencies closer to 1-5 Hertz
% must be damped to avoid catastrophic failures in storms and
% earthquakes.  The mathematics of this kind of engineering is
% well advanced, relying heavily on finite-element simulations
% on large computers.  One of the world's most famous damping
% devices must be the 660-ton steel pendulum suspended from
% the 92nd to the 87th floor of the skyscraper
% Taipei 101 in Taiwan.  This damper has become a tourist
% attraction and even has its own mascot, known as the Damper Baby.
% </latex>

%%
% It is interesting to consider the implications of linearity
% and nonlinearity in the three problems we have examined so far. The shape
% of the first
% solution, for the linear problem (4.3), is scale-independent:
% if the initial condition is doubled, the solution
% will double too.  The second
% solution, for the nonlinear problem (4.8), is scale-dependent:
% if the initial condition is doubled, the solution will
% change shape as well
% as amplitude, developing even sharper sawtooths.  The third solution, for
% the linear damped oscillator (4.9), is the most interesting from this
% point of view.  Again, linearity
% implies scale-independence, so that doubling
% the data will double the solution.  What is new is that since the
% wave is decaying as a function
% of $t$, the scale-independence has a substantive
% implication: as $t$ varies, the amplitude can change but not the shape.
% Thus it follows from linearity that the crests of
% the wave solution to (4.9) must be separated by fixed intervals, and
% hence that the decay must be exponential.

%%
% <latex>
% Let us consider now the mathematics of a general problem of this form,
% a {\em second-order constant-coefficient scalar linear homogeneous
% differential equation,}
% $$ y'' + \varepsilon y' + \omega^2\kern .3pt y = 0, \eqno (4.10) $$
% where $\varepsilon$ and $\omega$ are constants.  (Since no initial or
% boundary conditions have been specified, the classification
% of (4.10) is {\sf f\kern 1pt LASH}.)\ \ As mentioned in
% the introduction, this is one of the five most
% important equations of this book, and it will keep
% reappearing in various contexts.  In applications $\varepsilon$ will
% usually be nonnegative, although mathematically, $\varepsilon<0 $
% and $\varepsilon>0 $ are really
% the same story --- for negating $\varepsilon$
% is equivalent to negating $t$, that is, time-reversal.
% None of the formulas or theorems care whether $\varepsilon$ is
% positive or negative.
% </latex>

%%
% <latex>
% Motivated by (4.7), let us look for a solution to (4.10) of the form
% $$ y(t) = e^{rt} $$
% for a constant $r$.  Inserting this trial solution
% in (4.10) gives the condition
% $$ r^2 + \varepsilon\kern .3pt r + \omega^2 = 0 , $$
% which is a quadratic equation with roots
% $$ r = -{\varepsilon\over 2} \pm \sqrt{{\varepsilon^2\over 4} - \omega^2}.  \eqno (4.11) $$
% We call these the {\bf characteristic roots} of (4.10).
% Equivalently, we could write
% $$ r = -{\varepsilon\over 2} \pm i \sqrt{\omega^2-{\varepsilon^2\over 4}} . \eqno (4.12) $$
% The product of $r_-^{}$ and $r_+^{}$ is always
% $\omega^2$, independently of the value of $\varepsilon$.
% </latex>

%%
% <latex>
% For simplicity, let us suppose that $\varepsilon$ and $\omega$ are real.
% Then the formula (4.12) is the more useful one if
% $ \varepsilon^2/4< \omega^2$, i.e., $\varepsilon < 2\kern .5pt\omega$,
% which we call the {\em underdamped} or {\em subcritically damped} case.
% The factor $i$ reveals that solutions are oscillatory.
% In particular, if $\varepsilon=0$, (4.12) reproduces the solutions
% $y = \exp(\pm i\kern .5pt \omega\kern .3pt  t)$ of (4.7).  For $\varepsilon>0$,
% we get exponentially decaying solutions of the form
% $$ y(t) = e^{-\varepsilon t/2} \exp\left(\pm i \kern .5pt 
% t\sqrt{\omega^2-{\varepsilon^2\over 4}}\,\right). \eqno (4.13) $$
% This formula explains the decaying wave seen in the solution
% to (4.9).  With $\varepsilon = 0.1$, the decay rate is
% $e^{-t/20}$, as we can confirm by adding dashed lines to
% the figure.
% </latex>
ep = 0.1; envelope = chebfun(@(t) exp(-ep*t/2),[0 60]);
title(sprintf( ...
  'Fig.~4.4.~~Showing the exponential envelope from (4.13)',T),FS,11)
hold on, plot(envelope,'--r',-envelope,'--r'), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Formula (4.13) also explains the slight increase of
% the period $T$ above $2\pi$ we have noted.  Adapting (4.4), we
% see that the period of the decaying wave (or more precisely
% of its oscillatory factor) is now
% $$ T = 2\pi\kern -2pt  \left/ \kern -2pt \sqrt{\omega^2 - 
% \displaystyle{\varepsilon^2\over 4}}\right. .
% \eqno (4.14) $$
% Thus damping slows down the frequency of oscillation.
% It is the squaring of~$\varepsilon$ in this
% formula that explains why the modification
% to $T$ in our experiment has been so slight.
% To leading order, the addition of linear damping does not
% change the frequency of an oscillating system;
% the change of frequency is at the second order.\footnote{For this
% reason engineers often begin an analysis of
% frequency of vibrations of machines and structures
% by ignoring damping.  The same insensitivity of frequency
% to damping benefited watchmakers in the days of watches
% with balance springs,
% a design method going back to Hooke himself.}  An
% evaluation of (4.14)
% confirms the value of $T$ reported in title of Figure~4.3:
% </latex>
om = 1; T = 2*pi/sqrt(om^2 - ep^2/4)

%%
% <latex>
% Continuing with the assumption that
% $\varepsilon$ and $\omega^2$ in (4.10) are real,
% suppose on the other hand $\varepsilon^2/4 > \omega^2$, i.e.,
% $\varepsilon > 2\kern .5pt\omega$, the {\em overdamped} or {\em supercritically
% damped} case.
% Now (4.11) becomes more useful than (4.12), and we see that
% the solutions to (4.10) do not oscillate but grow or decay
% exponentially.  If $\varepsilon>0$ and $\omega^2>0$ they both decay, though
% at different rates.
% </latex>

%%
% <latex>
% Equation (4.10) is a second-order linear equation,
% and that is why it has two linearly independent solutions and
% requires two initial conditions to determine a unique
% solution.\footnote{Recall that in any vector space, 
% two vectors $\bf v$ and $\bf w$ are linearly independent
% if $A{\bf v} + B{\bf w} = {\bf 0}$ holds only
% with $A=B=0$.  Here, our vectors are functions $y(t)$.}
% If the two values $r_+^{}$ and $r_-^{}$ of (4.11) or (4.12)
% are distinct, then the general solution to (4.10) is
% $$ y(t) = Ae^{r_+^{} t} + B e^{r_-^{} t} . \eqno (4.15) $$
% If $\omega^2=\varepsilon^2/4$, however, or
% $\varepsilon = 2\kern .5pt\omega,$ then $r_+^{}$ and $r_-^{}$ are the same.
% This is the {\em critically damped} case,
% and we have to look elsewhere to find a second
% linearly independent solution
% of the ODE.\ \ A suitable choice is
% $y(t) = t\kern .5pt  e^{rt}$, as
% can be verified by substituting this function
% directly in (4.10).
% In this case we accordingly get the general solution
% $$ y(t) = (A+Bt) e^{rt}. \eqno (4.16) $$
% </latex>

%%
% Here is a theorem to summarize these observations.

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 4.1.  Solution of
% second-order linear scalar constant-coefficient ODE
% (\textsf{\textbf{f\kern 1pt LASH}}).}
% The equation
% $$ y''+ \varepsilon y' + \omega^2\kern .3pt y = 0  \eqno (4.17) $$
% has a two-dimensional vector space of solutions.
% For $\varepsilon < 2\kern .5pt\omega$ (subcritical damping) the general solution is
% $$ y(t) = e^{-\varepsilon\kern .5pt t/2}
% \Bigl[ A\sin\bigl(t \sqrt{\omega^2- \varepsilon^2/4} \,\bigr)
% + B \cos\bigl(t\sqrt{\omega^2- \varepsilon^2/4}\, \bigr)\Bigr],
% \eqno (4.18) $$
% for $\varepsilon > 2\kern .5pt\omega$ (supercritical damping) it is
% $$ y(t) =
% A\exp\bigl(-\varepsilon\kern .5pt
% t/2+t \sqrt{\varepsilon^2/4-\omega^2} \,\bigr)
% + B \exp\bigl(-\varepsilon\kern .5pt t/2 -
% t \sqrt{\varepsilon^2/4-\omega^2}\, \bigr),
% \eqno (4.19) $$
% and for $\varepsilon = 2\kern .5pt\omega$ (critical damping) it is
% $$ y(t) = 
% (A+Bt) \exp\bigl(-\varepsilon\kern .5pt t/2\bigr).
% \eqno (4.20) $$
% }
% </latex>

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  The solutions described are certainly linearly
% independent, so all that must be proved is that the
% dimension of the space of solutions is not greater than $2$.
% This can be derived from the uniqueness statement of
% Theorem~11.2 if the second-order scalar problem
% (4.17) is written as a first-order problem in two variables.  \qed
% \smallskip
% </latex>

%%
% Here are three images showing solutions of
% (4.17) with subcritical, critical, and supercritical damping.
% Critical damping is plainly the fastest.
L = chebop(0,20); L.lbc = [1;0]; om = 1; epep = [0.1 2 10];
for j = 1:3
  subplot(1,3,j); ep = epep(j);
  L.op = @(t,y) diff(y,2) + ep*diff(y) + om^2*y;
  y = L\0; plot(y,CO,ivp), axis([0 20 -1.2 1.2])
  title(['$\varepsilon = ' num2str(ep) '$, $\omega = ' ...
     num2str(om) '$'],FS,11,IN,LT)
  r1 = -ep/2 + sqrt(ep^2/4-om^2); r2 = -ep/2 - sqrt(ep^2/4-om^2);
  text(1.3,-.91,['$r_+^{} = ' num2str(r1) '$'],FS,8,IN,LT)
  text(1.3,-1.10,['$r_-^{} = ' num2str(r2) '$'],FS,8)
  ss = {'subcritical','critical','supercritical'};
  text(18,.9,ss{j},FS,10,HA,RT)
end
%%
% <latex> \vskip 1.02em </latex>


%%
% <latex>
% One of the most important and interesting aspects of
% oscillatory systems is their response to external
% inputs --- the inhomogeneous variant of (4.17).
% This will be the subject of Chapter 8.
% </latex>

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

%%
% <latex>
% Here is an example adapted from D. B. Meade, ``ODE models for
% the parachute problem,'' {\em SIAM Review} 40 (1998), pp.~327--332.
% </latex>

%%
% If there were no air resistance,
% a falling skydiver would experience only
% the force of gravity: $h''(t)=-g$, where $h$ is height above the earth
% and $g=9.8$. (We use 
% MKS units, i.e., meters-kilograms-seconds.)
% The speed of fall would
% increase linearly and the height would decrease quadratically.

%%
% <latex>
% However, the air resistance quickly becomes significant as
% the skydiver's speed
% increases. The atmosphere creates a force called {\em drag } that opposes gravity.
% For reasons of fluid mechanics that we will not go into, this force
% can be reasonably modeled as $F_D = k v^2$, where $v=h'$ is the velocity
% and $k$ (in units kg/m) depends on properties of the air
% and the skydiver.\footnote{Thus the
% drag force here is different from that of (4.17), which
% depends just linearly on the velocity.  As a rule one sees linear
% drag for bodies moving slowly through a fluid (laminar flow) and
% quadratic drag at higher speeds (turbulent flow).}
% The gravity and drag forces
% are in balance when they have equal magnitudes, giving the condition
% $mg=kv^2$, where $m$ is the mass of the skydiver. Since there is no net
% force on the skydiver under this condition, the velocity will not change,
% and this value $v_T^{}=-\sqrt{mg/k}$ is known
% as the {\em terminal velocity.}
% </latex>

%%
% If we express Newton's law for the skydiver in terms of the velocity, we
% have 
% $$ v'=\frac{kv^2}{m} - g = g \left(\frac{v^2}{v_T^2} - 1\right), $$ 
% which is separable.  Integration leads to 
% $$ v(t) = -v_T^{} \kern .5pt 
% {\rm tanh}\kern -1pt \left( \frac{gt}{v_T^{}} + C \right) $$
% for an arbitrary constant $C$. Since this is an explicit function of $t$,
% $v$ can be integrated again to get $h(t)$, but we will omit that
% result here and move to numerical solutions for simplicity.
% A reasonable approximation is
% $v_T^{}\approx 55$ m/s for the standard ``belly-down'' position
% in free fall. This enables us to calculate height and velocity for a dive
% starting at 4 km altitude.
clf, g = 9.8; vT = 55; h0 = 4000;
L = chebop(@(t,h) diff(h,2)-g*(diff(h)^2/vT^2 - 1),[0 75]);
L.lbc = [h0;0]; hfree = L\0; plot(diff(hfree),CO,ivpnl)
xlabel('time (s)',FS,9), ylabel('velocity (m/s)',FS,9)
title('Fig.~4.6.~~Skydiver velocity during free fall',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Terminal velocity is reached in about 15 seconds.
% </latex>

%%
% The reason skydivers need parachutes is that without them,
% their terminal velocities would be terminal!  The function of a 
% parachute is to greatly increase the drag and thereby 
% induce a much smaller terminal velocity, which we may
% take to be, say, 10 m/s.
% A typical altitude at
% which to open the parachute is $h=800{\rm m}$, which occurs 
% for our model at the time
tp = roots(hfree-800)

%%
% <latex>
% \noindent
% So our skydiver can enjoy a whole minute
% of free fall.  Then, with the parachute open, the
% second phase of her dive is given by
% </latex>
vT = 10; h0 = 800;
L = chebop(@(t,h) diff(h,2)-g*(diff(h)^2/vT^2 - 1),tp+[0 90]);
L.lbc = [hfree(tp);deriv(hfree,tp)]; hchute = L\0;  
plot(diff(hchute),CO,ivpnl)
xlabel('time (s)',FS,9), ylabel('velocity (m/s)',FS,9)
title('Fig.~4.7.~~Skydiver velocity after parachute opens',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The dive ends with $h=0$, at time
% </latex>
t0 = roots(hchute)

%%
% <latex>
% \noindent
% We can plot the entire course of the dive by concatenating the two
% phases.
% </latex>
h = join(hfree{0,tp},hchute); plot(h,CO,ivpnl)
xlabel('time (s)',FS,9), ylabel('height (m)',FS,9)
title('Fig.~4.8.~~The whole skydive',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Taking the derivative gives the velocity as a function of time.
% </latex>
plot(diff(h),CO,ivpnl)
xlabel('time (s)',FS,9), ylabel('velocity',FS,9)
title('Fig.~4.9.~~Skydiver velocity',FS,11), ylim([-70 10])
text(4,-8,'jump from the plane'), text(tp+2,-47,'parachute opens')
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% {\sc History.}  Much of this chapter goes back to
% the brilliant 17th century rivals Robert Hooke, Christiaan 
% Huygens, and
% Isaac Newton.  The names of Hooke and Newton are associated
% with the two most
% famous force laws $F$ in the equation $F(y) = my''$.  Hooke's
% law puts $F(y) $ proportional 
% to $y$, leading to simple harmonic motion and elasticity.
% Huygens studied such oscillations too and was the inventor of
% the pendulum clock.
% Newton's law for gravity puts $F(y)$ proportional to $y^{-2}$, 
% leading to elliptical orbits and cosmology.  When $y$ is effectively
% constant, as with gravity at the surface of the earth,
% so is the force $F$ --- the situation experienced
% by our skydiver and also by an apple falling from Newton's tree.
% \smallskip
% </latex>

%%
% <latex>
% {\sc Our favorite reference.}  For a salutary view of
% how ODE\kern .5pt s meet reality, see J. C. Dixon, {\em The Shock
% Absorber Handbook,} 2nd ed., Wiley, 2007.  Chapter 2, on 
% vibration theory, presents the material we have discussed
% in the vivid context of behavior of motor vehicles.  To optimize
% comfort, passenger cars are usually
% designed in the underdamped regime, whereas to
% optimize handling, race cars are closer to critically damped.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 4.}
% The scalar, linear, autonomous second-order ODE\/ $y'' +
% \varepsilon y'+\omega^2\kern .5pt y=0$
% has a two-dimensional vector space of solutions.  If\/~$\omega$ is
% real and $\varepsilon>0,$ all solutions decay exponentially.
% For $\varepsilon <2\kern .5pt\omega$ (subcritical damping),
% the solution is an oscillatory function times
% an exponential envelope, and for $\varepsilon >2\kern .5pt\omega$ (supercritical
% damping) it is a superposition of two exponentials with
% different time constants. \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em \underline{Exercise $4.1$}.  Hooke and Newton elliptical orbits.}
% Taking $y(t)$ to be complex as usual for tracking orbits in a plane,
% we can write $y'' = -y|y|^{p-1}$ for the classic problem of 
% a point orbiting in a central field.
% Produce seven plots showing the orbits 
% over $t\in [\kern .3pt 0,40\kern .3pt ]$ emanating from
% the initial conditions $y(0) = 1,$ $y'(0) = 0.5i$ for
% $p = -2, -1.5, -1, \dots, 1$; in each case use
% \verb|axis([-1 1 -1 1]),| \verb|axis square|.
% The case $p=1$ corresponds to Hooke's law, and the orbit is an ellipse with
% center $y=0$.  The case $p=-2$ corresponds to Newton's law of
% gravitation, and the orbit is an ellipse with one focus at $y=0$. 
% (V.~I. Arnol'd uses
% the terms {\em Hooke ellipse\/} and {\em Newton ellipse.}  If the ellipses
% are viewed as sets in the complex plane, then the square of
% a Hooke ellipse is a Newton ellipse.)
% \par
% {\em \underline{Exercise $4.2$}. Kepler's equal-area law.}
% For a particle in an elliptical
% orbit in an inverse-square field, the area swept out per 
% unit time is constant.  
% {\em (a)} If $y(t)$ is an orbit computed as in
% Exercise~4.1, show mathematically why the command
% \verb|A = cumsum(imag(diff(y)*conj(y)))| computes a chebfun representing
% the area $A(t)$ swept out as a function of time $t$
% (\verb|cumsum| is Chebfun's indefinite integral operator).
% {\em (\kern .7pt b)} Perform this calculation and plot
% the result (with $p=-2$ for the inverse-square
% force).  What is $A(40)\kern .5pt?$
% \par
% {\em Exercise $4.3$.  Exact equations.}  Sometimes an ODE is a derivative
% of an equation of lower order.  For example, the equation
% $y'' + t\kern .3pt y' + y = 0$ can be written $(y'+t\kern .3pt y)' = 0$.
% Use this observation to find the analytical solution to this ODE
% {\em (a)} in general, and {\em (\kern .7pt b)} with $y(0)=1$, $y'(0) = 1$.
% \par
% {\em Exercise $4.4$.  Falling to earth.}  Suppose
% Earth is a uniform sphere of mass $M$
% and radius~$R$, and a small body
% of mass $m\ll M$ is initially motionless
% at a height $h(0) = h_0$ above the surface
% and attracted downward by the gravitational force
% $F = GmM/(R+h)^2$, where $G$ is the universal
% gravitational constant.
% {\em (a)} Write down an ODE IVP whose solution gives
% $h(t)$ as a function of $t$.
% {\em (\kern .7pt b)} Solve it analytically to get a formula
% for~$t_c^{}$, the time at which the body hits the surface.
% {\em (c)} Rewrite this formula in terms of explicit numbers
% using the approximate values $M = 6\times 10^{24}$ kg,
% $R = 6 \times 10^6$ m, and
% $G = 6\times 10^{-11} \hbox{m}^3 \hbox{kg}^{-1} \hbox{s}^{-2}$.
% {\em (d)} What height $h_0$ leads to $t_c^{} = 1 \hbox{ minute}\kern 1pt$?
% \par
% {\em Exercise $4.5$.  Second solution for critical damping.}
% Confirm analytically that $y(t) = t e^{rt}$ is a solution of (4.10)
% in the critically damped case $\varepsilon = 2\kern .5pt\omega$.
% \par
% {\em \underline{Exercise $4.6$}.  Critical damping example.}
% {\em (a)} What value $\omega=\omega_c^{}$ in
% (4.9) gives critical damping?
% {\em (\kern .7pt b)} Fix $\omega$ at this critical value
% and compute solutions with the $\varepsilon$ parameter $0.1$
% of (4.9) taking values $0.09$, $0.1$, and $0.11$.  Superimpose
% semilogy plots of $|y(t)|$ against $t$ for these three choices.
% Comment on why one of the three curves
% has a shape different from the others, and on why two of them
% have nearly the same slope.
% \par
% {\em \underline{Exercise $4.7$}.  A weak nonlinear spring.}
% Suppose a particle with rest position $y=0$ is subject to an
% asymmetric restoring force $y'' = - y\exp(-y)$. 
% {\em (a)} Plot $y(t)$ for $t\in [\kern .3pt 0,200\kern.3pt]$
% for initial conditions $y(0) = 0$
% and $y'(0) = 0.5$, $1$, and $1.5$.  
% {\em (\kern .7pt b)} Define the idea of an ``escape velocity'' for this system for
% trajectories beginning at $y=0$.
% Based on the experiments of {\em (a)}, and other experiments if you
% wish, what can you say about its numerical value?
% {\em (c)} By thinking in terms of kinetic and potential energy, determine
% the escape velocity analytically.  (The potential energy can be
% defined as the work required to bring the particle in from
% $\infty$, where work is defined as an integral of force times
% distance.)
% \par
% {\em \underline{Exercise $4.8$}.  Acceleration of the skydiver.}
% Repeat the calculation for the skydiver application and
% plot the acceleration as a function of time.  What is the maximum
% acceleration in MKS units, and as a multiple of the
% gravitational acceleration $g\kern .5pt$?  This number is
% unrealistically large, and would probably kill the skydiver.
% Explain in a general way what kind of effects you think
% have been omitted from the model that in practice would
% make the acceleration not so extreme.
% \par
% {\em \underline{Exercise $4.9$}.  Hard nonlinearity}.
% Most springs get stiffer
% at high displacements. To see the effect on the
% period $T$, compute and plot $T$ for the equation $y'' + y + y^3=0$
% for solutions with initial velocity $y'(0)=0$ and 
% initial amplitudes $y(0) = 0.0, 0.1, \dots, 5$.
% Explain why the period changes as $y(0)$ increases.
% \par
% {\em \underline{Exercise $4.10$}.  Friction and sliding of a penny.}
% A penny of mass 1
% sits at position $y(t)$, starting motionless at $y(0) = 0$, on a wooden
% board inclined at angle~$\theta$ from the horizontal.  The board
% is slowly tilted up according to $\theta = 0.1t$ for $t\in [\kern 0.3pt
% 0,10\kern .3pt ]$.  The gravitational force downward on the penny
% can be resolved into $f_t^{} = \sin(\theta)$ in the $y$ direction along the board
% and $f_n^{} = \cos(\theta)$ normal to the board. Meanwhile there is a 
% frictional force $f_f^{} =- \mu_k^{} f_n$ on the penny,
% where $\mu_k^{}=0.5$ is a coefficient
% of kinetic friction, so long as the penny is sliding.  On the other
% hand if
% $f_t^{}\le \mu_s^{} f_n$, where $\mu_s^{}=1$ is a coefficient
% of static friction, then the penny is stationary, with
% $f_f^{} = -f_t^{}$.
% {\em (a)} Formulate this problem as an IVP in Chebfun and
% compute and plot the solution $y(t)$.  Where is the penny at $t=10\kern .5pt$?
% How many times continuously differentiable is the trajectory $y(t)$?
% {\em (\kern .7pt b)} Repeat both parts of {\em (a)} under the assumption 
% $\mu_s^{} = \mu_k^{} = 0.5$.
% \par
% {\em Exercise $4.11$.  Nondimensionalization.}  Show how
% by introducing a variable $u$ equal to a suitable multiple
% of $y$, the equation $y'' + ay' + by = 0$ can be reduced to
% $u'' + \varepsilon u' + u = 0$.  What is $\varepsilon$ in
% terms of $a$ and $b\kern .5pt$?
% </latex>
