CS 310 Team Lab 9: Advanced Functions: fzero, quad, and ode45

SOLUTION - Spring 2013

Contents

Problem 1: Volume of a Cut Hemisphere

For function definitions, see

1.a: Define (and test) areaCutHemisphere()

clear
area_0 = areaCutHemisphere(0);
disp(['areaCutHemisphere(0) = ', num2str(area_0), ...
      '  ( expected result: ', num2str(0.5*pi*9) ' )'])
area_3 = areaCutHemisphere(3);
disp(['areaCutHemisphere(3) = ', num2str(area_3), ...
      '  ( expected result: 0 )'])

areaCutHemisphere(0) = 14.1372  ( expected result: 14.1372 )
areaCutHemisphere(3) = 0  ( expected result: 0 )

1.b & 1.c: Compute volume of portions of the hemisphere

vol_entire = quad(@areaCutHemisphere, -3, 3);  % vol of entire hemisphere

% 1.c: Compute volume of half of hemisphere
vol_half = quad(@areaCutHemisphere, -3, 0);   % vol of half the hemisphere

% Since the volume of an entire sphere is (4/3)pi r^3, we can check our
% results. The sphere has radius 3, so the entire sphere has volumen 36pi.
vol_sphere = (4/3)*pi*3^3;
disp(['volume of entire hemisphere = ', num2str(vol_entire), ...
      '  ( expected result: ', num2str(vol_sphere/2), ' )'])
disp(['volume of half hemisphere = ', num2str(vol_half), ...
      '  ( expected result: ', num2str(vol_sphere/4), ' )'])

% Define volumeCutHemisphere() and compare its results to those above
disp(['volumeCutHemisphere(3) = ', num2str(volumeCutHemisphere(3))])
disp(['volumeCutHemisphere(0) = ', num2str(volumeCutHemisphere(0))])

volume of entire hemisphere = 56.5487  ( expected result: 56.5487 )
volume of half hemisphere = 28.2743  ( expected result: 28.2743 )
volumeCutHemisphere(3) = 56.5487
volumeCutHemisphere(0) = 28.2743

1.d: Find value of t so volume of cut-off hemisphere is 20

t_cut = fzero(@toBeZeroCutHemisphere20, 0);
disp(['Cut hemisphere at t = ', num2str(t_cut)])

% check result from fzero is correct by plugging into volumeCutHemisphere()
vol_t_cut = volumeCutHemisphere(t_cut);
disp(['Volume at cut = ', num2str(vol_t_cut)])

Cut hemisphere at t = -0.59301
Volume at cut = 20

Problem 2: Jumping Out of an Airplane (First-Order ODE problem)

See also parachuteODE.m

clear
timespan = [0 20];
v0 = 0;
[t_soln, w_soln] = ode45('parachuteODE', timespan, v0);
plot(t_soln, w_soln)

2.c: Plot the time (in order to see the step size)

plot(t_soln, 'o')

Problem 3: How fast AND how far is the parachutist falling?

See also parachuteODE2.m

clear; clf;
timespan = [0 20];
w0 = [0 0];
[t_soln, w_soln] = ode45('parachuteODE2', timespan, w0);

% w_soln contains 2 columns, one for v and one for y
v_soln = w_soln(:, 1);
y_soln = w_soln(:, 2);

% plot velocity (v) using a blue dotted line and distance fallen (y) using
% a green solid line
plot(t_soln, v_soln, 'b:', t_soln, y_soln, 'g')

Questions

Q: How far will the parachutist fall in 20 seconds?

A: Since the timespan is 0 to 20 seconds, the distance fallen in 20 seconds will be the last value in the y_soln vector.

Q: How fast is the parachutist falling after 20 seconds?

A: Since the timespan is 0 to 20 seconds, the velocity after 20 seconds will be the last value in the v_soln vector.

Q: What is the value of the terminal velocity?

A: The terminal velocity is 68.6 m/s We can get this by solving the original parachute ODE for longer timespans and seeing what the last value in the v_soln vector is.

disp(['Distance fallen in 20 seconds: ', num2str(y_soln(length(y_soln)))])
disp(['Velocity after 20 seconds: ', num2str(v_soln(length(v_soln)))])

for t_end = [50 : 25 : 200]
    [t_soln, v_soln] = ode45('parachuteODE', [0 t_end], 0);
    disp(['At time ', num2str(t_end), ', velocity = ', ...
         num2str(v_soln(length(v_soln)))])
end

Distance fallen in 20 seconds: 919.3792
Velocity after 20 seconds: 64.6601
At time 50, velocity = 68.5457
At time 75, velocity = 68.5985
At time 100, velocity = 68.6
At time 125, velocity = 68.6
At time 150, velocity = 68.6
At time 175, velocity = 68.6
At time 200, velocity = 68.5999

Published with MATLAB® 7.14