SOLVING SECOND ORDER ODE TO SIMULATE TRANSIENT BEHAVIOR OF SIMPLE PENDULUM USING MATLAB
AIM: Simulation of simple Pendulum motion by solving 2nd order ODE.
OBJECTIVES:
1) Write a program in MATLAB to solve second order ODE.
2) Write a program in MATLAB to simulate Simple Pendulum.
3) Plotting of Results.
4) Create an animation file of Simulation.
THEORY:
The simple Pendulum is an idealized mathematical model of a pendulum. This is a weight (or bob) on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it will swing back and forth at a constant amplitude. a model of this type of pendulum is shown in the figure below.
The motion of these pendulum is a function of 2nd Order ODE and can be solved by the equation as given below:
`(d^2theta)/dt^2 + (b/m)(d theta)/dt + (g/L)sintheta`
In the above equation,
g = gravity constant in m/s^2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
In order to solve the equation it must be converted in the first order ODE, which is shown in the following formulas,
BODY OF CONTENT:
PROGRAM:
%Program to solve a 2nd order ODE and to simulate a Pendulum
clear all
close all
clc
%Inputs
%Gravity constant
g = 9.81 %m/s^2
%Length of the Pendulum
L = 1 %m
%Mass of the ball
m = 1 %kg
%Damping co-efficient
b = 0.05
%Initial condition
theta_0 = [0,3]; %[Angular Displacement, Angular Velocity]
%Time span
t_span = linspace(0,20,500);
%Calling Function to solve ODE
[t, outputs] = ode45(@(t, theta) ode_func(t,theta,b,g,L,m), t_span, theta_0);
%Plotting of Transient behavior of Pendulum
plot(t, outputs(:,1))
hold on
plot(t, outputs(:,2))
xlabel(\'Time\')
ylabel(\'Angular Displacement & Angular Velocity\')
title(\'Transient Behavior of Pendulum\')
legend(\'Angular Displacement\',\'Angular Velocity\')
hold off
% Creating Animation of the Pendulum
ct = 1;
for i = 1:length(t_span)
%Fixed end of the Pendilum
x0 = 0;
y0 = 0;
%Free end of the Pendulum
x1 = L*sin(outputs(i,1));
y1 = -L*cos(outputs(i,1));
%Plotting
plot([-1 1],[0 0],\'linewidth\',5,\'color\',\'k\')
hold on
plot([0 0],[0 -1.5],\'linestyle\',\'--\',\'color\',\'k\')
plot([x0 x1],[y0 y1],\'linewidth\',3,\'color\',\'r\')
%Creating Sphere using mark command
plot(x1,y1,\'marker\',\'o\',\'markersize\',30,\'markerFacecolor\',\'b\')
axis([-1.5 1.5 -2 0.5])
pause(0.0005)
M(ct) = getframe(gcf);
ct = ct+1;
hold off
end
movie(M)
videofile = VideoWriter(\'Pendulum Simulation.avi\', \'Uncompressed AVI\');
open(videofile)
writeVideo(videofile, M)
close(videofile)Explaination of code to plot transient behavior of pendulum:
> In this code initially all the inputs are inserted, after that the vector of angular displacement and angular velocity for 0 & 3 is defined, then the time period is defined for 20 seconds with 500 spans.
> The New Function is created solve 2nd order ODE as \"ode_func\" and the function is coded as below:
%Program to create function to solve 2nd order ODE for Pendulum
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> In this function (ode45) the equation for the motion of Pendulum is defined and executed.
> From the First Plotting the output are as below:
> Now the Plotting for the animation has to be coded in the programm which is also shown in first programm syntax.
> for making an animation file for loop is used and the pendulum is plotted. the fixed position and free end of the pendulum is defined in the program.
> Pivot point is coded and the string and the bob also plotted. for plotting of bob the mark command is used.
> After the programmimg the movie command used to make a n animation of the outputs.
Output of all positions of the pendulum:
OUTPUT:
The output of the programm:
The link for the output file or the Animation file is given below: