Skydive2Demo.m
Contents
Overview
This script illustrates the solution of a first-order IVP
for finding the terminal velocity of a skydiver. The analytical solution is given by
Here stands for time, for velocity, and
- : gravitational constant
- : mass
- : seond-order drag coefficient
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.