Week 3 - Solving second order ODEs

Aim :-

• To write a program in MATLAB to solve 2^nd order Ordinary Differential Equation of a Simple Pendulum with damping.
• Simulate the motion of this pendulum between t=0 to 20 seconds,for angular displacement = 0 & angular velocity = 3 rad/s at t = 0.

Governing Equation :-

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

where g = acceleration due to gravity (in m/s^2)

l = length of the pendulum (in m)

b = damping coefficient

m = mass of the bob (in kg)

• The above equation represents the 2^nd order ODE which describes the motion of a simple pendulum with damping.
• Upon solving the above equation we get the position of the pendulum as function of time i.e.,.
• Now,we’ll use this data to animate the motion of the pendulum .

MATLAB Code :-

clear all
close all
clc


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


% Initial Conditions [Displacement,Velocity] 
theta_0 = [0;3];

% Time Limits (For obtaining solution to my ODE by integration we're defining time limits)    
t_span = linspace(0,20,20);

% Solving ODE (Calling ODE function)
[t,results] = ode45(@(t,theta)  ode_func(t,theta,b,g,l,m),t_span,theta_0);

displacement = results(:,1);
velocity = results(:,2);
s=1;

% Simulation 1 : 
for i = 1:length(t_span) % 20 values of time in total.
   
      x0 = l*sin(displacement(i));
      y0 = l*cos(displacement(i));

      x1 = 0 ;
      y1 = 5 ;


plot([x0 x1],[y0 y1],'-o','LineWidth',2,'MarkerSize',10,'MarkerEdgeColor','r','MarkerFaceColor',[1 0 0],'MarkerIndices',1:length(y0))
axis([-1 1 0 5])
xlabel('x-axis')
ylabel('y-axis')
title('Simple Pendulum')

M(s) = getframe(gcf);
s = s+1;
end



movie(M) % Matlab is going to create a movie with all the frames in M
% Creating avi files to dump out all the frames
% Creating videofile using 'VideoWriter' command
videofile = VideoWriter("Simple Pendulum Simulator",'Uncompressed AVI');
open(videofile) % Opening this videofile
writeVideo(videofile,M) % Writing in the frames that we captured in the videofile 
close(videofile)  % closing the videofile

% Pause command is not used,since it pauses the code which isn't very efficient.
 

















 

Function Code :-

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];

end

Workflow :-

• Allocate the values of b,m,l & g as given in the problem.
• Define initial conditions for displacement and velocity as 0 rad & 3 rad/s respectively.
• Define the time values in t_span variable using linspace command for 0 to 20 sec for integrating the ODE.
• Using ode45 solver to solve the ODE and saving the displacement and velocity in results variable column wise.
• Now,accessing these column vector values in displacement and velocity variable.
• Using ‘s’ variable for counting the no. of frames and setting its value initially to 1.
• Looping is done to get each values of displacement wrt time values and plotting them ultimately using plot command.
• Bob of the pendulum is made using respective ‘Marker’ commands.
• Fixixng the reference axis so that the pendulum appears in the middle of the x-axis.
• Getting in all the frames in M using getframe.
• Incrementing the frame counter variable ‘s’ to s+1 after every loop execution.
• With all the frames in variable ‘M(s)’ matlab is going to create a movie.
• Creating videofile using ‘VideoWriter’ command and also creating avi files to dump out all the frames.
• Opening the videofile using open command.
• Using writeVideo we’re writing in all the frames that we captured.in the videofile.
• close command is used to close this videofile.

Syntax used for ode45 solver :-

• [t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], 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.

Youtube link :- https://youtu.be/YAAoS1mdiqE