Week 3 - Solving second order ODEs

Aim: To write a program to solve a 2nd order ODE that represents the motion of a simple pendulum & to animate its motion.

Inputs from the question:

g = 9.81 m/s2, gravity

L = 1 metre, length of the pendulum

m = 1 kg, mass of the ball

B = 0.05, damping coefficient

Objective:

To Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Function Program:

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/1)*sin(theta1)
dtheta_dt = [dtheta1_dt ; dtheta2_dt]
end

Main Program:

% Program to simulate motion of a simple pendulum

clear all
close all
clc

% Inputs
b = 0.05;  % damping coefficient
g = 9.81;  % gravity in m/(s^2)
l = 1;     % length of the pendulum in metre
m = 1;     % mass of ball in kg

% Initial conditions  [Initial position & Initial velocity]
theta_0 = [0;3];

% Time
t_span = linspace(0,20,500);

% Solving ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

% Plots
figure(1)
plot(t,results(:,1),'color', 'b');
hold on;
plot(t,results(:,2),'color','r');
xlabel('time')
ylabel('angular velocity and angular displacement')
grid on

ct = 1;
% Pendulum Motion Animation
figure(2)
for i=1:length(t_span)
x0 = 0;
y0 = 0;
x1 = 1*sind(results(i,1));
y1 = -1*cosd(results(i,1));
plot([x0 x1], [y0 y1], [-0.5 0.5], [0 0], 'linewidth', 3);
hold on
plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','r','markeredgecolor','k')
axis(gca,'equal');
grid on
pause(0.001)
hold off

% Gif
   k = 0.5;
   l_span = [1:5];
     printf('creating GIF- Progress... %d/%dn', k, length(l_span));
     img = print ('-RGBImage');
     imwrite(img,'animation.gif','Delaytime',0.005,'Compression','bzip','WriteMode','Append');
   
ct = ct +1;
end

Results:

Plot

Pendulum position as per question

Harmonic Oscillation of Pendulum: GIF