Simulate Transient Behaviour of simple oscillation pendulum with damping using MATLAB.
Aim : To simulate the transient behaviour of a simple oscillation pendulum using MATLAB.
Governing Equation : In Engineering, ODE is used to describe the transient behavior of a system. This ODE represents the equation of motion of a simple pendulum with damping.
`(d^2 θ)/(dt^(2 ) )+b/m d(θ/dt)+ g/l sinθ=0`
In the above equation,
g = gravity in m/s^2
L = length of the pendulum in meter
m = mass of the ball in kg
b = damping coefficient
θ = angular displacement in radians.
t = time in sec.
Fig 1. Simple Oscillation 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" with respect to time.
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.
Objective :
- Write a program to simulate the transient behavior of a simple pendulum and to create an animation of it's motion.
- Work flow/steps :
Step 1: Hence we are using MATLAB as coading language it solves only single order ODE equations, since the given ODE equation is converted into two single order ODE equations.
`(d^2 θ)/(dt^(2 ) )+b/m dθ/dt+ g/l sinθ=0`
Let us consider, `θ=θ_1`
`(d^2 θ_1)/(dt^(2 ) )+b/m (dθ_1)/dt+ g/l sinθ_1=0`
Put `( dθ_1)/dt=θ_2` -------------- Equation (1)
`(dθ_2)/(dt)+b/m θ_2+ g/l sinθ_1=0`
`(dθ_2)/(dt)=- b/m θ_2- g/l sinθ_1` --------------- Equation (2)
In matrix form,
`d/(dt) [ ■(θ_1@θ_2 )]=[■(θ_2@(- b/m θ_2- g/l sinθ_1 ) )]`
Step 2 : Writing function for above ODE using MATLAB.
Note : The function name(i.e.,odefunc) should match with file name to run a code.
function [dthetadt] = odefunc(t,theta,b,g,l,m)
theta1 = theta(1)
theta2 = theta(2)
dtheta1dt = theta2
dtheta2dt = -(b/m)*theta2 - (g/l)*sin(theta1)
dthetadt = [dtheta1dt;dtheta2dt]
end
Step 3: Write down pre-defined data & initial conditions from problem statement.
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.
Step 4: Write a program to plot a figure.
% Topic : Simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.
% Let us clear & close any previous plots or code by using these cammands.
close all
clear all
clear
% In Engineering, ODE is used to describe the transient behavior of a system.
% Now our objective is to solve 2nd Order ODE which represents a transient
% behaviour of a simple pendulum.
% As we know that the 2nd Order ODE equation for Simple pendulum with
% damping is represented by,
% Equation : (dtheta/dt)^2 + (b/m)(dtheta/dt) + (g/l)sin(theta) = 0
% 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, theta = angular displacement is in radians.
% Define terms & Note down given useful data from problem statement.
%Damping coefficient(b):
b = 0.05
%gravity(g)
g = 9.81
%length(l) of the pendulum
l = 1
% mass of the ball
m = 1
% Initial conditions:
theta_0 = [0;3]; %[angular displacement;velocity]
% time points:
t_span = linspace(0,20,500)
% Now we have to define a function to solve given ODE using new file.
% Note your function name should match with file name.
% Solve ODE: Here we are calling function
[t,results] = ode45(@(t,theta) odefunc(t,theta,b,g,l,m),t_span,theta_0);
% Subplot actually used to divide a plot axis to our convinient.
subplot(5,4,[13,14,17,18])
% Plot a figure
plot(t,results(:,1),'color','b')
hold on
plot(t,results(:,2),'color','r')
hold off
subplot(5,4,[15,16,19,20])
plot(results(:,1),results(:,2))
xlabel('position-Radians')
ylabel('angular velocity')
Step 5: Write a program to animate using getframe and movie command.
% To make animation we capture all frames using this step:
ct = 1
for i = 1:length(results(:,1))
% Initial position
x0 = 0
y0 = 0
% Position at resp angular position
x1 = l*sin(results(i,l));
y1 = -l*cos(results(i,l));
% subplot used divide a figure
subplot(5,4,[1:12]);
% For fixed Line.
plot([-1 1],[0 0],'linewidth',4,'color','k')
hold on
% For Pendulum string.
line([x0 x1],[y0 y1],'linewidth',2,'color','b')
hold on
% For ball or mass
% Use markers for shape of ball i.e, '-o'
plot(x1,y1,'-o','markers',20,'markerfacecolor','r')
hold on
grid on
axis([-1.5 1.5 -1.5 0.5])
pause(0.03) % Execution of Program Stops this interval of time.
hold off
% we are capturing figure(i.e plot) as a movie frame using this syntax.
% capturing your frames using getframe command
S(ct) = getframe(gcf)
ct = ct +1
end
% Recording our movie frame.
movie(S)
% Creating file
videofile = VideoWriter('Osc_Pend.avi','Uncompressed AVI')
% open save file
open(videofile)
% Writing/feeding movie frame in saved file.
writeVideo(videofile,S)
% closing the opened file(It is important to close any opened file)
close(videofile)
Note : You find your animated movie in your program saved file.
Body of the content :
- Using MATLAB as a programming Language.
- 01.Defining ODE function :
% Note : any 2nd Order ODE is converted into two single order ODE.
% let us consider theta = theta1, dtheta1/dt = theta2
function [dthetadt] = odefunc(t,theta,b,g,l,m)
theta1 = theta(1)
theta2 = theta(2)
dtheta1dt = theta2
dtheta2dt = -(b/m)*theta2 - (g/l)*sin(theta1)
dthetadt = [dtheta1dt;dtheta2dt]
end02. Main Program :
% Topic : Simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.
% Let us clear & close any previous plots or code by using these cammands.
close all
clear all
clear
% In Engineering, ODE is used to describe the transient behavior of a system.
% Now our objective is to solve 2nd Order ODE which represents a transient
% behaviour of a simple pendulum.
% As we know that the 2nd Order ODE equation for Simple pendulum with
% damping is represented by,
% Equation : (dtheta/dt)^2 + (b/m)(dtheta/dt) + (g/l)sin(theta) = 0
% 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, theta = angular displacement is in radians.
% Define terms & Note down given useful data from problem statement.
%Damping coefficient(b):
b = 0.05
%gravity(g)
g = 9.81
%length(l) of the pendulum
l = 1
% mass of the ball
m = 1
% Initial conditions:
theta_0 = [0;3]; %[angular displacement;velocity]
% time points:
t_span = linspace(0,20,500)
% Now we have to define a function to solve given ODE using new file.
% Note your function name should match with file name.
% Solve ODE: Here we are calling function
[t,results] = ode45(@(t,theta) odefunc(t,theta,b,g,l,m),t_span,theta_0);
% Subplot actually used to divide a plot axis to our convinient.
subplot(5,4,[13,14,17,18])
% Plot a figure
plot(t,results(:,1),'color','b')
hold on
plot(t,results(:,2),'color','r')
xlabel('Intervals of time - sec')
ylabel('Angular Velocity & displacement')
legend('Angular Displacement','Angular velocity','location','southeast')
hold off
subplot(5,4,[15,16,19,20])
plot(results(:,1),results(:,2))
xlabel('position-Radians')
ylabel('angular velocity')
% To make animation we capture all frames using this step:
ct = 1
for i = 1:length(results(:,1))
% Initial position
x0 = 0
y0 = 0
% Position at resp angular position
x1 = l*sin(results(i,l));
y1 = -l*cos(results(i,l));
% subplot used divide a figure
subplot(5,4,[1:12]);
% For fixed Line.
plot([-1 1],[0 0],'linewidth',4,'color','k')
hold on
% For Pendulum string.
line([x0 x1],[y0 y1],'linewidth',2,'color','b')
hold on
% For ball or mass
% Use markers for shape of ball i.e, '-o'
plot(x1,y1,'-o','markers',20,'markerfacecolor','r')
hold on
grid on
axis([-1.5 1.5 -1.5 0.5])
pause(0.03) % Execution of Program Stops this interval of time.
hold off
% we are capturing figure(i.e plot) as a movie frame using this syntax.
% capturing frames using getframe command
S(ct) = getframe(gcf)
ct = ct +1
end
% Recording our movie frame.
movie(S)
% Creating file
videofile = VideoWriter('Osc_Pend.avi','Uncompressed AVI')
% open save file
open(videofile)
% Writing/feeding movie frame in saved file.
writeVideo(videofile,S)
% closing the opened file(It is important to close any opened file)
close(videofile)Outputs :
- Plot shows a relation between Angular Displacement - Velocity & time.
From above plot as time passes the angular displacement – velocity goes on decreasing.
2. Plot Shows a relation between Angular Velocity & position-Radians.
From above plot we can establish a relation between angular velocity & position.
3. Plotting a animation of Transient motion of simple pendulum.
4. Output at command window.
Errors:
- Entered comment line infront of getframe command.
Action : I think this error is because of writing a comment line after a command code , I just removed a comment line and write over the command line.
Before : S(ct) = getframe(gcf) : % capturing your frames using getframe command
After : % capturing your frames using getframe command
S(ct) = getframe(gcf)
2. I found a new error but it does not effect my result anyway.
Action : How to rid over this type of error?? Ans : Actually I alsa don’t Know, But one thing to note I copied my all codes and pasted in new sheet and running once again will coming over this type of error. But I thought it is a temporary solution. Please let me know if any one know this.
Conclusion :
- This is all about how to simulate transient behavior of simple pendulum using MATLAB.
- During coading I found some errors as shown above and difficulty in using subplot command(as a beginners you also face a same problem but don’t loose hope just practice).\
Referances :
1. https://en.wikipedia.org/wiki/Pendulum
2. Skill Lync.