Stimulation of simple pendulum with damping using Mat-lab

Objective:

The aim of the project is to stimulate the movement of the simple pendulum(0-20s) which moves according to the below second order differential equation(ODE) and to create an animation of its movement.

`(d^2 θ)/(dt^2 )+b/m*(dθ/dt)+g/l (sin⁡θ )=0`

where:

b - damping coefficient

l - length of the string

g - acceleration due to gravity

m - the mass of the bob

Procedure:

The equation of simple pendulum using damping is a second-order differential equation(ODE). So, it must be divided into two lower order differential equations as

`(dθ_1)/dt=θ_2`

`(dθ_2)/dt=-(b/m)*θ_2-(g/l)*sinθ_1`

However, some assumptions are made for converting the second order ode to first order ode.

they are

`θ=θ_1`

`dθ/dt=(dθ_1)/dt=θ_2`

`(d^2 θ)/(dt^2 )=d/dt ((dθ_1)/dt)=(dθ_2)/dt`

As the material is fixed the damping coefficient \'b\' is also fixed ( because it is a material property). Length of the string and mass of the bob are also specified.

%inputs
b=0.05; % damping coeffient
g=9.81; % acceleration due to the gravity
l=1; %length of the string
m=1; %mass of the bob

The initial displacement and velocity are taken as \'0\' and\'3\' respectively. The stimulation motion takes place for 0-20 sec as shown below.

%intial conditiond
theta_0=[0;3];
ct=1 %counter variable
%time Interval
t_span=linspace(0,20,500);

Then, slove the ode by calling the ode function which is defined and plot time(t ) and results that are obtained, the row specify the displacement and the second row specifies the velocity.

%calling ode function
[t,results]=ode45(@(t,theta) ode_func(t,theta),t_span,theta_0)

aa = plot(t,results(:,1)) % displacement vs time graph
hold on 
bb = plot(t,results(:,2)) %velocity vs time graph
legend([aa,bb],\'displacement\',\'velocity\');
hold off
plot(results(:,1),results(:,2))
xlabel(\'time(s)\')
ylabel(\'plot\')

The output plot is obtained as shown the figure(1). The ode function code is written below, then loop for the animation is created where the coordinates of the bob(x1,y1) and the fixed point(x0,y0)

output = results(:,1); %dispalcement of the pendulum
output1 = results(:,2); % velocity
% loop for animation
for i= 1:length(output)
    OUTPUT=output(i)
    x0 = 0;
    y0 = 0;
    x1 = -(l*sind(OUTPUT));
    y1 = -(l*(cosd(OUTPUT)));

after, the plot for these points and use marker command to circle bob and then create a movie file for the animation of the movement of a simple pendulum.

We can alternatively change the values of b or l or m values to get different animations and to study the nature of the pendulum

The full code of the program is given below:

ODE Function code:

function [dtheta_dt]= ode_func(t,theta)

b=0.05; % damping coeffient
g=9.81; % acceleration due to the gravity
l=1; %length of the string
m=1; %mass of the bob


theta1=theta(1);
theta2=theta(2)
%first ode
dtheta1_dt=theta2;
%second ode
dtheta2_dt=-((b/m)*(theta2))-((g/l)*sin(theta1));
dtheta_dt=[dtheta1_dt;dtheta2_dt]
end

ODE Main code:

%inputs
b=0.05; % damping coeffient
g=9.81; % acceleration due to the gravity
l=1; %length of the string
m=1; %mass of the bob
%intial conditiond
theta_0=[0;3];
ct=1 %counter variable
%time Interval
t_span=linspace(0,20,500);

%calling ode function
[t,results]=ode45(@(t,theta) ode_func(t,theta),t_span,theta_0)

aa = plot(t,results(:,1)) % displacement vs time graph
hold on 
bb = plot(t,results(:,2)) %velocity vs time graph
legend([aa,bb],\'displacement\',\'velocity\'); % for represinting the velocity and displacement curves
hold off
plot(results(:,1),results(:,2))
xlabel(\'time(s)\')
ylabel(\'plot\')


output = results(:,1); %dispalcement of the pendulum
output1 = results(:,2); % velocity
% loop for animation
for i= 1:length(output)
    OUTPUT=output(i)
    x0 = 0;
    y0 = 0;
    x1 = -(l*sind(OUTPUT));
    y1 = -(l*(cosd(OUTPUT)));
    %plotting
    plot([-0.01 0.01],[0 0],\'linewidth\',3,\'color\',\'c\') %plotting the base of the pendulum
    hold on
    plot([x0 x1],[y0 y1],\'linewidth\',2,\'color\',\'b\')
    hold on
    grid on
    plot(x0,y0,\'marker\',\'o\',\'markersize\',5,\'markerFaceColor\',\'r\'); %plotting the origin/fixed point of pendulum using marker command
    plot(x1,y1,\'marker\',\'o\',\'markersize\',20,\'markerEdgeColor\',\'k\'); %plotting the bob using the marker command
    axis([-0.02 0.02 -1.5 0.5])
    pause(0.003)
    hold off
   
    
    
    
    M(ct)=getframe(gcf)
    ct=ct+1
end

movie(M)
file=VideoWriter(\'simplependulum.avi\', \'uncompressed AVI\')
open(file)
writeVideo(file,M)
close(file)

Error:

• Initially, the motion of the pendulum is inverted(upside motion). it is rectified by changing the axis values to [-0.02 0.02 -1.5 0.5]
• Another error is case sensitive error which is rectified by writing the code correctly.

Result:

The plot (figure(1)) below represents the variation of the velocity and displacement for different amplitude and time.

Conclusion:

The above program stimulates the simple pendulum motion(second order differential equation) and animation of the movement of a simple pendulum is created as shown below