Physics
The Art of Free Fall

Things fall. Most everyone is familiar with the legend of the apple that inspired Newton to formulate the classical laws of physics. Modern scientists, in physics and in other fields, get their inspiration not only from direct physical observation, as in the case of Newton and his apple, but also from the observations of 'virtual' systems, where computers are used to model the behaviour of physical systems. In this lesson we will use Mathematica to build computational models that will represent the falling of objects due to gravity, both with and without the influence of air resistance.

We can begin by setting some initial constant values, representing the initial state of some object, (in honor of Newton let's call it an apple), which will be initially under gravitational acceleration (Grav=-9.81, measured here in metric units: meters per second per second), have no initial velocity, (Vel0=0.0), and have an initial position of 6 meters 'above' some relative point (Pos0=6.0).  All of these will be associated with an initial time,Tim0, set at zero. We ask Mathematica to confirm our output as a list of variables.

Grav =  -9.81;
Vel0 =  0.0;
Pos0 =  6.0;
Tim0 =  0.0;
{Tim0, Pos0, Vel0, Grav}

While the starting conditions (represented with the appended 'naught') will remain constant, as reminders of our starting point, we also establish state variables for time, position,  velocity and acceleration that we will expect to change as we simulate the falling of the apple. These states all start out equal to the initial states.

Acc  =  Grav;
Vel  =  Vel0;
Pos  =  Pos0;
Tim  =  Tim0;
{Tim, Pos, Vel, Acc}

In order to observe an object moving, time must pass by.  To represent the passage of time in our model, we choose a length of time, known as a time step, and then use the information from where the object was at the beginning of the time step to figure out where it will be at the end of the step. 
	In science, change is often represented by the greek letter delta, (Δ), so we will refer to our small change in time as 'delta time' through the variable Deltime, measure in seconds.

Deltime  =  0.05;

Now it's time for a change!  We anticipate adding some Drag on the apple due to air resistance, but for now we will neglect the drag value by setting it to zero. First we adjust the acceleration to balance the effects of gravity and air resistance. Acceleration is the change in velocity over time, so we multiply by our change in time to determine the velocity change and add this value to the old Velocity to get a new value. Velocity has a similar relationship and can be similarly added to the old Position value. We also add a change in time to adjust the current state in Time.

Deltime  =  0.05;
Drag  =  0.0;
Acc  =  Grav + Drag;
Vel  =  Vel + Acc*Deltime;
Pos   =  Pos + Vel*Deltime;
Tim  =  Tim + Deltime;
{Tim, Pos, Vel, Acc}

In many programming languages there is a shorthand for adding a value to a running total, a process also known as incrementing. Assigning a value using += tells Mathematica to increment the value on the left with the value on the right. Using this notation, we can repeat the process above to see the state of our apple after two time steps, or at .1 seconds after it started falling.

Drag  =  0.0;
Acc  =  Grav + Drag;
Vel  +=  Acc*Deltime;
Pos  +=  Vel*Deltime;
Tim  +=  Deltime;
{Tim, Pos, Vel, Acc}

If want to see how the apple continues to fall, we should repeat this process a number of times. We let Mathematica do this work for us by using a loop function. One common loop function takes the form Do[expr,{i,imin,imax}]. The loop uses a counter, in this case represented by the variable i, to count from imin  to imax, each time evaluating expr.  To repeat our time step procedure a set number of times we simply place the command inside a Do loop. We have already done the process twice, so eighteen more cycles should show where the apple is after one second.

Do[
(Drag  =  0.0; 
Acc  =  
          Grav + Drag; 
Vel  +=  
          Acc*Deltime; 
Pos  +=  
          Vel*Deltime; 
Tim  +=  Deltime;), {i, 1, 
        18}]
{Tim, Pos, Vel, Acc}

So, after one second, the apple is about .85 meters above our measuring level (probably the ground). Want if we wanted to know how far the apple had fallen in half a second, or at another earlier time?  We would need to repeat the calculation, because even though our procedure may have calculated the value, it has long since overwritten it with new values. One way to keep track of that information is to use the Table command, which acts like the Do command as far as repeating an expression, but actually stores the results in a list, which in this case we will name as StateList. Notice we reset the state values before generating this list (Why?).

Vel =  Vel0;
Pos  =  Pos0;
Tim  =  Tim0;
Deltime  =  0.05;
Steps  =  20;
StateList  = 
      Table[(
Drag  =  
            0.0; 
Acc  =  
            Grav + Drag; 
Vel  +=  
            Acc*Deltime; 
Pos  +=  
            Vel*Deltime; 
Tim  +=  
            Deltime; 
{Tim, Pos, Vel, Acc}), {i, 1, 
          Steps}]

Having generated the list, we might notice that our list does not include our initial state, which might be useful.  We use one of Mathematica's list capabilities to Prepend the initial state to our StateList.  Thinking ahead, we also see where it might be easier to set the end time of a simulation as opposed to setting the number of steps, which can be easily calculated from the start and endtimes and the Deltime time change.

Vel =  Vel0;
Pos  =  Pos0;
Tim  =  Tim0;
Deltime  =  0.05;
Endtime  =  1.0;
Steps  =  (Endtime -  Tim0)/
          Deltime;
StateList  = 
        Table[(
Drag  =  
              0.0; 
Acc  =  
              Grav + Drag; 
Vel  +=  
              Acc*Deltime; 
Pos  +=  
              Vel*Deltime; 
Tim  +=  
              Deltime; 
{Tim, Pos, Vel, Acc}), {i, 1, 
            Steps}];
StateList  =  
      Prepend[StateList, {Tim0, Pos0, Vel0, Grav}]

So now we have a list of states corresponding to each time.  We can look at individual parts of the list either using the Part command, or the shorthand double-bracket form, both of which require an index, which is a number representing a location within the list. If we wish to access our initial state, for example, we would index the first position in the list.

Part[StateList, 1]

StateList[[1]]

Because each part of the list is also a list, we can use two indices to access particular parts. Noting that position is the second element in our list of states, we can determine the position at the beginning of the 5th time step via

StateList[[5, 2]]

Requesting the positions for All time steps is relatively straightforward:

StateList[[All, 2]]

ListPlot provides the means for us to view the data graphically.  If we supply a single list of numbers, Mathematica plots them on the y-axis, using their position in the list as the x-axis value. The option AxesLabel provides the means for labelling each axis.

ListPlot[StateList[[All, 2]], 
        AxesLabel -> {"Time Step", "Position"}];
ListPlot[StateList[[All, 3]], 
        AxesLabel -> {"Time Step", "Velocity"}];
ListPlot[StateList[[All, 4]], 
        AxesLabel -> {"Time Step", "Acceleration"}];

Whenever we find ourselves repeating a calculation or a process, it often makes sense to define a reusable program or function in Mathematica or another language to make things easier, particularly if the process has well defined input and output. Here we take the sequence of initial inputs and use them as parameters for a function we name Trajectory, which will follow the same process as before to determine and output a list of states that correspond to whatever values are given for the parameters. The := assignment is used to tell Mathematica "don't evaluate this now, evaluate it when it is called for later".

Trajectory[Grav_, Vel0_, Pos0_, Tim0_, Endtime_, 
          Deltime_]  := 
(Vel  =  
            Vel0; 
Pos  =  
            Pos0; 
Tim  =  
            Tim0; 
Steps  =  (Endtime - Tim0)/
              Deltime; 
Traj  = 
            Table[
(Drag  =  
                  0.0; 
Acc  =  
                  Grav  +  Drag; 
Vel  +=  
                  Acc*Deltime; 
Pos  +=  
                  Vel*Deltime; 
Tim  +=  
                  Deltime; 
{Tim, Pos, Vel, Acc}), {i, 1, 
                Steps}]; 
 
          Prepend[Traj, {Tim0, Pos0, Vel0, Grav}]);

Having defined a function, we can now generate a list of states by supplying parameters, either in number form, or using defined variables.  (Can you tell which numbers represent what in the parameter list?)

GravEng  =  -32.2;
StateList  =  
        Trajectory[GravEng, 35.0, 0.0, 0.0, 2.0, 0.1];

As before, we can view our list graphically using ListPlot. If instead of a single list, we supply a list of two elements, ListPlot will plot the first on the x-axis and the second on the y-axis. Note how we use a list of indices to select particular parts of StateList.  We can store an entire graph result for later use by assigning a variable name, as in the case with PPlot1, below.

PPlot1 = 
ListPlot[StateList[[All, {1, 2}]], 
          AxesLabel -> {"Time", "Position"}];
ListPlot[StateList[[All, {1, 3}]], 
        AxesLabel -> {"Time", "Velocity"}];
ListPlot[StateList[[All, {1, 4}]], 
        AxesLabel -> {"Time", "Acceleration"}];

To introduce the idea of terminal velocity, we will need to add air resistance to our model. A simplified air-resistance law usually presumes that drag is proportional to some power of the velocity (usually (v^2)or (v^3)), and acting in opposition to this velocity. In order to get our signs correct, we will need and absolute value function to isolate the magnitude of the velocity (speed).

Abs[-3.0]

...and we will use the Sign function to determine the direction of the velocity.

Sign[-3.0]

We define a new trajectory function by modifying the previous framework - adding two parameters representing a drag coefficient and an exponent (for the velocity), and setting up the calculation of the acceleration due to the drag of air resistance. Notice that our new function, TrajwDrag, actually makes the old function obsolete, as in order to get a no-drag trajectory, we can simply use a drag coefficient of zero!

\(\(TrajwDrag[Grav_, Vel0_, Pos0_, 
		Tim0_, 
          Endtime_, Deltime_, 
		DragCo_, 
          DragExp_]\  := 
\((Vel\  = \ 
            Vel0; 
Pos\  = \ 
            Pos0; 
Tim\  = \ 
            Tim0; 
Steps\  = \ \((Endtime - Tim0)\)/
              Deltime; 
Traj\  = 
            Table[
\((Drag\  = \ \(-\ Sign[Vel]\)*DragCo*
                    Abs[Vel]^DragExp; 
Acc\  = \ 
                  Grav + Drag; 
Vel\  += \ 
                  Acc*Deltime; 
Pos\  += \ 
                  Vel*Deltime; 
Tim\  += \ 
                  Deltime; 
{Tim, Pos, Vel, 
                  Acc})\), {Steps}]; 
          Prepend[Traj, {Tim0, Pos0, Vel0, Grav}])\);\)\)

We create and plot a new simulation, using the plotting option PlotRange to ensure that we can see the entire range of data, (from  Minimum and Maximum values) as well as the x-axis at zero.

Grav  =  -9.81;
StateList  =  
        TrajwDrag[Grav, 0.0, 200.0, 0.0, 10.0, 0.1, 0.02, 
          2.0];
PosList  =  StateList[[All, 2]];
VelList  =  StateList[[All, 3]];
AccList  =  StateList[[All, 4]];
PPlot2 = 
ListPlot[StateList[[All, {1, 2}]], 
          AxesLabel -> {"Time", "Position"}, 
          PlotRange -> {0, Max[PosList]}];
ListPlot[StateList[[All, {1, 3}]], 
        AxesLabel -> {"Time", "Velocity"}, 
        PlotRange -> {Min[VelList], 0}];
ListPlot[StateList[[All, {1, 4}]], 
        AxesLabel -> {"Time", "Acceleration"}, 
        PlotRange -> {Min[AccList], 0}];

Look at the graphs above.  How do they differ from the drag-free simulations? Which graph gives the best indication of when the object reaches 'terminal velocity'? What is that terminal velocity?
	Newton's equations are well known, and solvable with calculus (which was actually developed to solve them!). Mathematica has the ability to solve equations, both algebraic (using Solve) and differential (using DSolve). Both cases require a list of equations (established using ==). For differential equations, the apostrophe ' is used to designate derivatives, and the equations must include initial and boundary conditions. DSolve takes a list of equations, a dependent variable to solve for, and the independent variable (that the derivatives rely upon), and returns a substitution rule with the solution of the equation - you should recognize the form (1/2) a (t^2)  +  (Subscript[v, 0]) t  +  (Subscript[x, 0]).

Eqns  = {pos''[t]  ==  GravEng, pos'[0] ==  35.0, 
          pos[0]  ==  0.0};
Subrule = DSolve[Eqns, pos[t], t]

We use the ReplaceAll function, or the shorthand counterpart, /. , to create a function that solves the equations

Soltn[t_]  =  ReplaceAll[pos[t], Subrule]
Soltn[t_]  =  pos[t] /.  Subrule
Soltn[2.0]

Plotting a function differs from plotting a list. For this case we use Plot, which requires a function and a list specifying a domain variable and a range for that variable. We save the plot for the analytical solution by assigning a variable name.

PPlot1an  = Plot[Soltn[t], {t, 0, 2}];

We can then compare our numerical simulation with the mathematical solution by plotting both results on a numerical graph.

Show[PPlot1, PPlot1an];;

Notice that our numerical simulation is only an estimate of the true solution.  How can we make the estimate better? The only parameter we could change would be our time step, Deltime. Let's compare plots using a series of three different time steps in the Trajectory simulation.

Test1  =  
        Trajectory[GravEng, 35.0, 0.0, 0.0, 2.0, 
          0.3];
Test2  =  
        Trajectory[GravEng, 35.0, 0.0, 0.0, 2.0, 
          0.1];
Test3  =  
        Trajectory[GravEng, 35.0, 0.0, 0.0, 2.0, 
          0.02];
PlTest1  =  
        ListPlot[Test1[[All, {1, 2}]]];
PlTest2  =  
        ListPlot[Test2[[All, {1, 2}]]];
PlTest3  =  
        ListPlot[Test3[[All, {1, 2}]]];
Show[PlTest1, PlTest2, PlTest3, PPlot1an];

Can we do the same sort of analysis including the presence of air resistance? Let's set up the equations including drag:

EqnswDrag  = {pos''[t]  ==  
            Grav - Sign[pos'[t]]*0.02*Abs[pos'[t]]^2, pos'[
              0] ==  0.0, pos[0]  ==  200.0};

