Intro to Scientific Computing

PHYS 27/193

Physics Department
University of the Pacific

A Worked Example: Theory meets Experiment

In Homework 2 you visited Vertical Visions, a BASE jumping video, photography, and stunt company who have data collected from a long freefall (here) You should all have a file with this data.

I decided to check out how good this data is.

Also, the distance data is in feet and speed data in miles/hour. It would be easier to convert this to meters and m/s. We can do this by taking advantage of the using feature when plotting data.

To find how to convert, I will use the Unix units program, a wonderful little text based program that gives you conversion factors for one set of units to MANY others.

In an xterm, you should be able to type "units" at the prompt and get: 

hadron[~]|>units

2438 units, 71 prefixes, 32 nonlinear units

You have: feet
You want: meters
        * 0.3048
        / 3.2808399
You have: mph
You want: m/s
       * 0.44704
        / 2.2369363


Now, I can set two conversion factor constants in gnuplot that I will use later. 

gnuplot> f2m = 0.3048
gnuplot> mph2mps = 0.44704


Theoretical Expectations from Physics

You probably know that for an object in free fall in the Earth's atmosphere, it begins to drop as if there is no air resistance, so the distance from the starting point increases as

where g = 9.8 m/s2.

At first, before the object is moving very fast, this will be a good approximation. So let's check it.

Do: 

gnuplot>plot "MYfreefall.dat" u 1:(f2m*$3), 0.5 * 9.8 * x**2

At the beginning of the graph (for small t), does the function match the data? Do you understand how I converted it to meters by simply using the using command and the constant f2m? Make sure you do.

