Week 3 - Solving second order ODEs

AIM: To simulate the transient behavior of the simple pendulum. (motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.)

EQUATIONS:

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

CODE:

% Simple Pendulum's oscillations
clear all
close all
clc
% Inputs
l = 1; %length of the pendulum in m
g = 9.81;% accelaration due to gravity in m/s^2
m = 1; % mass of the pendulum bob
b = 0.05; % damping coefficient
% Simulation of the motion in 0-20sec, theta = 0; d_theta_0/dt = 3rad/s at t=0
theta_0 = [0;3];
% time points
t_span = linspace(0,20,500);
% solve ODE
[t, results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

% Variation of Position of pendulum w.r.t time
x0 = 0; % Initial position
y0 = 0;
x1 = l*sin(results(:,1)); % position at time 't'
y1 = -l*cos(results(:,1));

ct = 1;
for i = 1:length(x1)
    plot([x0,x1(i)],[y0,y1(i)],'-o','MarkerFaceColor','r')
    text(1,1,['Time =  ',num2str(t(i,1))]);
    axis([-2 2 -2 2])
    pause(0.1)
    nk(ct) = getframe(gcf);
    ct = ct+1;
end

movie(nk)
video_file = VideoWriter('simplepend.avi','Uncompressed AVI');
open(video_file);
writeVideo(video_file,nk);
close(video_file);

EXPLANATION:

Prior to writing the program a function is written inorder to break down the 2nd order ODE into two first order ODEs as given below.

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

In the main program, we assign the values for b,g,l&m; specify the inital conditions i.e., Initial displacement = 0 & Initial Velocity = 3m/s

Time points are defined (500) for 0 to 20sec.

Then we call the function to compute the result for the corresponding times.

The results are then plotted.

LINK for animation:

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