SOLVING SECOND ORDER ODE OF SIMPLE PENDULUM
AIM: -
Program for second order ODE of the simple pendulum by using the MATLAB.
To simulate and create an animation of the motion of the simple pendulum.
GOVERNING EQUATIONS: -
(d^2 θ)/(dt^2) + b/m dθ/dt + g/L sinθ = 0
theta1= theta(1)
theta2= theta(2)
dtheta1/dt = theta2;
dtheta2/dt = -(b/m)*theta2-(g/l)*sin(theta1);
dtheta/dt = [dtheta1/dt; dtheta2/dt]
x1 = l*sin(THETA1)
y1 = -l*cos(THETA1)
Objective: -
To write a Program for second order ODE of the simple pendulum by using the MATLAB.
To simulate and create an animation of the motion of the simple pendulum.
Body of the content: -
By using the MATLAB, writing program for the second order ODE of the simple pendulum and plotting the results and to create a movie of an animation for the motion of the simple pendulum
The code is: -
clear all
close all
clc
b = 0.05; % b is the damping coefficient
g = 9.81; % g is the gravity (m/s^2)
l = 1; % l is the length of the pendulum (meters).
m = 1; % m is the mass (kg).
% here we've two first order ODEs ( velocity and displacement )
%initial condition
theta_0 = [0;3]; % initial displacement of pendulum is 0 rad. & initial velocity is 3 rad/s.
%defining value for time
t_span = linspace(0,20,250); % t_span is the time array i.e., 0 to 20 seconds.
% calling our ODE function
% here ode45 is an inbuilt solver in matlab, our ode depends upon time and theta
% ode_func is a function created with the input arguments
% solve ODE
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);
figure(1)
hold on
plot(t, results(:,1));
plot(t, results(:,2));
xlabel('time')
ylabel('plot')
ct=1
for i= 1:length(results(:,1))%length of results(:,1) means the no. of values in results(:,1).
THETA1=results(i,1);% here i in the for loop is a loop counter that it tells how many times the for loop going to executed.
% here the benefit with for loop is it's not a array, everytime it executes only one value upto i times.
x0 = 0;% coordinates for the fixed end of the link.
y0 = 0;
% Note:- for each for loop the loop counter variable should be different.
x1 = l*sin(THETA1);% here the angle is in radians. so, we use "cos" instead of "cosd" which takes the values in degrees.
y1 = -l*cos(THETA1);% coordinates for movable end of the link.
% plotting the link and weight.
figure(2)
plot([x0 x1],[y0 y1], [-0.6 0.6],[0 0], 'linewidth',6)
hold on
line([0 0], [0.35 -1.35], 'linewidth',1, 'linestyle','- . ')
plot(x1,y1, 'marker', 'o','markerfacecolor','m','markersize',30)
hold off
title('SIMPLE PENDULUM')
axis([-2 2 -1.5 0.5])
pause(0.1)
M(ct) = getframe(gcf); % ct represents each and every frame.In the getframe command calling the 'get current figgure (gcf)' and it'll be saved in an array called 'M'
ct=ct+1;
end
movie(M)% it means matlab generate a movie with all the individual frames.
videofile = VideoWriter('ode.avi', 'uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
ODE function
Note: - in the above code there called an ODE function
The code is :-
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
% as we know that "theta" is our solution array
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
Errors encountered while running the code: -
1.Here I forgot to define the hold off command and it became error
Solution: - It’s rectified by adding hold off command there.
Plot Results: -
As per the objective the plot got executed and this is the picture of that plot.
Movie link: - This is the movie of the motion of the simple pendulum animation file, google drive link below.
https://drive.google.com/file/d/1Oj0hIxINsAiANjHxOSzAWqdZD9TO08bo/view?usp=sharing
Conclusion: -
Program was written for the second order ODE of the simple pendulum,
Successfully created a movie of its motion animation and
Most of the errors I found were the typo errors and got solutions to them due to these things I realized and learnt how some new commands are existed.