Second Order ODE Solving For Pendulum Motion

Aim: - A program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

Sub Objectives:-

• Write a program to solve the ODE
• Plot the displacement v/s time & velocity v/s time graphs
• Create an animation of the motion of the pendulum

Theory:-

In Mathematics, an ordinary differential equation(ODE) is a differential equation containing one or more functions of one indepndent variable and the deriatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.

Here, we are using the ODE that represents the motion of a simple pendulum with damping. The general equation for this is:

`(frac{d^2theta}{dt^2}) + (b/m)frac{d theta}{dt} + (g/L)sintheta = 0`............(1)

We are going to use this 2nd order ODE equation to find out the velocity and position of the pendulum ball with respect to time.

b = damping coefficient

m = mass of the body

g = gravitational force

L = length of the pendulum

Solving Methodology

We say `theta = theta1`

`frac{d theta}{dt} = frac{d theta1}{dt} = theta2`` .............................(2)

and

`frac{d^2 theta}{dt^2} = frac{d^2 theta1}{dt^2} = frac{d}{dt}(frac{d theta1}{dt}) = frac{d}{dt}(theta2)` .............................(3)

Substituting eq 2 and 3 in 1, we have

`frac{d theta2}{dt} + b/m (theta2) + g/L(sintheta1) = 0`

Converting the equations in matrix form

`((frac{d theta1}{dt}),(frac{d theta2}{dt})) = ((theta2),(-(b/m)theta2-(g/L)sintheta1))`

Code

ODE Function Program

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

Program to simulate pendulum motion and creating animation

% Program to solve the ODE and simulate pendulum motion
clear all
close all
clc

%Inputs
b = 0.05;  %damping coefficient value
g = 9.81; %gravitational value in m/s^2
l = 1;    %length of the pendulum in m
m = 1;   %mass of the body in kg

theta_0 = [0;3];  %Initial condition for displacement and velocity
t_span = linspace(0,20,400);  %Time span points

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

%Plotting
figure(1)
plot(t,results(:,1))
hold on
grid on
plot(t,results(:,2))
ylabel('Plot')
xlabel('Time')
legend('Displacement','Velocity')

figure(2)
ct = 1;
for i = 1:length(results(:,1))
    THETA_1 = results(i,1);
    %Initial coordinates
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(THETA_1);
    y1 = -l*cos(THETA_1);
    
    %Plot
    plot([-0.5 0.5],[0,0],'color','g')
    line([x0 x1],[y0 y1],'color','r')
    hold on
    plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','b')
    grid on
    axis([-1.5 1.5 -1.5 1])
    hold off
    pause(0.0001)
    M(ct) = getframe(gcf);
    ct = ct + 1;
end
movie(M)
videofile = VideoWriter('Second_ODE_Pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Results

The displacement v/s time and velocity v/s time graphs were plotted.

The animation video for the pendulum motion was created and uploaded in youtube.

Link for video:- https://youtu.be/Z90RHyQIrpU