To simulate the transient behavior of a simple pendulum and to create an animation of its motion
Aim: To simulate the transient behaviour of a simple pendulum and to create an animation of its motion
Introduction:
In Engineering, ODE is used to describe the transient behaviour of a system. A simple example is a pendulum. 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.
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.
The aim of this study to write a program in Matlab that will simulate the pendulum motion.
Use:
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, jangular velocity=3 rad/sec at time t=0.
To solve the above second order equation, we break in single order ODE equations and then we can solve them in Matlab. i.e.
theta1 = theta(1)
theta2 = theta(2)
dtheta1/dt = theta2 ……………………………………..(i)
dtheta2/dt = (-b/m)*theta2 - (g/l)*sin(theta1) ………….(ii)
Program 1: function used in ode45 is ode_func, code is written below…
% Functions used as \'ode_func\' as code 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
Program 2: whole code
clear all
close all
clc
% Initial input values
b = 0.05; % Damping cofficient
g = 9.81; % gravity in m/s^2
l = 1; % Length of pendulum in m.
m = 1; % Mass of the ball in Kg
% Initial conditions used.
% Taking Angular displacement [rad] is \'0\' & Angular velocity [rad/s] is \'3\' rad/s.
theta0 = [0;3];
% Time span between 0 to 20 seconds
time = linspace(0,20,500);
% Calculation of Angular displacement and angular velocity by using ode45
[T, Results] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), time, theta0);
% Plottting of Angular displacement and angular velocity w.r.t. Time.
figure(1)
plot(T, Results(:,1),\'linewidth\',0.75,\'color\',\'r\') % plot of angular displacement vs time
hold on
plot(T, Results(:,2),\'linewidth\',0.75,\'color\',\'b\') % plot of angular velocity vs time
xlabel(\'Time [sec]\');
ylabel(\'Angular displacement [rad] & Angular velocity [rad/s]\');
axis([0 21 -3 3]);
hold on
% Animation of pendulum
% To create the animation we used the \'for\' loop
% here \'ct\' is a counter veriable
ct = 1;
Theta = Results(:,1);
for i = 1: length(Theta)
theta1 = Theta(i)
x0 = 0;
y0 = 0;
x1 = l*sin(theta1);
y1 = -l*cos(theta1);
figure(2)
plot([-1 1],[0 0], \'linewidth\',5, \'color\', \'g\');
hold on
plot([x0 x1],[y0 y1],\'linewidth\',3, \'color\',\'b\');
hold on
% Now we one pin joint and ball of pendulum with the help of \'markers\' tool
plot( x0, y0, \'H\', \'markersize\', 10,\'MarkerEdgeColor\', \'r\', \'MarkerFaceColor\',\'r\');
plot( x1, y1, \'o\', \'markersize\', 30,\'MarkerEdgeColor\', \'r\', \'MarkerFaceColor\',\'r\');
hold off
axis([-1.5 1.5 -1.5 0.5])
pause(0.003)
M(ct) = getframe(gcf);
ct= ct+1;
end
% To make an animation by collecting each frame following commands are used
movie(M)
videofile = VideoWriter(\'pendulum.avi\',\'uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Code explanations:
- Before simulation of pendulum motion we need some initial inputs that is the value of damping coefficient (b), length of pendulum (l), mass of pendulum ball (m), and gravity (g) etc.
- Initial conditions used as Angular displacement [rad] is \'0\' & Angular velocity [rad/s] is \'3\' rad/s. and time span of motion is 0 to 20 seconds.
- For Calculation of Angular displacement and angular velocity with respect to time, we used ‘ode45’ in build function program.
- Then we plot the angular displacement and angular velocity w.r.t.time.
- We used for loop to simulate the motion of pendulum with each value of angular displacement and plotting the bar and length of pendulum . And using ‘marker’ tool to plot pendulum ball in animation.
- Then we save each frame and using ‘movie’ command to save the motion animation of pendulum.
Output:
Animation: