Week 3 - Solving second order ODEs

AIM -

Aim of this Challange is to animate the movement of pendulam by solving the equation of motion of a simple pendulum with damping by using MATLAB.

OBJECTIVE -

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

Write a program in MATLAB that will simulate the pendulum motion, then use the information to animate the motion of the pendulum.

THEORY -

Simple Pendulum -

A simple pendulum device is represented as the point mass attached to a light inextensible string and suspended from a fixed support. A simple pendulum shows periodic motion, and it occurs in the vertical plane and is mainly driven by the gravitational force.

It is interesting to note that the time period of the real simple pendulum remains constant even if the amplitude is changed but if the acceleration due to gravity changes the time period of the simple pendulum also changes. This property is also used to determine or verify acceleration due to gravity. The time period also depends on the length of the string to which the Bob of the simple pendulum is attached. 

MATLAB CODE -

1. Function for Solution

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

Explanation -

Creating a function to get the value for `(d theta) /dt` by solving first and secound ODEs. Using the given values for

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

2. Final CODE -

clear all
close all
clc

%Given Inputs 
b = 0.05;  %c is the damping coefficient
g = 9.81;  %acceleration due to gravity
l = 1;     %length of the pendulum
m = 0.1;   %mass of the bob

%intial condition for displacement
theta_0 =[0;3];

%time variation in the motion of pendulum
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);
plot(t,results(:,1),'color','r')
hold on
plot(t,results(:,2),'color','b')
hold on
subplot(5,4,[15,16,19,20]);
plot(results(:,1),results(:,2))
axis([-1 1 -4 4])
xlabel('position(rad)');
ylabel('Angular Velocity');

ct = 1 ;

%running loop for different values of theta
for i= 1:length(results(:,1))
  x0=0;
  y0=0;

  x1 = l*sin(results(i,1));
  y1 = -1*cos(results(i,1));
  
  subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
  plot([1 1],[0 0],'LineWidth',4,'color','g');
  hold on
  line([x0 x1] , [y0 y1],'linewidth',2,'color','b');
  hold on
  %marker is used to create bob
  plot(x1,y1,'o','markers',20,'markerfacecolor','g');
  
  grid on;
  axis([-1.5 1.5 -1.5 0.5])
  pause(0.03);
  hold off
  M(ct) = getframe(gcf);
  ct = ct+1;
end
movie(M)
videofile=VideoWriter('Pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Explanation -

• First giving the inputs for g= gravity in m/s2, l = length of the pendulum in m, m = mass of the ball in kg, b=damping coefficient .Intializing condition for displacement, time variation in the motion of pendulum
• Then solving the ODE and ploting the time plot
• Running a loop for different values for `theta`
• Inside the loop describing the position for wire and bob
• At the end making the animation for pendulum

Results -

In the graphical plot and the animation we can see that the simulayion oscilation motion is in between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

We see that the Angular displacement is damped slightly due to the lower value of damping coefficient (0.05).

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