Mathematical Models of a Falling Rock

1. Introduction and Background
1.1 Variables and Parameters
1.2 Notation
1.3 Governing Physical Principles

2. Mathematical Models 
2.1 A Discrete in Time Model
2.1 A Continuous in Time Model

3. Algorithm for a Discrete in Time Model
3.1 Pseudocode
3.2 C Code
3.3 Fortran Code

4. Computational Experiments and Model Refinements
4.1 Computational Experiments
4.2 Model Refinements

1. Introduction and Background

A mathematical model can be described as a set of equations or expressions which describe some physical phenomenon. Mathematical models are frequently used to describe the state of a physical system as the system evolves in time. As a concrete example of a physical system, consider a rock dropped with initial velocity V0 (meters/second) from a height of Z0 (meters). So that we can formulate a very simple mathematical model describing this physical system, we make the simplifying assumptions that (1) air resistance is zero and (2) the gravitational force is constant.

This document has two main goals: (1) present mathematical models of a simple physical phenomenon and (2) describe some mathematical modeling terminology.

1.1 Variables and Parameters

In our example the physical system is a rock dropped with initial velocity V0 from a height of Z0. In a mathematical model, the state variables are quantities that characterize the behavior of the physical system. In the case of a falling rock, the state of the physical system is completely characterized by the velocity of the rock and the vertical distance that the rock has traveled from its release point. The state of the system will evolve with time, i.e., the rock will fall. As a consequence the state variables in our model, the rock height and the rock velocity, are time dependent.

The input variables, output variables, and parameters for our mathematical model of the falling rock are

The distinction between input variables and parameters is somewhat arbitrary. One could argue that the gravitational constant should be an input variable. The point of view often adopted is that input variables characterize a single physical problem and parameters determine the context or setting of the physical problem. For example two falling rock problems with different input variables are: (1) drop the rock from 1000 meters; and (2) throw the rock downward with velocity 2 meters/second from a height of 500 meters. The reason for classifying the gravitational constant as a parameter is that it defines the context of the problem. That is, the gravitational constant does not change from one falling rock problem to another as long as we are near the earth's surface. If we were to move our falling rock experiment to the moon, then the value of the gravitational constant would change.

1.2 Notation

It is common practice to use shorthand notation for the variables and parameters. A convenient way of introducing the notation to be used in a mathematical model is the following tabular format.
 

Quantity

Symbol

Units

time t seconds
acceleration due to gravity g (g=9.8)  meters/second/second
velocity of the rock V or V(t) meters/second
rock height (distance above ground)  Z or Z(t) meters
initial velocity of the rock V0 or V(0)   meters/second
initial height of the rock Z0 or Z(0) meters

1.3 Governing Physical Principles

A mathematical model for a rock falling with constant acceleration can be developed from the definitions of velocity and acceleration.

Velocity: Velocity equals the rate of change of position.

Acceleration: Acceleration equals the rate of change of velocity.

These last two sentences are the physical principles from which we will develop two mathematical models for a falling rock.

2. Mathematical Models

In this section we will develop discrete in time and a continuous in time mathematical models for a falling rock.
 

2.1 A Discrete in Time Model

To develop a discrete in time mathematical model, we mark off a time line as indicated in Figure 1.

  <-----  h  ----->


--x-------|-------x-------|-------x-------|-------x-- (t-axis)

  t0      t1/2     t1      t3/2    t2    . . .
Figure 1. A time line with grid points marked by "x".

The points on the t-axis marked with the symbol " x " are called grid points or mesh points. They represent the discrete times at which the mathematical model will approximate the height of the rock. The points marked with the symbol " | " are the midpoints of the line segments joining adjacent time grid points. These midpoints will prove to be convenient time values for approximating rock velocities.

The grid we have defined is called a uniform grid or a regular grid since each pair of adjacent grid points are separated by a constant time interval. This time interval between grid points, h, is called the time step, the grid spacing or the mesh spacing.

Following the notation introduced in Figure 1, we can represent an arbitrary grid point on the t-axis by

where we have assumed that t = 0.  Note that ti+1- ti = h, which is just a restatement of the fact that successive grid points are separated by a time interval of h.

