Week 3 - Solving second order ODEs

% This is the code for solving the second order ode equation of forced
% damped vibrations and the equation is
% (d^2 theta/dt^2)+((b/m)*(dtheta/dt))+((g/l)*sin(theta)), for solving the
% second order ordinary differrential equation in matlab we would first convert the
% equation into first order equation by assuming dtheta/dt = theta2
% which would then be subsituted into the original equation and hence the
% second order equation would be converted to first equation which is as
% follows dtheta2/dt + ((b/m)*theta2)+((g/l)*sin(theta)) and then array of
% two vectors would be created wherein theta2 = dtheta/dt and dtheta2/dt =
% -((b/m)*theta2)+((g/l)*sin(theta)) which is used to solve the equation
% using ode45 function of matlab which integrate the seperated arrays
% and results so obtained would be scaler quantities of angular displacement and
% angular velocity, here the array of angular displacement is used to animate the
% motion of pendulum


b = 0.05;   %coefficient of damping in Nm/s
g = 9.81;   %gravitational acceleration in m/s2
m = 1;      %mass of the object in kg      
l = 1;      %length of the object in m
r = 1;      %length of the string of the pendulum in m

%tspan is the limit used by ode45 to calculate the steps for the
%runge-kutta method or it can also be called the limits for integration
% here 0 to 20 sec is the time limit in sec as per the challenge and it is 
% furthur divided into 200 steps for better calculation and representation
tspan = linspace(0,20,200); 

%y0 is the initial condition provided for the runge-kutta method or in the
%integration of the two functions first row represents the initial  
%angular displacement and second row is for the initial velocity of
%pendulum
y0 = [0; 3]; 

%ode45 is used to solve the equation of forced damped vibrations here t is
%for the time steps provided to obtain the results at that time and results
%would contain the solution of the equation at that given time, here odefunc is
%the function created to convert second order equation to first order
%equation and b, m, g, l are the constant of the equation the t used after
%ode45 syntax shows that the equation is the function of time and theta is
%dependent on time and time is independent variable. tspan and y0 are
%limits and initial conditions respectively
[t, results] = ode45(@(t,theta) odefunc(t, theta, b, g, l, m),tspan,y0);

%for loop is used for creatig the animation of pendulum using the results

ct = 1; %counter for capturing each frame in for loop
for i = 1:length(results(:,1)) 
    
    %initial coordinates for of the string of pendulum
    x0 = 0;
    y0 = 0;
    
    %here x and y component of the pendlum are reversed and negative sign
    %is used to animate pendulum as per challenge
    
    %length of string is multiplied with cosine of the angular displacement
    y1 = -r*cos(results(i,1));

    %length of string is multiplied with sine of the angular displacement
    x1 = r*sin(results(i,1));

    %plotting
    plot([x0 x1],[y0 y1],'*',[x0 x1],[y0 y1])
    axis([-1.1 1.1 -1.1 0])
    
    pause (0.05)    %pause for the better animation
    M(ct) = getframe(gcf);  %capturing plot of each position of pendulum for animation
    ct = ct+1;  %increasing counter value for storing next frame
end

%Animatting
movie(M)    %calling the arrays of the frame


videofile = VideoWriter('pendulum.avi','Uncompressed AVI'); %creating videofile
open(videofile);                %opening the file
writeVideo(videofile,M);        %adding the frames to the file
close(videofile);               %closing and storing the file


%function for converting second order ode equation to first order ode
function [dtheta_dt] = odefunc(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