Project 7: Skydiver with air resistance

Wednesday, April 23, at 4:00 PM: No pseudocode. Instead, please write detailed equations showing Heun's method. Write down on paper exactly how you will go from current values of height and velocity to the values at the next timestep. Take it in two pieces. First, write down the equations you will use to get the new velocity. Then, write down the equations you will use to get the new height. Use Heun's method in each case.
Monday, April 28, at 9:00 AM: The finished MATLAB code and analysis

Your job in this project: calculate the height and speed of a man who jumps out of an airplane, just like last week. This time, however, you are take into account two complicating factors:

air resistance slows him down
the density of air decreases with altitude

You will also have to compare Euler's method to Heun's method.

The situation

As before, an airplane moves in level flight at a height of H = 10,000 meters above the ground. Joe opens the door and steps out at t = 0. What are his

height y(t)
velocity v(t)

as a function of time?

The complicating factors are two: first, Joe feels air resistance as he falls. The force of air resistance is

                                                    2
         F(air)  =  C * (density) * (area) * (speed)

where

         F(air)          is the force of air resistance (Newtons)

         C               is the drag coeff, C = 0.6 for a person

         density         is the density of the air at current altitude 
                               (see below)

         area            is Joe's cross-section, 0.5 square meters

         speed           is Joe's current speed (meters/sec)

You will also need to know Joe's mass:

         M               is Joe's mass, 80 kilograms

The second complicating factor is the change of the density of air in the Earth's atmosphere with altitude. We don't ordinarily notice it, because the scale height of the change is much larger than the altitudes we typically cover in a short time. To a decent approximation, the density of air is

         rho(height)  =   (1.21 kg/m^3) * exp(-height/8000 m)

It might be very useful to write a very short MATLAB function which calculates the density of air given a height above the ground. That might make your main program shorter and easier to understand, since some of the details would be "hidden" in the other function.

The assignment

You should write a MATLAB routine which calculates position and velocity as a function of time. It should look like this:

  function skydiver2(start_height, timestep, method, output_file)

where

     start_height            is the starting height, in meters

     timestep                is the size of the timestep (t1 - t0) to use
	                           in calculations, in seconds

     method                  is 'Euler' to use Euler's method
                                'Heun'  to use Heun's method

     output_file             is the name of a file in which you will
                                   write values of time, height, and speed

As you go through a loop, calculating height and speed at each time, you should write the values into a text file, for later use.

The output file should contain data just like last week's project; it should be three columns of numbers, something like this:

  0.0          10000.0         -0.0
  1.0           9995.1         -9.8
  2.0           9980.4        -19.6

where the first column is time (in seconds), the second column is height above the ground (in meters), and the third column is velocity (in meters per second, upwards positive and downwards negative).

Here's what you need to do:

First, compare the two methods when there is no air resistance (just set the C coefficient to zero). Calculate Joe's height and velocity at a time of 10 seconds after he jumps, using both methods and timesteps of 1.0, 0.5, and 0.25 seconds. Calculate Joe's height and velocity at this time analytically as well. Determine the fractional relative error in height and velocity for each method, for each choice of timestep.
Discuss how each method improved as the timestep was decreased. Be quantitative.
Now run the calculations with air resistance. Use both methods, and timesteps of 1.0 and 0.5, and 0.25 seconds.
Consider the height and velocity of the skydiver after exactly 20 seconds have passed. If a method works perfectly, it should yield the same values no matter what timestep is used. How much do the values for Euler's method change when one cuts the timestep in half? Be quantitative:
The difference in height at 20 seconds for a 1-second timestep vs. a 0.5-second timestep is 20.2 meters. That's a fractional change of 0.005 in height. When we decrease timestep from 0.5 seconds to 0.25 seconds, the height at 20 seconds changes by 4 meters. That's a fractional change of 0.001.
What about Heun's method? Be quantitative.
Again, discuss the improvement in each method as the timestep was decreased.
When we cut the timestep in half, the fractional change decreases by a factor of five ...

Bells and Whistles

Add one more complication: at a time of exactly 10 seconds after he jumps, the skydiver pops open a parachute. The parachute's cross-section area grows linearly with time, like so:
```
           t < 10 seconds      area = 0
      10 < t < 12 seconds      area = 30 sq.m * ((t - 10)/2)
      12 < t      seconds      area = 30 sq.m
  
```
In other words, the parachute opens gradually from 0 square meters to 30 square meters over 2 seconds. Modify your program so that it calculates Joe's motion during the entire fall, both before, during, and after the parachute opens. Look carefully at the velocity and height just after Joe opens the parachute. Does Joe's height ever increase? Should it? Think carefully, and explain.
Write a routine which uses the fourth-order Runge-Kutta technique to integrate the equations of motion. Compare the results to those using Euler's method.

This page maintained by Michael Richmond. Last modified Apr 29, 2003.