Week 3 - Solving second order ODEs

Damped Oscillations of Pendulum:

The general differential equation which governs the behaviour of system is given as:

d^2 (θ)/d(t^2)+ (b/m)*(dθ/dt)+ (g/l)*sinθ=0

In above equaton, b is daming coefficient, m is the mass of pendulum and l is length of string while θ is angular dispalcement.

It is clearly observed that the above equation is 2nd order differential equation.

First we have to convert it into 2 sets of 1st order differential equations.

Let θ=θ1;

Hence dθ/dt=d(θ1)/dt =θ2............... 1st equation

And d^2(θ)/dt=d^2(θ1)/dt=d(θ2)/dt

Hence d(θ2)/dt=-(b/m)*θ2-(g/l)*sinθ1 ............. 2nd equation

d(θ1)/dt= θ2

d(θ2)/dt= -(b/m)*θ2-(g/l)*sinθ1

In matrix form, the above system of equation can be written as,

d/dt(θ')=[θ2;-(b/m)*θ2-(g/l)*sinθ1]

The code is given below:

clc
clear all
close all

b=0.05;               % dampin coefficient in Kgm/s
m=1;                  % mass of pendulum ball in Kg
l=1;                  % legth of pendulum string in m
g=9.81;               % acceleration due to gravity in m/s^2

%initial conditions
theta_0=[0;3]         % initial displacement and velocity of pendulum

% time array
t_span=linspace(0,20,200)

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

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('angular displacement and velocity')
xlabel('time')
legend('Angular displacemnt in radians','Angular velocity in Rad/s')
title('Angular Displacement and Angular Velocity of Damped Oscillations')
hold off

% Animation to show motion of pendulum:
counter=1;
for i=1:length(t)
    x0=0;            % x0 is coordinate of fixed end of pendulum string
    y0=0;            % y0 is coordinate of fixed end of pendulum string
    l=1;             % length of pendulum string in metre
    x1=l*sin(results(i,1));
    y1=(-1)*l*cos(results(i,1));
    
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',2,'color','y')     % Draw the line which represents pendulum string
    hold on
    plot([2 -2],[0 0],'color','k')                      % Draw axes and quadrant lines
    plot(x1,y1,'o','MarkerSize',35,'color','r')         % Show the pendulum ball using marker
    axis([-1 1 -1 1])
    hold off
    
    M(counter)=getframe(gcf);
    counter=counter+1;
end
movie(M)
videofile=VideoWriter('Pendulum_Animation.avi','uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)
    

function[dtheta_dt]=ode_func(t,theta,b,m,l,g)
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

For animation, please visit following link:

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