A program to simulate the transient behavior of a simple pendulum using ODE
Objectitve
A program to stimulate transient behavior of simple pendulum usnig matlab program with the help of ODE and create a animation movie showing the motion of a simple pendulum
Introduction
fig - Simple Pendulum
A simple example of ODE is simple pendulum which is shown above in figure. It is pivoted at some point with bob attached to it by string. When this pendulum is provided with some external force it oscilate making some angle callled as θ.
Let,
g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg
b=damping coefficient.
We have ODE for second order as
`d^2θ/dt^2+b/m⋅dθ/dt+g/lsinθ=0`
To solve this differential equation we need to break it down into single ODE.
Let say that `θ = θ1`
Differentiating it to 1st order we get
`( dθ)/dt=(dθ1)/dt=θ2`------`eq (1)`
Differentiating again to 2nd order we get
`(d^2θ)/dt^2=(d^2θ1)/dt^2=d/dt((dθ1)/dt)=(dθ2)/dt`-----`eq (2)`
Substituting the values of eq(1) & (2) in original ODE
`(d^2θ)/dt^2+b/m⋅(dθ)/dt+g/lsinθ1=0`
`(d^2θ1)/dt^2+b/m⋅(dθ1)/dt+g/lsinθ=0`
`(d^2θ2)/dt^2+b/m⋅(dθ2)/dt+g/lsinθ=0`
`(dθ1)/dt=θ2` -----`ODE(1)`
`(dθ2)/dt=-b/m*(θ2)+g/l*sinθ1` ------`ODE(2)`
`d/dt[θ]=[θ1:θ2]`
Now this is a solution vector because it contains two quantity `θ2` & `b/m(θ2) - g/lsinθ`. If we solve this two ODE using matlab we can get the result for angular displacement and angular velocity.
Step wise Description of Matlab Programming
- First we define the input parameter as
b = 0.05 =Damping coefficient in m
g = 9.81 = Gravity in m/s^2
L = 1 = Length in m
m = 1 = Mass in Kg
- Initial conditions in term of theta_0 is taken as
Initial Displacement = 0 m
Initial angular velocity is 3 rad/s
- Time Span is 20 sec
- As we know to solve first order ODE, matlab has inbuilt solver which solves nonstiff differential equations — medium order method with two input variables.
[time,results] = ode45(@(time,theta) ode_func(time,theta)b,g,l,m),t_span,theta_0), where t_span[0 20], integrates the system of differential equations theta = f(time,results) from 0 to 20 with initial conditions theta_0.
@
Uses:
Function handle construction and reference
Description: The @ symbol forms a handle to either the named function that follows the anonymous function that follows the @ sign.
- To solve the differential equal we need to call the function which buit in matlab as
ode_func(time,theta,b,g,L,m),t_span,theta_0)
- After computing the results we have to plot the graph for angular displacement and angular velocity for the reults using plot command.
- To make animation in matlab we need to use 'for' loop and with the help of counter option we can create annimation.
- Explanation for for loop -
For i=1:length(results(:,1)) This loop will run for the length of result i.e, all row and 1 column from the result of displacement]
Displacement =results(:,1) Displacement is defined as the result of all row and 1 column
Position_of_pendulum = Displacement(i) Position of pendulum is defined by Displacement of result.
- To define the position of the bob we take
x0 = 0
y0 = 0
x1 = l*sin(Position_of_pendulum)
y1 = -l*cos(Position_of_pendulum)
x1 & y1 are defined by resolving the pendulum as below
The displacement is the arc length s. We see from figure that the net force on the bob is tangent to the arc and equals -mgsinθ. The weight mg has the components mgcosθ along the string and mgsinθ tangent to the arc. Tension in the string exactly cancels the components mgcosθ parallel to the string. This leaves a net restoring force back towards the equilibrium position at θ = 0.
The displacement 's' is directly proportional to θ. When θ is expressed in radians, the arc length is a cicrcle is to its radius(L in the instance) by this
s = Lθ
Here s is nothing but the x and y component
- Now to plot the animantion we use 'Plot' command and create horizontal line first. After it we create vertical line (string) using second plot in same plot with the help of 'hold on' command.
- Now to make the bob, 'Marker' command is used by which we can represent the mass attached to string and compute 'hold off' to close the plot window.
- movie(M) - This is used to make movie of M
videofile = VideoWriter('Simple_Pendulum_Damped.avi','Uncompressed AVI') - To create a videofile we call videofile command.
open(videofile) - We use this command to open video file to write the video
writeVideo(videofile,M) - To wiite the video file of movie name M
close(videofile) -Finally to close video file after writing the video this command is used
Matlab programming
1) Program when damping coefficient (b) = 0.05
%A program to simulate the transient behaviour of a simple pendulum
clear all
close all
clc
%Defining inputs parameters
b = 0.05; %Damping coefficient in m
g = 9.81; %Gravity in m/s^2
L = 1; %Length in m
m = 1; %Mass in Kg
%Initial conditions
theta_0 = [0;3] %Initial Displacement = 0 m & Initial angular velocity is 3 rad/s
%Time points
t_span = linspace(0,20,500); %time span is 0 to 20 sec and lispace will divide the time into 500 iterations
%Solve ODE
[t, results] = ode45(@(time, theta) ode_func(time,theta,b,g,L,m),t_span,theta_0);
%ode45 is inbuit solver in matlab that will solve ode with 2 variable only and it will integrate the differential eqaution i.e, theta with respect to t
%The @ symbol forms a handle to either the named function that follows the @ sign
figure(1)
plot(t,results(:,1),'linewidth',2,'color','blue')
hold on %hold for plotting another plot in same plot
plot(t,results(:,2),'linewidth',2,'color','red')
ylabel ('Plot') %Label y axis
xlabel('Time') %Label x axis
legend('Angular Displacement','Angular velocity')
ct=1
for i=1:length(results(:,1)) %This loop will run for the length of result i.e, all row and 1 column from the result
Displacement =results(:,1)
Position_of_pendulum = Displacement(i) %Position of pendulum is nothing but theta
x0 = 0;
y0 = 0;
x1 = L*sin(Position_of_pendulum);
y1 = -L*cos(Position_of_pendulum);
figure(2)
plot([-1 1],[0 0],'linewidth',5,'color','r') %Plotting for horizontal line in the axis plane where the pendulum will be attached
hold on %To hold the current for computing another plot in the same plot
plot([x0 x1],[y0 y1],'linewidth',3,'color','b') %Plotting for string
hold on %To hold the current for computing another plot in the same plot
plot(x1,y1,'-o','MarkerFaceColor','g','MarkerSize',20) %Plotting for the bob attached to the string
axis([-1 1 -1.5 0]) %Axis in which this annimation will run
M(ct) = getframe(gcf); %To capture all the frame in single animation
ct=ct+1; %When the first frame get capture in M(ct) counter will get increase by 1 and will run till it complete i = length(t)
hold off %To hold off the plot
end
movie(M)
videofile = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI') %Command to make videofile with file name Simple_Pendulum with extension .avi
open(videofile) %open video file to write the video
writeVideo(videofile,M) %Wiite the video file of movie name M
close(videofile) %close video file after writing the videoWe get the output reuslt as Angular Displacement and Angular Velocity which are not damped as show in figure at damping coefficient (b) = 0.05
fig - When damping coefficient (b) = 0.05
Youtube link for annimation of video -
Nature of pendulum when the damping coefficient (b ) = 0.05
2) Program when damping coefficient (b) = 0.08
%A program to simulate the transient behaviour of a simple pendulum
clear all
close all
clc
%Defining inputs parameters
b = 0.8; %Damping coefficient in m
g = 9.81; %Gravity in m/s^2
L = 1; %Length in m
m = 1; %Mass in Kg
%Initial conditions
theta_0 = [0;3] %Initial Displacement = 0 m & Initial angular velocity is 3 rad/s
%Time points
t_span = linspace(0,20,500); %time span is 0 to 20 sec and lispace will divide the time into 500 iterations
%Solve ODE
[t, results] = ode45(@(time, theta) ode_func(time,theta,b,g,L,m),t_span,theta_0);
%ode45 is inbuit solver in matlab that will solve ode with 2 variable only and it will integrate the differential eqaution i.e, theta with respect to t
%The @ symbol forms a handle to either the named function that follows the @ sign
figure(1)
plot(t,results(:,1),'linewidth',2,'color','blue')
hold on %hold for plotting another plot in same plot
plot(t,results(:,2),'linewidth',2,'color','red')
ylabel ('Plot') %Label y axis
xlabel('Time') %Label x axis
legend('Angular Displacement','Angular velocity')
ct=1
for i=1:length(results(:,1)) %This loop will run for the length of result i.e, all row and 1 column from the result
Displacement =results(:,1)
Position_of_pendulum = Displacement(i) %Position of pendulum is nothing but theta
x0 = 0;
y0 = 0;
x1 = L*sin(Position_of_pendulum);
y1 = -L*cos(Position_of_pendulum);
figure(2)
plot([-1 1],[0 0],'linewidth',5,'color','r') %Plotting for horizontal line in the axis plane where the pendulum will be attached
hold on %To hold the current for computing another plot in the same plot
plot([x0 x1],[y0 y1],'linewidth',3,'color','b') %Plotting for string
hold on %To hold the current for computing another plot in the same plot
plot(x1,y1,'-o','MarkerFaceColor','g','MarkerSize',20) %Plotting for the bob attached to the string
axis([-1 1 -1.5 0]) %Axis in which this annimation will run
M(ct) = getframe(gcf); %To capture all the frame in single animation
ct=ct+1; %When the first frame get capture in M(ct) counter will get increase by 1 and will run till it complete i = length(t)
hold off %To hold off the plot
end
movie(M)
videofile = VideoWriter('Simple_Pendulum_Damped.avi','Uncompressed AVI') %Command to make videofile with file name Simple_Pendulum with extension .avi
open(videofile) %open video file to write the video
writeVideo(videofile,M) %Wiite the video file of movie name M
close(videofile) %close video file after writing the videoWe get the output reuslt as Angular Displacement and Angular Velocity which are damped after increasing th value of daming coefficient as show in figure for damping coefficient (b) = 0.8
fig - When damping coefficient (b) = 0.8
Youtube link for annimation of video -
Nature of pendulum when the damping coefficient (b ) = 0.8
Error faced during running program
1) Parenthesses Error
fig - Parenthesses Error
Forgoted to add close the bracket. ) added to resolve the error.
2) Variable 'l' not d=defined
fig - Variable 'l' not d=defined
Conclusion
By executing the matlab program for simple pendulum we can say that damping coefficient is the factor which control the damping of pendulum as shown in the plot. As soon as we increase the value of damping coefficient from 0.5 to 0.8 the time taken to damped the motion of pendulum decreases which can be seen in the youtube video