\input{labpreamble}
\begin{document}

\subsubsection*{PHSX 221L\hfill Computer Simulation of a Falling Object\hfill Name:\underline{\hspace*{2in}}}

% 
% \underline{Title:} Computer Simulation of a Falling Object

\underline{Object:} (a) To run a pre-written program which uses 
numerical methods to compute the position and speed of a 
freely-falling object dropped from rest at any given time.  (b) To 
run the same simulation with a spreadsheet application.  (c) To 
modify this simulation to use average values of velocity, and to take 
into account air resistance.  (d) To compare each of the above 
solutions with the exact solution obtained from mathematical methods.

\underline{Theory:} Newton's 2nd law,
$\Sigma\boldsymbol{F}=m\boldsymbol{a}$, is the governing equation for
any accelerating object.  For the one-dimensional motion considered
here, we drop the vectors.

\underline{For free-fall} the only force is gravity ($F_{g}=mg$), so 
$mg=ma$ and acceleration is a constant $g = 9.80$~m/s$^2$.  We also 
assume the objects are dropped from rest.

From the definitions of velocity and acceleration we can write two 
coupled differential equations:
\begin{equation}\label{consta}
a= {dv\over dt} = g = 9.80
\end{equation}
\begin{equation}
v ={ dy\over dt}
\end{equation}

Converting these to finite difference equations by changing the 
differential operator $d$ to $\Delta t$ we have
\begin{equation}
v_{i+1}= v_{i} + a_{i}\Delta t
\end{equation}
\begin{equation}\label{simpleY}
y_{i+1}=y_{i} +v_{i}\Delta t
\end{equation}
where $a_{i}$ is always $g = 9.80$, and $v_{0} = y_{0} = 0$.  This is 
Euler's Method (see section 6.5 in Serway).

These equations produce computations which can be tabulated.  The 
$(i+1)$th values are computed from the $i$th values, \emph{i.e.,} each new 
row in the table is computed from the previous row.
% \vspace{-.2in}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}\hline
   $i$ &   $t$       &   $a$    &   $v$   &   $y$    \\ \hline
    0  &   0         &   9.80   &    0    &    0     \\
	1  &  $\Delta t$ &   9.80   & $v_{1}$ & $y_{1}$  \\
	2  & 2$\Delta t$ &   9.80   & $v_{2}$ & $y_{2}$  \\
\vdots &  \vdots     &  \vdots  & \vdots  & \vdots   \\
\end{tabular}
\end{center}

Notice that the missing $\frac{1}{2}g(\Delta t)^{2}$ term will cause
error in the last column.  

\pagebreak

Of course equation \ref{consta} can be integrated (twice) analytically to 
obtain the exact solution
\begin{equation}
v=gt
\end{equation}
\begin{equation}\label{analyticsimpleY}
y=\frac{1}{2}gt^{2}
\end{equation}
where both constants of integration are zero ({\em i.e.} at $t=0$, 
$v_{0}=y_{0}=0$).  Then the values in the generated table can be 
compared to the exact solutions.

The solution by numerical methods is in general only an approximation
to the exact solution.  Theoretically, the numerical solution should
converge to the exact solution in the limit as 
$\Delta t\rightarrow 0$, but there are also round-off and
truncation errors.  The basic Euler method doesn't obscure the physics 
as more sophisticated methods (such as the Runge-Kutta) might; but it 
is seldom used because the error is too large.  However, a simple 
modification to the Euler method produces significantly better 
results; it entails using the average velocity during a time interval 
instead of just the velocity at the beginning of the interval to 
compute the next value of $y$.

