Introductory Programming
Fall 2006

1) We started out using ode45 to solve a single, first-order equation:

ode45(@f, [0, 1], 0)

where the input variables are the name of the function file, the time
range, and the initial value of y.

In f.m, we wrote a function that computes the time derivative of y.

function dydt = f(t, y)
    dydt = t + y
end


2) Then we solved a system of three first-order equations:

[T, M] = ode45(@lorenz, [0,30], [3, 2, 23]);

The initial condition is a vector with the values x, y and z.

function dVdt = f(t, V)
    x = V(1);
    y = V(2);
    z = V(3);
 
    sigma = 10;
    b = 8/3;
    r = 28;
 
    dxdt = sigma * (y-x);
    dydt = x * (r-z) - y;
    dzdt = x*y - b*z;
 
    dVdt = [dxdt; dydt; dzdt];    % pack the results in a column vector
end


The result from ode45 is a vector T that contains the time
values, and a matrix M that contains one column for each
variable x, y and z.

To plot the components:  plot3(M(:,1), M(:,2), M(:,2))


3) Then we figured out how to rewrite a second order equation as
a system of first order equations:

ode45(@freefall, [0, 30], [4000, 0])

The initial condition is a vector that (we understand) contains
the initial position and velocity.

function dxdt = f(t, V)
    r = V(1);      % the first component is position
    v = V(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


4) The next step is to generalize this to a system of second-order
equations.  For example, the equations that govern a projectile can be
written with two second-order equations (as long as the projectile
stays in a plane, two dimensions will suffice)

R -> [rx; ry]
V -> [vx; vy]

dRdt = V
dVdt = A = [ax; ay]

ax = 0
ay = -g

Translating freefall to deal with more than one dimension is
remarkably simple if we use vectors:

function dWdt = proj(t, W)
    R = W(1:2);
    V = W(3:4);

    dRdt = V;                          
    dVdt = acceleration(t, R, V);

    dWdt = [dRdt; dVdt];    % pack the results in a column vector
end

function dVdt = acceleration(t, R, V)
    g = 9.8;      % acceleration of gravity in m/s^2
    dVdt = [0; -g];
end


All I did was to change the names of the variables that
are now vectors (and I changed V to W to avoid confusion)

W -> the vector we get from ode45 that contains [R; V]
V -> the currect velocity vector


* Put the previous code in a file named proj.m and use it to compute
the distance a projectile would travel in the air, given an initial
velocity V = [10; 10] m/s.

* Write a function named pdist that takes two input variables,
speed and angle, and computes the distance a projectile with
the given speed and launch angle would travel.

* What is the launch angle that maximizes the distance travelled?


4b) Of course, a projectile in a vacuum is boring.  Wind
resistance makes things more interesting!

The vector form of drag force looks like this:

Fd = -1/2 C rho A v V

where

Fd -> the drag force in vector form [Fdx; Fdy]
C  -> the drag coefficient (0.3 for a baseball)
rho-> the density of air (1.3 kg / m^3)
A  -> cross sectional area of the object (0.0042 m^2 for a baseball)
v  -> the magnitude of the current velocity
V  -> the current velocity in vector form [Vx, Vy]

* Write a function named drag_force that takes the current velocity V
as an input variable and returns Fd as an output variable.

(You might want to show me this function before you go on.)

* Now modify acceleration to compute and add in acceleration due to
drag.  The mass of the baseball is 0.145 kg.

* Use your code to compute the trajectory of a baseball that leaves
home plate at an inital height of 1m, initial speed of 40 m/s,
and launch angle of 45 degrees from horizontal.

You might find the following code helpful.

function distance = baseball(velocity, angle)
    options = odeset('Events', @events);
    theta = angle * pi / 180;
    vx = velocity * cos(theta);
    vy = velocity * sin(theta);
    [T, M] = ode45(@proj, [0, 100], [0, 1, vx, vy], options);
    plot(M(:,1), M(:,2))
    distance = M(end,1);
end

function [value,isterminal,direction] = events(t,W)
    % Locate the time when height passes through zero in a 
    % decreasing direction and stop integration.
    value = W(2);        % Detect height = 0
    isterminal = 1;      % Stop the integration
    direction = -1;      % Negative direction only
end

You can download it from http://wb/ip/code/baseball.m

* With wind resistance, is the optimal launch angle greater or less
than 45 degrees?

* How much farther would the ball travel at Coors Field, home of
the Colorado Rockies, which is at an altitude of 5200 ft, where
the density of air is only 0.73 kg / m^3?

* The Green Monster in Fenway Park is about 12m high and about 97 m
from home plate along the left field line.  What is the minimum
speed a ball must leave the bat in order to clear the monster
(assuming it goes off at the optimal angle)?


Regarding the drag coefficient of a baseball:

From http://carini.physics.indiana.edu/E105/texture.html

From Figure 2.1 in Adair's book, you can see estimates of how the
drag coefficient varies with speed for several types of balls. For a
baseball, the drag coefficient begins to drop below 0.5 (the smooth
sphere value) for speeds around 55-60 mph (25-27 m/s), dropping to
about 0.3 for speeds of 90 mph (41 m/s). As a result, the terminal
speed for a baseball is significantly greater than it would be for a
smooth ball (95 mph or 43 m/s, instead of 70 mph or 32 m/s). Since
most pitches ( other hard throws) and most hits have speeds greater
than 25 m/s, the seams play an important role in reducing drag on a
baseball as we currently play it. For example, it would be difficult
to hit a smooth baseball more than 300 ft (95 m). Instead, for initial
speeds of 120 mph (55 m/s) the drag coefficient of a baseball is about
half as large and hits of 450 - 500 feet are not uncommon for strong
hitters.


* How hard would it be to adapt your program so that the coefficient
of drag depended on velocity?