Week 3 - Solving second order ODEs

Here is the code for plotting Time vs Angulardisplacement and Angularvelocity

All Inputs b=0.05(damping coefficient), g=9.81m/s2(gravity), l=1m(length of the pendulum), m=1kg(mass).

close all
clear all
clc
b = .05 %Damping coefficient
g = 9.81 %Gravity(m/s2)
l = 1 %Length of the pendulum(m)
m = 1 %Mass of the pendulum(kg)

% Initial condition
theta_0 = [0;3]; %Initial displacement = 0, Initial velocity = 3

% Time Points
t_span = linspace(0,20,500); %Time array with 500 values from 0 to 20 sec 

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


plot(t,results(:,1)) %Time(sec) vs Angular displacement(rad)
hold on
plot(t,results(:,2)) %Time(sec) vs Angular velocity(rad/s)
xlabel('time')
ylabel('plot')
legend('Angular displacement','Angular velocity') % Blue = Ang.disp, red = Ang.vel

figure(2)
x0 = 0; % (x0,y0) is one end of the pendulum which is fixed (ORIGIN)
y0 = 0;
ct = 1;
for i = 1:length(t_span) % total no.of values depends on t_span length
    x1 = l*sin(results(i,1)); % x-coordinate of the other end 
    y1 = -l*cos(results(i,1)); % y-coordinate of the other end
    plot([-1 1],[0 0],'linewidth',6,'Color','g'); % Support rod(Fixed line)
    axis([-2 2 -2 2]); % X-axis & Y-axis limits
    line([x0 x1],[y0 y1],'linewidth',2,'color','b'); % Pendulum length
    hold on;
    plot(x0,y0,'O','MarkerSize',3,'markerfacecolor','y'); % Fixed point
    plot(x1,y1,'O','MarkerSize',20,'markerfacecolor','r'); % Fixed point
    title('Simple Pendulum')
    legend('support rod','length of pendulum','fixed point','body of pendulum');
    grid on;
    hold off;
    M(ct) = getframe(gcf); % It captures all figures one by one
    ct = ct+1;
end
movie(M);
videofile = VideoWriter('Simple Pendulum','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile) 

x1 = L times sin(theta) and y1 = -L times cos(theta)

Legend command is used for labelling.

In every figure there is line between (0,0) and (x1,y1),these are all joined together by getframe command.

Here is the ode function where t,theta,b,g,l,m are inputs.

theta1 is first element of theta and theta2 is second element of theta array.

dtheta1_dt = theta2.

dtheta2_dt is defined from the second order differential eqation given in the question.

dtheta_dt = (dtheta1_dt dtheta2_dt)'

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

Here is the youtube link for animation video

https://youtu.be/EGKRrVzdC5Q