%% 4. Second-order equations and damped oscillations %% % % \setcounter{page}{35} % %% % % 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)$. % %% % % 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$. % %% % % 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$. % 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) %% % \vskip 1.02em %% % 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) %% % \vskip 1.02em %% % % \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$. % %% % % 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$. % 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) %% % \vskip 1.02em %% % % \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. % %% % % 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. % %% % 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. %% % % 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. % %% % % 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$. % %% % % 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. % 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 %% % \vskip 1.02em %% % % \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: % om = 1; T = 2*pi/sqrt(om^2 - ep^2/4) %% % % 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. % %% % % 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)$$ % %% % Here is a theorem to summarize these observations. %% % % \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)$$ % } % %% % % \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 % %% % 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 %% % \vskip 1.02em %% % % 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. % %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: skydiver}\\[-3pt] % \hrulefill % \end{center} % %% % % Here is an example adapted from D. B. Meade, ODE models for % the parachute problem,'' {\em SIAM Review} 40 (1998), pp.~327--332. % %% % 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. %% % % 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.} % %% % 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) %% % \vskip 1.02em %% % % \noindent % Terminal velocity is reached in about 15 seconds. % %% % 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) %% % % \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 % 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) %% % \vskip 1.02em %% % % \noindent % The dive ends with $h=0$, at time % t0 = roots(hchute) %% % % \noindent % We can plot the entire course of the dive by concatenating the two % phases. % 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) %% % \vskip 1.02em %% % % \noindent % Taking the derivative gives the velocity as a function of time. % 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') %% % \vskip 1.02em %% % % {\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 % %% % % {\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 % %% % % \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} % %% % % \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$? %