Week 3 - Solving second order ODEs

AIM: To write a program in MATLAB that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

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.

INPUTS

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

time = 0-20 sec,

for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Procedure:

• Since the given problem represents a second order differential equation, converting this 2nd order ODE into 1st order ODE by following below steps

theta = theta1

`(d theta)/ dt` = `(d theta_1)/dt`

`(d theta_1^2)/dt^2 = (d theta_2)/dt`

CODE:

clear all
close all
clc

%INPUTS
b = 0.05; %damping coefficient
g = 9.81; %gravity
l = 1; %length of pendulum
m = 1;  %mass of pendulum

%intial condition
theta_t = [0;3];

%time points
t_span = linspace(0,20,100);
%---------------------------------------------------------------------------------------------------------
%solving ode
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_t);

plot(t,results(:,1),'b')
hold on
plot(t,results(:,2))
ylabel ('plot')
xlabel('time')
title('2nd order ode')

%--------------------------------------------------------------------------------------------------------------
%animation
ct =1;
theta_1 = results(:,1);
for i =1:length(theta_1)
   %starting co-ordinates
    x0 = 0 ;
    y0 = 0;
    
    %end co-ordinates
    theta_2 = theta_1(i);
    x1 = l*sin(theta_2);
    y1 = -1*cos(theta_2);
   
    %------------------------------------------------------------------------------------------------------
    figure(2)
   %Pendulum bob
    plot([x0 x1], [y0 y1], 'o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','b','linewidth',5)
    hold on
    
    
   %defining axis
    axis([-2 2 -2 2])
    
    
   %string of pendulum
    line([x0 x1], [y0 y1], 'linewidth',4,'color','r')
    hold on
    
    
    %Support for pendulum / base for pendulum
    plot([-1 1], [0 0], 'linewidth',5,'color','[0 0.3 0]')
    title ('Pendulum')
    M(ct) = getframe(gcf)
    hold off
    ct =ct+1;
    
end
%-------------------------------------------------------------------------------------------------------------
%Creating movie
movie(M)
video_file= VideoWriter('Pendulum_challenge.avi', 'uncompressed AVI');
open(video_file)
writeVideo(video_file,M)
close(video_file)

• Here we make use of inbuilt ode45 function to convert the 2nd order ODE into 1st otder ode

Program for ODE function:

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

RESULTS:

• The animation video is uploaded in google drive and link to acces it is attached below.

https://drive.google.com/file/d/15MeQ9VDthU0M0rfZ-pPG2HzUd6DEFQEG/view?usp=sharing

ERRORS FACED:

1. Pendulum was moving in opposite direction, this error was corrected by including negative sign in y1
2. Forgot to add the string this error was corrected by including command plot for string of pendulum
3. Pendulum was continuosly swinging , this error was corrected by including hold off command.

Conclusion:

• MATLAB code to solve the ODE is succesfully compiled.
• This ODE represents the equation of motion of a simple pendulum with damping and the animation is succesfully created