Week 3 - Solving second order ODEs

%Main program
%A program to solve the second order Ordinary differential equation of a
%simple pendulum and to animate its motion.
%The following values are considered
b=0.05; %damping co-efficient
g=9.81;%acceleration due to gravity (m/s^2)
l=1; %length of the string of the pendulum(m)
m=1; %mass of the ball attached to pendulum (kg)

%the initial conditions are considered as (th_0 = theta values initially)
th_0 = [0;3];%the initial values of angular displacement(radians) and angular velocity(rad/s) respectively
%time span considered for solving is
t_span = linspace(0,20,100);%for time = 0sec to 20 sec with a 100 number of values in between

%solve the ODE which governs the motion of the simple pendulum using
%ode45 solver with function ode_pendulum defined in a separate program and
%is called .
%Here,the results in output are two columsns , first being displacement and
%other being velocity
[t,results] = ode45(@(t,th) ode_pendulum(t,th,b,l,g,m),t_span,th_0);

%Plot between time Vs angular displacement as well as angular velocity to
%understand the effect of damping on displacement and velocity based on
%their amplitude( the heights and troughs of the plots)
figure(1)
plot(t,results(:,1),'color','r')
hold on
plot(t,results(:,2),'color','g')
xlabel('time (s)')
ylabel('angular_ displacement(rad)& angular_ velocity(rad/s)')
legend('angular displacement','angular velocity')

%To plot the animation of the simple pendulum,'for loop' is used
ang_disp = results(:,1);%defining the angle which gives displacement and position of the pendulum
ct=1; % to have increment starting from 1 for the 'for loop'

for i= 1:length(ang_disp) %the length of the values considered for plotting
    T1 = ang_disp(i); %the angle used for plotting in every frame
   %the co-ordinates which generates the line and orients it as per the
   %angle to orient at the position
    x0=0;
    y0=0;

    x1= l*cos(((3*pi)/2)+T1); %the string is considered to be aligned at 270 degrees initially
    y1= l*sin(((3*pi)/2)+T1); % for convenience of understanding

   %to make the animation by plotting different frames as per the loop
   figure(2)

   plot([x0 x1],[y0 y1],'linewidth',0.5,'color','b');
   hold on
   plot(x1,y1,'o','markers',5,'MarkerFaceColor','r');

   plot([-1 1],[0 0],'linewidth',0.5,'color','r');
   hold off
   axis([-2 2 -2 2])

   title('simple pendulum animation')
   pause(0.01)
   grid on;
   M(ct)=getframe(gcf); %too capture the frames of the plots
   ct=ct+1; %for increment into next iteration
end

%To create the videoo file of the animation and to save it.
movie(M)
videofile=VideoWriter('simple_pendulum animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)



%Function program
%A program to define the ode function which governs the motion of the
%simple pendulum,it is a second order differential equation .
%the output id the differentiation of theta with respect to time and the
%inputs depend upon time(t) ,theta = th (angular displacement),damping
%co-efficient(b),length of pendulum (l),acceleration due to
%gravity(g),mass(m)

function [dth_dt] = ode_pendulum(t,th,b,l,g,m)
th1= th(1);
th2= th(2); %th(2) = theta(2) =d(theta(1))/dt
dth1_dt = th2 ;
dth2_dt = -(b/m)*(th2) - (g/l)*sin(th1);%d(theta(2))/dt is the given by that formula
dth_dt = [dth1_dt ; dth2_dt];%the output results , in two columns
end