% SIMIODETemplate-TeacherORStudent.tex and simiode.cls for LaTeX2e
% Version 29 July 2015.
%
% This class (simiode.cls) is a modification of that used by
% the journal PRIMUS - Problems, Resources, and Issues in
% Mathematics Undergraduate Studies.
%
% See http://www.tandfonline.com/loi/upri .
%
% The class was developed by  Marc Frantz with suggestions of the
% PRIMUS Editorial Board, especially John Thoo.
%
% 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.
%
\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.
%
% Basically the 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 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.
% This option was created  by Eric Sullivan, Carroll College, Helena MT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{hyperref}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will
% look in final form, so that you can better plan your
% presentation.  To see your manuscript formatted very
% similarly to the on-line published format of  SIMIODE,
% comment the following line out by inserting a percent sign % (%) at its beginning.

\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\\SHUTTLECOCK FALL}



\markboth{Modeling Buoyancy}{Modeling Buoyancy}




\author{Brian Winkel,
Director SIMIODE\\
Chardon OH 44024 USA\\
BrianWinkel@simiode.org
}


\makeStitlePDFLaTex



\begin{abstract}
We are given data on the time and position of a shuttlecock as it falls to the ground from a set height. We attempt to model the falling object and we confront the different resistance terms and models. Finally, we introduce a new way of comparing models, the Akaike Information Criterion, and apply it to our models. 
\end{abstract}




\section*{SCENARIO DESCRIPTION}

\begin{figure}[h]
\centering
\includegraphics[width=5.3cm]{3-19-shuttlecock.eps}

\vspace*{-2.2 true cm}
\caption{\label{fig:3-19-shuttlecock}  Shuttlecock, the flying object used in the game of badminton.}
\end{figure}

We are given data on time and distance traveled below the starting height of 2 m of a free falling shuttlecock (see Figure~\ref{fig:3-19-shuttlecock})   in Table~\ref{table:data}. We will make a mathematical model of this phenomenon and estimate the parameters in the model. We shall also compare the model's predictive power with the actual data. 





\begin{table}[h]
\begin{center}\footnotesize 
\begin{tabular}{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}\hline
Time (s) &\hspace*{-.20  true cm} 0 \hspace*{-.15 true cm} &\hspace*{-.20  true cm}0.347\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.47\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.519\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.582\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.65\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.674\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.717\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.766\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.823\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.87\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.031\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.193\hspace*{-.20  true cm} &\hspace*{-.20  true cm}  1.354\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.501\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.726\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.873\hspace*{-.20  true cm} \\ \hline
Distance (m) &\hspace*{-.32  true cm} 0\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 0.61\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.22\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.52\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 1.83\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 2.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm}2.13\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 2.44\hspace*{-.20  true cm} &\hspace*{-.20  true cm}  2.74\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 3.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 4.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 5.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 6.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 7.00\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 8.50\hspace*{-.20  true cm} &\hspace*{-.20  true cm} 9.50\hspace*{-.20  true cm} \\ \hline
\end{tabular}
\parbox{8.6 true cm}{\caption{\label{table:data} Time and distance traveled data on a free falling shuttlecock which is dropped from a height of 2 m at rest. }}
\end{center}
\end{table}

\begin{enumerate}
  \item Draw a Free Body Diagram with all the forces acting upon the shuttlecock as it falls.
  \item Use Newton's Second Law of Motion which states that the sum of the forces acting on a body is equal to the product of the body's mass and acceleration to produce an equation in which one side is $m\, y''(t)$, where $y(t)$ is the positive distance fallen (in m) at time $t$ (in sec) and $m$ is the mass of the shuttlecock. Comment on the reasonableness of your model.
  \item Knowing the fact that this experiment was done in Villanova PA USA find a way to estimate all parameters in your model.
  \item With these determined parameters compare your model with the data.
  \item Check with others or perhaps look up some on-line efforts with modeling falling bodies (with and without resistance).
  \item Compare your different models and comment on them.
 \end{enumerate}





