Solving ODE for a Simple Pendulum

TRANSIENT BEHAVIOUR OF A PENDULUM

AIM: A MATLAB Program to solve an ordinary differential equation of motion of a pendulum

OBJECTIVES:

1. Solve an ODE which represents the equation of motion of a simple pendulum with damping.
2. Understand the transient behaviour of a pendulum
3. Factors affecting the displacement and velocity of the given motion

EQUATIONS USED:

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

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

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

`d/dt [[theta_1], [theta_2]] = [[theta_2], [-b/m theta_2 - g/lsintheta]]`

FUNCTIONS:

• ``(∂θ_1)/(∂t)=θ_2``
• `(d^2theta)/dt^2 + b/m (d theta)/dt + g/l sintheta = 0`

DETAILS:

Programming Language used: MATLAB v2020a

Concept Used: ordinary diffrential equations and pendulum motion

Loops used: For loop and If-else (nested)

Errors Faced: "Error Using Movie: Movie contains uninitialized frame"

Methods to Overcome Error: increment of counter should be decleared after the getframe cmd

Problems Faced: multiple plots are created if we use "line' cmd instea of 'plot' thus making it difficult to analyze

THEORY:

• A simple pendulum consists of a metal bob hung by a string connected to a pivot point. The bob oscillates in to and fro motion because of its mass and momentum.
• The movement from mean to one extreme and then to the other extreme and again back to the mean is called a complete oscillation.
• The number of oscillation the pendulum completed in one second is its frequency and the time required for one complete oscillation is called the time period.

CODE USED AND STEP-WISE EXPLAINATION:

clear all
close all
clc
% inputs
b = 0.05;           % damping co-efficient
g = 9.81;           % gravitational const. in m/s^2
l = 1;              % gives length of string in m
m = 0.1;            % gives mass og the bob in kg
% initial condition
theta_0= [0;3];     % giving initial values of the angle 
% time_plot
t_span = linspace(0,20,500);  % creates a array from 0 - 20 of 500 numbers
% solving ode by calling the inbuilt function and a madeup function 'ode_func'
[t, results] = ode45(@(t,theta) ode_func(t, theta,b,l,g,m),t_span,theta_0);   

% plotting graph for angular volocity and angular displacement
figure(1)               % to plot disp vs velocity in figure 1
plot(t,results(:,1))    % plotting time vs all the rows of first column 
hold on                 % the initial graph is held
plot(t,results(:,2))    % plots time vs all the the rows of second column results
legend('angular disp.', 'angular vel.')     % assigns legend to the graph
axis([0 20 -4 4])       % defines the limit of the x and y axis
xlabel('time')          % names the x axis   
ylabel('plot')          % names the y axis
title('Curve characterictics of simple pendulum')      % gives title to plot
hold off                % releasing the plot

% plotting morion of the pendulum
ct = 1;                      % initial value of counter
figure(2);                   % creating a figure for animation
for i = 1:length(t_span)     % this loop iterates through each value of t-span
    % initial cordinates
    x0 = 0;     
    y0 = 0;
    % final cordinates
    x1 = l*sin(results(i,1));      % gives final cordinates wrt to x axis
    y1 = -l*cos(results(i,1));     % gives final cordinates wrt to y axis
    
    plot([-1 1],[0 0],'k','LineWidth',10);  % plots a graph of given specifications
    hold on
    plot(0,0,'wo');                         % initial base line for the support
    plot([0 0 ],[-1.3 0.15],'g--','LineWidth',2);% represnt final position og pendulum
    % pendulum string
    plot([x0 x1], [y0 y1],'k','LineWidth',3);  % creates a string
    % pendulum bob
    plot(x1,y1,'ro','MarkerFaceColor','b','MarkerSize',20); %creates bob of given size
    title('A Simple Pendulum');             % gives title to plot
    text(-0.4 ,-1.5,['time =', num2str(t_span(i)),'Sec' ],'fontsize',16);
    text(-1 ,-1.8,['Aungular Vel =', num2str(results(i,2)),'rad/Sec' ],'fontsize',16);
    axis([-1.5 1.5 -2.3 0.5 ]);         % gives limit to the axis
    if i == 1        
        pause(2);                      % sets an if loop for pausing the animation
    else
        pause(0.0001);                 % to get convinient view to the user
    end
    M(ct) = getframe(gcf);             % stores the frames in an array 'M'
    hold off;
    ct = ct+1;                         % increments the counter by one
end

% making a movie from frames
movie(M);                              
% create a video file and to save the same
videofile = VideoWriter('Pendulum2.avi','Uncompressed AVI');
open(videofile);
% writing the movie in videofile
writeVideo(videofile,M);
close(videofile);

Function Code:

function[dtheta_dt] = ode_func(t,theta,b,l,g,m)  % inputs of the function
    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

% Refer the function used under equations section to get the IDEA

RESULTS:

The graph of angular displacement and angular velocity vs time is obtained as:

The animation of Simple pendulum :

https://youtu.be/p9fKzbyyuJA

ERRORS FACED:

SOLUTION: the output dtheta_dt must be a row vector not a column vector:

dtheta_dt = [dtheta1_dt ; dtheta2_dt];

THE rest of the errors were typos and hence corrected using manual debugging

CONCLUSION:

1. From the above program we could analyse the motion od a simple pendulum using ODE
2. We could also study the factors affecting displacement and velocity of the pendulum
3. We also understand the variation in displacement and time with respect to the amplitude of the bob