Solving second order ODE using MATLAB

AIM :

1. To solve the second-order ODE of a simple pendulum

2. To plot the graph between angular displacement vs time and angular velocity vs time

3. To make the animation of a simple pendulum

GOVERNING EQUATIONS :

1.Second order differential equation of simple pendulum is

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

To convert this into first order ODE replace

`θ= θ_2`

`(d θ_1)/dt = θ_2`

Our Euation will transform into

`(dθ_2)/dt^2 +b/(ml)θ _2 +g/l sinθ _1 =0`

where

θ = angular displacement rad/sec

m = mass of the pendulum bob kg

g= accleration due to gravity in `m/s^2`

l = length of the string in meters

t = time in seconds

b = damping coefficient

2.NEWTON'S SECOND LAW

`F=ma`

F = force in newtons

m = mass of the body in kg

a =accleration in `m/s^2`

OBJECTIVES OF THE PROJECT: The main object of this project is to solve the Second-order ODE of a simple pendulum and make use of it to make an animation of a simple pendulum.

I used Matlab online for this project.

I wrote a function called ode_function .m

So that I can call this function whenever I want to sole ODE.

Here is how it looks

%ODE45 Solve non-stiff differential equations, medium order method.

% Simple pendulum ODE Equation
% To formulate the ODE equation in the form of matrices 
function [dtheta_dt] = ode_function(t,theta,b, g,l,m)
    
    %convert the 2nd order ODE to 1st order ODE
    % d(θ_2)/dt + (b/m*l) *θ_2 +(g/l)*sin( θ_1) = 0
    
    theta1 = theta(1);
    theta2 = theta(2);
    dtheta1_dt = theta2;
    dtheta2_dt = -(b/m*l)*(theta2) - (g/l)*sin(theta1);
    dtheta_dt = [ dtheta1_dt ;  dtheta2_dt];
     
end

the above function takes the input arguments t, theta,b,g,l,m and returns the output argument dtheta_dt

To create an animation we should first create many figures and store them under one variable

The position of a simple pendulum is a function of time(t) and theta( θ)

We have to make use of the functions like plot, movie, ode45 and other basic function to build a simple pendulum figure using for loop to change theta

The code for the animation of a simple pendulum is simple_pendulum_ODE.m

clear all 
close all
clc

%To plot 2nd order ODE of Simple Pendulum

%Consider a simple pendulum having this intial requirements
% mass  = m in kg
m=1;
% length = L in meters
L=1;
% gravitaional constant = g in m/s^2
g= 9.81;
%damping coefficient = b
b=0.05;

%intitial condition 
theta_0=[0;3];

%time points
t_span=linspace(0,20,200);

%solving ODE
[ t, results ] = ode45(@(t,theta) ode_function(t,theta,b,g,L,m),t_span,theta_0);

figure(1)
%plotting the displacement and velocity

plot(t,results(:,1))
hold on 
plot(t,results(:,2))
ylabel('angular dispacement & angular velocity')
xlabel('time')
legend('angular displacement','angular velocity')
title('ANGULAR DISPALCEMENT, ANGULAR VELOCITY  vs TIME')

%To create an image that represents the position of pendulum
%Consider a fixed point (x0,y0) aboout which pendulum swings
x0=0;
y0=0;

temp=results(:,1);
%Position is determined by function of theta and time t
%We are assuming the pendulum is in max postion towards right side 
%let the angle made by string with verical is  alpha α
ct=1;
for i= 1: length(temp)
      alpha=temp(i);
       x1 = L*sin(alpha);
       y1= -L*cos(alpha);
       
       %plotting simple pendulum
       figure(2)
       plot([x0 x1],[y0 y1],'color','r','LineWidth',2)
       axis([-2, 2, -1.5, 0.3 ])
       hold on
       
   plot(x1,y1,'o','MarkerSize',14,'MarkerFaceColor','g')
      hold on 
      plot([-1,1],[0,0],'linewidth',5,'color','b')
      hold on
      plot([0,0],[0,-15],'linestyle','- -')
      
      title('Animation of Simple Pendulum')
      pause(0.3)
     hold off
      
      
 M(ct) = getframe(gcf);
 ct=ct+1;
 
end
movie(M)
videofile=VideoWriter('animation_of_simple_pendulum.avi','Uncompressed AVI');

open(videofile)
writeVideo(videofile,M)
close(videofile)

the link for the animation of the simple pendulum is provided here

https://youtu.be/8XWOs25LrR

OUTPUT: The graph between angular displacement vs time and angular velocity vs time is as shown in following figure

ERRORS need to be looked up at:

1. Make sure that you are using hold on command to get the plots under one figure

2. Don't forget to initialise the input arguments before passing them to a function

3. Use the debugging tool to correct the errors , it is always useful.

4. As I used Matlab online I faced some problem to create an animation file, I suggest Matlab downloaded version.