Week 3 - Solving second order ODEs

Aim :

To solve the ODE which represents the equatio of motion of a simple pendulum with damping. To write a program that solves the ODE which represnts the equatio of motion of a simple pendulum with damping and create an animated video of output obtained by solving this ODE.

Theory :

Here,

• g = gravity in m/s^2.
• L = Length of the pendulum in meters.
• b = damping coefficient.
• m = mass of the ball in kg.

A simple pendulum consists of a mass 'm' hanging with a string of length 'L' and fixed at a pivot point. When displaced to an initial angle 'θ' and released, the pendulum will swing back and forth with periodic motion.

given inputs,

• L = 1 meter,
• m = 1kg,
• b = 0.05,
• g = 9.81m/s^2.

Initial conditions,

• Angular displacement = 0,
• Angular velocity = 3rad/sec,
• time t = 0,
• The motion is simulated between 0-20 sec.

MATLAB code :

To solve ODE a function will be called :

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

Main program to solve ODE of a pendulum :

%solving second order ODE of a pendulum

%Inputs
b = 0.05
g = 9.81
l = 1
m =0.1

theta_0 = [0;5]

t_span = linspace(0,20,500)

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting
figure(1)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
title('Damped harmonic motion of pendulum')
legend('Angular displacement','Angular velocity')
xlabel('time')
ylabel('plot')

%Animation of the pendulum
ct = 1;
for i =1:length(results(:,1))
    x_start = 0
    y_start = 0
    x_end = 1*sin(results(i,1))
    y_end = -1*cos(results(i,1))
    
    figure(2)
    plot([-1 1],[0 0],'linewidth',5)
    hold on
    plot([x_start x_end],[y_start y_end],'linewidth',5)
    hold on
    
    plot(x_end,y_end,'o','markersize',20,'markerFacecolor','y')
    title('Animation of transient behaviour of a pendulum')
    axis([-1.5 1.5 -1.5 1])
    pause(0.03)
    hold off
    
    M(ct) = getframe(gcf)
    ct = ct+1
end

movie(M)
videofile=VideoWriter('Pendulum Animation.avi','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Figure :

Conclusion :

The x_start, y_start and x_ end, y_end values are

Animation of the pendulum :

https://www.youtube.com/watch?v=ktJtWi2vTEE