Week 3 - Solving second order ODEs

clear all
close all
clc

b=0.05;   %damping coefficient
l = 1;    %length of the pendulum in m,
m = 0.1;  %mass of the ball in kg,
g=9.81;   %gravity in m/s^2,


% Initial condition
theta_0 = [0;3];  %Angular displacement in rad and angular velocity in rad/s

%time points
t_span = linspace(0,10,20);

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

% To defined pivot point of pendulum

x0 = 0; %initial point
y0 = 0; %initial point
plot(x0,y0,'Marker','o','MarkerSize',5)

% to extract the first coloum of a result

y = results(:,1); %theta/angilar displacement

% to create images for each time value and stored in it array 'M'
ct = 1;

for i = 1:length(t_span);  % i = 1,2,3,....,20
    x1 = l*sin(y(i));
    y1 = -l*cos(y(i));
    plot([x0 x1],[y0 y1],'Marker','o');
    hold on
    plot(x1 ,y1,'Marker','o','MarkerSize',25,'MarkerFaceColor','r');
    hold off
    title({'Simple Harmonic motion';'of the Pendulum'},'color','g');
    axis([-1 1 -2 0]) % axis range -1 1 are on x-axis and -2 0 0n y-axis
    pause(0.3)
    M(ct) = getframe(gcf);
    ct = ct+1;

end 
   
  movie(M);
  videofile = VideoWriter('Harmonic_motion_of_pendulum.avi','Uncomppressed AVI');
  open(videofile);
  writeVideo(videofile,M);
 close(videofile);



  %Sub-program for ode_fun

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