Week 3 - Solving second order ODEs

• AIM: Writing a program that can simulate the simple pendulum motion.
• A pendulum is a weight suspended from a pivot so that it can swing freely.
• When pendulum is displaced from its equilibrium position, due to gravitational force it tends to achieve the equilibrium position again. And oscillating around its equilibrium position.
• Due to damping the amplitude of the oscillation gradually decreases.

• Motion of the simple pendulum can be represented by following Ordinary Differential Equation

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/L sinθ=0`

Where, b = Damping coefficient

m = Mass

g = Gravitational Acceleration

θ = Angular displacement

• Program

clear all
close all
clc

%Program for simulating pendulum motion
%Input
%Damping Coefficient
b = 0.05;

%Gravitationa Acceleration(mm/sec^2)
g = 9.81;

%Length of pendulum(metre)
l = 1;

%Weight attached(kg)
m = 1;
 
%Initial condition
%theta = [angular displacement(radian);angular velocity(rad/sec)]
theta = [0;3];

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

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
title('Velocity & Angular Displacement Plot')
ylabel('plot')
xlabel('time')
hold off

k=1;
for i=1:length(results)
   x1 = 0;
   y1 = 0;
   
   %Rotating axis by 90 degrees for simulation
   THETA = results(i,1) - (pi/2);
   
   %Calculating coordinated for position of pendulum
   x2 = l*cos(THETA);
   y2 = l*sin(THETA);
   
  
   figure(2)
   %Plotting the string
   plot([x1 x2],[y1 y2],'linewidth',2)
   
   hold on
   
   %Plotting the fixed wall
   plot([-1 1],[0 0],'linewidth',3)
   
   %Plotting the mass
   plot(x2,y2 ,'-o','markers',20)
   hold off
   
   axis([-1 1 -1.5 0.25]);
   pause(0.03)
   
   %Storing frame for simulation
   F(k) = getframe(gcf);
   k = k + 1;
   
end

%Creating and saving the video file
movie(F)
videofile = VideoWriter('Pendulum_Sim.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,F)
close(videofile)

Link for the simulation..

https://www.youtube.com/watch?v=Vou1ylSDOUA

• Problem Faced
• For calculating the position of the pendulum used following code

x2 = l*cosd(THETA);

y2 = l*sind(THETA);

• In above code the value of THETA is considered in degree
• But the in put angular velocity is given in rad/sec.
• Hence, changed above code to following code in which the value of THETA is considered in radian.

x2 = l*cos(THETA);

y2 = l*sin(THETA);