Week 3 - Solving second order ODEs

Program to solve the Ordinary differential equation and to simulate simple pendulum motion.

Program code:

clear all
close all
clc

b = 0.05; - Damping Coefficient
g = 9.81; - Gravity
l = 1; - Length
m = 1; - Mass

% initial condition
theta_0 = [0,3]; - Initial displacement and Angular velocity

% time points
t_span = linspace(0,20,300); - Travelling time of a Pendulum

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

angular_displacement = results(:,1);

ct = 1; - Counter
for i = 1:length(angular_displacement)
THETA1 = angular_displacement(i);

x0 = 0;
y0 = 0;

x1 = l*sind(THETA1);
y1 = -l*cosd(THETA1);

figure(2)

%plotting
plot([x0 x1],[y0 y1],'LineWidth',3,'color','r');
hold on
plot(x1,y1,'marker','.','markersize',100)
axis([-0.1 0.1 -1.2 0.2]);
%pause(0.03);
M(ct)=getframe(gcf);
ct = ct +1;
hold off


end

movie(M)
videofile = VideoWriter('ordinary_differential_equation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Functional code:

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

Results:

The output of the program below for your reference,

https://drive.google.com/file/d/1i9YR0snRkRfwgPBpLNP889y4uOze2d-1/view?usp=sharing