MOTION OF SIMPLE PENDULUM AND SOLVING SECOND ORDER ODE
AIM
To write solves the ODE. This ODE represents the equation of motion of a simple pendulum with damping.
OBJECTIVE
To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping and create a animated video of output obtain by solving this ODE.
THEORY
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on Newton's 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. Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.
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.
From the above equation, we create two ODEs’ i.e.,
Now we write a function code in which we are going to use some parameters that are going to use in the main program. We basically call this function in our program to solve our problem.
When these equations are solved, they are converted into two first-order ODEs which are then converted to matrices. These matrices can be solved in MATLAB to obtain the velocity and displacement values with respect to time.
In matrix
by solving these we will get the θ1,θ1 and θ2,θ2 values of the simple pendulum.
MAIN CODE
% Program to describe and animate transient behaviour of a simple pendulum system using ODE at Matlab.
% inputs of simple pendulum
b=0.05;
g=9.81;
l=1;
m=1;
% initial co-ordinate of the pendulum
theta_0 = [0;3];
% To solve ODE, integrate the systen ODE to uptain solution,
t_span = linspace(0,20,500);
% Solve ODE
[t,result] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
% ploting displacement and angular velocity with flow of time
figure(1)
plot(t,result(:,1),'color','r')
hold on
plot(t,result(:,2),'color','b')
xlabel('Time in seconds')
ylabel('Displacement(rad/sec)/Angular velocity(rad/sec)')
legend('Displavement','Angular velocity')
title('Simple pendulum')
% animation of simple pendulum
ct = 1;
disp_theta = result(:,1);
for i = 1:length(disp_theta)
% initial value
x0=0;
y0=0;
theta_displacement = disp_theta(i);
x1 = l*sin(theta_displacement);
y1 = l*(-cos(theta_displacement));
% ploting and animation
figure(2)
plot([-1 1],[0 0],'linewidth',3,'color','r')
hold on
line([x0 x1],[y0 y1],'linewidth',2,'color','b')
plot(x0,y0,'marker','o','markersize',5,'markerfacecolor','r')
plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','r')
axis([-1.5 1.5 -1.5 1])
pause(0.003)
hold off
M(ct) = getframe(gcf);
ct =ct+1;
end
movie(M)
videofile=VideoWriter('Simple_Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)MAIN CODE EXPLAINATION WITH CODE
% Program to describe and animate transient behaviour of a simple pendulum system using ODE at Matlab.
% 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.
% Equation= d^2theta/dt^2 + b/m*dtheta/dt + g/l*sin(theta)=0
% inputs of simple pendulum
b=0.05; % damper coefficient value of pendulum
g=9.81; % gravity force of pendulum(m/s^2)
l=1; % length of the pendulum (m)
m=1; % mass of the pendulum ball(kg)
% initial co-ordinate of the pendulum
theta_0 = [0;3]; %(rad/sec) we are working on @nd order ODE, the 1st one is displacement,2nd one is Angular velocity
% To solve ODE, integrate the systen ODE to uptain solution,
% input value of time array for simulation
t_span = linspace(0,20,500);
% Solve ODE
% [t,result]-output of time,output result having displacement,velocity array value depends time span.
% ode45- inbuilt matlab ode solver
% @(t,theta) - t,theta is main variable, ODE functions based on t, theta. Means its dependent,independent variable.
% ode_func- calling function solving equation and input values
% t_span- Time array which how far integration need to done.
% theta0 - Initial condition of pendulum corresponding to time
[t,result] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
% ploting displacement and angular velocity with flow of time
figure(1) % plot of frist figure for flow of velocity and displacement
plot(t,result(:,1),'color','r') % line flow plot of displacement depend upon time
hold on % hold all line plot in same figure
plot(t,result(:,2),'color','b') % line flow plot of Angular velocity depend upon time
xlabel('Time in seconds') % name of x axis
ylabel('Displacement(rad/sec)/Angular velocity(rad/sec)') % name of y axis
legend('Displavement','Angular velocity') % Indication plot line or to show the name of the line
title('Simple pendulum') % title of the plot
% animation of simple pendulum flow using computed value of time,displacement, angular velocity(solution vector)
% using loop to create animation
ct = 1; % ct is frist counter program, its create loop connection.
disp_theta = result(:,1);% frist row of theta is displacement theta
% loop program for pendulum animation, helps to display ball movement.
% array cant create animation, but loop can create array values animation.
for i = 1:length(disp_theta)
% initial value of pendulum, which is fixed value of one end of the pendulum
x0=0; % fixed point of one end represent x axis
y0=0; % fixed point of one end represent y axis
theta_displacement = disp_theta(i);
% compute co-orrdinate of pendulum ball point, Maybe it can be start point or other end place point of pendulum
x1 = l*sin(theta_displacement); % pendulum ball point represent by x-axis
y1 = l*(-cos(theta_displacement)); % pendulum ball point represent by y axis
% ploting and animation
figure(2) %% plot of second figure plot
plot([-1 1],[0 0],'linewidth',3,'color','r')% pendulum two end point connection line plot
hold on % hold all line plot in same figure
line([x0 x1],[y0 y1],'linewidth',2,'color','b')% line of pendulum connection between ball and fixed point
plot(x0,y0,'marker','o','markersize',5,'markerfacecolor','r') % indicating or marking fixed point marking
plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','r') % indicating or marking ball point
axis([-1.5 1.5 -1.5 1]) % axis limits
pause(0.003) % video running speed command
hold off % cancel all pendulum position in single figure
M(ct) = getframe(gcf); % access the frist figure to start animation
ct =ct+1; % adding figure one by one with frist counter program figure to create hole animation process
end % end the program
% create movie and animation code
movie(M)
videofile=VideoWriter('Simple_Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
FUNCTION CODE
% function program for ODE input and equation solver, compute solution vector
function[dtheta_dt]=ode_func(t,theta,b,g,l,m);
% theta is solution array, which have theta1,theta2 breaking equation
theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];FUNCTION CODE EXPLAINATION WITH CODE
% function program for ODE input and equation solver, compute solution vector
function[dtheta_dt]=ode_func(t,theta,b,g,l,m);% pendulum inputs
% theta is solution array, which have theta1,theta2 breaking equation
theta1 = theta(1); % value of theta1 is displacement of pendulum
theta2 = theta(2); % value of theta2 is velocity of pendulum
dtheta1_dt = theta2; % derivation breaking stand to compute solution vector
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); % breaking stand to compute solution vector
dtheta_dt = [dtheta1_dt;dtheta2_dt]; % solution vector to compute displacement and velocity array
end
OUTPUT SCREENSHOT
OUTPUT PLOT
OUTPUT PLOT ALL POINT MOVEMENT PENDULUM
OUTPUT PLOT PENDULUM BALL MOVEMENT
OUTPUT ANIMATION YOUTUBE LINK
CONCULUTION
The simulation of the transient behaviour of a simple pendulum is succesfuly created by using program in MATLAB.