\teacher{
\section*{NOTES FOR TEACHERS} 

\subsection*{Approach used in source\cite{PeastrelEtAl1980}} 

There are many articles, principally in the physics teaching literature, about modeling falling bodies, e.g., coffee filters, balloons, shuttlecocks, for they illustrate the concepts of resistance and terminal velocity. Often physicists will do things quite differently than one might do in a mathematics class. 


For example, in modeling the free fall of a shuttlecock the authors in our source article \cite{PeastrelEtAl1980} use the basic quadratic model, i.e. resistance to the falling shuttlecock is proportional to the square of its velocity and in opposition to the direction of the fall. The only other force considered is from the acceleration due to gravity. Hence the model, presupposed is 


\begin{equation}\label{eq:ModelPre}
m v'(t)= -K v(t)^2 +m g\, , \quad  v(0)=0\, ,
\end{equation}
which when we divide out by $m$ and identify $k = K/m$ becomes 

\begin{equation}\label{eq:ModelFine}
v'(t)= -k v(t)^2 + g \, , \quad  v(0)=0\, .
\end{equation}
 After determining the terminal velocity from setting $0 = m v'(t)= -K v(t)^2 +m g$, i.e. no change 
in velocity, and solving for $v_T$ which we substitute for $v(t)$ in the  steady state, we find that $v_T = \frac{\sqrt{g}}{\sqrt{k}}$. 
 


Upon substituting this value for $v_T$ into (\ref{eq:ModelFine}), using $g =9.8024\, m/s^2$ from the local area value obtained for Villanova PA (actually Philadelphia PA, which is nearby) which is found at the web site {\tt http://www.physics.montana.edu/demonstrations/video/1\_mechanics/demos/localgravity} \linebreak {\tt  chart.html}, and integrating to obtain a distance traveled function for the shuttlecock, $y(t)$, we find that the solution for the distance traveled, $y(t)$, is 
  

\begin{equation}\label{eq:SourceModelGen}
y(t) =  0.102016 (v_T)^2 \log\left(\cosh\left(\frac{9.8024 t}{v_T}\right)\right) \, . 
\end{equation}

From (\ref{eq:SourceModelGen}), using a least squares approach, i.e. minimizing the sum of the differences between the model and the data squared as a function of $v_T$, the authors find that $v_T =6.80\, m/s^2$ and they have their completed model: 

\begin{equation}\label{eq:SourceModel}
y(t)=4.72067\log(\cosh(1.441t)). 
\end{equation}

We can see from the plot in Figure~\ref{fig:3-19-SourceModelPlotWithData} of this model over the data that the model does an exceptional job of predicting the data. 


\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-SourceModelPlotWithData.eps}


\caption{\label{fig:3-19-SourceModelPlotWithData} Plot of the model from (\cite{PeastrelEtAl1980}) and the original data. }
\end{figure}



\subsection*{Other Approaches and Models }

We offer several additional models and complete analyses in a Mathematica notebook (along with the above source model approach), 3-19-Mma-T-ShuttleCockFall-TeacherVersion.nb, and a pdf version of that file for those who do not use Mathematica. The file is fully annotated so that the efforts can be replicated in another computer algebra system such as SAGE or Maple. 

\subsection*{Case of no resistance $m*v'(t)= m*g$} 
Here there are no parameters, just the fixed model $m*v'(t)= m*g$. We encourage students to pursue this model because it is easy to solve and they have some recollection about Galileo saying 
something like objects fall at the same rate (alas, they forgot the caveat - in a vacuum!) However, the model and the data do not agree as can be seen in Figure~\ref{fig:3-3-19-NoResistanceModelPlotWithData} and students push on to consider resistance in their Free Body Diagram if they have not already done so. 



\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-NoResistanceModelPlotWithData.eps}


\parbox{7 true cm}{\caption{\label{fig:3-3-19-NoResistanceModelPlotWithData}  Plot of the model with no resistance, $m*v'(t)= m*g$, and the original data.}}
\end{figure}

 
\subsection*{Case of $r =1$ in $v'(t)= -k*v(t)^r + g = -k*v(t)+g$} 

Here we solve the differential equation for $v(t)$, integrate it to get a distance traveled function, $y(t)$, using an initial condition of $y(0) = 0$, i.e. displacement from original drop position of 2 m high down to 0 m, and then proceed to do a minimization of sum of square errors on the parameter $a$ in the form of the model $v'(t)= -a*v(t)^r +g$. We obtain the model (\ref{eq:r1Resist}) for the distance traveled, $y(t)$, 

\begin{equation}\label{eq:r1Resist}
y(t)=9.23238t+8.69552e^{ -1.06174t} -8.69552.  
\end{equation}

\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-r1ResistanceModelPlotWithData.eps}


\parbox{8 true cm}{\caption{\label{fig:3-19-r1ResistanceModelPlotWithData}  Plot of the model with resistance using $r = 1$, $v'(t)= -k*v(t)^1 + g$, and the original data. }}
\end{figure}


 

This model looks reasonable (see Figure \ref{fig:3-19-r1ResistanceModelPlotWithData}) except the bows in the data and the model are conspicuously different. Indeed, we should expect the residuals (i.e. the differences between model and data to be random, not as regular as they are in Figure 5 and they are not. Therefore, this model is not a good one. 



\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-r1ResidualsResistanceModelPlot.eps}


\parbox{11.6 true cm}{\caption{\label{fig:3-19-r1ResidualsResistanceModelPlot}  Plot of the residuals (differences between model and data) from the model with resistance using $r = 1$, $v'(t)= -k*v(t)^1 +g$, and the original data. }}
\end{figure}



 
Thus we move on to the next case. 

\subsection*{Case of $r =2$ in $v'(t)= -k* v(t)^r + g = -k* v(t)^2 +g$} 

This is the case with which the authors in [3] worked. However, here we shall solve the differential equation and apply least squares technique to determine the parameter $k$ in the model, $v'(t)= -k* v(t)^2 +g$. When we do this we find that our model for the distance traveled, $y(t)$, is 

\begin{equation}\label{eq:Caser2} y(t)=4.72067\log(\cosh(1.441t))\, . 
\end{equation} 

We plot the model over the data in Figure~\ref{fig:3-19-r2ResistanceModelPlotWithData} and say this is spot on! 


\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-r2ResistanceModelPlotWithData.eps}


\parbox{8 true cm}{\caption{\label{fig:3-19-r2ResistanceModelPlotWithData}  Plot of the model with resistance using $r = 2$, $v'(t)= -k * v(t)^2 +g$, and the original data. }}
\end{figure}



Moreover, the residuals look random as well (see Figure~\ref{fig:3-19-r2ResidualsResistanceModelPlot}), i.e. no pattern, indicating errors are normally distributed about the model. 




\begin{figure}[h]
\centering
\includegraphics[width=7.3cm]{3-19-r2ResidualsResistanceModelPlot.eps}

\parbox{8 true cm}{\caption{\label{fig:3-19-r2ResidualsResistanceModelPlot}  Plot of the residuals(differences between model and data) from the model with resistance using $r = 2$, $v'(t)= -k * v(t)^2 +g$, and the original data.}}
\end{figure}


\subsection*{General Case of $a$ and $r$ in $v'(t)= -a * v(t)^r + g$} 

This is a tougher case, not only because it does not presume $r$, but because in that case the model cannot be solved analytically for a general $r$ value. Thus what we basically do is set up a grid over a reasonable region of $a$ and $r$, actually $0 \leq  a \leq 0.4$ and $1.5. \leq r \leq 2.5$ with a mesh of 0.001.  Then for each pair of parameters, $(a, r)$, we numerically solve the corresponding differential equation model, $v'(t)= -a * v(t)^r + g$, with $v(0)= 0$; integrate that solution to obtain a solution for $y(t)$, the distance traveled by the shuttlecock;  compute the corresponding sum of square errors between the model and the data; and store the data point, $(a, r, {\rm SSE}(a,r))$. We find the data point which has the smallest ${\rm SSE}(a,b)$ and use that as our optimal set of parameters for our final model in this case.    On our computer(an XPS Dell computer with Intel(R)Core(TM) i7-3770 CPU with 16GB of memory and running at 3.40 Ghz) it takes about 45 minutes to run through this grid. When done, the best fitting parameters are $a =0.205$ and $r =2.019$. These give a plot and sum of square errors which is practically identical to that rendered in our interpretation of the source article and in the case $r =2$. 

\subsection*{Comparing Models Using the Akaike Information Criterion(AIC)}
 
The Akaike Information Criterion (AIC) \cite{Akaike1974} is a way to compare models in terms of their accuracy  using sum of square errors and number of parameters. We know that if we have $m$ data points, ${\rm data} = \{(t_i,y_i)\vert i=1,2,3, \ldots , m\}$, then we can put a polynomial of degree $m -1$ exactly through the data with a sum of square errors of 0. However, this is costly with $m$ unknowns, the coefficients of the fitting polynomial of degree $m - 1$, and even though it gives a perfect fit ($SSE = 0$) it uses a large number of parameters. Moreover, this polynomial explains nothing about the underlying phenomenon we are modeling.  The AIC compensates for each of these notions - SSE and number of parameters. 

The Akaike Information Criterion, where $k$ is number of parameters to fit, $n$ is the number of observations, and RSS is the minimum residual sum of square errors (or just sum of square errors, SSE) is given in (\ref{eq:AIC}): 

\begin{equation}\label{eq:AIC}
AIC =2(1+k)+n \log\left(\frac{\rm RSS}{n}\right) \,.
\end{equation}
NB: When SSE is small then the log function is negative. 
Smaller (even negative) values of AIC represent a higher level of statistical support for the corresponding model. A difference of 2 in the AIC for two models is needed to suggest a significant difference between these two models.
 
\begin{table}
\centering
\begin{tabular}{|l|c|c|c|c|c|} \hline
Model & $r$ & $k = a$ & SSE & AIC   \\ \hline
Source & 2 && 0.00920356 & {\bf -123.811}  \\ \hline 
No Resistance && 0 & 128.369& 36.3689 \\ \hline  
Case &$r = 1$ &  1.06174 & 0.523775  &-55.1584  \\ \hline 
Case &$r = 2$ &  0.211834 & 0.00920356 & {\bf -123.863}  \\ \hline 
General $r$ & 2.019 & .0205 & 0.00919924 &  {\bf -121.871} \\ \hline  
\end{tabular}

\parbox{6.8 true cm}{\caption{\label{table:AICSummary} 
Information we need about our models to decide which model is best is presented.}}
\vspace*{-.7 true cm}
\end{table} 

 
Thus if AIC's are close to each other then the models are comparable. The smallest AIC means best model - in terms of balancing ACCURACY(RSS) and COMPLEXITY(k).
 
So after we complete all model analyses we will see if one (or several) of the models is best. We offer the parameters, sum of square errors, and AIC values for the best parameters in our models in Table~\ref{table:AICSummary}, where the boldface AIC values are close to each other and far below the other AIC's so these three models have about the same value. Indeed, they are all of the quadratic, $r =2$ or $r \approx 2$, model with respect to the power of velocity in the model for resistance. 

This all makes sense since the visuals on all three models are almost identical and much better than other models offered. 

\section*{CONCLUSION}
 

This is a very nice activity as it requires modeling from some basic physics principles (Newton's Second Law of Motion); it uses real data - albeit collected by someone else, but students could very easily collect this data with a motion sensor device and a shuttlecock; has an option for analytic solutions, either by hand or by computer algebra system; and offers several modeling paths which can then be compared. 

\begin{thebibliography}{99}
\bibitem{Akaike1974} Akaike, H. 1974. A new look at the statistical model identifcation. {\it IEEE Transactions on Automatic Control\/}. 19(6): 716-723.  

\bibitem{Ledder2013} Ledder, G. 2013. {\it Mathematics for the Life Sciences\/}. New York: Springer.

\bibitem{PeastrelEtAl1980} Peastrel, M., R.Lynch, and A. Armenti, Jr. 1980. Terminal velocity of a shuttlecock in vertical fall. {\it American Journal of Physics\/}. 48(7): 511-513. 


\end{thebibliography}
}

\end{document}
 