Fairly quickly, the force of air resistance becomes apparent and once the force of gravity and the force of air resistance are equal, the object will no longer accelerate (right? F = ma = Fgrav - Fair for those who haven't had Physics in a while).

Once the acceleration is zero, the speed will be constant. This speed is called terminal velocity. We can see that on our plot, after about 8 seconds--the position increases linearly. We can estimate the terminal velocity by plotting a line along with the data and varying its slope until we get it to match the linear part of the data, after 8 sec. It helps to offset the line with a negative b constant (a*x - b).

Try to estimate the terminal velocity.


Ok. The position data in column 3 seems to make sense: quadratic at the beginnig, and linear at the end, as expected. We'll deal with the velocity data later.

Now for the theory. The governing equation here is Newton's

F = ma

We all know that one, right? However, What it really means is this.

where F is the total force Fgrav - Fair, m is the mass, and the acceleration is the second derivative of position x(t) with respect to time. I am being pedantic here for those who have not had calculus for some time. I will shortly get to differential equations--if you have not had DiffEQs, don't worry--this is just an example. Nonetheless, follow as much as you can.

Recall that the force on a falling object due to gravity is constant

Now, it turns out that except for rather slow speeds, the force of air resistance can be approximated by a force which grows as the square of velocity

The constant k depends on a number of things

Let's simplify the differential equation:

In the last differential equation, I defined b = k/m. This is the one we have to solve for the function x(t).

I will skip the solution of this equation (it is done in many books on Differential Equations, as well as many jr. level Physics texts. There is also a good webpage review that I found: Auburn Engineering Course).

The function x(t) is

where g is the acceleration of gravity, and b = k/m is as defined above. Cosh() is the hyperbolic cosine. You should check that this satisfies the differential equation in the last line above!



Theory meets Experiment

If you visit the webpage above to look at the method of solution, or even if you just verify that the solution satisfies the DiffEQ, you can see that this was a non-trivial result. How well does it match the experimental data?

All we have to do to check is to fit the solution to the data. Remember that gnuplot (almost) always uses x for the horizontal axis even though we are calling it t in our theoretical discussion. Similarly, we have called the position x(t) even though it's been plotted on the y axis. ...Just wanted to be clear. 

gnuplot>f(x) = 1/b*log(cosh(sqrt(b*g)*x))
gnuplot> fit f(x) "MYfreefall.dat" u 1:(f2m*$3) via g, b
...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================

g               = 10.0539          +/- 0.0284       (0.2824%)
b               = 0.00338874       +/- 2.414e-05    (0.7123%)
...

gnuplot>plot "MYfreefall.dat" u 1:(f2m*$3), f(x)

BANG!

The theoretical function is really good fit to the data. Also the fitted value for g is rather close to the known value (9.8): 

gnuplot> print (g-9.8)/9.8 * 100
2.59118938657317

within about 2.6%.

Notice that in the fit report above, the error on the value for g is 0.284 = 2.8%. Thus the data is consistent with the known value of g=9.8 m/s.

That's not too bad for pulling some data from an extreme sports webpage and applying the laws of physics!

Velocity

Next lets look at the velocity. Given the position, x(t), the velocity is just the time derivative

Check this. You may have to consult your math book to recall the derivative of cosh(x).

Now we can check the velocity data against the data. Recall that the velocity data is in the stupid units: mph = miles per hour. However, our constant mph2mps takes care of that. So, define the function v(x), and plot the data and theory. 

gnuplot>v(x) = sqrt(g/b)*sinh(sqrt(b*g)*x)/cosh(sqrt(b*g)*x)
gnuplot>plot "MYfreefall.dat" u 1:(mph2mps*$2), v(x)

The match between theory and experiment is not good.

The key thing to realize here is that velocities are usually measured by taking the initial and final positions over some time interval dt, then dividing the change in position dt by dt to get the average speed over that interval. Thus the velocity really represents the speed closer to the middle of the time interval. Thus the velocity recorded at, say, 3 seconds really should be reported as the speed at 2.5 seconds, since the last measurement of position was made at 2 seconds. Does this make sense to you? In other words, the velocity data should really be shifted back (left) by 0.5 second in time.

That's easy to do with using! 

gnuplot>plot [0:20] "MYfreefall.dat" u ($1-0.5):(mph2mps*$2), v(x)

Much better. I also added an extended plot range to the graph so we could see the approach to terminal velocity, and how long it takes to get to a constant speed.

We also see that the "Terminal Speed" reported in the original data is not quite correct.
The falling object has not yet obtained constant speed.

You might want to plot thus using a "using" trick like:

u 1:($1>7 ? $2 : 1.0/0)

which will block plotting any points after 7 seconds. Do you see how it works?

A crazy thing to try

Let's see if we can figure out the mass of the jumper!

From our fit, we know that b, the coefficient of air resistance, is

b = 0.0033

and

Hence

We can estimate the parameters in this equation. The density of air, , is about 1 kg/m3. The "Shape Coefficient", C, for a prone falling human is around 1 (it's probably slightly higher but to our level of accuracy, 1 will do; see the article on Drag Coefficients at Wikipedia).

Hint: Don't Trust Wikipedia There's lots of great information on Wikipedia. However, there is unfortunately lots of just wrong stuff mixed in, and there very little way to tell the difference without some work. The best thing is to find facts that have references (from reputable sources) and check them.

The surface of a falling person is also hard to estimate, but let's say a 6 foot person (~1.8 m) times a width of half a meter. That's 0.9 m2. That means that the mass of our falling person is about

136 kg = 300 lbs.

That seems a bit high, but it could easily be




Arrows and Labels

The graph at the top of this page captures my theoretical curves and data for the position and velocity from the extreme free fall experiment.

You can see in the plot that I added Labels and Arrows in the graph, which help annotate it to make it cleared. This is an example of presentation graphics--making a picture that you will publish to help someone understand your work.

To add these things is easy.

Arrows

The syntax of adding an arrow to a graph is: 

gnuplot>set arrow from 1,0.3 to 2,0.5

This puts an arrow starting at position (x=1, y=0.3) on the graph, to the position (2,0.5), putting an arrowhead pointing toward the last point. If you don't want arrowheads (meaning that you are just drawing lines!), you can add the option nohead at the end of the line.

You don't have to do anything in the plot command. Once you have set arrows, they will appear in any plot you make. You can "erase" all arrows by doing 

gnuplot>unset arrow

If you remember the order in which you set the arrows, or if you explicitly gave them numbers (as shown below) you can just do

unset arrow n

where n is the number; first arrow is 1, next is 2, etc.

You can also do: show arrow

(it's arrow not arrows)

I used this method when I made the plot of the energy levels in a quantum dot on the solution page of Homework 3. Have a look at the figure at the bottom of that page. The actual commands are here: 

set arrow 1 from 0, 0.28 to 1, 0.28 nohead
set arrow 2 from 0, 2.502 to 1, 2.502 nohead
set arrow 3 from 0, 6.976 to 1, 6.976 nohead
set arrow 4 from 0, 13.643 to 1, 13.643 nohead
set arrow 5 from 0, 22.455 to 1, 22.455 nohead
set arrow 6 from 0, 33.307 to 1, 33.307 nohead
set arrow 7 from 0, 45.811 to 1, 45.811 nohead
f(x) = (x<0 || x>1) ? 50 : 0
plot [-1:2][-10:70]  f(x) w line 3

The y values of the arrow positions are the energy values of the solutions, so this gives a nice graphical way to represent them.

Labels

Adding some text in the plot is just as easy. You simply do: 

set label 1 "Distance" at 10.3, 145
set label 2 "Velocity" at 1.2, 450

which put the text strings in ""s at the locations specified. You can left (the default) or right justify the test by adding left or right (or center or rotate, etc. See help set label) after the location.


Continuing Lines

On a plot command line, you often have too much to type. If you need to continue on the next line, just add \ (a single "backslash") at the end of the line. This will give you a new (shorter) prompt to continue the line. Here's an example: 

gnuplot> plot [0:100][-2:2] "data.dat" u 1:($3/200) w l 2,\
> "data.dat" u 1:4 w l 4

When you hit return, gnuplot will act as if you had typed all the above on one line. You can continue adding new lines with \ as many times as you want. Try it.

Plot cos(x), sin(x), and tan(x) by putting each on a separate line. 

gnuplot> plot cos(x),\
> sin(x),\
> tan(x)