Skydive2Demo.m

Contents

Overview

This script illustrates the solution of a first-order IVP

$$ v'(t) = g - \frac{c}{m} v^2(t), \quad v(0) = 0 $$

for finding the terminal velocity of a skydiver. The analytical solution is given by

$$ v(t) = \sqrt{\frac{g m}{c}} \tanh \left( t \sqrt{\frac{g c}{m}}\right). $$

Here $t$ stands for time, $v$ for velocity, and

The right-hand side function is coded in Skydive2.m

Initialization

clear all
clc
close all
g = 9.81; m = 68.1; c = 0.225;              % specific constants
v = @(t) sqrt(g*m/c)*tanh(t*sqrt(g*c/m));   % analytical solution
f = @Skydive2;                              % function handle for right-hand side function
t0 = 0; y0 = 0;                             % initial time and initial velocity
tend = 15;                                  % end time
N = 100;                                    % number of time steps
h = (tend-t0)/N;                            % stepsize

Call Euler

We will discuss this function later

[t,y] = Euler(t0,y0,f,h,N);

Plot

hold on
set(gca,'Fontsize',14)
xlabel('t','FontSize',14);
ylabel('Velocity','FontSize',14);
plot(t,y,'r','LineWidth',2);                % Euler solution

pause
tt=linspace(t0,tend,201);
plot(tt,v(tt),'g','LineWidth',2);           % analytical solution
legend('Euler solution','Analytical solution','Location','NorthWest');
hold off
disp(sprintf('Terminal velocity: %f m/s.', sqrt(g*m/c)))

Terminal velocity: 54.489999 m/s.


© Greg Fasshauer, IIT. Published with MATLAB® 7.11