Motion of pendulum using matlab

Objective:

1. To write a programme to solve the ODE which is of second order using matlab.

2. To create a video showing the motion of pendulum.

Description:

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

The below ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

t = 0-20 sec,

Initial conditions:

angular displacement=0,

angular velocity=3 rad/sec2

Matlab code to get displacements and velocity plots:

close all
clear all
clc

%given variable conditions
l=1; m=1; b=0.05; g=9.81;

% range of time
t=[0 20];

%creating the ode using anonymous function
f = @(t,x) [x(2);-(b*x(2)/m)-(g*sin(x(1))/l)];

initial_conditions=[0 3];
[t,x] = ode45(f,t,initial_conditions);

displacement = x(:,1)
velocity = x(:,2);
acceleration = diff(velocity)./diff(t);


plot(t,displacement,'linewidth',1.5)
hold on
plot (t,velocity,'linewidth',1.5)
xlabel('time')
ylabel('Velocity,Displacement')
legend('Displacement','Velocity')
grid on
hold off

Matlab code to create a video of pendulum motion:

figure(2) %for animation

O = [0,0]; %to set the location of pivot point
axis(gca,'equal'); %aspect ratio of plot
axis ([-1 1 -1.5 0.1]); %limits of plot so that axis is independent of motion
grid on;


for i=1:length(t) %to get loop
 
   P= l*[sin(x(i,1)) -cos(x(i,1))]; %setting the position of bob at the tip of the string using sin and cos components
   
   pivot= viscircles(O,0.01); %circle in origin O as reference point with 0.01 radius
   
   string=line([O(1) P(1)],[O(2) P(2)]); %to create a line between origin and bob
 
   bob=viscircles(P,0.05) %circle at P as reference point with 0.05 radius
   
   pause(0.001); %time interval to update the plot
   
   M(i)=getframe(gcf); %to fetch frame for each plot at i
   
 if i<length(t)  %delete previous objects till the final loop
       delete(string);
       delete(bob);
       delete(pivot)
 end
end

%to open and create a video using video writer
video=VideoWriter('pendulum.avi','uncompressed AVI');
video.FrameRate=30;  %to set the frame rate for video
open(video);         %to open the video
writeVideo(video,M); %to attach fetched frames together
close(video);        % to close the process

Output:

As per the intitial conditions, angular displacement started from 0 and angular velocity started from 3.

From the above plot we can observe that, angular displacement and angular velocity is goes on deminishing still it aquires the stable position.

The motion of the pendulum is as shown in below video

The entire motion of the pendulum in one plot is shown below