In the notation table of Section 2.3, we have chosen to denote the physical height and velocity of the rock by Z(t) and V(t), respectively. A common notational convention in mathematical modeling is to use variables of an opposite case to denote the mathematical model approximations of physical quantities. Following this convention, we define zi to be a numerical approximation to Z(ti) and vi to be a numerical approximation to V(ti).

To continue our development of a mathematical model of a falling rock, recall that if an object moves along a line at a constant rate R for a time T, then the object moves a distance D. In words, "distance traveled (D) equals rate (R) multiplied by time (T)", or

D = R*T

Put another way, if the initial position of the object is PI and the new position of the object T seconds later is PN, then

PN = PI + R*T    (1)

Similarly, if an object moves along a line at an initial rate of RI and accelerates at a constant rate A for a time interval T then its new velocity is given by

VN = VI + A*T    (2)

Equations (1) and (2) are special cases of the governing physical principles stated in Section 2.4.

Recall that V0 is the initial velocity of the rock. Let v1/2 denote the numerical model velocity of the rock at t1/2 = h/2, the first point marked " | ". Then, since g, the acceleration due to gravity, is assumed constant, we conclude from Equation (2) that:

v1/2 = V0 - g*(h/2)   (3a)

Similar reasoning for successive points marked " | " yields

v3/2 = v1/2 - g*h ,   v5/2 = v3/2 - g*h, ...    (3b)

Now we turn to the part of the mathematical model which deals with estimating the position of the rock as it falls. Recall that Z0 is the initial height of the rock. If the initial velocity is zero or negative (downward), then during the time interval t = 0 to t = h, the velocity of the rock increases. To use Equation (1) to develop a discrete mathematical model of the falling rock, we must decide on a constant fall velocity to assign to the rock for the period from t = 0 to t = h. If we choose V0 to represent the fall velocity over this period, we will underestimate the distance traveled. If we choose v1 = V0 - g*h to represent the fall velocity from t = 0 to t = h, we overestimate the distance traveled. Motivated by these last two observations, we formulate our mathematical model of the height of the falling rock using the velocity v1/2 = V0 -g*(h/2) to represent the fall velocity during the period from t = 0 to t = h. Using the constant v1/2 to represent the changing velocity over this period is a modeling assumption.

Now that we have assumed a constant fall rate during t = 0 to t = h, we can apply Equaton (1) to estimate the vertical distance from the ground to the rock. In Equation (1) let PI = Z0 (the initial height), R = v1/2, T = h and PN = z1 to conclude

z1 = Z0 + v1/2*g   (4a)

Similar reasoning for successive grid points marked " x " in Figure 1 yields

z2 = z1+ v3/2*g,  z3 = z2 + v5/2*g, ...    (4b)

Equations (3a,b) and (4a,b) are the core of our discrete in time mathematical model for a falling rock. We will postpone the discussion of how to implement this mathematical model until Section 3.

2.2 A Continuous in Time Mathematical Model

For our falling rock example, a continuous in time mathematical model can be formulated very simply in terms of derivatives as follows.

V'(t) = -g,   V(0) = V0

Z'(t) = V(t),   Z(0) = Z0

where the prime ' denotes differentiation with respect to time.

Those who have studied calculus will recognize that the exact solution to the above mathematical model is given by

V(t) = -g*t + V0    and    Z(t) = - 0.5*g*t2 + V0*t + Z0

Later we will compare the mathematical model predictions from the discrete in time method with the continuous in time method.

3. Algorithm for the Discrete Mathematical Model

In this section we will proceed with a computer implementation of the discrete mathematical model of a falling rock. Recall that the discrete mathematical model is based on Equations (3a,b) and (4a,b).

3.1 Pseudocode

Before setting out to write C or Fortran code to implement a mathematical model, it is a good idea to make sure the coding notation and coding logic are clear. A first step is the creation of language-independent pseudocode which we later can implement in a particular programming language. In this section, we develop such a pseudocode for the discrete-in-time mathematical model of a falling rock. There are many formatting styles for pseudocode, one of which is illustrated below.


Algorithm: Discrete in Time Mathematical Model of a Falling Rock

Step 1. Document

This algorithm uses the discrete in time mathematical model defined by Equations (3a,b) and (4a,b) to approximate the velocity and height of a rock dropped with initial velocity V0 from an initial height of Z0. The model is based on the assumptions that (1) air resistance is negligible and (2) acceleration due to gravity is constant.

INPUT

