MOTION SIMULATION OF PENDULUM USING SECOND ORDER ODE

Aim:

To write a program to define the pendulum movement with given input parameters angular displacement and angular velocity with time 0-20 secs.

Introduction to ODE:

In Engineering, an Ordinary Differential Equation (ODE) is used to describe the transient behavior of a system. A simple example is a pendulum.

The way the pendulum moves depends on Newton's second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

The input parameters given are:

g= 9.81 m/s^2

l = 1m (length of the pendulum in metre)

m = 1 kg; (mass of the ball in kg)

b =0.05 (damping coefficient in Ns/m)

Governing equation:

To write the second ODE program, the below governing equation is used:

(d^2 θ)/(dt^2 )+b/m*dθ/dt+g/l sinθ=0

lets say that,

θ= θ_1

dθ/dt= θ_2

(d^2 θ)/dt^2 =(d^2 θ_1)/dt^2 = d/dt ((dθ_1)/dt)= (d(θ_2))/dt

(d^2 θ_1)/dt^2 + b/m*(dθ_1)/dt+ g/l*sinθ_1

(dθ_2)/dt= -b/m*θ_2- g/l sinθ_1

d/dt [■(θ_1@θ_2 )]=[■(θ_2@-b/m θ_2-g/l sinθ_1 )]

Explanation:

the ode45 is an inbuilt solver in the Matlab.usually, they are used to Solve nonstiff differential equations.

syntax: [t,y] = ode45(odefun,tspan,y0,options)

'odefun' is used to define or call the function for solving the ordinary differential equation.

this solver uses the integration settings defined by options.

the y0 takes up the initial conditions, here we have conditions as angular displacement = 0 and angular velocity = 3rad/sec at time t = 0.

we take t_span within the range of 0-20 secs with 500 equally spaced values.

since 'y' is an array of 500 rows and 2 columns, graph is plotted with time and y values and the plot is saved using 'saveas' command in JPG format.

a loop is executed to extract values from the 'y' array to measure the position of the ball for time 't'. Marker command is used to define a ball and size is kept as 15.

the 'getframe(gcf)' command is executed to store the frame or plot for making the movie.

last using 'VideoWriter' command, the videofile is created on the name 'second_ode_pendulum.mp4' in MPEG format.

Body of the Code:

clear all
close all
clc

b = 0.05;                                         %Damping coeeficient(Ns/m)
g = 9.81;                                         %Acceleration due to gravity(m/s^2)
l = 1;                                            %length of pendulum(m)
m = 0.1;                                          %Mass of the ball(kg)

%initial conditions
theta_0 = [0;3];                       %angular displacement = 0; angular velocity = 3

t_span = linspace(0,20,500);           %time is ranged btw 0-20

%calling function ODE
[t, y] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span,theta_0);   %calling ODE function
figure(1)
p1 = plot(t, y(:,1),'linewidth',2)
hold on
p2 = plot(t, y(:,2),'linewidth',2)
xlabel('Time')
ylabel('Angular displacement and Angular velocity')
legend('Angular displacement','Angular velocity')
title('Pendulum motion')
saveas(figure(1),'Pendulum motion','jpg')
%counter for movie making
ct=1;
%creating a for loop for oscillating the pendulum

for i=1:length(y)
    THETA = y(i,1);
    
    %initial boundary conditon
    x0=0;
    y0=0;

    %Defining x and y coordinates
    x1=  l*sin(THETA);
    y1= -l*cos(THETA);
    
    figure(2)
    
    hold on
    plot([x0 x1],[y0 y1],'linewidth',2);
    plot(x1,y1,'-o','markersize',15,'markeredgecolor','r','markerfacecolor',[0.6 0.6 0.6]);
    plot(x0,y0,'-*','markersize',1,'markeredgecolor','b','markerfacecolor',[0.6 0.6 0.6]);
    title('Pendulum motion')
    axis([- 1 1 -1 0]);
    pause(0.03);
    
    M(ct)=getframe(gcf)
    ct=ct+1;
    clf()                     %clf() code is  used for clearing the current graphic image
end
%%creating the animation and movie(MPEG-4)

movie(M);
videofile = VideoWriter('second_ode_pendulum.mp4','MPEG-4');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function for Ordinary differential equation:

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

link for the motion of pendulum:

https://youtu.be/bg_UH7o0AIQ

Plot for the Angular velocity/ Angular displacement vs time:

Results:

thus, the motion of a pendulum has been executed with Matlab programming with 'ode45' function. also, the movie is created for the same in MPEG format.

from the plot, we can see that the angular velocity and angular displacement becomes narrow after some time which means the system is damped.