% 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 
	%\footnote{DELETE THIS FOOTNOTE BEFORE SUBMISSION.  In the \TeX \, source code there is a toggle/option to produce a Teacher %Version or a Student Version by ``carving'' out with  $\backslash$teacher\{\} that material which should only appear in the %Teacher version. See the source code.} 
	\else STUDENT VERSION \fi\\Euler and How Far Can You Throw a Ball?}



\markboth{Euler and How Far Can You Throw a Ball?}{Euler and How Far Can You Throw a Ball?}

\author{Chris McCarthy\\
BMCC City University of New York\\
New York NY 10007 USA\\
cmccarthy@bmcc.cuny.edu
}

\makeStitlePDFLaTex


\begin{abstract}
If a tennis ball is thrown through the air it will eventually hit the ground due to gravity. Using Euler's method, write a short script (Python, Matlab, R, etc.) to find the trajectory of the ball which will maximize the distance the ball lands from the thrower taking into account air resistance (drag).



 This scenario can be used to introduce students to Euler's method for higher order ordinary differential equations. To do this project students should be comfortable with the concept of second order ODE's, vectors, and  $F = ma$ from elementary physics.  A little programming knowledge is a plus. But beyond that, not much is needed. The student version contains the necessary physics background, including a brief
 discussion of drag, as well as a quick introduction to Euler's method for higher order ODE's. A Python script to solve the problem, if we neglect drag, is provided to the students. Their job is to modify the script to take into  account drag. A full script, which solves the problem and includes drag is provided for teachers. 
\end{abstract}


\section*{SCENARIO DESCRIPTION}

 If a tennis ball of mass $m$  is thrown through the air it will eventually hit the ground due to gravity.  If you can 
 throw a tennis ball 12 meters/second (about 26.8 mph) how far can you throw it; meaning how far away from you can you make it land? Assume that when the ball leaves your hand it is at a height of 2 meters (about 6 feet).  With the help of Euler's Method, write a short script (Python, Matlab, R, etc.) to find the ideal launch angle $\theta$ to throw the ball,  so that it will result in the ball landing as far as possible from you. Also find that  maximal distance and plot that trajectory. Make sure to include air resistance (drag) in your model. Don't worry about rotational effects.

\medskip

