Week 3 - Solving second order ODEs

Aim:

Write a program in MATLAB to solve 2nd order ODE for a simple pendulum and to plot angular displacement and angular velocity graph to analyze its motion to observe transient behavior of a pendulum.

Solution:

Assume a simple pendulum with a pendulum bob which has a mass of m and suspended with a string.

Now, for initial condition, θ = θ1

Hence,

dθ/dt =dθ1/dt =θ2

Further,

d2θ/dt2 =d2θ1/dt2=d/dt(dθ1/dt)=d/dt(θ2)

2nd order ODE for mass-pendulum damping system,

d2θ/dt2 +b/m ⋅dθ/dt +g/l*sin(θ1)=0

So,

d2θ1/dt2 +b/m ⋅dθ1dt +g/l*sin(θ1)=0

d2θ2/dt2 +b/m ⋅θ2 +g/l*sin(θ1)=0

1st ODE

dθ1/dt =θ2

2nd ODE

dθ2/dt =−b/mθ2−g/l*sin(θ1)

Nomenclature:

b = damping coefficient

l = length of string, m

m = mass of bob, kg

g = acceleration due to gravity, m/s^2

Θ = angular displacement

Program:

clear all 
close all
clc 

b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%Initial Condition
theta_0 = [0;5];

%time points 
t_span = linspace(0,10,500);

%solve ODE
[t,results] = ode45(@(t,theta) ode_fun(t,theta,b,g,l,m),t_span,theta_0);


plot(t,results(:,1))
hold on
plot(t,results(:,2))
hold on 
ylabel('plot')
xlabel('time')


% creating an animtion 

theta = results(:,1)
ct = 1

for i = 1 : length(theta)
    Theta = theta(i)

    x0 = 0
    y0 = 0

    x1 = l*sin(Theta)
    y1 = -l*cos(Theta)

    plot([x0,x1],[y0,y1],'LineWidth',3)
    hold on
    plot(x1,y1,'o','MarkerSize',15,'MarkerFaceColor','r')
    hold off

    axis([-1 1 -1.5 0])
    grid on 
    pause(0.3)
    M(ct) = getframe(gcf)
    ct = ct + 1
end

movie(M);
videofile = VideoWriter('pandulam.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

function[dtheta_dt] = ode_fun(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

Output:

Conclusion:

The effect of damping can be clearly seen from the graph as the graph of angular displacement and angular velocity displacement depicts that the motion of pendulum is gradually decreasing with time. Hence, any 2nd order ODE can be solved with MATLAB

Link to animation:

https://drive.google.com/file/d/1gafCQ9ZKRdbp1M1EsnLoIv4x2IlLlWYR/view?usp=sharing