Week 3 - Solving second order ODEs

%program for simple pendulum using ode45.
clear all
close all
clc

%input variables for the secound order ode of a simple pendulum
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

theta_0 = [0; 3];

t_span = linspace(0,20,500);

%ode 45 function to solve second order ode.
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span, theta_0 );

%plotting of angular velocity and displacement
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
grid on
title('motion of pendulum')
xlabel('time in sec')
ylabel('angular velocity(rad/s)/displacement(rad)')
legend('displacement','angular velocity')

angle = results(:,1);

%count command for iterations.
ct = 1;

for i = 1 : length(angle)
    omega = angle(i);
    figure(2)
    title('motion of pendulum')
    
    %fixed corordinates of one end of simple pendulum.
    x0 = 0;
    y0 = 0;
    
    %moving coordinates of the bob.
    x1 = l*cos((3*pi/2)+omega);
    y1 = l*sin((3*pi/2)+omega);
    
    %plotting of the bob motion.
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g');
    hold on;
    plot(x1,y1,'o','markers',15,'MarkerFaceColor','k');
    hold on;
    plot([-3 3],[0 0],'linewidth',2,'color','b');
    hold off
    axis([-3 3 -2 2])
    grid on 
    pause(0.1)
    M(ct) = getframe(gcf);
    ct = ct+1;
end

%to create a animation of the simple pendulum
movie(M)
videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ode function program:

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end