% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{hyperref}
\usepackage{url}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\Projectile Motion with Altitude-Dependent Gravity}



\markboth{Projectile Motion with Altitude-Dependent Gravity}{Projectile Motion with Altitude-Dependent Gravity}

\author{Jakob Kotas\\
Menlo College\\
Atherton CA 94027 USA\\
jakob.kotas@menlo.edu
}

\makeStitlePDFLaTex


\begin{abstract}
Basic projectile motion without air resistance typically assumes gravity is constant. In reality, the acceleration due to gravity is proportional to the inverse-square of the distance between the centers of mass of the Earth and the projectile. When projectiles are near to Earth's surface, this deviation is negligible, but it becomes important when dealing with projectiles at high altitudes. We lay out a technique narrative for such a case which is a second-order, linear, constant-coefficient, and non-homogeneous ordinary differential equation (ODE). Through the course of the scenario, students will solve the ODE analytically and investigate its behavior numerically. They are also walked through the idea of linearizing a nonlinear ODE and discover the pitfalls in solving the approximate version. Along the way, they are introduced to the hyperbolic sine and cosine functions.
\end{abstract}



\section*{INTRODUCTION AND SCENARIO DESCRIPTION}
In the most basic form of projectile motion, we assume that there is no air resistance, and that gravity is a constant $g$. Then, the ordinary differential equation (ODE) that describes the resultant vertical motion is:

\begin{equation*}
    y'' = g
\end{equation*}

\noindent where $y(t)$ is the height of the projectile at time $t$. Let us use the convention that $y>0$ are heights above the surface and correspondingly $g<0$. Let us also assume initial conditions:

\begin{equation*}
    y(0) = 0, \hspace{1cm} y'(0) = v_0
\end{equation*}

\noindent where the object is launched from the surface $y=0$, and that $v_0>0$ is a known upward initial velocity. Let us write $y'$ as $v$, so that $v$ is the vertical velocity. We then have:

\begin{equation*}
    v' = g
\end{equation*}

{\bf 1.} Integrate both sides of this equation with respect to dt, and apply one initial condition to solve for the unknown constant. Then change $v$ back to $y'$, integrate again, and apply the other initial condition.

\teacher{
\textcolor{red}{
\noindent Integrate once:
\begin{equation*}
    v = gt + k_1
\end{equation*}
\noindent Apply one initial condition:
\begin{equation*}
    v(0) = k_1 = v_0
\end{equation*}\begin{equation*}
    v = gt + v_0
\end{equation*}\begin{equation*}
    y' = gt + v_0
\end{equation*}
\noindent Integrate again:
\begin{equation*}
    y = \frac{1}{2} gt^2 + v_0t + k_2
\end{equation*}
\noindent Apply other initial condition:
\begin{equation*}
    y(0) = k_2 = 0
\end{equation*}\begin{equation*}
    y = \frac{1}{2} gt^2 + v_0t
\end{equation*}
}}

To check your answer for part {\bf 1}, you should have gotten:

\begin{equation*}
    y(t) = \frac{1}{2}gt^2 + v_0 t,
\end{equation*}

\noindent which is a well-known result in physics. Among the many modeling assumptions that goes into this ODE is that gravity $g$ is a constant. In reality, we know that gravity is a function of height above the surface. In particular, as per Newton's law of universal gravitation, it is proportional to the inverse-square of the distance $r$ between the center of mass of the Earth (or whatever heavenly body we happen to be on), and the center of mass of our projectile.

Let us call the gravity our projectile experiences $\gamma = \gamma(r)$.

\begin{equation*}
    \gamma \propto \frac{1}{r^2}
\end{equation*}

Since we are calling the surface of the Earth $y=0$, we can shift coordinates by arguing that $r=y+R$ where $R$ is the radius of the Earth and $y>0$. That gives us:

\begin{equation*}
    \gamma \propto \frac{1}{(y+R)^2}
\end{equation*}

\noindent and the constant of proportionality can be found since we know $\gamma = g$ when $y = 0$ (or equivalently, when $r = R$).

\begin{equation*}
    \gamma = \frac{gR^2}{(y+R)^2}
\end{equation*}

If we go back to our original ODE and replace $g$ with $\gamma$, we have:

\begin{equation*}
    y'' = \frac{gR^2}{(y+R)^2}.
\end{equation*}

