Solving second order ODEs using Matlab

AIM-Write a program in matlab that represents the equation of motion of a simple pendulum with damping.

PROGRAMMING LANGUAGE USED IS MATLAB

OBJECTIVE OF THE PROJECT

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time. 

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

LET US FIRST TRY TO UNDERSTAND HOW THE DIFFERENTIAL EQUATION  IS DERIVED USING NEWTONS SECOND LAW 

THE PICTURES ABOVE SHOW THE DERIVATION OF THE ODE EQUATION FOR A SIMPLE PENDULUM WITH DAMPING

THE ODE EQUATION IS GIVEN BY

In the above equation

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

INITIAL CONDITIONS SPECIFIED

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

We are Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

THE CODE IS BEING DIVIDED INTO TWO PARTS 

1)In the first part we are using the ode 45 function to solve the differential equation and pendulum animation is generated using loops.

2)The ode is being defined using functions.

FIRST PART OF THE CODE

clear all
close all
clc
%Giving the parameters
%Giving the damping coefficient
b=0.05;
%Giving the acceleration due to gravity
g=9.81; %S.I unit is m/s^2
%length of the pendulum
%Actually we are considering the spring mass system but the differential
%equation for a spring mass system is same as that of the pendulum so we
%use the equation of pendulum according to our convinience.
l=1; %Lenght of the pendulum in m.
m=1;%mass of the bob in kg.

%Giving the initial values
%Giving the value for theta
theta_th=[0;3];
%Time points
t_span=linspace(0,20,500);%It will be used to perform integration .
%It is being divided into 500 parts to plot the graph between time and
%dispacement and time and velocity
%Using the ODE solver to solve the ODE
[t,results]=ode45(@(t,theta)odesmd_func(t,theta,b,g,l,m),t_span,theta_th);
plot(t,results(:,1))
hold on
plot(t,results(:,2))

xlabel('time');
ylabel('displacement,velocity')

THETA_1=results(:,1);
ct=1;
%%Initialising the loop for Animation
for i=1:length(t)
r=THETA_1(i);
t_k=t(i);
s_k=results(:,1);
%Providing the initial state of the pendlum
x0=0;
y0=0;
%Prviding the final state of the pendulum
x1=1*sin(r); %The coordinate axes chossed is the 4th coordinate.
y1=-1*cos(r);%The coordinate axes choosed is the 4th coordinate -ve of y axis is being taken
%That is why we are using - sign here
figure(2) %If we donot give figure 2 then matlab will plot the pendulum in the same figure
plot([-1,1],[0,0],'linewidth',2,'color','yellow');%Base line

line([x0,x1],[y0,y1],'linewidth',1.5,'color','green');%%pendulum
xlabel('time');
ylabel('angular displacement');
title(['Particle at t= ',num2str(t_k),'seconds']);

hold on;

%%pendulum ball
line([0,0],[0,-1],'linewidth',0.5,'color','red','Linestyle','--')%plotting the mean position line

plot(x1,y1,'o','markersize',30,'markerfacecolor','r','MarkerEdgecolor','b');%ball
axis([-2 2 -2 1]);
%Creating a Rectangle
rectangle('Position',[-1.5 -1.5 3 2],'Curvature',0.2)



pause(0.02);
hold off;
M(ct)=getframe(gcf);%Capturing the frame
ct=ct+1;%%Iterating the frame
end

%%Creating Movie
movie(M)
videomifile=VideoWriter('Second_order_ode.avi','uncompressed AVI');
open(videomifile)
writeVideo(videomifile,M)
close(videomifile)

SECOND PART OF THE CODE

function[dtheta_dt]=odesmd_func(t,theta,b,g,l,m)
theta1=theta(1)%Here theta 1 is the dispacement
theta2=theta(2)%And the velocity is given by theta 2

dtheta1_dt=theta2;
dtheta2_dt=(-b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

CODE EXPLANATION

Matlab has a pre defined function known as ODE45 it automatically solves differential equations.

The initial conditions used here are:-

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

We are using for loops to create an animation of the damped pendulum by plotting its angula displacement vs its time

we are also plotting the the graph of displacement ,velocity of the pendulum vs time.

WE HAVE TO SOLVE 2 ODE 

Using the concepts of functions we are solving providing the ode to the ode solver and the plotting the graphs

We are using the concepts of loops to create the animation of the ODE

OUTPUT OBTAINED

The differential equation is solved using the ode solver in matlab and the results obtained are shown below

THE ANIMATION 

We are Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

https://youtu.be/ppa2op7TkSw-THE YOUTUBE LINK FOR THE ABOVE VIDEO 

THE GRAPH FOR THE DISPLACEMENT VS TIME AND VELOCTY VS TIME IS GIVEN BY THE PICTURE BELOW

 ERRORS FACED WHILE PROGRAMMING

----//Errors faced while programming//----

This is not a complicated code so the only error that can occur is while creating the functions and only minute issues and can be solved by the following methods

Error 1-Saving the name of the file different as that given in matlabsolution--The file name should be the same as that of the function name other wise it may cause error

Error 2- writing the equations wrong ,always write the correct equations and cross check with the value obtained mathematically with the values obtained through the code.

Error 3- While using the plot function always check to give the right parameters

Error4 --Matlab is Case Sensitive so take care of that

Error5-- Make sure to provide % sign before command statements . Use help command to get information of any function or go to the mathworks website to know more about the function