To simulate the transient behaviour of a simple pendulum using MatLab

Aim
Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.

Objective
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" with respect to time. Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.
Our objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

Introduction

Pendulum

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. And we are going to use Differential Equation to solve.

where b = damping force
m = mass
g = gravity
l = length
here the equation are in second order Ordinary Differential Equation in Matlab it will not accept the values in that, so we need to do some changes in equation.Matlab will do calculation only with matrix manipulations.

Equation Derive

`theta = theta1`

`d_theta/dt = d_(theta1)/dt=theta2`

`(d^2theta)/dt^2 = (d^2theta)/dt^2=d/dt((d_theta1)/dt)=d/dt(theta2)`

`(d^2theta)/dt^2+b/m d_theta/dt+g/lsintheta=0`

`(d^2theta^1)/dt^2 + b/m(d_theta1)/dt + g/lsintheta1=0``

`(d^2theta^2)/dt^2 + b/m(d_theta2)/dt + g/lsintheta1=0`

`(d_theta1)/dt = theta2`

`(d_theta2)/dt =(-b)/mtheta2-g/lsintheta1`

`d/dt[(theta1),(theta2)]=[(theta2),(-b/mtheta2-g/lsintheta1)]`

where `(theta1) & (theta2)` are taken as angular displacement because the action of moving a part or object from its position or from orgin.

where in right side matrix will be a angular velocity with rate of change of position with respect to time.

Angular displacement

A body is the angle in radians (degrees, revolutions) through which a point revolves around a centre or line has been rotated in a specified sense about a specified axis. When a body rotates about its axis, the motion cannot simply be analyzed as a particle, as in circular motion it undergoes a changing velocity and acceleration at any time (t). When dealing with the rotation of a body, it becomes simpler to consider the body itself rigid. A body is generally considered rigid when the separations between all the particles remains constant throughout the body's motion, so for example parts of its mass are not flying off. In a realistic sense, all things can be deformable, however this impact is minimal and negligible. Thus the rotation of a rigid body over a fixed axis is referred to as rotational motion.

Angular Velocity
In physics, angular velocity refers to how fast an object rotates or revolves relative to another point, i.e. how fast the angular position or orientation of an object changes with time. \
There are two types of angular velocity:

1. orbital angular velocity
2. spin angular velocity

Orbital Angular Velocity

Orbital angular velocity refers to how fast a point object revolves about a fixed origin, i.e. the time rate of change of its angular position relative to the origin. orbital angular velocity which depends on the choice of origin.

Spin Angular Velocity

Spin angular velocity refers to how fast a rigid body rotates with respect to its centre of rotation. Spin angular velocity is independent of the choice of origin.

In our case it is Orbital angular velocity for simple pendulum.

Input the value in matlab.

My Inputs

L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2.
angular displacement=0
angular velocity=3 rad/sec
t=0

Program

Calling Function

In calling function we have programed to run the formula which

`d/dt[(theta1),(theta2)]=[(theta2),(-b/mtheta2-g/lsintheta1)]`

left side equation are in left and similarly for other equation.

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

close all
clear all
clc

%Input
l = 1;
m = 1;
b = 0.05;
g = 9.81;

% at initial condition
theta_0 = [0;3];

%time value
t_span = linspace(0,70,500);

%solver ODE
[t, result] = ode45(@(t,theta) function_pendulum(t, theta, b,g,l,m),t_span,theta_0);
%plotting
figure(1)
plot(t,result(:,1))
hold on
plot(t,result(:,2))

ct = 1
for P = 1:length(result(:,1))
    Length = result(:,1);
  
    x0 = 0;
    y0 = 0;
    x1 = l*sin(result(P,1));
    y1 = -l*cos(result(P,1));
    
    figure(2)
    plot([-1.5 1.5],[0,0],'linewidth',5,'color','r')
    grid on
    hold on
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b')
    plot(x1,y1,'o','MarkerSize',10,'MarkerFaceColor','g')
    hold off
    axis([-1.5 1.5 -1.5 0.4])
    pause (0.03)
    
   M(ct) = getframe(gcf);
 
    ct = ct + 1;
end

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

At first we are defining with initial conditions,then defining with time interval from 0 to 70 with eqaual 500 range in between that using

Linspace(0,90,500);

Then calling the function to solve value with ode45 which is pre defined in Matlab for solving intergration values.

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

ode45

Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs (t,y) even if one of the inputs is not used. but we have value more than two input so we define our formula inside the ode45 with respect to change in time and value of `theta`.

ploting that value using plot command,

Animating

First we need to use counter before for loop, then following with for loop we need to define result value (Angular displacement) there only we will notice change in position so we have taken that value.

Then ploting that with change in time with respect to position.

Caputring each position using

M(ct) = getframe(gcf);

then incrementing the counter values till n time.

Using movie(m) command creating a movie with all the frame.

Result

The graph shown that show down of its angular displacement with respect to number of time it swings.

Animated Video Link

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

Conclusion

Program has been done to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.