AIM: Solving the second order ODE in matlab.

let's solve second order ODE by range Kutta method.

our aim is to simulate the unsteady state system of the pendulum.

lets solve second order ode in Matlab. the equation is below

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Matlab only understands first order ode, therefore we have to convert this second order ode into the first-order ode.

let's assume

Θ=Θ(1)

(dΘ/dt)=(dΘ(1)/dt) Θ(2) %% assuming the first derivative of Θ(1)=Θ(2)
%%placing the value of Θ

(d^2Θ/dt^2)=(d^2Θ(1)/dt^2) =dΘ(2)/dt %%%placing the value of Θ

putting the values in the main equation and rewrite the equation

let's convert this equation in matrix form.

deriving this equation in MatLab by function method

for function

%pendulam equation d^2theta/dt^2+(b/m)(dtheta/dt)+(g/l)sin(theta)=0
%g = gravity in m/s2,
%L = length of the pendulum in m,
%m = mass of the ball in kg,
%b=damping coefficient.

function [thetadot]=solving_ode_for_pendulum(t,theta,g,l,m,b)
theta1=theta(1);
theta2=theta(2);
%angular displacement
thetadot(1)=theta(2);

%anguar velocity;
thetadot(2)=(-g/l)*sin(theta1)-(b/m)*(theta2);

%thetadot is 2 row 1 column matrix
thetadot=[thetadot(1);thetadot(2)];

end

some basic trigonometric for x1 and y1

solving the equation by ode 45 method.

% defining the variables and solving the equations by runge kutta method
clc 
close all
clear all

% defining the values
g=9.81;%in m/s2,
l=1;% in m,
m=1;%in kg
b=0.05;
 % this is second order ode we need two conditions to solve
 %initial conditions
 theta0=[0;3];
 %time span
 t0=linspace(0,20,500) 
 % solving ode 
 [t,theta]=ode45(@(t,theta)solving_ode_for_pendulum(t,theta,g,l,m,b),t0,theta0)
 figure(1)
 plot(t,theta(:,1),'c')
 hold on
 plot(t,theta(:,2),'r')
 xlabel('time')
 ylabel('angular displacement and velocity')
 legend('angular displacement','angular velocity')
 hold off
 figure(2)
 plot(theta(:,1),theta(:,2))
  xlabel('angular displacement')
  ylabel('angular velocity')
        
 %creating the animation
for i=1:length(theta(:,1))
    x0=0;
    y0=0;
    x1=1*sin(theta(i,1));
    y1=-1*cos(theta(i,1));
    
    figure(3)
    %drawing the fix line
    plot([-5 1],[0 0],'linewidth',5,'color','r')
    hold on;
    %string for pendulum
    plot([x0 x1],[y0 y1])
   hold on;
   %pendulum mass
   plot(x1,y1,'-o','markers',30,'markerfacecolor','y')
   axis([-2,2,-2,1])
   grid on
   pause(0.002)
   hold off
   animation(i)=getframe(gcf);
end

movie(animation)

videofile=VideoWriter('pendulum_animation.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,animation)
close(videofile)
 

plot of angular velocity and displacement with respect to time.

plot of angular displacement and angular velocity

• Here we can conclude that at t=0 sec and displacement and velocity are taking place and it is continuously decreasing with respect to time. at one particular time(at the end), it will gain the original position. the rate of decreasing of displacement and velocity are different.

youtube link for an animation of a pendulum :