Week 3 - Solving second order ODEs

Aim :- To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

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.

• Write a program in Matlab\Octave that will simulate the pendulum motion for L=1 metre,m=1 kg, b=0.05, g=9.81 m/s2, Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.

Introduction:- In Engineering, Ordinary Differential Equation is used to describe the transient behavior of a system. A simple example is a pendulum.The way the pendulum moves depends on the Newtons second law. When Newtons second law is applied, we get a second order Ordinary Differential Equation that describes the position of the "ball" with respect to time.

The transient behaviour of the pendulum can be solved using the second order differntial equation.

Script Explanation:-

• Defining given parameter

• intial condition at t=0 sec

• Defining time limit

• solving ordinary diffrential equation

• plotting the time Vs angular velocity & Displacement

• counter for frame calculation

• Defining initial position of ball

• Defining final position of ball

• Creation of frame for animation

• code to save the video file

################ Matlab code for solving Ordinary Diffrential Equation:- ################

%function for solving the ode
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
theta_1 = theta(1);  %displacement
theta_2 = theta(2); %velocity
dtheta_1_dt = theta_2; %angular velocity
dtheta_2_dt = -(b/m)*theta_2 - (g/l)*sin(theta_1); %acceleration
dtheta_dt = [dtheta_1_dt; dtheta_2_dt];
end

clear all
close all
clc
%inputs
b=0.05; % damping cofficient
length=1; %lenth of thread(m)
gravity=9.81;% m/s^2
mass=1; % mass of ball (kg)
%inital condtion
theta_0=[0;3]; 
%time points
time_span=linspace(0,20,300);
%solve ode
[t,results]= ode45(@(t,theta) ode_func(t, theta,b,gravity,length,mass),time_span,theta_0);
%plotting the result
plot(t, results(:,1),'color','g')
hold on
plot(t,results(:,2),'color','black')
xlabel('time')
ylabel('angular dispacement & velocity')
legend('angular dispacement','angular velocity')
hold off
%Animating the transient motion of thesimple pendulum
ct=1;
for i=1:length(results(:,1))
 %initial position
 x0=0;
 y0=0;
 %final postion
 x1=length*sin(results(i,1));
 y1=-length*cos(results(i,1));

 figure(2)
 plot([0 0],[-1.5 0.5])
 hold on
 plot([-1 1],[0 0],'linewidth',2); %creating a beam
 plot([x0 x1],[y0 y1],'linewidth',2,'color','b') %creating the line from initial value to final value
 axis([-1.5 1.5 -1.5 0.1])
 plot(x1, y1,'o','MarkerSize',20,'color','r','markerfacecolor','black','markeredgecolor','black') %create the ball attached at end of the rope
 hold off
 pause(0.001)
 M(ct)=getframe(gcf)

 ct=ct+1;

end
%making the video file
movie(M)
videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Output of the code :-

Animation : - After converting .Avi format file to .gif file

Conclusion:- Here, we can conclude that at t=0 sec and displacement and velocity are taking place and it is continuously decreasing with respect to time, at one particular time(at the end), it will gain the original position and the rate of decreasing of angular displacement and velocity are different.