\underline{When an object moves through a fluid}, such as air, the
drag force $F_{D}$ is $\frac{1}{2}DA\rho v^{2}$ where $D$ is the drag
coefficient, $A$ is the cross sectional area of the object normal to
the its path, and $\rho$ is the density of the fluid.  For situations
of slow motion through a viscous fluid (low Reynolds number) $D$
depends on the other variables in a way that leaves the drag force
proportional to $v$.  Stoke's Law, $F_{D}=6\pi\eta r v$ (where $\eta$
is the viscosity), then applies.  For situations of faster motion
through less viscous fluids (moderately high Reynolds number) the drag
force is proportional to the square of the speed, $F_{D}\propto v^{2}$,
because $D$ is almost constant.  We examine the second case here, as
is generally done for macroscopic objects in air.  (See also section
6.4 in Serway.)

For falling through air, there are two forces on the left 
side of Newton's 2nd law.
\begin{equation}
\Sigma F= mg -kv^{2}
\end{equation}
where we use $k$ for $\frac{1}{2}DA\rho$.  Newton's law then gives
\begin{equation}
mg-kv^{2}=ma.
\end{equation}
Acceleration is decreasing ($a=g-(k/m)v^{2}$) and the speed 
approaches a constant value called terminal speed, $v_{t}$.

 
From the definitions of velocity and acceleration we can write two 
coupled differential equations:
\begin{equation}\label{vara}
a= {dv\over dt} = g -\frac{k}{m}v^{2}
\end{equation}
\begin{equation}
v ={ dy\over dt}
\end{equation}

\newpage

Converting these to finite difference equations by changing the 
differential operator $d$ to $\Delta t$ we have
\begin{equation}\label{aair}
a_{i+1}=g-\frac{k}{m}v_{i}^{2}
\end{equation}
\begin{equation}\label{vair}
v_{i+1}= v_{i} + a_{i}\Delta t
\end{equation}
\begin{equation}\label{yair}
y_{i+1}=y_{i} +v_{i}\Delta t
\end{equation}
where $a_{0} = g$, and $v_{0} = y_{0} = 0$.

% \pagebreak
These equations produce computations which can be tabulated.  In the 
following example $g=10.0$, $\Delta t = .10$, $m=5.0$, and $k=0.050$ in S.I. units.
\begin{table}[ht]
	\centering
\begin{tabular}{|c|c|c|c|c|}\hline
   $i$ &   $t$ (s) &$a$ (m/s$^2$)& $v$ (m/s)& $y$ (m)   \\ \hline
    0  &   0.00    &    10.000   &    0.000 &   0.000   \\
	1  &   0.10    &    10.000   &    1.000 &   0.000   \\
	2  &   0.20    &     9.990   &    2.000 &   0.100   \\
	3  &   0.30    &     9.960   &    2.999 &   0.300   \\
	4  &   0.40    &     9.910   &    3.995 &   0.600   \\
\vdots &  \vdots   &   \vdots    &  \vdots  &  \vdots   \\
\end{tabular}	
  	\caption{With Air Resistance.}\label{tbl:airresisttable}
\end{table}

% This is tabular (like the first one) instead of the table above.
% \begin{center}
% \begin{tabular}{|c|c|c|c|c|}\hline
%    $i$ &   $t$ (s) &$a$ (m/s$^2$)& $v$ (m/s)& $y$ (m)   \\ \hline
%     0  &   0.00    &    10.000   &    0.000 &   0.000   \\
% 	1  &   0.10    &    10.000   &    1.000 &   0.000   \\
% 	2  &   0.20    &     9.990   &    2.000 &   0.100   \\
% 	3  &   0.30    &     9.960   &    2.999 &   0.300   \\
% 	4  &   0.40    &     9.910   &    3.995 &   0.600   \\
% \vdots &  \vdots   &  \vdots    & \vdots  & \vdots   \\
% \end{tabular}
% \end{center}

\goodbreak

Though analytic integration of differential equations is often not possible, 
equation \ref{vara} can also be integrated (twice) analytically to obtain 
the exact solutions
\begin{equation}
v = v_{t}\tanh({t/\tau})
\end{equation}
\begin{equation}
y=v_{t}\,\tau \ln\cosh({t/\tau})
\end{equation}
where $\tau = \sqrt{m/kg}$ and is called the characteristic time, and 
$v_{t} = \sqrt{mg/k}$; again both constants of 
integration are zero.  The result can be compared with that obtained 
from the numerical procedure.

