Chebfun Logo
Oxford University
Mathematical Institute

A Bouncing Ball

Filomena Di Tommaso, February 26th, 2013

(Chebfun example ode/BouncingBall.m)
[Tags: #linearODE, #Ball, #physics]

This example simulates a bouncing ball. The trajectory of the ball is a piecewise quadratic chebfun.

Initialise:

clc, clear all, close all
trajectory = chebfun;

Parameters:

g = 9.81;                   % Force of gravity
x0 = 0;                     % Initial point
y0 = 0;
v0 = 30;                    % Initial velocity
theta(1) = pi/4;            % Starting angle
k = 0.9;                    % Coefficient of rebound
r = 1;                      % Radius of the ball

i = 1;                      % Counter of the number of parabolas
v_0x = v0*cos(theta);       % Initial velocity in the x direction
v_0y(i) = v0*sin(theta);    % Initial velocity in the y direction

The main loop:

while ( theta > 0.15 )

    % Time at which the ball impacts the ground:
    t_ground(i) = (v_0y(i)+sqrt(v_0y(i)^2+2*g*y0))/g;

    % x coordinate of the impact point:
    x_ground(i) = x0+v_0x*t_ground(i);

    % Motion law;
    x = chebfun('x', [x0, x_ground(i)]);
    y = y0 + (-g/(2*v_0x^2))*(x-x0).^2 + (v_0y(i)/v_0x)*(x-x0);

    % Piecewise chebfun:
    trajectory = [trajectory ; y];

    % Velocity at the impact point:
    v_fy=v_0y(i)-g*t_ground(i);

    % The y component of the velocity decrease according to the coefficient of
    % rebound k:
    v_0y(i+1)=-k*v_fy;

    x0 = x_ground(i);
    y0 = 0;

    % Corner rebound ball:
    theta(i) = atan(v_0y(i+1)/v_0x);

    i = i+1;
end

The final part of the trajectory is horizontal:

final_part = chebfun(@(y) 0*y);
trajectory = [trajectory ; final_part];

I shift the trajectory vertically by a quantity equal to the radius r in order to visualize the impact with the ground:

trajectory = trajectory + r;

Parameter to plot the circle:

n = 50;
t = linspace(0, 2*pi, n);
t_ground = [0, t_ground, 5];
x_ground = [0, x_ground, 900];
% pause(2)

c = 1;
ks = 120;   % Coefficient of elasticity of a ball
m = 0.313;  % Weight of the ball in gr

Plot the result nicely:

axis([0, 900 , 0, 25])

% Since I can't set axis equal to see the ball I must increase the eccentricity.
a = 25;
b = 1;
M = [];
for j = 1:i
    % Temporal interval
    time = linspace(0, t_ground(j+1), n);
    xt = x_ground(j) + v_0x*time;
    ballx = xt(1) + a*r*cos(t);
    bally = trajectory(xt(1)) + b*r*sin(t);
    h = patch(ballx, bally, 'g'); hold on
    g = plot(ballx(1), bally(1), '*r'); % Point on the ball
    drawnow
%     M = [M getframe()];
    for k = 2:n
        if ( b < 1 )
            b = b + ((1-b)/(i)*j);
        else
            b = 1;
        end
        ballx = xt(k) + a*r*cos(t);
        bally = trajectory(xt(k)) + b*r*sin(t);
        set(h, 'XData', ballx, 'YData', bally)
        set(g, 'XData', ballx(n-k+1), 'YData', bally(n-k+1))
        drawnow
%         M = [M getframe()];
    end
    b = r - (1/2*m*v_0y(j)^2)/ks;
    ballx = xt(k) + a*r*cos(t);
    bally = trajectory(xt(k)) + b*r*sin(t);
    set(h, 'XData', ballx, 'YData', bally)
    set(g, 'XData', ballx(n-k+1), 'YData', bally(n-k+1))
    drawnow
%     M = [M getframe()];

    set(h, 'XData', ballx, 'YData', bally, 'Visible', 'off')
    set(g, 'XData', ballx(n-k+1), 'YData', bally(n-k+1), 'Visible', 'off')


end

We now write the images to an animated .gif!

% for j = 1:5:numel(M)
%     im = frame2im(M(j));
%     [imind,cm] = rgb2ind(im,16);
%     if j == 1;
%         imwrite(imind, cm, 'html/BouncingBall.gif', 'gif', ...
%             'Loopcount', inf, 'DelayTime', 1e-5);
%     else
%         imwrite(imind, cm, 'html/BouncingBall.gif', 'gif', ...
%             'WriteMode', 'append', 'DelayTime', 1e-5);
%     end
% end

HTML Comment Box is loading comments...

Please contact us with any questions and comments.
Copyright © 2013, The University of Oxford & The Chebfun Team.