Ordinary Differential Equation Of A Pendulum (Matlab/Simulink)

Aim

1. Write a Matlab code that will simulate the pendulum motion using second order ODE .
2. Represent the above code using simulink .

Introduction -

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum. 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.

Objective -

1 .

• To simulate second order ODE of Pendulum .
• To create animation video for Pendulum motion .

2.

• To simulate the differential equation in simulink by calling a matlab script .

Governing equations -

The equation we got from newtons secnd law of rotation .

In the above equation , we take the following values .

gravity in m/s2 (g) = 9.81

length of the pendulum in m (L) = 1

mass of the ball in kg (m) = 2

damping coefficient (b) = 0.9

Now we solve the second order ODE to first order ,

θ= θ₁

Now to find the position of the pendulum ,

Fist point be (x0,y0) = (0,0) .

Second point - (x1,y1) = ( l*sin(θ) , - l*cos(θ) ).

MATLAB-

Code -

clear
close all
clc

% inputs 
g = 9.81 ;
b = 0;
l  = 1 ;
m = 2 ;


A = [0;3]; % initial position and velocity 

% time span 
t_span = linspace(0,50,500);

% ODE solver used to solve first oder differential equation 
[ t, results,T_1] = ode45(@(t,A) ODE_Function_1(t,A,b,g,l,m) ,t_span,A);


% time period 
%T_1 = 2*pi*sqrt(l/g);

%plot 
figure;
subplot(2,1,1); % used to create different plots in a figure .
plot(t,results(:,1), 'linewidth',2); %  angular displacement 
hold on ;
plot(t,results(:,2),'linewidth',2); % velocity 
xlabel('Time');
ylabel('plot');
title(' plot V/S Time');
grid on ;
legend('angular displacement','velocity ');
subplot(2,1,2);
plot(results(:,1),results(:,2),'linewidth',2,'color','g');
xlabel('Angular displacement');
ylabel('velocity');
title('Angular displacement  V/S Velocity ') ;

grid on ;


% animation  for simulation of pendulum 
ct =1;
 for i = 1:length(t_span)
    u= results(i,1);
    
    x0 =0;     % initial position 
    y0 = 0;
    
    x1 = l*sin(u);      % final position of the string 
    y1 = -l*cos(u);
   
    % plot
    figure (3)
     
    % plottinng wall 
    
    plot([-1,1] ,[0,0],'linewidth',3,'color','g');
    hold on ;
    plot([-1,1] ,[0.01,0.01],'linewidth',3,'color','g');
    plot([-1,1] ,[0.02,0.02],'linewidth',3,'color','g');
    hold on ;
    % plottinng string
    
     plot([x0 x1] , [y0 y1] ,'linewidth',1,'color','r');
    
    % marker
    
    plot(0,-0.01,'o','MarkerSize',7,'markeredgecolor','k');
    
    %plotting ball 
    
    plot(x1,y1,'o','MarkerSize',15,'MarkerFaceColor','b');
    grid on ;
    xlabel('x-axis');
    ylabel('y-axix');
    title(' Pendulum Simulation ') ;
   
    axis([-2.5 2.5 -2.5  2.5  ]);
    hold off ;
    pause(0.01);
    
     M(ct)=getframe();
     ct=ct+1 ;
    
    
end 


%mywriter = VideoWriter('Pendulum_1.avi','Uncompressed Avi');
%open(mywriter);             % to open the movie file 
%writeVideo(mywriter,M);      % write to the file and storing the images  
%close(mywriter);             % close the file when done .
 

Explanation -

1. The first three lines of code is to clear the screen and start fresh .
2. The inputs taken for the pendulum is defined also the time span is defined .
3. The initial position and velocity is taken as a matrix in a variable "A" .
4. The ODE solver "ode45" is used to solve the first oder differential equation which we got by solving second order ODE of pendulum .The ode45 is an in-built matlab solver .

The syntax is -

[t,results] = ode45(@(t,A) ODE_Function_1(t,A,b,g,l,m) ,t_span,A) ,

here the solver is told that the dependent and independent variables are "A" and"t" in the first order ODE

then a function is defined where the different inputs are passed ,which calculates the values using the

equations mentioned above .Then the initial conditions are passed which are used by the solver .

5. The plot command are then used to plot the graph between angular position and velocity with time .
6. A " subplot" command is used to plot - syntax - rows, column ,position .

function [dtheta_dt , T ] =  ODE_Function_1(t,A,b,g,l,m)
theta1 = A(1)
theta2 = A(2)
dtheta1_dt = theta2;
T = m*l^2*(theta2)^2 + m*g*cos(theta1) ;
dtheta2_dt =  -b/(m*l^2)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [ dtheta1_dt ; dtheta2_dt ];


end

7. The Function "ODE_Function_1" is used to calculate the values . It is defined as follows .

8 .After plotting the animation of the simulation is created .

9.The for loop is varied till the time span of simulation , the angular position obtained from ODE solver is used to calculate the final position of the string where the ball is atached .

10. A "mywriter" variable is defined to which the frames are written to ,this stores the frame

and forms the animation video for the pendulum motion .

result -

When the damping co-efficient is varied we can see different results -

1. (b=0 ) Undamped

2. (b=0.4) Lighty Damped

3 .(b>1) Critically Damped

Errors -

In the line "53" a negative sign was missing , which was the error , due to which the the pendullum animation came wrong .

Conclusion -

1. The main objective of simulating the second oder ODE of a pendulum is completed In matlab .
2. The animation vedieo of the motion has been made .

Animation Video -

Link -https://www.youtube.com/watch?v=IR9AJnqEBYw

SIMULINK -

It is a block diagram environment for Model-Based Design.The following ODE can be solved in simulink using block diagrams .

Block Diagram -

Explanation -

• The blocks which have been used are present in the simulink library .
• The block diagram was strted from considering the value for "theta dot" -( which is the second derivative of angular position -theta)
• it is then passed through an integrator block_1 and again the velocity which is obtained is integrated through block_2 ,which gives angular position .
• At the start the angular position is taken as zero . The valuse obtained from the integrator are multiplied using a gain operator .
• It then passes through a subtract block which at each position subtracts the value and maintains the flow of code .

Now the simulink code can be called through a matlab script where we can define the input values .

clear
close all
clc

g = 9.81 ;
b = 0.9;
l = 2;
m = 2 ;

X = linspace(0,200,500);

for i = 1: length(X)   % loop to allocate the time_span to the array "T" .
    T = X(i);
end

% Load model 
open_system('pendulum_simulink');

% Simulate model
sim('pendulum_simulink');

.The script will run and will create the variables , which are necesary to run the simulink model .

. The above values are used to run the simulink .

Results -

1. Three plots have been shown . The "red" plot is angular position , "blue" is velocity , " angular acceleration " all with respect to time .

Conclusion

1. The perfect simulink model for the second oder ODE is obtained .
2. A graph representing position ,velocity and acceleration is obtained .

Thank You .