To solving the second order ordinary differential equation for making the pendulum animation and ploting the graph by using MATLAB
AIM
To simulate the transient behaviour of simple pendulum, plotting the angular velocity and displacement in function of time by using MATLAB.
INTRODUCTION
A pendulum is a weight suspended from a pivot so that it can swing freely. When the pendulum is displayed in sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate back toward the equilibrium position. When released, the restoring force acting on pendulum's mass cause it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.
SIMPLE PENDULUM DIAGRAM
DAMPING OSCILLATIONS WITH RESPECT TO TIME
Critically damped - A system moves as quickly as possible toward equilibrium without oscillating about the equlibrium
Undamped - A system oscillates at constant frequency without experiencing decay of its amplitude.
Underdamped - A system moves quickly to equlibrium, but will oscillate about the equilibrium point as it does so.
OBJECTIVE
1) To solving the ordinary differential equation in MATLAB for pendulum motion.
2) Simulating the pendulum oscillation to type appropriate coding in MATLAB.
3) Learn how to break it down from higher order to lower order differential equation.
4) Plotting the graph for angular velocity and angular displcement with respect to time.
5) To visualize how the string will move for each iteration by the command of for loop.
ODE EQUATION
1) The abbreviation of ode is ordinary differential equation is a differential equation containing one or more functions of one independent variable and the derivates of those functions.
2) For transient motion at any system, the second differential equation is help us to finding the unknown variables at instantenous point in the system trajectory.
3) For finding pendulum motion, this second order differencial equation used and it's written below.
`(d^2theta)/(dt^2)+b/m (d theta)/(dt)+g/l sintheta = 0` ------- 1
Pendulum physical specification
Length, l = 1 m
Mass, m = 1 kg
Damping coefficient, b = 0.05 Ns/m
Acceleration force , g = 9.81 m/sec^2
4) In order to finding the values, convert the equation from second order to first order in mathematically.
`theta = theta1`
`(d theta)/dt = (d theta1)/dt = theta2` ----- 2
`(d^2 theta)/dt^2 = (d^2 theta1)/dt^2 = d/dt((d theta1)/dt) = d/dt(theta2)` ----- 3
Subsitute the 2 & 3 value in equation 1
`(d^2 theta)/dt^2+b/m (d theta)/dt + g/l sintheta = 0`
`(d^2 theta1)/dt^2 + b/m (d theta1)/dt + g/l sintheta1 = 0`
`(d theta2)/dt +b/m theta2 + g/l sintheta1 =0`
`(d theta1)/dt = theta2`
`(d theta2)/dt = -b/m theta2 - g/lsintheta1`
`d/dt[ theta1; theta2] = [theta2;-b/m theta2 - g/lsintheta1]`
PROGRAMMING FOR SIMULATING THE PENDULUM
clear all % Remove all the variables and corresponding values from the workspace tab
close all % Closing the all figures
clc % Clean the data from command window
% Inputs value for simple pendulum
% Length of the string in m
l = 1;
% Mass of the ball in kg
m = 1;
% damping coefficient in Ns.m^-1
b = 0.05;
% Gravitational force in m.s^-2
g = 9.81;
% Initial condition for pendulum oscillation
theta_0 = [0; 3]; % it's a column vector, the initial condition of displacement and angular velocity is 0 & 3
% Time intervals
t_span = linspace(0,20,200); % linspace command is use to divide in linearly spaced values from starting value to ending value
% Solving ordinary differential equation
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure(1) % Making the figure
plot(t,results(:,1)) % plotting the time intervals in x axis and displacement in y axis, matrix describes in all row values & first column values
hold on % Hold on command is help to plotting overlap it previous one
plot (t,results(:,2)) % Plotting the time intervals in x axis and velocity in y axis,matrix describes in all row & second column result values
xlabel('Time in sec') % x axis in time
ylabel('Angular displacement & Angular Velocity') % y axis in angular velocity & displacement
legend('Angular dispalcement','Angular velocity')
title('Transient behaviour of a Simple pendulum')
figure(2) % Making the new figure
ct = 1 % Counter is count the loop
for i = 1:length(results(:,1))
x0 = 0; % pivot point at x coordinate
y0 = 0; % pivot point at y coordinate
x1 = (1*sin(results(i,1))); % Length of the string projection on x axis
y1 = (-l*cos(results(i,1))); % Length of the string projection on y axis
plot([-1 1],[0 0],'linewidth',10,'color','b') % to plot the thick line with given bounded distance
line([x0 x1],[y0 y1],'linewidth',5,'color','g') % line will plotting refer to x &y coordinates
hold on % Hold on is to allow the overlapping
plot(x1,y1,'o','markerfacecolor','r','markersize',30) % Plot the sphere at point of x1 & y1 with given dimension
axis([-2 2 -2 1]) % to fix the plot graph maximum & minimum distance of both x & y coordinates
grid on % To showing the grid on the graph
title('simple pendulum animation') % to create the title
pause(0.001) % pause between current figure to next loop figure
M(ct) = getframe(gcf); % M is a variable,getframe is save the current figure window
hold off % Hold off is stop the overlapping
ct = ct+1 % Increment to next counter
end % end the for loop
movie(M) % stitching the figures to make movie
videofile = VideoWriter('Simple pendulum.avi','Uncompressed AVI') % making a videofile with file name and file type
open(videofile) % open the videofile for writing the videofile
writeVideo(videofile , M) % writing into the videofile from collection all M figure
close(videofile) % Closing the videofile after completion of writing.EXPLANATION OF CODING LOGIC
1) Give the physical specification parameter for simple pendulum such as length of t string, mass of the ball, damping coefficient and gravitational force acting on a ball.
2) To set initial condition for simple pendulum such as initial angular displacement and initial angular velocity.
3) Set time how far system will oscillate with respect to time and make it in time array with linearly spaced values.
4) Make a function for solving lower order ordinary differential equation and that has definition of funtion in ode45.
5) We have using the ode45 command has been calculating the displacement and velocity with respect individual time array.
6) Results are printed in matrix form 200 rows, 2 columns. In row velocity and displacement printed with respect to each time intervals.
7) Plot the graph x axis is time and y axis is angular dispalcement, the same thing i do in matlab script window texting appropriate code.
8) Overlapping another graph x axis is time and y axis is angular velocity, same as in previous plot manner.
9) Both the plots are plotted in single figure.
SYNTAX OF ODE45
[t,results] = ode45(@(t,theta) ode_func(t,theta,bg,l,m),t_span,theta_0)
1) ode45 is a inbuild solver in matlab it's specially for solving the oridinary differential equation.
2) t and theta is main inputs, ode is function of t,theta and it is dependent and independent variables in matlab.
3) t_span is a time array it should provide a initial condition,provide integral value in each iteration from the time array.
4) Ode_func is definition of ode45, it instruct to how ode function is define.
5) [t,results] is a calling function in ode45, t is the time and results are angular displacement & velocity both are unknown it denoted in dtheta_dt for manual referance.
6) In each time array value is provided to the t_span in only one value, simultaneously two values provided in theta_0.
FOR LOOPING
1) for loop ( i ) is stored in only one value at current iteration. In next iteration, for loop will providing the next value in set of array belongs to commands.
2) We only takened as displacement value to subsituted in for looping condition
3) Counter have count the iteration after completion of one iteration or one looping , it will increment to second iteration in subsequently.
4) length(result(:,1)) colon indicates all the rows has taken , 1 indicates first column only. The result(:,1) values are 200 and it's 200*1 matrix.
5) That ( i ) value is passes to coordinates, it make sense to calculate the point for x1 & y1.
6) To create a code for subsequently calculating the coordinates points and plot them on a graph.
7) Pause it on each iteration, that plot will have been saved in M by command of getframe(gcf), it basically saved in current figure window.
8) Increment to next iteration.
9) To covert the figure to movie by using movie command.
PROGRAM FOR FUNCTION
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];
endEXPLANATION OF FUNCTION CODE
1) In previously, we convert the second order differential equation to first order differential equation in mathematically. That transform equation is perform into function
2) Function is basically reuseable code, in this ode_func is a template.
3) ode45 is to give two theta value, first value of theta is theta1 and second value of theta is theta2. That theta1 and theta2 value is subsitute into the dtheta1_dt and dthet2_dt as well.
4) Finally, dtheta_dt is calculated it has column vector, dtheta_dt is stored in name of results in ode45 command.
5) This is operate function of time in seconds.
RESULTS
PLOTTING FOR UNDERDAMP OSCILLATIONS IN PENDULUM
1) The angular displacement is initially 0 and angular velocity is initially 3, over the period of time oscillation had get damped out. So that oscillation displacement and velocity slightly falling down to the equilibrium position.
2) This oscillation is underdamping condition.
Reasons for falling down the oscillation time
1) Due air resistance against the ball or bob travel.
ANIMATION VIDEO
You tube link : https://youtu.be/2DrpL1BFAIw
ERRORS FACED WHILE PROGRAMMING
1) I using the comma in between the angular dispalcement & Angular Velocity it shows the error in line 28. After that removed the comma and quote in between both angular displacement and velocity.
1) I forget the closing parenthesis in line 41. After that i put the closing parenthesis in correctly.
1) I saving the function in different folder, command window showing the error i can't run this function.
2)After i saving in appropriate folder, matlab gets run smoothly.
CONCLUSION
Finally, to simulated the transient behaviour of simple pendulum and plot the angular displacement and velocity in graph for its coded by MATLAB.