AIM:

To write a MATLAB code for simulating the transient behavior of a simple pendulum and to create an animation of its motion.

SIMPLE PENDULUM:

A simple pendulum consists of a ball (point-mass) m hanging from a (massless) string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with the periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained 

`tau = I alpha qquad Rightarrow qquad -mg Lsintheta = mL^2 frac{d^2theta}{dt^2}`

and rearranged as

`frac{d^2theta}{dt^2} + frac{g}{L}sintheta = 0`

If the amplitude of angular displacement is small enough, so the small-angle approximation  holds true, then the equation of motion reduces to the equation of simple harmonic motion

 `frac{d^2theta}{dt^2} + frac{g}{L}theta = 0`

The simple harmonic solution is

`theta(t) = theta_o cos(omega t) `

where `theta_o`is the initial angular displacement, and `omega = sqrt{g/L}`the natural frequency of the motion. The period of this system (time for one oscillation) is

`T = frac{2pi}{omega} = 2pisqrt{frac{L}{g}} .`

The period of a pendulum does not depend on the mass of the ball, but only on the length of the string. Two pendula with different masses but the same length will have the same period. Two pendula with different lengths will different periods; the pendulum with the longer string will have a longer period.

With the assumption of small angles, the frequency and period of the pendulum are independent of the initial angular displacement amplitude. A pendulum will have the same period regardless of its initial angle

REAL NONLINEAR SIMPLE PENDULUM:

When the angular displacement amplitude of the pendulum is large enough that the small-angle approximation no longer holds, then the equation of motion must remain in its nonlinear form

`frac{d^2theta}{dt^2} + frac{g}{L}sintheta = 0`

This differential equation does not have a closed-form solution but instead must be solved numerically using a computer. MATLAB numerically solves this differential equation very easily with the built-in function ode45.

DAMPED OSCILLATION OF SIMPLE PENDULUM:

Dampers are used in oscillating systems for primarily dissipating the energy of the system to 

• have a controlled motion or prevent oscillation
• to alter the frequency of oscillation which helps in avoiding resonance
• to reduce the dynamic stresses induced which boosts fatigue life.

By applying Newton's second law for rotational systems

 `Ialpha+bl^2omega+mglsintheta=0`

`ml^2(d^2theta)/dt^2+bl^2(d theta)/dt+mglsintheta=0` 

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

This is a differential equation of `2^(nd)` order. The inbuilt function used in this code 'ode45' cannot process a `2^(nd)` order ODE. Hence we convert this equation into two `1^(st)` order ODEs, and then feed them into a user-defined function 'ode_func' which is further served as an input to the inbuilt function ode45.

Breaking down into `1^(st)`order ODEs,

Assume `theta=theta_1` be a solution of the `2^(nd)` order ODE and,

`(d theta_1)/dt=theta_2`.........................................................................1

Now we convert the `2^(nd)` order ODE into  two `1^(st)` order ODEs with variables `theta_1 and theta_2`.

`implies``(d^2theta)/dt^2=(d^2theta_1)/dt^2=(d theta_2)/dt^`

substituting the above in,

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

`(d^2theta_1)/dt^2+b/m(d theta_1)/dt+g/lsintheta_1=0`

`(d theta_2)/dt+b/m(d theta_1)/dt+g/lsintheta_1=0`...................................2

LHS array of the below equation is the output of the user-defined function 'ode_func', which is inturn fed an as input to the inbuilt function 'ode45',

`[(d theta_1)/dt`       `(d theta_2)/dt]`=`[theta_2`        `(-g/lsintheta_1- b/mtheta1)]`

 The 'ode45' along with above matrix takes in time interval 't_0' array and the initial conditions of  angular displacement(in radians) and and angular velocity (in radians/second) in form of arrays and gives out a (lengthX2) matrix where length is the length of t_0 array.And the 2 colums contain angular displacement and angular velocity corresponding to the time interval element of 't_0'.

PROGRAM:

 MAIN PROGRAM:

clear all
close all
clc

%input parameters
l=1;
m=1;
b=0.05;
g=9.81;

%time interval taken as an array
t_0=linspace(0,20,300);

%initial conditions that contains (displacement in radians)
%and angular velocity in (rad/s)
x_0=[0;3];

% using the inbuilt function ode45 that uses time interval, the initial condition
%and user-defined function as its input.

%we have done function calling for the user define function while feeding
%it as an input to the ode45
%'result' is the final array, that ode45 gives us back which contains two columns
%first one has displacement and the second one has angular velocity after every
%0.06 seconds interval
[t,result]=ode45(@(t,x) ode_func(t,x,l,m,b,g),t_0,x_0);

% to ensure the velocity and displacement curves w.r.t time 
%comes  in the same figure we have used
%hold on
figure(1)

%plot the first column of the result array w.r.t time
plot(t,result(:,1));


hold on

%plot the second column of the result array w.r.t time 
plot(t,result(:,2));
legend('Angular displacement','Angular Velocity');
xlabel('TIME')
ylabel('AMPLITUDE')

%defining this variable 'gh' to access individual elements of array M 
gh=0;

%using for loop to access each element in the 'displacement column' of the array,'result'
%and assigning it to a variable theta
%theta is further used to locate the exact position of the pendulum at the
%end of each interval of 0.06 sec

    for j=1:length(result(:,1));
        theta=result(j,1);
        x=l*sin(theta);
        y=-l*cos(theta);
        figure(2)
         
         plot([-1 1],[0 0],'color','r','linewidth',3)
         hold on
         plot([0 x],[0 y],'linewidth',1)
         plot(x,y,'-o','Markersize',5,'Markerfacecolor',[0.1 0 0])
         
         
         axis([-2 2 -2 2])

         %using this pause to ensure that the displacement of the pendulum is in sync 
         %with their specific time intervals
         %there are 299 such time intervals
         %so a new plot is plotted for each time interval according to the
         %corresponding value in the 'result' arrays displacement column  
         pause(0.06)
           hold off
       gh=gh+1;
       % writing each frame of the plot to an array called M 
       M(gh)=getframe(gcf);
        
    end
    
%generating a movie with that array 'M'
movie(M)

%OPENING A NEW FILE NAMED Simple_Pendulum of FILE type Uncompressed AVI
f=VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');
open(f)

%writing the contents of array M into the NEW FILE
writeVideo(f,M)
close(f)

 USER DEFINED FUNCTION USED AS AN INPUT FOR ode45:

function [dx_dt]=ode_func(t,x,l,m,b,g)
x1=x(1);
x2=x(2);
dx1_dt=x2;
dx2_dt=(-b/m)*x2-(g/l)*sin(x1);
dx_dt=[dx1_dt;dx2_dt];
end

 RESULT:

The amplitude of the Angular displacement ( in radians ) and Angular Velocity(in radians/seconds) is plotted with respect to Time( in seconds).

The simulation of the damped oscillation is performed for the first 20 seconds.

YOUTUBE LINK:https://youtu.be/LnLj6WFym-s