\noindent which is all well and good. But this ODE is very hard to solve, because it is nonlinear!

{\bf 2.} Why is this ODE nonlinear?

\teacher{
\textcolor{red}{
\noindent It cannot be written in the form
\begin{equation*}
    a(t)y'' + b(t)y' + c(t)y = d(t)
\end{equation*}
because the $y$ on the right-hand side is squared and in the denominator of the fraction.
}}

Let us linearize the ODE to make it easier to solve. The solution we get then will be an approximate solution, but we expect it to still be reasonable as long as $y\approx 0$. The benefit is that the linear version is much easier to solve.

{\bf 3.} Linearize the right-hand side of our ODE, which we will call $f(y)$. Our linearization is the first-order Taylor series, expanding around $y=0$.

\begin{equation*}
    f(y) \approx f(0) + f'(0)y
\end{equation*}

\noindent within some interval of convergence around $y=0$. Remember that $y$ is our independent variable and $g$ and $R$ are constants.

\teacher{
\textcolor{red}{
\begin{equation*}
    f(y) \approx g - \frac{2gy}{R}
\end{equation*}}
}

{\bf 4.} Use this new $f(y)$ as our right hand side; this is the equation we'll be using for the rest of this exercise. Characterize this equation: What is its order? Is it linear or nonlinear? Homogeneous or non-homogeneous? 

\teacher{
\textcolor{red}{
\begin{equation*}
    y'' + \frac{2g}{R} y = g
\end{equation*}
Second-order, linear (by design), non-homogeneous.}
}

{\bf 5.} Use your knowledge of this type of ODE to find the general solution. Your answer should be in terms of the constants $g$ and $R$.

\teacher{
\textcolor{red}{
\noindent We use the method of undetermined coefficients.
\noindent Homogeneous equation:
\begin{equation*}
    y'' + \frac{2g}{R} y = 0
\end{equation*}
Characteristic equation: let $y_h=e^{\lambda t}$
\begin{equation*}
    \lambda^2 + \frac{2g}{R} = 0
\end{equation*}\begin{equation*}
    \lambda = \pm \sqrt{\frac{-2g}{R}}
\end{equation*}
Recall that $g<0$ and $R>0$ therefore both $\lambda$ are pure real. Homogeneous solution:
\begin{equation*}
    y_h = c_1 e^{\left(\sqrt{-2g/R}\right) t} + c_2 e^{-\left(\sqrt{-2g/R}\right) t}
\end{equation*}
Particular solution: let $y_p=k$, an unknown constant.
\begin{equation*}
    0 + \frac{2g}{R} k = g
\end{equation*}\begin{equation*}
    k = \frac{R}{2}
\end{equation*}
Let $\omega = \sqrt{\frac{-2g}{R}}$. Solution:
\begin{equation*}
    y = y_h + y_p = c_1 e^{\omega t} + c_2 e^{-\omega t} + \frac{R}{2}
\end{equation*}
}}

Next, we can rewrite the answer to part {\bf 5} which will make the algebra a bit easier later on. We will make use of the hyperbolic trig functions ``hyperbolic sine" (sinh) and ``hyperbolic cosine" (cosh), which are defined as follows:

\begin{equation*}
    \sinh\left( x \right) = \frac{e^x - e^{-x}}{2}
\end{equation*}
\begin{equation*}
    \cosh\left( x \right) = \frac{e^x + e^{-x}}{2}
\end{equation*}

Using these functions, we can rewrite something which looks like:
\begin{equation*}
    c_1 e^x + c_2 e^{-x}
\end{equation*}

Split each term into two halves:
\begin{equation*}
    \left( \frac{c_1}{2} e^x + \frac{c_1}{2} e^x \right) + \left( \frac{c_2}{2} e^{-x} + \frac{c_2}{2} e^{-x} \right)
\end{equation*}

Include some terms which add up to zero:
\begin{equation*}
    \left( \frac{c_1}{2} e^x + \frac{c_1}{2} e^x \right) + \left( \frac{c_2}{2} e^{-x} + \frac{c_2}{2} e^{-x} \right) + \underbrace{\left(\frac{c_2}{2} e^x - \frac{c_2}{2} e^x \right)}_{0} + \underbrace{\left(\frac{c_1}{2} e^{-x} - \frac{c_1}{2} e^{-x} \right)}_{0}
\end{equation*}

