Week 3 - Solving second order ODEs

Aim:-

To simulate the motion of simple pendulum whose length is 1 meter , mass of the bob is 1 kg , damping ration is 0.05 , g is 9.81m/s^2. Simulate the motion between 0-20sec, for angular displacement = 0 , angular velocity =3 rad/sec at a time t=0. Use this information to animate the motion of the pendulum.

Equations used in the program:-

`(d^2 θ)/(dt^2 )+b/m*dθ/dt+g/l*sinθ .... (1.1)`

This ODE represents the equation of motion of a simple pendulum with damping.

`s=l*sinθ ..........(1.2)`

This is the equation for the angular displacement from the equilibrium. `sinθ~θ(radians)`,when the displacement from the equilibrium is small. 

And equation for   angular velocity is given by

`ω=(ds)/(dθ) .................(1.3)`

Since we were given the values at the equilibrium line we are using this formula.

Objectives of the project:

Writing a program that solves equation 1.1, by using the data in the question write a code that can simulate the motion of a pendulum.

Body of the content:

We have used MATLAB to write the code to simulate the forward kinematics of a 2R robotic arm, MATLAB uses C++ in its interface.

The code to simulate the motion of a simple pendulum:-

%program to simulate pedulum motion.
clear all
close all
clc

%inputs 
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition 
%here the value Zero represents angular displacement and the value 3 represents angular velocity which are given in the question.
theta_0=[0;3] 

%time intervals 
% here we are creating an time period array to calculate the position of the pendulum for each value of time.
t_span= linspace(0,20,150);

%solve ode
%ode45 is an inbuilt function to solve ode's in matlab ; here our ode is a function of t, theta. here we are calling a fuction called ode_func which we wrote. 
[t,theta] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

%ploting for tracing of angular velocity and displacement
%here we are using this command to show two graphs in single plot .
subplot(2,1,1)

%here we are ploting between t and theta(all values in the row ,1 st column) to get a graph between angular displacement and time. here samll b represents blue colour.
plot(t,theta(:,1),'b','LineWidth',2)

%this command is used to namme the y axis 
ylabel('angular displacement')

%this command is used to name the x axis
xlabel('time in s')

% this command is used to show the grid in the plot. 
grid on
subplot(2,1,2)
hold on 
%here we are ploting between t and theta(all values in the row ,2nd column) to get a graph between angular velocity and time.here k represents black colour. 
plot(t,theta(:,2),'k','LineWidth',2)
ylabel('angular velocity ')
xlabel('time in s')
grid on
% creating an animation out of the plot 
s0= 0;%angular diplacement
v0= 0;%angular velocity
ct= 1;%counter represents eacch and every frame.

for i=1:length(t)
    % using the formula to calculate angular displacement we are calculating the position of the pendulum at that particular theta value.
    s1=l*sin(theta(i,1));

    % using the formula to calculate angular velocity we are calculating the velocity of the pendulum at that particular theta value.
    v1=-l*cos(theta(i,1));
    
    figure(2)
    %here we are ploting this to represent support rod of the pendulum 
    plot([-1 1],[0 0],'LineWidth',6,'color','g');

    % we are creating an axis of (-2 to 2 in x direction ) and (-2 to 2 in y direction )
    axis([-2 2 -2 2]);

    % we are creating a plot of pendulum rod's (length of the pendulum) different positions. 
    line([s0,s1],[v0,v1],'linewidth',2,'color','b');

    %we are using a hold on command so that we dont get overlap of one plot onto other.
    hold on

    % we are using a command called 'o' as a marker and by using markersize command we are able to customize the size ofthe marker. we are creating a fixed point. 
    plot(s0,v0,'O','MarkerSize',3,'MarkerFaceColor','y');

    % we are creating bob of the pendulum
    plot(s1,v1,'O','MarkerSize',20,'MarkerFaceColor','r');

    % we are assigning a title to the plot.
    title('Damping of Simple Pendulum');

    %we are assigning a legend to the plot.
    legend('Support rod','length of the pendulum','Fixed point','Bob of pendulum');
    grid on;% to show grid lines we are using grid on command 

    % since we have to get seperate plots or frames  in order to make a video instead of a animation (gif image). we swithcing off hold command. 
    hold off

    %we are creating an array of frames using the plot we generated. 
    M(ct)= getframe(gcf);
    %the conter variable is being incremented
    ct=ct+1;

