Motion of a Damped Simple Pendulum using ode45 Solver in MaATLAB
AIM: To write a program in Matlab to simulate the transient behavior of a simple pendulum and to create an animation of its motion.
THEORY:
The Ordinary Differential Equation (ODE) which represents the equation of motion of a Simple Pendulum with damping is given below :
where,
b=damping coefficient,
m = mass of the ball (in kg),
g = acceleration due to gravity (in m/s^2),
L = length of the pendulum (in m).
Assuming,
l=1 metre; m=1 kg; b=0.05; g=9.81 m/s^2 and Simulating the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.
PROCEDURE:
- In order to solve this 2nd order Order Differential Equation(ODE's), we first convert this into two 1st order ODE's. The results are stored in 2 columns, the angular displacement and the angular velocity values are stored in the 1st and the 2nd column respectively.
- We call this function using the Matlab predefined command "ode45".
- We then input the initial data such as time, angular displacement and angular velocity along with the parameters l,g,m, and b.
- We specify the time interval using the linspace command.
- We can generate the graph separately for both the parameters with respect to time using plot command.
- We use the 'for' loop to plot every position of the pendulum for the specified time interval. Whenever the 'for' loop increments the value of the pointer, we obtain the successive displacement value of the pendulum.
- We plot [x y] and assign a marker as per requirements to represent the bob of the pendulum. This plot is stored in an array of the counter, and then the counter is incremented and the loop is closed.
- We capture the frames as captured using a getframe command and stored in a variable called M(ct). We can take the output of these frames from M using the movie command and videowriter command. We also specify the required video format and denote the file name.
CODE:
1. Main Code
close all
clear all
clc
%Inputs
l = 1; %Length
b = 0.05; %Damping coefficient
g = 9.81; %Acceleration due to gravity
m = 1; %Mass of the pendulum
%Counter Variable
ct = 1;
%Time Interval
t_span = linspace(0,20,500);
%Initial Conditions
theta_0 = [0,3];
%Calling the ode function
[t,Results]=ode45(@(t,theta) ode_pendulum(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
hold on;
plot(t,Results(:,1),'color','red');
plot(t,Results(:,2),'color','blue');
hold off;
xlabel('time(s)');
ylabel('Plot');
legend('Angular Displacement','Angular Velocity')
title('Angular Displacement and Angular Velocity')
axis([0 20 -5 5])
%Animation
%Creating Plots which shows position of Pendulum at different times
figure(2)
theta_m = Results(:,1);
for i = 1:length(theta_m)
theta_y = theta_m(i);
x0 = 1;
y0 = 3;
x1 = x0 + 3*sin(theta_y);
y1 = y0 - 3*cos(theta_y);
%plotting the pendulum
plot([-1 3],[3 3],'color','r','linewidth',4);
hold on
plot([x0 x1],[y0 y1],'color','black','linewidth',3);
hold on
plot(x1,y1,'o','markers',20,'markerfacecolor','b')
hold off
grid on
axis([-3 5 -3 4])
title('Simple Pendulum1')
pause(0.0012)
M(ct) = getframe(gcf);
ct = ct + 1;
end
movie(M)
videofile = VideoWriter('Simple Pendulum1.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
2. ODE FUNCTION
function [dtheta_dt] = ode_pendulum(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
RESULTS:
(1)
(2) Animation
Youtube Link:
https://www.youtube.com/watch?v=KA4k8WsKgkQ&feature=youtu.be
ERRORS:
Error in calling the "ode45" function. We have to make sure we use the same function as to how we saved in the first function.
Error in plotting the pendulum. The mistake was plot([0 3],[0 3],'color','r','linewidth',4).
Correction: plot([-1 3],[3 3],'color','r','linewidth',4);