Factor:
\begin{equation*}
    = \frac{c_1 + c_2}{2} \left(e^x - e^{-x}\right) + \frac{c_1-c_2}{2}\left(e^x+e^{-x}\right)
\end{equation*}

Rename:
\begin{equation*}
    \tilde{c_1} = c_1+c_2, \hspace{1cm} \tilde{c_2} = c_1-c_2
\end{equation*}

And finally we have:
\begin{equation*}
    = \tilde{c_1} \cosh(x) + \tilde{c_2} \sinh(x)
\end{equation*}

{\bf 6.} Use the relationship above to rewrite your answer to part {\bf 5} in terms of sinh and cosh.

\teacher{
\textcolor{red}{
\begin{equation*}
    y = \tilde{c_1} \cosh\left( \omega t \right) + \tilde{c_2} \sinh \left( \omega t \right) + \frac{R}{2}
\end{equation*}
}}

{\bf 7.} Apply the initial conditions $y(0)=0$, $y'(0)=v_0$ to solve for the unknown constants in terms of $v_0$. Then rewrite the solution under these initial conditions. (Hint: some possibly useful facts about hyperbolic trig functions are given below.)

\begin{equation*}
    \sinh(0) = 0 \hspace{2cm} \cosh(0) = 1
\end{equation*}
\begin{equation*}
    \frac{d}{dx} \sinh(x) = \cosh(x) \hspace{2cm} 
    \frac{d}{dx} \cosh(x) = \sinh(x)
\end{equation*}

\teacher{
\textcolor{red}{
\begin{equation*}\displaystyle
    \left\{
        \begin{array}{rcccl}
            y(0)&=&\tilde{c_1} + \frac{R}{2}&=&0\\
            y'(0)&=&\omega \tilde{c_2}&=&v_0
        \end{array}
    \right.
\end{equation*}
\begin{equation*}\displaystyle
    \left\{
        \begin{array}{rcl}
            \tilde{c_1}&=&-\frac{R}{2}\\
            \tilde{c_2}&=&\frac{v_0}{\omega}\\
        \end{array}
    \right.
\end{equation*}
General solution:
\begin{equation*}
    y = \frac{v_0}{\omega} \sinh \left( \omega t \right) + \frac{R}{2} \Big( 1 - \cosh \left( \omega t \right) \Big)
\end{equation*}
}}

Congratulations, you've solved the equations of motion! Now let's investigate how the graph looks visually and how it compares to the equations of motion when gravity was a constant. For starters, you can take the following values for the parameters:

\begin{center}
$R=6.38\times 10^6$ m\\
$g=-9.8$ m/s$^2$\\
$v_0=1000$ m/s
\end{center}

{\bf 8a.} Based on the graph, what is the highest altitude the projectile reaches? Compare this to our old equations of motion with constant gravity, given below. What was the highest altitude in that case?

\begin{equation*}
    y = \frac{1}{2}g t^2 + v_0 t
\end{equation*}

\teacher{
\textcolor{red}{
Figure below shows trajectory with constant gravity (blue, dashed) and variable gravity (red, solid). Horizontal axis is time in seconds. Vertical axis is height in km.
\begin{center}
\includegraphics[width=0.4\textwidth]{5a.png}
\end{center}
With constant gravity: max altitude = 51.020 km\\
With variable gravity: max altitude = 51.435 km\\
We see little difference between the two trajectories.
}}

{\bf b.} Answer the same questions as part {\bf 8a} but with $v_0=5000$ m/s instead. What do you notice?

\teacher{
\textcolor{red}{
Figure below shows trajectory with constant gravity (blue, dashed) and variable gravity (red, solid). Horizontal axis is time in seconds. Vertical axis is height in km.
\begin{center}
\includegraphics[width=0.4\textwidth]{5b.png}
\end{center}
With constant gravity: max altitude = 1275.510 km\\
With variable gravity: max altitude = 1762.294 km\\
There is a larger discrepancy between the two graphs.
}}

{\bf c.} Answer the same questions as {\bf 6a} but with $v_0=6000$ m/s instead. What do you notice here? The answer should be quite different than parts {\bf 8a} and {\bf b}! Why does this happen?

\teacher{
\textcolor{red}{
Figure below shows trajectory with constant gravity (blue, dashed) and variable gravity (red, solid). Horizontal axis is time in seconds. Vertical axis is height in km.
\begin{center}
\includegraphics[width=0.4\textwidth]{5c.png}
\end{center}
With constant gravity: max altitude = 1836.735 km\\
With variable gravity: max altitude unbounded\\
We are violating the assumption that $y\approx 0$ that we made in linearizing the ODE. In our linearized version, the sign of gravity switches and the object gets pulled upward if it surpasses a certain height. In particular that height is when:
\begin{equation*}
    f(y) = g - \frac{2gy}{R} = 0
\end{equation*}\begin{equation*}
    y = \frac{R}{2}.
\end{equation*}
Therefore we only expect our approximate solution to be valid when $y \ll R/2$. To fix this issue we would need to consider the full nonlinear ODE given these initial conditions to get qualitatively correct behavior.
}}

{\bf d.} (Bonus!) What is the condition on $v_0$ (in terms of $g$ and $R$) where the behavior observed in {\bf 6c} will not occur?

\teacher{
\textcolor{red}{
We stated earlier that the ODE solution is:
\begin{equation*}
    y = \frac{v_0}{\omega} \sinh \left( \omega t \right) + \frac{R}{2} \Big( 1 - \cosh \left( \omega t \right) \Big)    
\end{equation*}
Rewriting this in terms of exponentials and grouping like terms gives:
\begin{equation*}
    y = \left( \frac{v_0}{2\omega} - \frac{R}{4} \right) e^{\omega t} - \left( \frac{v_0}{2\omega} + \frac{R}{4} \right) e^{-\omega t} + \frac{R}{2}
\end{equation*}
We want $\displaystyle \lim_{t\rightarrow\infty} y = -\infty$, and we have that $v_0$, $R$, $\omega$ are all strictly positive. The $e^{-\omega t}$ term becomes negligible as $t\rightarrow\infty$, and $R/2$ is a finite constant. Thus the only term we must worry about is the $e^{\omega t}$ term. If the coefficient of $e^{\omega t}$ is positive, then $\displaystyle \lim_{t\rightarrow\infty} y = \infty$, and if the coefficient of $e^{\omega t}$ is zero, then $\displaystyle \lim_{t\rightarrow\infty} y = R/2$. We only have $\displaystyle \lim_{t\rightarrow\infty} y = -\infty$ if the coefficient of $e^{\omega t}$ is negative. Thus our condition is:
\begin{equation*}
    \frac{v_0}{2\omega} - \frac{R}{4} < 0
\end{equation*}\begin{equation*}
    v_0 < \sqrt{\frac{-gR}{2}}
\end{equation*}
Using the values of $g$ and $R$ we chose, that gives:
\begin{equation*}
    v_0 \lessapprox 5591.2 \textrm{ m/s}
\end{equation*}
The quantity $\sqrt{-gR/2}$ \textit{resembles} an ``escape velocity" of sorts for the system.
}}




\teacher{
\section*{NOTES FOR TEACHERS}
This activity is helpful in reinforcing the solution technique of undetermined coefficients for second-order, linear, non-homogeneous equations. It would fit best immediately after discussion of such equations.

Solving in terms of parameters rather than numerically may be a challenge for some students. The fact that we have defined $g<0$ leads to the quantity $-2g/R$ being positive and therefore $\sqrt{-2g/R}$ is real; this notation reinforces the need to consider parameters' sign when a value is not given. (Similarly, the real solutions are exponentials or sinh and cosh, as opposed to sin and cos had $g$ been positive.)

Our decision to rename $\sqrt{-2g/R}$ as $\omega$ in the teacher version leads to cleaner step-by-step solutions. This is not pointed out in the student version. The instructor may wish to guide students toward this idea if they get bogged down in algebra.

This activity could be self-guided or completed in small groups. Anticipate spending one class day ($\approx 60$ minutes) on it. Students will need a graphing utility for parts {\bf 8a-c}.



%\begin{thebibliography}{99}




%\bibitem{Gause1934} Gause, G.~F. 1971. {\it The Struggle for Existence}. New York:  Dover Publications, Inc. First published in 1934 by The Williams \& Wilkins Company and available completely on the world wide web at \url{http://www.ggause.com/Contgau.htm}. Accessed  31 March 2008.



%\bibitem{Hutchinson1978a} Hutchinson, G. Evelyn. 1978. {\it Introduction to Population Ecology\/}. New Haven CT:  Yale University Press.

%\end{thebibliography}
}





\end{document}