\goodbreak

\underline{Procedure:} 
\vspace{-.1in}
\begin{enumerate}

\item Type the Euler's Method free-fall program into a computer using
any language you like (choose one that is common in your major, or
consider BASIC if you haven't ever done any programming before).  Also
do it on a spreadsheet for comparison of these two tools.  Run it
using a total time of 10 seconds and 25 intervals (\emph{i.e.} $\Delta
t = 0.40$~s) and compute the distance fallen.  Then reduce $\Delta t$
to 0.10~s and find the distance fallen after 10 seconds using 100
intervals.  Compare each of these with the exact solution given by
equation \ref{analyticsimpleY}.  Which run shows the best agreement? 
Why?  How close are they to the exact answers?

\item Modify the program and the spreadsheet to use the Improved Euler
Method by averaging velocities over each interval.  That is, since you
have already computed $v_{i+1}$ before you compute $y_{i+1}$, use the
average of $v_{i+1}$ and $v_{i}$ instead of just $v_{i}$ in equation
\ref{simpleY}.  Now run the program again with 25 time intervals and
compare with the exact solution.  Is the agreement better?  Why?  How 
good is it?

\item Now add air resistance and assume the falling body has a mass of
4.0~kg and has a shape, size, and speed which cause the air
resistance to be proportional to the square of the velocity.  If the
coefficient of resistance $k$ is 0.035~kg/m, then $F_{D} = 0.035
v^{2}$ if $F$ is in newtons and $v$ is in m/s.  Write a new version of
your program and your spreadsheet that will calculate the distance
fallen in 10~seconds using equations \ref{aair}--\ref{yair}.  (You can
verify that your program is working correctly by trying the sample
values above first and comparing with table \ref{tbl:airresisttable}. 
The results in the table are obtained using the unmodified Euler
Method.)  How does this result agree with an exact solution obtained
by integration?

\item Plot $a$, $v$, and $y$ vs.~$t$ on sets of vertically stacked axes 
(with the horizontal $t$ axes parallel) for both the 
free-fall case and (on another set of axes) the case with air resistance.

\end{enumerate}


\underline{Conclusions:} 
\vspace{-.1in}
\begin{enumerate}

\item How did you taylor this lab to fit your major?

\item What did you learn from this lab?

\item Compare the appropriateness of the two tools: a program in a 
programming language and a spreadsheet application.

\end{enumerate}


\newpage 

\vspace{-.25in}

\begin{verbatim}
10 rem this BASIC program performs numerical integration to calculate the
20 rem speed and distance of fall as a function of time for a
30 rem freely-falling body
40 print "Enter the total time of fall and the number of intervals"
50 input "Total time = ";tt
60 input "Number of intervals = ";n
70 dt = tt/n : rem dt is delta t
80 i = 0 : t = 0 : a = 9.8 : v = 0 : y = 0 : rem next initialize other vars
90 rem i is counter, t is time, a is accel, v is velocity, y is distance.
100 print "   i      t          a              v              y"
110 print "___________________________________________________________"
120 for i = 0 to n
130 print using "####";i;" ";
140 print using "####.###";t,a,v,y
150 t = t+dt
160 v = v+a*dt
170 y = y+v*dt
180 next i
\end{verbatim}

% \vspace{-.25in}
% 
% \begin{verbatim}
% /* FreeFall.c  Example */
% /* This C program performs numerical integration to calculate the speed*/
% /* and distance of fall as a function of time for a free-falling body */
% #include <stdio.h>
% #include <conio.h>     /* for clrscr(), getch(), etc.*/
% main()
% {   int N,i,;
%     float t,a,v,y,T,dt;
%     clrscr();
%     printf("\n\n This program will calculate the speed and distance ");
%     printf("\n of fall as a function of time for a freely-falling body.");
%     printf("\n\n Enter the total time of fall, and the number of intervals.");
%     printf("\n     Total time =                ");    scanf("%f", &T);
%     printf("\n     Number of intervals = ");    scanf("%d", &N);
%     dt = T/N;
%     i = t = y = v = 0;       a = 9.80;
%     printf("\n\n    i       t        a        v        y");
%     printf("\n\_______________________________");
%     printf("\n\n%4d %9.4f $12.6f $12.6f", i, t, a, v, y);
%     for (i = 1; i <=N; ++i)
%       {  t = t + dt;
%          v = v + a*dt;
%          y = y + v*dt;
%          printf("\n\n%4d %9.4f $12.6f $12.6f", i, t, a, v, y);
%          getch(); }
%     return 0; }
% \end{verbatim}

\newpage

\begin{verbatim}
/* FreeFall.java  Example */
/* This Java program performs numerical integration to calculate the speed*/
/* and distance of fall as a function of time for a free-falling body */
/* Input the total time and number of intervals on the command line. */
class FreeFall
{
   public static void main(String[] args)
   {
      int N = 25, i = 1;
      double t = 0.0, a = 9.8, v = 0.0, y = 0.0, T = 10.0, dt;
        
      if (args.length > 0)
      {
         T = Integer.parseInt( args[0] );
         N = Integer.parseInt( args[1] );
      }   
      dt = T/N;
        
      System.out.println("   i      t       a       v       y");
      System.out.println("______________________________________");
      for (i = 0; i <= N; ++i)
      {
         System.out.println( "   " + i + "    " + t + "    " 
                             + a + "    " + v + "    " + y);
         t = t + dt;
         y = y + v*dt;
         v = v + a*dt;
      }
   }
}
\end{verbatim}

\newpage

\begin{verbatim}
/*
Purpose:	Contains a C++ program that performs numerical integration to
			calculate the speed and distance of fall as a function of
			time for a free-falling body
File:		FreeFall.cpp
Author:		Garth Sorenson
Date:		28 Sep 2000
*/


#include <iostream>
using std::cin;
using std::cout;
using std::endl;
using std::ios;

#include <iomanip>
using std::setw;
using std::setprecision;
using std::setiosflags;


/*
Purpose:	Controls free-fall program
Input:		total time as floatin point, number of intervals as an integer
Output:		interval, time, acceleration, velocity, distance for each inteval
*/
int main()
{
	// local data
	int N,		// number of intervals
		i;		// loop control variable
	double T,	// total time
		   t,	// current time
		   a,	// acceleration
		   v,	// velocity
		   y,	// distance
		   dt;	// delta time

	// greeting
	cout << "This program will calculate the speed and distance\n"
		 << "of fall as a function of time for a freely-falling body.\n\n"
		 << "Enter the total time of fall, and the number of intervals.\n"
		 << endl;


	// get total time and number of intervals from user
	cout << "Total Time = ";
	cin >> T;
	cout << "Number of Intervals = ";
	cin >> N;
	cin.ignore(80, '\n');

	// initialize delta time, current time, 
	// distance, velocity, and acceleration
	dt = T / static_cast<double>(N);
	i = t = y = v = 0;
	a = 9.80;

	// display table header
	cout << setw(3) << "i" << setw(8) << "t" << setw(10) << "a"
         << setw(10) << "v" << setw(10) << "y" << endl
		 << "-------------------------------------------" << endl;

	// set formatting
	cout << setiosflags(ios::fixed | ios::showpoint)
		 << setprecision(2);

	// display interval values
	for (i = 0; i <= N; i++)
	{
		cout << setw(3) << i << setw(10) << t << setw(10) << a
			 << setw(10) << v << setw(10) << y << endl;
		t += dt;
		v += a * dt;
		y += v * dt;
		cin.get();
	}

	return 0;
}
\end{verbatim}

\end{document}