end 

movie(M);%we are merging all the frames into one using movie command.
videofile=VideoWriter('Simple Pendulum','Uncompressed AVI');%we are creating a videofile and assigning its name and type of video fime.
open(videofile);%we are opening the video file.
writeVideo(videofile,M);%we are writing the data from M to video file. 
close(videofile);%we are closing the vidoe file.

Code for the function used inside the program:-

function [dtheta_dt]= ode_func(t,theta,b,g,l,m)

% here we are assigning theta1 as variable theta's 1st element  
theta1 = theta(1)
%here we are assigning theta2 as theta2 as variable theta's 2nd element  
theta2 = theta(2)
% here we are assuming dtheta /dt = theta2 to transform 2nd order ode to 1st order. 
dtheta1_dt = theta2;
%this is the converted 1st order ode after substituting theta2 in place of dtheta/dt in the given ODE.
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1);
% we are creating an array of 2 columns to represent it in the form of AX=B here X=dtheta_dt. 
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

Explanation of the function:-

• We are giving inputs inside the brackets of the function while we are writing a function so that we use them.
• here we are assigning theta1 as variable theta's 1st element
• we are assigning theta2 as theta2 as variable theta's 2nd element
• we are assuming dtheta /dt = theta2 to transform 2nd order ode to 1st order.
• we are creating an array of 2 columns to represent it in the form of AX=B here X=dtheta_dt.

Explanation of the code ;

• We used a “ % ” to write a comment and we used “ clear all ”  command to make sure that the work space is clear and we used “ close all ” command to make sure to close all plots and figures from previous run. And “ clc ” command to make sure the command window is clean.
• We have given inputs to the program from the given data.
• We have taken a variable theta_0 in which one column is angular displacement and other is angular velocity
• Next we are creating an array of time interval of 0-20 secs with 150 divisions
• We are solving the ODE to get values of t and theta (position of  the pendulum at that particular time).
• We are solving the ODE using a inbuilt function called ode45. Our ODE depends on t, theta. And to be able to differentiate, ode45 integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.
• After calculating the values of t and theta , we plot a graph between t and theta.
• We use a inbuilt command called subplot to plot two plots separately in one plot.
• In 1^st plot we are plotting between t and theta (with values in all rows in 1^st column) to visualize data between time and angular displacement.
• In 2^nd  plot we are plotting between t and theta (with values in all rows in 2^nd column) to visualize data between time and angular velocity.
• We used xlabel and ylabel commands to give names to xaxis and yaxis of the plot.
• To create an animation we have 1^st taken initial time, position and angular velocity of the pendulum.
• Then using a for loop we have calculated all the values at different positions and velocities at different time intervals.
• We have used equations 1.2 and 1.3 to get the positions and velocities at different positions.
• To create an animation we have to 1^st create a support rod so we are plotting a horizontal line of length 2 units at the origin with linewidth 6.
• Then we are creating a axis to be able to visualize full figure.
• Now we are plotting the line or the pendulum rod (length of the pendulum) within the loop to get different plots which we can use to create a video file.
• To create a clean animation we are using hold on command.
• Now while plotting we are specifying the marker symbol and markersize to be able to create a bob for the pendulum and fixed pin at the origin to make the animation look clean.
• We are giving a tittle to the plot by using the tittle command.
• We are showing the legends using legend command.
• Next we are creating an array M and saving all the plots as frames and we are using a counter variable and updating it in the loop.
• Then we are merging all the frames of array M into one using movie command.
• Then we are creating a video file with tittle simple pendulum and type of the video is uncompressed avi.
• Then we are opening the video file and writing all the data of array M to the video file and closing the video file.

Results:

Are shared in the attachment and for animation link is shared below.

https://drive.google.com/file/d/1b8W8GpAvt_wnxcBPyReZDUgvDMjYa4WP/view?usp=sharing

As we can see from the above graphs we can observe a sine graph in plot ‘a’ which and cosine graph in plot ‘b’ from equation 1.2 and 1.3 we can say that out result is pretty accurate.

Conclusion :-

• The angular displacement of the pendulum depends on the length of the pendulum. If the length of the pendulum is more then the pendulum will have more number of rotations.
• Using subplot we can plot multiple plots in one single plot.
• Ode45 function integrates the system of differential equation from an initial t_0 to t_f with initial condition.
• We were able to create an animation from solving the differential equation.