V0 -- initial velocity (double precision)
Z0 -- initial height (double precision)
g -- gravitational constant of acceleration (double precision)
h -- time step (double precision)
imax -- number of time steps to perform (integer)

OUTPUT

for i = 0 to imax
if height is positive print:
ti -- t-value a time step i (double precision)
zi -- approximate height at time ti (double precision)
vi -- approximate velocity at time ti (double precision)

AUXILIARY VARIABLES

vi_left -- approximate velocity at time ti - .5*h (double precision)
vi_right -- approximate velocity at time ti + .5*h (double precision)
i -- loop counter (integer)

Step 2. Initialize numerical solution at t = 0

Set ti = 0
      vi = V0
      zi = Z0
Print ti, vi, zi

Step 3. Calculate velocity at t = .5*h

Set vi_left = V0 - 0.5*g*h

Step 4. Begin time stepping

For i = 1,2,...,imax
    Do steps 5-7

Step 5. Advance solution one time step

Set ti = ti + h
      zi = zi + vi_left*h
      vi_right = vi_left - g*h
      vi = .5*(vi_left + vi_right)

Step 6. Output numerical solution

If zi < 0, print "rock is on the ground" and terminate calculation
else
Output ti, zi, vi

Step 7. Prepare for next time step

Set vi_left = vi_right


3.2 C Code

This C source code for implementing the pseudocode of Section 3.1 was developed and tested using the gnu C compiler. The format string " %lf " for dealing with variables of type double may have to be changed to read " %f " when the code is compiled by other compilers.

3.3 Fortran Code

Here is Fortran source code implementation of the pseudocode of Section 3.1.

4. Computational Experiments and Model Refinements

4.1 Computational Experiments

To become familiar with the behavior of the discrete in time mathematical model for a falling rock, perform the following experiments.

Experiment 1. Compile the C or Fortran code. Then make several runs with Z0 = 100 meters, V0 = 0 meters/second and in different runs choose different values for the time step h. Choose imax sufficiently large so that the rock hits the ground. Compare the predicted "time to impact" values for the different time steps. How does the model behave if the user should accidentally enter a negative value of Z0 and/or time step h?

Experiment 2. Compile the C or Fortran code. Choose Z0 = 100 meters and, based on your observations from Experiment 1, choose a fixed value of time step h. Then make several runs with different V0 values and try to determine the value of V0 that will result in a rock trajectory that rises to a maximum height of 150 meters before it begins to fall back to the ground. If you have studied calculus, also try to use the continuous in time model to answer this question.

Experiment 3. Modify the C or Fortran code so that the modified code prints out the difference between the continuous in time model and the discrete in time model at each grid point. Run the modified code for various choices of Z0, V0 and h. What do you conclude about the accuracy of the discrete in time model?

4.2 Model Refinements

Model Refinement 1. Suppose it has been observed that the falling rock does not continue to accelerate, but instead reaches a terminal velocity of VMAX (VMAX < 0).  In the continuous in time model, if we make the modeling assumption that

V'(t) = - g(1 - V(t)/VMAX)

would this produce the correct qualitative behavior? Can you find the exact solution to a continuous in time model with the above expression for V'(t)? Can you support such a modeling assumption based on a force balance between the gravitational force and the frictional (due to air resistance) force acting on the rock?

Many other similar modeling assumptions will produce mathematical models in which a terminal velocity is attainable. For example,

V'(t) = - g[1 - (V(t)/VMAX)2]

will result in a falling rock mathematical model with a terminal velocity. How does a modeler know which mathematical model is a better description of the physical problem?

Model Refinement 2. Suppose it has been observed that the falling rock does not continue to accelerate, but instead reaches a terminal velocity of VMAX. In the discrete in time model, if we make the modeling assumption that

vi_right = vi_left - g*(1 - vi_left/VMAX)    provided vi_left > VMAX

and

vi_right = vi_left    if vi_left <= VMAX

would this produce the correct qualitative behavior? Can you modify one of the source codes to incorporate the above modeling assumption?

Model Refinement 3. After reading the spring-mass modeling project, can you reformulate the falling rock model based on Newton's 2nd Law and the definition of velocity? If you assume that air resistance increases in proportion to fall velocity, can you give a force balance argument that leads to a discrete in time model in which the rock reaches a terminal velocity? Repeat for the continuous in time model, if you have studied calculus.