Simulation of Transient Behaviour of Simple Pendulum using Matlab

Aim:

To simulate the transient behaviour of a simple pendulum using second order ODE's and to create an animation of it's motion using MATLAB software.

Theory:

Simple pendulum

A simple pendulum has a small-diameter bob and a string that has a very small mass but is strong adequate not to stretch significantly.

The period of a simple pendulum is given by T = 2π √l/g

Where, l = Length of a simple pendulum; g = acceleration due to gravity.

From the equation, we can write the relation between the time period of a pendulum and acceleration due to gravity as T is inversely proportional to √g. This relation is referred to as the fourth law of the simple pendulum. In ideal theory, the suspending material will have no mass but in realism that is not feasible.

Differential Equation of Oscillations

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.

Here values taken as below

L=1 metre,

m=1 kg,

b=0.05,

g=9.81 m/s2.

ode45

Solve nonstiff differential equations - medium order method
Syntax
[t,y] = ode45(odefun,tspan,y0)

Program steps:

clear all 
close all
clc

%input values 

b = 0.05;
l=1;
m=1;
g=9.81;
ct=1;

%initial conditions

theta_a=[0;3];   %initial condition of angular displacement= 0 and angular velocity is 3 m/s

t_span = linspace(0,20,500);  %time input

%solve ODE 
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,l,m,g),t_span,theta_a);


%Plotting
figure(1)
plot(t,results(:,1))     %getting the value of displacement
hold on
plot(t,results(:,2))     %getting the value of velocity
ylabel('plot')
xlabel('time')
hold off

for i=1:length(t)          %it will the length of t         
    tilt = results(:,1);   %the variable named tilt inputed the all results in row1 
    move = tilt(i);        %so upto the loop times it will run 
  
 % points of pendulum 
    
    x0=0;                               
    y0=0;
    
    x1=l*sin(move);
    y1=l*cos(move);
    
   figure(2)
plot([-0.5 0.5],[0 0],'linewidth',2,'color','b')
hold on
plot ([x0 x1], [y0 -y1], 'linewidth',2,'color','r')
axis([-2 1 -2 1]) 
pause (0.0003)
hold on
p = plot(x1,-y1,'O','MarkerFaceColor','red');
hold off
M(ct) = getframe(gcf);  %getting the frame and storing in M
ct=ct+1;    
end
movie(M)
videofile = VideoWriter('pendulum oscillation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Equation solving

The differential equation is converted into two first order differential equation

Note: If the differential equation is of 3rd order then we need to convert into 3 first order equation for solving

`(d^2 theta)/dt^2 + (b/m)(d theta)/dt+g/l(sintheta)=0`

`(d theta_2)/dt = -b(theta_2)/m-g/l(sintheta1)`

`(d theta_1)/dt = theta_2`

Steps done in Program :

• Used the 'ode45' special solver in Matlab
• 'ode45' solution solver is used to solve systems of equation that involve a mass matrix
• Create a function with the given variable to get the result of calculating 1st and 2nd equations in matrix form
• Plot a Displacement-Velocity graph with respect to time
• Create a 'for' loop to plot the graph with repetitive result recorded at different angular displacement with respect to time.
• Use getframe and counter command inorder to create a movie
• The above graph represents the angular displacement and velocity with respect to 20 sec
• Damping Coefficient (b) is a material property that indicates whether it will bounce back or return energy to a system
• Create a movie using 'VideoWriter' command after ending the loop

Function format:

function [dtheta_dt] = ode_func(t,theta,b,l,m,g)

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

Steps :

• Create a function with all the given variables
• Function is created here inorder to solve the 1st and 2nd equation in matrix form
• Calling this function in ode45 solver

Output

Animation of transient behaviour of a simple pendulum is attached below:

https://www.youtube.com/watch?v=QlyZHfZ0-7U

Errors faced while executed the program

Modified damping coefficient (b) to get a result as mentioned below:

L=1 metre,

m=1 kg,

b=1,

g=9.81 m/s2.

clear all 
close all
clc

%input values 

b = 1;
l= 0.05;
m= 1;
g= 9.81;
ct= 1;

%initial conditions

theta_a=[0;3];   %initial condition of angular displacement= 0 and angular velocity is 3 m/s

t_span = linspace(0,20,500);  %time input

%solve ODE 
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,l,m,g),t_span,theta_a);


%Plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
hold off

for i=1:length(t)
    tilt = results(:,1);
    move = tilt(i);

    x0=0;
    y0=0;
    
    x1=l*sin(move);
    y1=l*cos(move);
    
   figure(2)
plot([-0.5 0.5],[0 0],'linewidth',2,'color','b')
hold on
plot ([x0 x1], [y0 -y1], 'linewidth',2,'color','r')
axis([-2 1 -2 1]) 
pause (0.0003)
hold on
p = plot(x1,-y1,'O','MarkerFaceColor','red');
hold off
M(ct) = getframe(gcf);  %getting the frame and storing in M
ct=ct+1;    
end
movie(M)
videofile = VideoWriter('pendulum oscillation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Output (by modified the value of damping coefficient)

Conclusion:

The plot describes as the time increases the velocity and displacement decreases since this is due to damping factor. The simulation using MATLAB is also completed successfully.