Week 3 - Solving second order ODEs

AIM:- To simulate the transient behavior of a simple pendulum and to create an animation on its motion on MATLAB.

OBJECTIVE:-

• To write a program that solves the second-order ODE corresponding to the motion of a simple pendulum and to plot its angular displacement and angular velocity with respect to time.

Where, g= gravity (m/s)

l= length (m)

b= damping coefficient (N-m/s)

m= mass (kg)

• To create the animation of the motion of pedulm.

INTRODUCTION:-

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing, and a right swing is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

Function Code

function [dtheta_dt] = ode_function(t, theta, b,g,l,m)

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

Explanation

• In this function program of the second ODE is introduced by 2 elements theta1 and theta3.
• Then the theta1 and theta2 are given by dtheta1_dt and dtheta_dt to a column matrix.

Main code of 2nd ODE

%Second order ODE
%Input 
b= 0.05; %Damping coeff
g= 9.81; %Acceleration due to gravity
l= 1; %Length
m=0.1; %mass

% Initial condition 
theta_0 = [0;5]; 

% Time points 
t_span = linspace(0,20,200); 
% Solving ODE ................... 
[t, results] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m), t_span, theta_0); 

% Plotting ....................
figure(1) 
plot(t,results(:,1)) 
hold on 
plot(t,results(:,2)) 
xlabel('Time') 
ylabel('Velocity & Displacement') 

% Animation ................... 
ct = 1; 
for i = 1:length(results(:,1)) 
x0 = 0; 
y0 = 0; 
x1 = l*sin(results(i,1)+pi); 
y1 = l*cos(results(i,1)+pi); 
figure(2) 

plot([x0 x1],[y0 y1],'-o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','linewidth',4)
axis([-2 2 -2 0]) 
hold on 
plot([-1 1],[0 0],'linewidth',5,'color','[0 0.3 0]') 
M(ct) = getframe(gcf);
hold off 
ct = ct + 1; 
end 

% Movie & Videofile ........................ 
movie(M) 
videofile=VideoWriter('Pendulum_Challenge.avi','Uncompressed AVI') 
open(videofile) 
writeVideo(videofile,M) 
close(videofile)

Explanation of code

• First, all the parameters b,g,l,m are introduced.
• Then initial conditions are given and the time span is given by linspace command.
• Using ode45 command the ode function is introduced.
• Then using the plot command the velocity & displacement in the y-axis and time is the x-axis.
• Then using the movie command the animation is saved in avi format.

Plotting of graph:-

The link of Pendulum 2nd order ode animation is given bellow

Pendulum animation