Week 3 - Solving second order ODEs

https://youtu.be/lmRT6suBM-Y

% program to plot the variation in angular displacement and angular velocity
% of a simple pendulum w.r.t time for a given mass of bob, length of string 
% and damping coefficient.

clear all
close all
clc

% The ODE governing the motion of pendulum is given as
% (d^2.theta/dt^2)+(b/m)(d.theta/dt)+(g/l)sin(theta)=0
% here b, represents the damping coefficient,
% m, represent the mass of the bob
% g represent acc. due to gravity
% l represent the length of string to which the bob is attached

b=0.35;         % assuming a damping coefficient value

g=9.81;         % g=9.81 m/s^2

l=1;            % assuming length of string as 1 m

m=1;            % assuming mass of bob as 1 kg

% intitialising counter ct, for writing movie to file
ct=1;

% cordinates(h,k) represent the position of pivot about which the pendulum
% rotates about

h=5;
k=5;

% time span for which the motion under study is given in the matrix'z'

z=linspace(0,20,500);

% initial condition array 'f' represent the conditions at time,t=0
% the colum vector represents the intial condition that is passed on to the
% ode45 function, with first row representing angular displacemt(rad) and
% the second row representing angulat velocity(rad/s)

f=[0;3];   % assuming angular displcement=0 rad and angular velocity= 3 rad/s at t=0

% solving the ODE using the ode45 function, the independent and dependent
% variables are passed on followed by the ODE to be integrated followed by
% the span of the independent variable upto which the equation is to be
% integrated followed tby the initial condition matrix

figure(1)    % opening the figure window 1 to plot the graph of t vs angular displ.
              % and t vs angular velocity
% 'results' matrix is a column vector with first column giving the values of 
% angular. displ. and second column giving values of angular velocity

% the ODE is inserted to the function as another pre-defined function as ode_func
[t,results]=ode45(@(t,theta)ode_func(theta,b,g,l,m),z,f);

% plotting the graph of t vs angular displacement

plot(t,results(:,1));

hold on

% plotting the graph of t vs angular velocity
title('t Vs angular dispacemt and t Vs angular velocity');

plot(t,results(:,2));

xlabel('t,time(s)');

legend('angular displacement(rad)','angular velocity(rad/s)');

% Assigning a string s1 to display the time counter in the plot window

s1='time(s)= ' ;   

% configuring coordinates(x2,y2) for tracing the path of the bob position
x2=h+l*sin(results(:,1))
y2=k-l*cos(results(:,1))


 for i=1:length(z)
     figure(2)
     % coordinates (x,y) reresents the position of end of string for every
     % point in time.
     
     x=h+(l*sin(results(i)));
     y=k-(l*cos(results(i)));
     
     % string s2 is assigned to display the time counter for displaying
     % time value using the command num2str
     s2=num2str(t(i));     
     
     % plotting the mean position of the pendulum
     
     plot([5 5],[5.2 4],'linewidth',0.5,'color',[0 0 0],'linestyle','--')
     hold on
     
     % Plotting the trail of movement of bob
     plot(x2,y2,'-','Markersize',0.1,'markerfacecolor','r','markeredgecolor','r');
     
     % plotting the string of the pendulum with corrsponding points (h,k)
     % and the point(x,y)
        
     plot([h x],[k y],'linewidth',1.5,'color',[1 0.1 0.8]);
     
     
     plot(x,y,'o','Markersize',20,'markerfacecolor','g','markeredgecolor',[0 0 0]);
     
     % plotting the circular bob at the end of the pendulum with a,b
     % representing the centre point abou which the circle is plotted
      
     % plotting the fixed line about which the pendulum swings
     plot([4.8 5.2],[5 5],'linewidth',3,'color',[1 0 0])
    
     % displaying the time value obtained using string
     text(5.4,5,[s1,s2])
     
     
     
     axis([4 6 3.5 5.5])
     
     legend('mean position','path traced by the bob');
     
     title('Motion of a simple pendulum')
     
     hold off
          
     S(ct)=getframe(gcf);
     ct=ct+1;
    
 end

movie(S);
Videofile=VideoWriter('pendulum.AVI');
Videofile.FrameRate=10;
open(Videofile)
writeVideo(Videofile,S)
close(Videofile)

% definition of the function ode_func
 
% function[dtheta_dt]=ode_func(theta,b,g,l,m)
% u=theta(1);  %assigning all the first column values(ang. displacement) to matrix u
% v=theta(2);  %assigning all the second column values(ang. velocity) to matrix v
% 
% dt1=v;
% dt2=-(b/m)*v-(g/l)*sin(u);
% dtheta_dt=[dt1;dt2];
% end



% errors faced
% 1.Images where getting overlapped in all frame.Solved by placing hold off command at % the end of for loop.(images attached)