function ball(initial_velocity, timestep)
% Program to compute the position and velocity of a ball
%    which is thrown vertically upward from ground level.
%    We make simplifying assumptions:
%          - no air resistance
%          - constant force of gravity
%
% User supplies the initial velocity of the ball and
%    the value of the timestep for numerical integration.
%    We use simple Euler's method to compute the motion.
%
% Print out lines which include mechanical energy during each timestep.
%
% Usage: ball (initial_velocity, timestep)
%
%    where
%             initial_velocity       is initial speed of ball
%                                    vertically upward (m/s)
%
%             timestep               is size of steps in time (s)
%
% Returns: nothing
%
% MWR 11/28/2012


% initialize variables
%    current height above the ground (m)
y = 0;
%    new height above the ground (m)
y_new = 0;
%    current velocity (positive upward) (m/s)
vy = initial_velocity;
%    new velocity (m/s)
vy_new = 0;
%    time since start of simulation (s)
t = 0;
%    acceleration due to gravity (positive up) (m/s^2)
g = -9.8;
%    mass of ball (kg)
m = 1.0;


% compute initial mechanical energy
ke = 0.5*m*vy*vy;
gpe = m*abs(g)*y;
initial_energy = ke + gpe;
% fprintf(1, ' ke %9.3f  gpe %9.3f  tot %9.3f \n', ke, gpe, energy);

% loop over time; continue as long as ball is above ground
while (y >= 0)

     %  compute mechanical energy
     ke = 0.5*m*vy*vy;
     gpe = m*abs(g)*y;
     energy = ke + gpe;
     frac_error = (energy - initial_energy) / initial_energy;


     %  print current time, position velocity, energy, frac error in energy
     fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  E %12.6e  f %9.4e \n', t, y, vy, energy, frac_error);

     %  use const accel to compute new velocity
     vy_new = vy + g*timestep;

     %  use current velocity to compute new position
     y_new = y + vy*timestep;

     %  update position and velocity for next time through loop
     vy = vy_new;
     y = y_new;

     %  update time
     t = t + timestep;

% end of loop over time
end


%  print current time, position velocity, energy, frac error in energy
fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  E %12.6e  f %9.4e \n', t, y, vy, energy, frac_error);



% end of program
end

Published with MATLAB® 7.13