Week 4.1 - Solving second order ODEs

AIM : 

To  write a program that solves the second order ODE

GOVERNING EQUATIONS : 

where 

g = gravity in m/s2

L = Length of the pendulum in meters

m = mass of the ball in kg

b=damping coefficient

OBJECTIVE : 

 To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

To Create an image that represents the position of the pendulum for each value of time

BODY OF THE CONTENT : 

1)PROGRAMMING LANGUAGE USED :

HTML( default )

2) CODE :

clear all
close all
clc
% Inputs
l=1 ;  % in metres

m=1; % in Kg

b=0.05 ; % damping coefficient 

g=9.81 ;  % in m/s^2

% Initial Conditions
theta_0 = [0 ;3];
% time points
t_span=linspace(0,20,500);
% solving ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%Plotting

plot(t,results(:,1))
hold on
plot(t,results(:,2))

ylabel('Plot')
xlabel('Time(seconds)') 
title('Damped Harmonic Motion of a Pendulum')
legend('Angular displacement ', 'Angular velocity')


ct =1;
 for i = 1:length(results(:,1))
      x_start = 0  ;   
       y_start = 0 ;
        x_end = l*sin(results(i,1)) ;
         y_end = -l*cos(results(i,1)) ;

          % plotting the animation
           figure(2)
            plot([-1 1], [0 0], 'linewidth', 5)
             hold on
              plot([x_start x_end],[y_start y_end],'linewidth',3)
               hold on

                plot(x_end,y_end,'O','markersize',10,'markerFacecolor','g')
                 title('Animation of transient behaviour of a Pendulum')
                  axis([-3 3 -3 1])
                   pause(0.05)
                    hold off 

                     M(ct) = getframe(gcf)
                      ct = ct+1

 end


 % Animation of the pendulum
  movie(M)
  videofile=VideoWriter('Pendulum Animation.avi', 'uncompressed avi')
  open(videofile)
  writeVideo(videofile,M)
  close(videofile)

 SUPPLEMENTARY PROGRAM LIKE FUNCTION USED IN PROGRAMMING :

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

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 

 STEPWISE EXPLANATION OF CODE :

clear all

close all

clc

% Inputs

% defining the inputs as below

l=1 ;  % in metres          % length of the string = 1 meter

m=1; % in Kg.                 % mass of the ball=1 Kg

b=0.05 ; % damping coefficient .      % damping coefficients of damper system

g=9.81 ;  % in m/s^2.               % acceleration due to gravity 

% Initial Conditions

theta_0 = [0 ;3];.            %solution vector consists of values of initial displacement (here it is 0) and initial velocity (here it is 3)

% time points

t_span=linspace(0,20,500);.  %  creating time array from 0 to 200 with 500 values using linspace command

% solving ODE

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%ode45 is inbuilt command that solves  ordinary partial differential equations in matlab

@(t,theta).        %  main inputs where t and theta are dependent and independent variables of ODE system

%Plotting

plot(t,results(:,1))        % plots t vs displacement that is stored in results of 1st column

hold on        

plot(t,results(:,2))         % plots t vs velocity that stored in results of  2nd column

ylabel('Plot')        % labels  y axis as plot

xlabel('Time(seconds)') .    % time in second is taken as x axis

title('Damped Harmonic Motion of a Pendulum')

 % keeps the title of the plot 

legend('Angular displacement ', 'Angular velocity')

 % command that used to distinguish between angular displacement and angular velocity

ct =1;

 for i = 1:length(results(:,1)).         % used for loop 

to iterate the loop For each value in time

      x_start = 0  ;   defining starting x coordinate

       y_start = 0 ; defining the starting y cordinate

        x_end = l*sin(results(i,1)) ; % ending x cordinate

         y_end = -l*cos(results(i,1)) ;  % ending y cordinate

          % plotting the animation

           figure(2)

            plot([-1 1], [0 0], 'linewidth', 5)

             hold on

              plot([x_start x_end],[y_start y_end],'linewidth',3)

 % plotting using starting and ending cordinates providing first X coordinates and then the Y coordinates  of two points 

  hold on  % command that makes plots do not get over return

 % to close video file we used close command

       FUNCTION CODE EXPLAINATION :

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

% ode function is the name of the function and the format used is for function definition

t, theta ,b ,g , l, m are the input parameters that are passing

theta1= theta(1)   %defining theta (1) as theta1 

theta2= theta(2) % defining theta(2) as theta2 to make the second order differential equation into single order differential equation

dtheta1_dt = theta2  % differentiating theta 1 with respect to time gives theta 2

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1)

% after transforming second order differential equation into single order differential equation then dtheta2_dt becomes  { -(b/m)*theta2 - (g/l)*sin(theta1)}

dtheta_dt = [dtheta1_dt; dtheta2_dt]

% solution vector for theta is transformed into above when single order differential equation is written in matrix form

end 

 % end is a keyword that terminate function

PLOTS :

 ANIMATION LINK  :

https://youtu.be/PumPK8NTm6I

CONCLUSION :

Hence , solved the second order differential equation that represents the equation of motion of simple pendulum under damping .