Introductory Programming
Fall 2006


Today's instructions

1) create a file named freefall.m in the usual place you put scripts

2) type in the following code

function dXdt = freefall(t, X)
   r = X(1);      % the first component is position
   v = X(2);      % the second component is velocity

   drdt = v;                          
   dvdt = acceleration(t, r, v);

   dXdt = [drdt; dvdt];    % pack the results in a column vector
end

function a = acceleration(t, r, v)
   g = 9.8;      % acceleration of gravity in m/s^2
   a = -g;
end

X is a vector with two elements: position and velocity.
The function freefall unpacks X, computes drdt and dvdt, and
then packs the result into a vector.

The function acceleration computes the acceleration of a
hypothetical body under force of gravity.  Notice that
acceleration is (in general) a function of time, position
and velocity, but we are starting with the simple case
where acceleration is constant.

3) evaluate freefall from the command line using the
   values t=0 (seconds), r=4000 (meters), v=0 (m/s)

   Assuming that it works, what does the result mean?

4) use ode45 to compute the position and velocity of the
   body over a time range from 0 to 10 seconds, with the
   initial values t=0 (seconds), r=4000 (meters), v=0 (m/s)

   How long does it take the body to fall from 4000 meters
   to the ground?

5) Modify the acceleration function to compute the net
   acceleration of a 70 kg skydiver with drag constant c = 0.22.

   F_drag = c v^2

   What is the skydiver's terminal velocity?  How long does 
   he/she take to hit the ground?

6) Modify acceleration again so that after 30 seconds of free-fall
   the skydiver deploys a parachute, which (almost) instantly
   increases the drag constant to 2.7.

   What is the terminal velocity now?  How long (after deployment)
   does it take to reach the ground?


Second order systems
--------------------

The nice people in Natick have given us a tool that can
estimate the solution to any first order diff EQ; that
is, anything that can be expressed as

y' = f(t, y)

All we have to do is write a function that implements f

It takes t and y as input variables and produces the
value of y' at t and y as a result.


So what do we do with second order systems?

y'' = f(t, y, y')

ode45 can't handle them like that.

One solution is to replace y  with r
                      and  y' with v

Then we can write

a = f(t, y, y')

dv/dt = y'' = a

dr/dt = v


And then define a vector [r, v]. Now

dX/dt = [dr/dt, dv/dt] = [v, a]

And ode45 can handle that fine.

Punchline: you can write an nth-order diff EQ as a system of
	   n first-order diff EQs