Solving 2nd Order ODE to determine the Motion of a Simple Pendulum

PURPOSE :

To write a program to solve the provided ODE which decribes the transient behaviour of a simple pendulum and to create an animation to display the motion of the pendulum

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

INPUTS :

g - gravity in m/s2 = 9.81 m/s2
L - Length of pendulum = 1 m
m - mass of bob = 1 kg
b - Damping Coefficient = 0.05

Time Span for motion = 0-20 sec
Initial Angular Displacement = 0
Initial Angular Velocity = 3 rad/sec

KEY POINTS :

* To Solve the 2nd order ode by using ode45 solver in matlab, it is required to split the provided ODE into two 1st order ODEs

Assume

`theta = theta_1`

`(d theta_1) /dt = theta_2`

Taking this as equation (i)

Substituting in provided 2nd order ODE :

`rArr (d^2theta_1)/dt^2 + b/m(d theta_1)/dt + g/l Sin theta_1`

`rArr (d theta_2)/dt + b/m theta_2 + g/l sin theta_1`

`rArr (d theta_2)/dt = -b/m theta_2 - g/l sin theta_1`

Taking this as equation (ii)

Writing the above 2 equations in matrix form :

`(d)/dt [[theta_1],[theta_2]] = [[theta_2],[-b/m theta_2 - g/l sintheta_1]]`

from the above equations, it can be noted that, the first equation solves for angular displacement while the second equation solves for angular velocity

* To Determine the positions of the pendulum wrt change in time

NOTE : Angular displacement is obtained from the solved ODE

Substituting the values of angular displacement in the co-ordinate shown above and plotting the points obtained with respect to a fixed origin point will yield the desired animation showing the motion of the simple pendulum.

PROGRAM :

%Program to solve ODE representing Equation of Motion of Simple Pendulum
clear all
close all
clc

%Input Data to the ODE Function
b = 0.05; % Damping Coefficient
g = 9.81; % Acceleration due to gravity in m/s2
l = 1; % Length of Pendulum String in m
m = 1; % Mass of Ball in kg

% Initial Angular Displacement = 0 rad
% Initial Angular Velocity = 3 rad/sec
% Initial Boundary Condition of Pendulum
theta_0 = [0;3];
% Time Frame
t_span = linspace(0,20,300);

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

%Plotting the Graphs and Creating the animation of Pendulum Motion

figure(1) %Rate of change of Displacement and velocity with respect to time
plot (t,results(:,1))
hold on
plot (t,results(:,2))
xlabel(\'Time (sec)\')
ylabel(\'Plot\')
axis ([0 20 -3.8 3.8])
legend(\'Angular Displacement (rad)\',\'Angular Velocity (rad/sec)\')
title(\'Variation Of Initial Parameters with Respect to Time\')

figure(2) %Angular Displacement vs Angular Velocity
plot(results(:,1),results(:,2))
xlabel(\'Angular Displacement (rad)\')
ylabel(\'Angular Velocity (rad/sec)\')
title(\'Angular Displacement vs Angular Velocity\')

figure(3) %Animating the Pendulum Motion
ct = 1;
pink = [1,0.2,0.7];
for i = 1:length(t_span)
    Amgular_Displacement = results(i,1);
    % Co-ordinates of Pendulum ends
    x0 = 0;
    y0 = 0;
    x1 = l*sin(Amgular_Displacement);
    y1 = -l*cos(Amgular_Displacement);
    %Plotting Pendulum String
    plot([x0 x1],[y0 y1],\'linewidth\',5,\'color\',\'b\')
    title(\'Motion Of Pendulum\',\'fontsize\',16)
    hold on
    %Plotting Pendulum Bob
    plot ([-1 1],[0,0],\'-.\',\'linewidth\',3,\'color\',\'k\')
    plot(x1,y1,\'o\',\'MarkerSize\',20,\'MarkerFaceColor\',pink,\'MarkerEdgeColor\',pink)
    axis ([-1.5 1.5 -1.5 0.5])
    pause (0.005)
    hold off
    M(ct) = getframe(gcf);
    ct = ct+1;
end

%Saving the Animation File
movie(M)
videofile = VideoWriter(\'ODE_Pendulum.avi\',\'uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile)  

FUNCTION USED :

A function was used to define the 2nd order ODE as two 1st order ODEs to be used in ODE45 solver.

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

RESULTS :

* Plotting the Variation of Angular Displacement and Angular Velocity with respect to time

* Plotting Variation of Angular Displacement with respect to Angular Velocity

* Animation to show the motion of pendulum within the time span

LINK :

https://youtu.be/TpNJKELUWzg

EXPLANATIUON OF THE CODE :

* Initially to clear any previously defined values, close any open windows and clear the command window, commands are provided

* Following that, all input parameters which are known are provided, Here, to split the value of the Time span between a defined range into an array with equal intervals between the values, Linspace command is used.

* Also, since we know that the result matrix is a matrix of 2 columns in which the first column defines the angular displacement and the second column defines the angular velocity, the initial conditions defined under theta_o are defined in same manner, a matrix of order (1x2)

* To solve the ODE, an inbuilt solver \'ODE45\' is used which calls the user-defined function and solves it with respect to the inital boundary conditions and time span

* Once the result has been computed, three graphs are plotted using the plot command

* To make the plots easier to understand, commands like \'xlabel\', \'ylabel\', \'legend\', \'title\' are used to provide headings and labels

* The concept of counter \'ct\' is introduced to aid in segregation of frames and an array \'pink\' is created to be used in subsequent plots for displaying the color of the same name.

* \"For Loop\" is used to iterate the different positions resulting due to varying values of Angular Displacement

* The position of the two ends of the pendulum are defined as co-ordinates (x0,y0) and (x1,y1) and are plotted

* To plot the pendulum bob, the plot command is used along with a marker to define a bob at all points of (x1,y1). The size and color of the marker is also controlled with respective commands

* The refresh rate is provided with the pause command

* The hold command is used to overlap plots and is switched off before refreshing to avoid overlap while displaying varying positions.

* To create an animation, an array M is created to store the values of the current frame for varying counter \'ct\' values. Then the array M is made a movie

* This movie is now written to a videofile after creating it and opening it and once written is closed.

* The video file is created in the same location as the program.

ERRORS FACED :

* Use of \'Results\' instead of \'results\'

CORRECTION : identified by re-checking

* Use of \'pink\' instead of pink to define the color

CORRECTION : identified by re-checking, since pink was a defined color rather than in-built color, quotes should not be used