Week 3 - Solving second order ODEs

Aim - Write a Matlab code that will simulate the pendulum motion using second order ODE .

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 .

2.To create animation video for Pendulum motion .

Governing equations -

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(θ) ).

Code -

clear
close all
clc

% inputs 
g = 9.81 ;
b = 0.5;
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 .

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

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*theta2 - (g/l)*sin(theta1);
dtheta_dt = [ dtheta1_dt ; dtheta2_dt ];


end

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 frames

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 -

1. 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 .
2. The animation vedieo of the motion has been made .

Animation video -

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