Week 3 - Solving second order ODEs

clear all
close all
clc
% A programme to simulate the motion of a pendulum using 2nd order ode.

% inputs
b = 0.05; %..........damping coefficient
g = 9.81; %..........gravitational acceleration(m/s^2)
l = 1;    %..........length of string(metre)
m = 1;    %..........mass of bob(kg)

% initial conditions for angular displacement and angular velocity of-
% pendulum.
% initial angular displacement = 0 rad
% initial angular velocity = 3 rad/s
theta_0 = [0;3];

% time points
% using linspace command the time span arrey is between 0 to 20 sec with
% 500 time instants.
t_span = linspace(0,20,500);

% solving ode using 'ode45' inbuilt solver 
% the contents in 'ode45' parantheses are such that,
% '@(t,theta)' defines that our ode is a function of 't' and 'theta'
% 'ode_func' is a seperate function that solves two 1st order ode.  
% also shows how ode function is define.
% 't_span' provides initial conditions and time arrey for making -
% integration 
% 'theta_0' provides initial condition corrosponding to starting time 
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);


% plotting the results for displacement and velocity 

figure(1); %................assigning seperate figure window

plot(t,results(:,1),'linewidth',2); %.....plotting 't' vs 'displacement'

hold on %...................for ploting next result on same frame 

plot(t,results(:,2),'linewidth',2); %.....plotting 't' vs 'velocity'

title('displacement and velocity damping curves');%..labeling title for
                                                    % the plot
legend('angular displacement','angular velocity');
                               %.....indicating displacement & velocity
                                    % curves
xlabel('time(t) in sec') %..........indicating x axis as time scale


% Animation

x0 = 0; %...................co-ordinates of fixed point of pendulum
y0 = 0;

ct = 1; %...................counter represents each and every frame

for i = 1:length(t) %.........'i' is a loop counter that contains time 
                              % arrey with 500 values and the code inside
                              % the forloop will exicute 500 times
    
   x1 = l*sin(results(i,1)); %....creating x co-ordinet arrey for different 
                                 % positions of swinging pendulum
   
   y1 = -l*cos(results(i,1)); %...creating y co-ordinet arrey for different 
                                 % positions of swinging pendulum
   
   
   figure(2) %...................assigning new figure window
   
   plot([-1 1],[0 0],'linewidth',3,'color','g');
                                   %......plotting horizontal line as fix
                                          % support beam
   
   axis([-2 2 -2 2]);%....................providing x & y limits
   
   line([x0 x1],[y0 y1],'linewidth',2,'color','r');
                              %......plotting pendulum string
   
   hold on;%........................to plot next results on same plot
   
   plot(x0,y0,'O','markersize',5,'markerfacecolor','y');
                               %......plotting pivot point of fixed end
                                     % of string
   
   plot(x1,y1,'O','markersize',20,'markerfacecolor','b');
                               %......plotting swinging bob
   
   title('damping of simple pendulum');%.....labeling title for the plot
  

   hold off %..........restricting to plot previous position after plotting
                      % each position
  
   M(ct) = getframe(gcf);%.....using 'getframe(gcf)' command storing all 
                              % frames into an arrey M.
   ct = ct + 1;%...............incrementing the counter to get all frames
                              % in sequence.

end

movie(M)%..........making a movie of the arrey 'M' in which all frames are
                  % stored

videofile = VideoWriter('simple_pendulum','Uncompressed AVI');
                        %.....using 'VideoWriter' command the video title
                             % and video format is given to videofile

open(videofile);%............initiating video

writeVideo(videofile,M);%....playing and saving the video

close(videofile);%...........closes video file

% the following function 'ode_func' is about solving two 1st order ode  
% which were obtain after converting 2nd order ode into two 1st order
% ode.

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

theta1 = theta(1);%.......assigning first value of solution arrey[theta]
theta2 = theta(2);%.......assigning second value of solution arrey

dtheta1_dt = theta2;%.....solving 1st ode

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);%.....solving 2nd ode

dtheta_dt = [dtheta1_dt; dtheta2_dt];%......combining two ode into a 
                                           % single arrey
end

Errors occured while programing:

this error was about invalid expression in videowriter command due to missing inverted commas for the

term Uncompressed AVI.

after writing this term in inverted commas, error get terminate.

Youtube link for animation video of simulation of simple pendulum:

https://youtu.be/cuE-Oyioj9c