Second Order ODE for a simple Pendulum with Damping Using Matlab
Second Order ODE for a simple Pendulum with Damping Using Matlab
Second order ODE for a simple pendulum with damping
YouTube link:-
https://www.youtube.com/watch?v=XSoj8k_meC4
Second order Ordinary Differential Equation of the position of the "ball" w.r.t time
Report is attache (Eqution photo taken from Attached Report)
Main programm for second order ODE for a simple pendulum with damping
% Solve the second order ODE for a simple pendulum with damping
clear all % Clear all Variables from memory
close all % close all charts and figures
clc % clear command window
% gravity in m/s^2
g = 9.81;
% length of the pendulum in m
L = 1;
% mass of the ball in kg
m = 1;
% damping coefficient
b = 0.05;
% Initial condition (angular velocity=3 rad/sec from time t=0)
theta_0 = [0;3];
% Time points
t_span = linspace(0,120,1200);
% Solve ODE
% ode45(@(t,theta) - find out Angular Displacement with respect to time
% ode_function(t,theta,b,g,L,m) - call ode function
% t_span - time array total 2000 from 0 to 120
% theta_0 - angular velocity=3 rad/sec from time t=0
[t,results] = ode45(@(t,theta) ode_function(t,theta,b,g,L,m),t_span, theta_0);
% subplot - To plot multiple plot in single figure
% (10,10,[61,70,91,100]) 10,10 is 10x10 matrix
% [61,70,91,100] - position of plot in 10x10 matrix
subplot(10,10,[61,70,91,100]);
% plot pendulum position by using results matrix first column with respect to time
plot(t,results(:,1),'linewidth',1) % linewidth use to define width of line
hold on % hold on use to put next plot in same chart
% plot pendulum position by using results matrix second column with respect to time
plot(t,results(:,2),'linewidth',1) % linewidth use to define width of line
hold off % hold off use for turn of the hold on command
title('Pendulum position') % title command use for provide title to plot
xlabel('Time') % lable x-axis
ylabel('Angular Displacement') % lable y-axis
grid on % turn on the grid
% subplot - To plot multiple plot in single figure
% (10,10,[8,10,38,40] 10,10 is 10x10 matrix
% [8,10,38,40] - position of plot in 10x10 matrix
subplot(10,10,[8,10,38,40])
plot(results(:,1),results(:,2),'linewidth',1) % plot both result column in single plot
title('Angular velocity Diagram') % title command use for provide title to plot
xlabel('Position radious'); % lable x-axis
ylabel('Angular velocity'); % lable y-axis
grid on % turn on the grid
ct = 1; % initial value for the plot number; set it 1
for i = 1:length(results(:,1)) % i is the range of first column of results
x0 = 0; % fixed point of link in x-axis
y0 = -0.04; % Fixed point of link in y-axis
x1 = L*sin(results(i,1)); % angular movement in x-axix
y1 = -L*cos(results(i,1)); % angular movement in y-axix
x2 = (L/2)*sin(results(i,1)); % movemet for text 'Link' in x-axis
y2 = -(L/2)*cos(results(i,1)); % movemet for text 'Link' in y-axis
% subplot - To plot multiple plot in single figure
% (10,10,[1,6,41,46]) - 10,10 is 10x10 matrix
% [1,6,41,46] - position of plot in 10x10 matrix
subplot(10,10,[1,6,41,46])
% Plot pendulum support line
plot([-0.6 0.6],[0 0],'linewidth',8,'color','k');
hold on % hold on use to put next plot in same chart
% Plot pendulum support dot
plot (0,-0.03,'o','markersize',12,'markerEdgecolor','k','markerFacecolor','k');
% Provid Link with length=1
line([x0 x1],[y0 y1],'linewidth',2,'color','b');
% Name the link with L=1
t = text((x2+0.04),y2,'L=1');
% give font size
s = t.FontSize;
% Mention font size of 12
t.FontSize = 12;
hold on; % hold on use to put next plot in same chart
% plot pendulum circle
% with mrkers give size to pendulum
% with markerfacecolor give face color to the pendulum circle
plot(x1,y1,'-o','markers',20,'markerfacecolor','r')
% Text command to provide 'm' as a name to pendulum circle
t = text((x1-0.04),(y1),'m');
% give font size
s = t.FontSize;
% Mention font size of 12
t.FontSize = 12;
title('Pendulum Damped Animation') % title command use for provide title to plot
xlabel('x-axis') % labeling of x-axis
ylabel('y-axis') % labeling of y-axis
grid on % enable grid in graph
axis([-1.5 1.5 -1.5 0.5]) % Fix axis by setting x and y limit
pause(0.003) % Pause interval between each frame
hold off % hold off use for turn of the hold on command
% Set frame position and size
a0=20; % Initial x-axis
b0=20; % Initial y-axis
width=900; % width of the frame
height=655; % hight of the frame
set(gcf,'position',[a0,b0,width,height]) % set poition and side using set command
% save all frame in M named array
% in M(ct), ct represents each and every frame to store in M
% use getframe command to get current frame window in each loop
M(ct)= getframe(gcf);
% increment of number of frame in each loop
ct=ct+1;
end
% move command generates movie using frames in M array
movie(M)
% Create videofile using VideoWrite command to collect all frames in M array
videofile = VideoWriter('pendulum_damped.avi','Uncompressed AVI');
% open video to collet all frames from M array
open(videofile)
% writeVideo capture frames in 'M' into videofile
writeVideo(videofile,M)
% at the end close videofile
close(videofile)ODE function Programm
% ODE Function
% create ODE function
% Function of second order ode and its entities
function [dtheta_dt] = ode_function(t,theta,b,g,L,m)
theta1 = theta(1) % In every loop take new value from theta2
theta2 = theta(2) % In every loop take new value of theta_0
dtheta1_dt = theta2; % as per diff, equation put dtheta/dt as theta2
% After putting the theta2 in equation
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);
% Value of diff. Equation in matrix form
dtheta_dt = [dtheta1_dt; dtheta2_dt];
% End of function
endPlot
1) represents the position of the pendulum
2) Angular dispacement vs Time Diagram
3) Angular velocity vs position radius Diagram
4) pendulum animation screenshot
YouTube link:-
https://www.youtube.com/watch?v=XSoj8k_meC4