Week 4.1 - Solving second order ODEs
A pendulum is a body suspended with a string to move freely back and forth under the influence of the gravity. When a pendulum is displaced from it’s equilibrium position, it is subjected to a restoring force due to gravity, that will accelerate it back towards the equilibrium position.
When, we have the second order ODE. Then we have to break it down in to the frist order ODE.
- Displacement
- Velocity
Challenge Pendulum Animation
function [dthetha_dt] = pand(t,thetha,b,g,l,m)
% To solve the 2nd order ODE in to the 1st order ODE.
thetha1 = thetha(1);
thetha2 = thetha(2);
% Angular Velocity
dthetha1_dt = thetha2;
% Angular Displacement
dthetha2_dt = -(b/m)*thetha2 - (g/l)*sin(thetha1);
% For final vector solution.
dthetha_dt = [dthetha1_dt ; dthetha2_dt];
end
% Main Code
clear all
close all
clc
b = 0.05; % Length of the string Damping coefficient.
g = 9.81; % Gravity Constant
l = 1; % Length of the string
m = 1; % Mass of the ball.
thetha_0 = [0;5]; % Initial Velocity is ‘0’ and the Initial displacement ‘5’.
t_span = linspace(0,20,500); % Time span for Pendulum Oscillation.
% To get solution, we need to integrate.
% Integration is done with in certain limits. And these limits are defined.
% Integration is done by the ODE45. It is a inbuilt command of MATLAB.
% ‘t’ is the integration time. And t_span is the time array.
[t, result] = ode45(@(t,thetha) pand(t,thetha,b,g,l,m),t_span, thetha_0);
% dthetha_dt is obtained by the ‘pand’ function. In the function command, we has %two output. So, we have write two variables [t,results]. Result has two columns. %Both are velocity and Displacement, with respects to the time.
figure(1)
plot(t,result(:,1))
hold on
plot(t,result(:,2))
hold off
xlabel('plot')
ylabel('time')
legend('Displacement','Velovity')
% Plot command is used to plot graph of ‘angular displacement’ and ‘angular
% velocity’ with respect of the time.
x0 = 0;
y0 = 0;
mv = 1;
for i = 1:length(t_span);
x1 = l*sin(result(i,1));
y1 = -l*cos(result(i,1));
figure(2)
plot([x0,x1],[y0,y1],[-1,1],[0,0],'r-','linewidth',4)
hold on
plot(x1,y1,'marker','.','markerfacecolor','m','markersize',100)
xlabel('Angular Displacement')
ylabel('Angular Velocity')
legend('Pendulam')
title('Pendulam Animation')
axis([-1.5,2,-1.5,2])
hold off
pendulam_movie(mv) = getframe(gcf);
mv = mv+1
pause = 0.004;
end
% for loop command is used to plot different graphs at different time iteration.
% By using hold on/off command, we can plot different graphs at the same axes.
% We have used inbuilt ‘getframe(gcf)’ command to animate the movie_pendulam.
movie(pendulam_movie)
videofile = VideoWriter('animation_pendulam.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,pendulam_movie)
close(videofile)