Unfortunately, the presence of the sign and absolute value functions make this set of differential equations very difficult (perhaps not yet possible) to solve by normal analysis. Indeed, if you try to use DSolve, Mathematica will return an error message, after a considerable amount of time (> 15 minutes on a standard PC) 'thinking'. However, as we have seen, if we take very small time steps, we can make a very good estimate of a real system.  Mathematica uses a process similar to the one we have built to 'solve' differential equations numerically, using a function called NDSolve.

Subrule = 
        NDSolve[EqnswDrag, pos[t], {t, 0, 10}];
Soltn[t_]  =  ReplaceAll[pos[t], Subrule];
Soltn[t_]  =  pos[t] /.  Subrule;
PPlot2an  = Plot[Soltn[t], {t, 0, 10}];

Even though this solution isn't mathematically 'perfect', we can still see how our TragwDrag function gets closer and closer, or converges to the solution as we make the time step smaller.  This is actually one of the ways mathemeticians and scientists can gain insight on 'unsolvable' math problems - by building a numerical method to solve the problem and observing the functions for increasingly accurate simulations.

Test4  =  
        TrajwDrag[Grav, 0.0, 200.0, 0.0, 10.0, 1.0,  .02, 
          2.0];
Test5  =  
        TrajwDrag[Grav, 0.0, 200.0, 0.0, 10.0, 0.3,  .02, 
          2.0];
Test6  =  
        TrajwDrag[Grav, 0.0, 200.0, 0.0, 10.0, 0.1,  .02, 
          2.0];
PlTest4  =  
        ListPlot[Test4[[All, {1, 2}]]];
PlTest5  =  
        ListPlot[Test5[[All, {1, 2}]]];
PlTest6  =  
        ListPlot[Test6[[All, {1, 2}]]];
Show[PlTest4, PlTest5, PlTest6, PPlot2an];

This notebook created by Dr. Garrett R. Love for the Shodor Education Foundation. Copyright 2001