\textbf{Physics Background and Euler's Method Review.}   Let $x(t) $ and $y(t)$ represent a ball's position at time $t$, in the $x,y$ plane, with $y$ being height. Vectors will be denoted by a \textbf{bold} font.   The ball's velocity is: 
\begin{equation}\label{VelDef1}
\textbf{v}  =  \dfrac{d}{dt}\langle x, y \rangle =  \langle \dot{x}, \dot{y} \rangle =  \langle v_x, v_y \rangle.
\end{equation}
The ball's acceleration is:
\begin{equation}\label{AccelerationDef1}
\textbf{a}  =  \dfrac{d}{dt}\langle v_x, v_y \rangle  = \langle \ddot{x}, \ddot{y} \rangle = \langle \dot{v}_x, \dot{v}_y \rangle.
\end{equation}
From Newton's Second Law of Motion we know that Force  = mass $\times$ acceleration: 
\begin{equation}\label{Newton}
\textbf{F} = m\textbf{a}
\end{equation}
For our models we will approximate the earth as being flat with gravity pointing straight down (in the negative $y$ direction) and so we will  approximate the force of gravity acting on a mass $m$ as:
\begin{equation}\label{G}
\mathbf{F_{g}} = \langle 0, -mg \rangle
\end{equation}
with 
$g = 9.8\ \text{meters}/\text{second}^2$ in the MKS system.
%\footnote{MKS = Meters Kilogram Seconds}   

\medskip

\textbf{Simple Model for a Thrown Ball (no drag).} 
To start with, we will ignore drag. Once we understand this simple model, we can make it  more realistic by adding drag.

Combining  (\ref{AccelerationDef1}), (\ref{Newton}), and (\ref{G}) we get the following second order ODE (which neglects drag):
\begin{equation} \label{ODE0}
\mathbf{F_{g}} = \langle 0, -mg \rangle = m \langle \ddot{x}, \ddot{y} \rangle = m\textbf{a}
\end{equation}
Writing the $\ddot{x}$ and $\ddot{y}$  parts of (\ref{ODE0}) on separate lines we get the 2nd order system of ODE's:
\begin{align} 
m \ddot{x} &=   0    \nonumber  \\ 
m\ddot{y} &= -mg \label{ODE1}   
\end{align}
To apply Euler's method, we need to convert $(\ref{ODE1})$ to a 
first order system of ODE's. We do the usual substitution.  As in (\ref{VelDef1}) 
and (\ref{AccelerationDef1}),  we let 
$\dot{x} =  v_x$, so that $\ddot{x} =  \dot{v}_x$; 
we let
$\dot{y} =  v_y$, so that $\ddot{y} =  \dot{v}_y$;
and we divide by $m$. This gives us  the  first order version of (\ref{ODE1}):
\begin{align}
\dot{x} &=  v_x \nonumber \\
\dot{y} &=  v_y  \nonumber \\
\dot{v}_x &=   0  \nonumber \\ 
\dot{v}_x &=  -g  \label{ODE2}
\end{align}
We can rewrite (\ref{ODE2}) in vector form as:
\begin{equation} \label{ODE3}
\dfrac{d}{dt}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y
\end{pmatrix}
=
\begin{pmatrix}
v_x  \\
v_y \\
0 \\
-g
\end{pmatrix}
\end{equation}
%\begin{pmatrix} requires amslatex package
%	a & b \\
%	c & d
%\end{pmatrix}
We can linearly approximate the state of the system 
$\begin{pmatrix}
x  \\
y \\
v_x \\
v_y
\end{pmatrix}
$ at time $t + \Delta t$ using
information about the system at time $t$:
\begin{equation} \label{ODE4}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y
\end{pmatrix}_{\text{at} \ t + \Delta t}
\approx
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y
\end{pmatrix}_{\text{at} \ t}
+
\dfrac{d}{dt}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y
\end{pmatrix}_{\text{at} \ t}
\cdot \Delta t
\end{equation}
It is convenient\footnote{Explicitly including time  allows us to easily keep track of time while doing Euler's method calculations in vector form. It is especially convenient for when we have to deal with non-autonomous ODE's.
	The ODE model discussed in this paper is autonomous.}  to explicitly include time $t$ in the linear approximation. When
we explicitly include time $t$ (\ref{ODE4}) becomes: 
\begin{equation} \label{ODE5}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{\text{at} \ t + \Delta t}
\approx
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{\text{at} \ t}
+
\dfrac{d}{dt}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{\text{at} \ t}
\cdot \Delta t
\end{equation}
We rewrite (\ref{ODE5}) using (\ref{ODE3}):
\begin{equation} \label{ODE6}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{\text{at} \ t + \Delta t}
\approx
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{\text{at} \ t}
+
\begin{pmatrix}
v_x  \\
v_y \\
0 \\
-g \\
1
\end{pmatrix}_{\text{at} \ t}
\cdot \Delta t
\end{equation}
Note, the $``1"$ appearing in the bottom of the third vector in (\ref{ODE6}) is because $\dfrac{dt}{dt} = 1$. 

\noindent Recall Euler's method produces a sequence of Euler points (or vectors)
$$
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{0},\
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{1},\
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{2},\
\ldots
$$
which approximate the state of the system (and the solution) at times $t = 0, \ \Delta t, \ 2 \Delta t \ldots$. 

\noindent 
So,
for example, 
$\begin{pmatrix}
	x  \\
	y \\
	v_x \\
	v_y \\
	t
\end{pmatrix}_{2}
$
is the Euler approximation of the solution at $t = 2 \Delta t$.

\noindent \textbf{Euler's Method Example}. Suppose the ball is thrown at time $t = 0$ seconds   from 
$\langle x,y \rangle = \langle 5, 10 \rangle$ meters with
initial velocity $\langle v_x, v_y \rangle = \langle 1, 2 \rangle$ meters/second.
Using a step size of $\Delta t = 0.1$  seconds, approximate the location of
the ball at time $t = 0.2$ seconds using the Euler method. Ignore air resistance (drag).
\\
\textbf{Answer.} The $0^{th}$ Euler point is the initial conditions.
\begin{equation} \label{Euler0Example1}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{0}
=
\begin{pmatrix}
5  \\
10 \\
1 \\
2 \\
0
\end{pmatrix}
\end{equation}
The $1^{st}$ Euler point approximates the state of the system after 0.1 seconds. 
By (\ref{ODE6}) and (\ref{Euler0Example1}) it is given by:
\begin{align} \label{Euler1Example1}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{1}
&=
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{0}
+
\begin{pmatrix}
v_x  \\
v_y \\
0 \\
-g \\
1
\end{pmatrix}_{0}
\Delta t \nonumber \\
&=
\begin{pmatrix}
5  \\
10 \\
1 \\
2 \\
0
\end{pmatrix}
+
\begin{pmatrix}
1  \\
2 \\
0 \\
-9.8 \\
1
\end{pmatrix}
(.1)  \nonumber \\
&=
\begin{pmatrix}
5.1  \\
10.2 \\
1.0 \\
1.02 \\
.1
\end{pmatrix}
\end{align}
The $2^{nd}$ Euler point approximates the state of the system after 0.2 seconds. 
By (\ref{ODE6}) and (\ref{Euler1Example1}) it is given by:
\begin{align} \label{Euler2Example1}
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{2}
&=
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{1}
+
\begin{pmatrix}
v_x  \\
v_y \\
0 \\
-g \\
1
\end{pmatrix}_{1}
\Delta t \nonumber \\
&=
\begin{pmatrix}
5.1  \\
10.2 \\
1.0 \\
1.02 \\
.1
\end{pmatrix}
+
\begin{pmatrix}
1.0  \\
1.02 \\
0 \\
-9.8 \\
1
\end{pmatrix}
(.1)  \nonumber \\
&=
\begin{pmatrix}
5.2  \\
10.302 \\
1.0 \\
0.04 \\
.2
\end{pmatrix}
\end{align}
So, at time $t = 0.2$ seconds the ball's position is $\langle x(0.2),\ y(0.2) \rangle \approx  \langle 5.2,\ 10.302 \rangle$ meters and the ball's velocity is $\langle v_x(0.2),\ v_y(0.2) \rangle \approx  \langle 1.0,\ 0.04 \rangle$ meters/second.

\textbf{Including Air Resistance.} 
The standard formula %\footnote{See \url{http://www.physics.usyd.edu.au/~cross/TRAJECTORIES/Trajectories.html}}	
\cite{CrossJune2011Url, brody2002physics}
for the magnitude of the drag force $F_D$ on a ball is given by
\begin{equation} \label{DragForceScalar}
F_D = \dfrac{C_D\  d\  A\ ||\mathbf{v}||^2}{2}
\end{equation}
where $C_D$ is the drag coefficient, for a tennis ball $C_D = 0.5$; $d$ is the density of air, which 
at 20 degrees Centigrade ($68^\circ$ F) is 1.21 $\text{kg/meter}^3$; $A$ is the cross sectional area of the ball,
which for a ball of radius $r$ is  $\pi r^2$, a tennis ball's radius is  $r = 0.0265$ meters;
the velocity squared $||\mathbf{v}||^2 = v_{x}^{2} + v_{y}^{2}$. 
The mass of a tennis ball is  0.058 kg. 
However, it is important to note that drag force does not depend on the mass of the ball.

The drag force vector $\mathbf{F_D}$ points in the opposite direction of the ball's motion. So:
\begin{equation} \label{DragForceVector}
\mathbf{F_D} = - F_D \dfrac{\mathbf{v}}{||\mathbf{v}||} = -\dfrac{C_D\  d\  A\ ||\mathbf{v}||^2}{2}\ \dfrac{\mathbf{v}}{||\mathbf{v}||} = -\dfrac{C_D\  d\  \pi r^2\ ||\mathbf{v}||}{2}\ \mathbf{v} = 
-c\, ||\mathbf{v}||\  \mathbf{v}  
\end{equation}
is the drag force vector on the ball due to air resistance. Note, we let $c = \dfrac{C_D\  d\  \pi r^2}{2}$.

Forces are vectors so they add together. So, to include the drag force in your model
you just 
Equation (\ref{ODE0}):
\begin{equation*} 
\mathbf{F_{g}} = m\textbf{a}
\end{equation*}
to
\begin{equation} \label{GFforce}
\mathbf{F_{g}} + \mathbf{F_D} = m\textbf{a}
\end{equation}

At this point you have everything you need to solve the problem and find
the launch angle $\theta$ which will maximize how far the ball will land from
you.  

\medskip

\textbf{Python Script} 

Accompanying this modeling scenario is a Python script called 
\url{ThrownBallEulerStudent.py}. It solves the above problem, but assumes
there is no air resistance. You can ``easily'' modify this script
to include drag. 

\medskip

\textbf{Details.} The accompanying Python script, \texttt{ThrownBallEulerStudent.py}, computes the tennis ball's trajectories from launch to hitting the ground, for launch angles $\theta = 0, 1, 2, \ldots, 90$ degrees using  Euler's method. It stores the $x$ coordinate for
where each trajectory hits the ground, along with the corresponding  launch angle 
$\theta$, as well as outputting this to the console. It then applies Euler's method to that launch angle $\theta$ which gave the ``best (maximal)'' result to
recalculate the optimal (best) trajectory.  It outputs and saves a plot of
the best trajectory, see below, and prints to the console  a paragraph that gives some information about the scenario, including the optimal 
launch angle and maximum $x$ distance.  

\medskip

\textbf{From Output of \texttt{ThrownBallEulerStudent.py}:}

\begin{center}
	\includegraphics[width=.5\linewidth]{NoDragGraph.png}
\end{center}
\texttt{Assuming drag = 0. If when you throw a tennis ball you release it at a height of 2  meters and a speed of 12 meters/second, and you want it to land furthest from you, you should throw the ball at an angle of 42  degrees: it will land about 16.6 meters away.}

\medskip

\textbf{Notes about Python}  
\begin{enumerate}
	\item You don't have to download Python to run or edit Python scripts. There are free, online Python compilers. A good one is at	\url{https://repl.it/languages/python3}.  You just copy and paste your script
	into the online compiler window and hit RUN.
	\item The fastest and best way to run Python is to install it on you computer.
	A free and widely used Python IDE, especially in the sciences and engineering, is called Anaconda \cite{AnacondaUrl}. 
	\item \texttt{ThrownBallEulerStudent.py} is written for Python 3. There are two versions of Python: Python 2 and the newer Python 3. They are almost exactly the same. Unfortunately, they are not 100\% the same. So, sometimes a script written for one version of Python can't be compiled by the other version of Python.
	\item To do $5^2$ in Python use $5**2$.
	\item To do $\sqrt{7}$ in Python  use np.sqrt(7).
\end{enumerate}

\textbf{Possible Scenario Extensions} 
\begin{enumerate}
	\item Wind.
	\item After the ball hits the ground, model how it bounces.
	\item Instead of level ground, suppose the ground is an inclined plane \cite{winkel2015}.
\end{enumerate}

\medskip

\teacher{
\textbf{Solution for Teacher Version} 

\begin{align*} 
\mathbf{F_{g}} + \mathbf{F_D} &= m\textbf{a} \\
\langle 0, -mg \rangle -c\, ||\mathbf{v}||\  \langle v_x, v_y \rangle   &= m \langle \ddot{x}, \ddot{y} \rangle 
\end{align*}
Write the $\ddot{x}$ and the $\ddot{y}$ parts on different lines. 
\begin{align*} 
-c||\mathbf{v}||\ v_x  &= m \ddot{x} \\
-mg  - c||\mathbf{v}|| \ v_y   &= m \ddot{y}  
\end{align*}
Divide by $m$.
\begin{align*} 
-\dfrac{c}{m}\, ||\mathbf{v}||\ v_x  &= \ddot{x} \\
-g  - \dfrac{c}{m}\,||\mathbf{v}|| \ v_y   &= \ddot{y}  
\end{align*}
Convert to first order.
\begin{align*} 
\dot{x} &= v_x\\
\dot{y} &= v_y \\
\dot{v}_x &=  -\dfrac{c}{m}\, \sqrt{v_{x}^{2} + v_{y}^{2}}\ v_x   \\
\dot{v}_y &= -g  -\dfrac{c}{m}\, \sqrt{v_{x}^{2} + v_{y}^{2}}\ v_y  
\end{align*}
The initial position, at $t = 0$ seconds, is $\langle 0, 2 \rangle$ meters and we throw the ball at launch angle
of $\theta$, so the initial velocity is $12 \langle \cos \theta, \sin \theta \rangle$ meters/second.  We have the Euler iterations stop when the ball hits the ground, i.e. when $y$  goes from positive to negative (or zero). 
For each value of $\theta$ we record how far in the $x$ direction the ball travels before hitting the ground. 
We can start by throwing the ball in the horizontal direction, i.e. $\theta = 0 $ degrees, and then increase in one degree increments until $\theta =  90$ degrees. Here is the Euler iteration scheme for $\theta$.

$$
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{0}
=
\begin{pmatrix}
0  \\
2 \\
12 \cos \theta \\
12 \sin \theta \\
0
\end{pmatrix}
$$

$$
\begin{pmatrix}
	x  \\
	y \\
	v_x \\
	v_y \\
	t
\end{pmatrix}_{n+1}
=
\begin{pmatrix}
x  \\
y \\
v_x \\
v_y \\
t
\end{pmatrix}_{n}
+
\begin{pmatrix}
v_x  \\
v_y \\
-\dfrac{c}{m}\, \sqrt{v_{x}^{2} + v_{y}^{2}}\ v_x \\
-g  -\dfrac{c}{m}\, \sqrt{v_{x}^{2} + v_{y}^{2}}\ v_y \\
1
\end{pmatrix}_{n}
\cdot \Delta t
$$


The accompanying Python script, \texttt{ThrownBallEuler.py}, computes
the tennis ball's trajectories from launch to hitting the ground, for $\theta = 0, 1, 2, \ldots, 90$ degrees using the Euler method. It stores the $x$ coordinate for
where each trajectory hits the ground, along with the corresponding launch angle 
$\theta$, as well as outputting this to the console. It then applies the Euler method to that launch angle $\theta$ which gave the ``best (maximal)'' result to
recalculate the optimal (best) trajectory.  It outputs and saves a plot of
the best trajectory, see below, and prints to the console  a paragraph that gives some information about the scenario, including the optimal 
launch angle and maximum $x$ distance.  

\medskip

\textbf{From Output of \texttt{ThrownBallEuler.py}:}

\begin{center}
	\includegraphics[width=.5\linewidth]{graph.png}
\end{center}
\texttt{If when you throw a tennis ball you release it at a height of 2  meters and a speed of 12 meters/second, and you want it to land furthest from you, you should throw the ball at an angle of 40  degrees: it will land about 14.7 meters away.}

} % End teacher part



\begin{thebibliography}{98}


\bibitem{AnacondaUrl} Anaconda. 2019. {\it Anaconda Distribution\/}. \href{https://www.anaconda.com/distribution/}{https://www.anaconda.com/distribution/}.  Accessed 18 August 2019.

\bibitem{brody2002physics} Brody, H. and R.~Cross, and C.~Lindsey. 2002. {\it The Physics and Technology of Tennis\/}. Suwanee GA: Racquet Tech Publishing.

\bibitem{CrossJune2011Url} Cross, R. 2011. {\it Ball Trajectories\/}. \href{http://www.physics.usyd.edu.au/~cross/TRAJECTORIES/Trajectories.html}{http://www.physics.usyd.edu.au/~cross/TRAJECTORIES/Trajectories.html}. Accessed 18 August 2019.


 





\bibitem{winkel2015} Winkel, B. 2015. 3-045-T-BallBounce.


\end{thebibliography}




\end{document}
