Solving second order ODE of Simple Pendulum with Damping System
SECOND ORDER ODE’s – SIMPLE PENDULUM WITH DAMPING
AIM :-
To create a program for solving the ODE that represents the equation of motion of a simple pendulum with damping.
PLATFORM USED :-
- Software – MATLAB R2019a
- Programming Language – Python
SYSTEM SPECIFICATIONS :-
- Manufacturer - DELL
- Model – Vostro 3478
- Operating System - Windows 10
- Processor - Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz 1.80GHz
- RAM – 8 GB
GOVERNING EQUATIONS :-
- Equation of motion of a simple pendulum with damping,
`(d^2θ)/(dt^2) + \frac{b}{m} \(dθ)/(dt) + \frac{g}{L} sin\theta = 0`
Where,
b - Damping Co-efficient [No Unit]
m - Mass of the ball [kg]
g - Gravity [`m/s^2`]
L - Length of the pendulum [m]
DESCRIPTION :-
INTRODUCTION :-
Initially we need to discrete the order ODE into order ODE’s for easy calculation of the ODE using ODE45 solver. Then using the input parameters (damping coefficient, gravity, length of the pendulum and mass of the bob) calculate the displacement and velocity w.r.t angle . Next plot the graph between time period and velocity.
LAW’s ASSOCIATED :-
The laws associated with this order ODE are,
- Newton’s II Law
- Law of Length
- Law of Mass
- Law of Acceleration
- Law of Amplitude
NEWTON’S II LAW :-
It states that the force (F) is equal to mass (m) times acceleration (a).
`F = m*a`
LAW OF LENGTH :-
It states that the time period (T) of simple pendulum is directly proportional to square root of it’s length (L).
`T α \sqrt L`
LAW OF MASS :-
It states that the time period (T) of simple pendulum is independent of mass and material of the bob.
LAW OF ACCELERATION :-
It states that the time period (T) of simple pendulum varies inversely as the square root of acceleration due to gravity (g).
`T α \sqrt\frac{1}{g}`
LAW OF ISOCHROISM :-
The period of a simple pendulum does not depend on the amplitude of oscillations, if the amplitude is small.
SIMPLE PENDULUM AND SPRING MASS DAMPER SYSTEM :-
SIMPLE PENDULUM :-
The equation that is associated with motion of simple pendulum which is derived from Newton’s II law of motion is,
`\frac{d^2\theta}{{dt}^2} + \frac{g}{L} sin\theta - - - - - - 1`
SPRING MASS DAMPER SYSTEM (INSTANTANEOUS VELOCITY) :-
The equation associated with spring mass with damping system based on Newton’s II law of motion is,
`\ \frac{d^2x}(\dt^2)\ + \frac{b}{m}\ \frac{dx}{dt} + \frac{k}{m}x = 0 - - - - - - 2`
Where,
x(t) - Position of mass at time t
b - Damping Coeffiicient
k - Spring Constant
m - Mass
By combining both the equations 1 and 2 the ODE with simple pendulum damping system is derived,
`(d^2θ)/(dt^2) + \frac{b}{m} \(dθ)/(dt) + \frac{g}{L} sin\theta = 0 - - - - - - A`
DISCRETIZATION OF ODE :-
We need to convert the 2nd order ODE to two 1st order ODE’s for calculation.
Initial condition is `\theta\ =\ \theta1,`
`\frac{d\theta}{dt}\ =\ \frac{d\theta1}{dt}\ =\ \theta2 - - - - - - 3`
`\frac{d^2\theta}{dt^2}\ =\ \frac{d}{dt}\ (\frac{d\theta}(dt))\ =\ \frac{d}{dt}\ (\theta2)\ =\ \frac{d\theta2}{dt} - - - - - - 4`
Substitute equations 3 and 4 in A,
`\frac{d^2\theta}{dt^2} + \frac{b}{m} \frac{d\theta}{dt} + \frac{g}{L} sin\theta1 = 0`
`\frac{d\theta2}{dt} + \frac{b}{m} \theta2\ + \frac{g}{L} sin\theta1 = 0`
`\frac{d\theta2}{dt} = -\ \frac{b}{m} \theta2\ - \frac{g}{L} sin\theta`
1st order ODE’s,
`\frac{d\theta1}{dt}\ =\ \theta2 - - - - - - 5`
`\frac{d\theta2}{dt} = -\ \frac{b}{m} \theta2\ - \frac{g}{L} sin\theta1 - - - - - - 6`
ODE45 SOLVER :-
ODE45 is an inbuilt solver for ODE’s (Ordinary Differential Equation) in MATLAB. There are many methods to use ODE45 solver in MATLAB but for solving this order ODE we are using “Pass Extra Parameters to ODE” method.
General code for using ODE45,
[t,y] = ode45(odefun,tsapn,y0)
Where,
t - Output parameter for storing the t_span values of array
y - Output parameter for storing displacement & velocity values
odefun - Function that stores the order ODE’s
tspan - Time period of pendulum ( Integration of ODE’s)
y0 - Initial Condition for solving the order ODE’s
PASS EXTRA PARAMETERS TO ODE FUNCTION :-
The code for solving the ODE’s using ODE45 solver by passing extra parameters are,
[t,y] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), tspan, theta)
Where,
t - Output parameter for storing the t_span values of array
y - Output parameter for storing displacement & velocity values
odefun - Function that stores the order ODE’s
tspan - Time period of pendulum ( Integration of ODE’s)
theta - Main Inputs for the ODE derivatives, (`{\frac{d\theta}{dt} }` - dtheta_dt)
(b,g,l,m) - Extra Parameters for solving in ODE’s
INPUT AND OUTPUT PARAMETERS :-
INPUT PARAMETERS :-
- Gravity (g) = 9.81 (`m/s^2`)
- Length of the pendulum (l) = 1 m
- Mass of the bob (m) = 1 kg
- Damping Coefficient (b) = 0.05
- Time Span (t_span) = [0 20] seconds (s) – (Condition for Integration of ODE’s)
- Initial Condition (theta) = [0 ; 3] – [Angular Displacement ; Angular Velocity]
- Time (t) = 0 (For solving the ODE right-hand side value)
OUTPUT PARAMETERS :-
- Time array (t)
- Agular Displacement and Angular Velocity (y) w.r.t “theta”.
CODING :-
CODE FOR SOLVING THE ODE BY ODE45 SOLVER :-
CODE :-
% Clearing all the previous memory in the Worksapce
clear all
% Terminating all the previous plots
close all
% Clearing all the values in Command Window
clc
% Input Parameters for passing in the ODE's
b = 0.05; % Damping Co-efiicient
g = 9.81; % Gravity of Earth
l = 1; % Length of Pendulum String
m = 0.2; % Mass of the Pendulum Bob
% Initial Conditions for solving 1st order ODE's "(dtheta1_dt ; dtheta2_dt)"
theta = [0;3];
% Time Interval - Initial conditions for performing Integration
tspan = linspace(0,20,200);
% ODE45 Solver
% Output Varaibles [t,y]
% ode45(@Defining the values for "dtheta_dt" differentiation)
% Calling the function to solve the 1st order ODE's
% Inputting the values for performing Integration and intial conditions
[t,y] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),tspan,theta);
% Plotting
% Generating seperate (first) plot for pendulum motion
figure(1)
% Plotting time and angular displacement values
plot(t,y(:,1),'linewidth',1.5);
% Retaining the plots
hold on
% Plotting time and angular velocity values
plot(t,y(:,2),'linewidth',1.5);
% Labelling X and Y axes
xlabel('Time (s)')
ylabel('Angular Displcement (rad) & Angular Velocity (1/s)')
% Labelling the values for easy understanding
legend('Angular Displacement','Angular Velocity')
% Assignig the title for the plot
title('Motion of Pendulum w.r.t Time')
% Initializing a variable 'ct' for generating output figures seperately
ct=1;
% Looping Statement
% Looping statement for plotting the axes of pendulum
% Generating index from the output (1st column) and storing it in 'i'
for i = 1:length(y(:,1))
% Assigning the generated index numbers to a variable 'time'
time = y(i,1);
% INPUTS
% Initializing the value for X and Y axes co-ordinates which is fixed
x0 = 0;
y0 = 0;
% Calculating the co-ordinates for angular displacement of bob
x1 = l*sin(time);
y1 = -l*cos(time);
% Plotting
% Generating seperate (second) plot for pendulum Image
figure(2)
% Length of the string
plot([x0 x1],[y0 y1],'color','b','linewidth',2);
% Retaining the plots
hold on
% Fixed End of the Pendulum
plot([-0.5 0.5],[0 0],'linewidth',2,'color','black');
% Pendulum Bob
plot(x1,y1,'o','MarkerSize',20,'MarkerEdgeColor','r','MarkerFaceColor',[1 0.6 0.6]);
% Stop retaining the plots
hold off
% Fixing axes for the plots
axis([-1.5 1.5 -1.5 0.3])
% Assignig the title for the plot
title('Motion of Simple Pendulum')
% Capturing all the individually generated figure frames into an array ‘M’
M(ct) = getframe(gcf);
% Incrementing the plotted figures to generate in a separate avi file
ct=ct+1;
% End of for loop
end
% Generating a movie with all the captured frames
movie(M);
% Create a new AVI file to include all the captured frames (uncompressed AVI file)
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI');
% Opening the video file that is generated
open(videofile);
% Writing all the separately captures figures to the video file to generate a video
writeVideo(videofile,M);
% Closing the video file
close(videofile);
PROGRAM PROCESS :-
- The compiler starts the program.
- Then it checks for the input values for the parameters is in correct format.
- The time array is generated using linspace command between (0-20) seconds.
- Then it executes the ode45 solver line.
- It checks for the function where the 1st order ODE’s are present.
- Then the t_span and theta values are taken as input and then it solves the ODE’s.
- After solving the ODE’s the ode45 solver integrates the solutions generated.
- Then plotting the graph between angular displacement and time.
- Again plotting graph between angular velocity and time.
- Next the compiler executes the for loop statement.
- The loop counter ‘i’ starts generating values of angular displacement.
- Then the X and Y co-ordinates for plotting the pendulum is calculated.
- Now a separated graph is plotted for the string length and then the fixed end.
- Again the pendulum bob is being plotted in the same plots previously generated.
- The axis is defined for the final plots and then title for the plots are assigned.
- The generated plots are captured into an array ‘M’.
- Then the figures that are plotted gets incremented by value 1.
- The for loop ends.
- Generating new movie file of all the captured frames.
- The movie file name is assigned “simple_pendulum” and stored in “avi” format.
- The video file is generated and stored in the desired location.
- Compiler terminates the program.
ERROR ENCOUNTERED :-
- ERROR 1 – Included the ‘hold on’ command before plotting the pendulum and bob.
- ERROR 2 – Forgot to include a ‘-’ (negative) symbol for calculating the Y co-ordinate (y1) for the pendulum string length.
OUTPUT :-
- Graph were plotted between Angular Displacement and Angular Velocity vs Time array.
- Animation file of simple pendulum motion was generated.
ERROR 1 :-
ERROR 2 :-
GRPAH PLOTTED :-
CODE FOR THE ODE FUNCTION :-
CODE :-
% Code to calculate the values for the 1st Order ODE's
% Creating a function for ODE's
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
% Assigning the variables theta1 and theta2
theta1 = theta(1);
theta2 = theta(2);
% First ODE
dtheta1_dt = theta2;
% Second ODE
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
% Assigning both the ODE's variables to 'dtheta_dt' in array format
dtheta_dt = [dtheta1_dt ; dtheta2_dt];
% Ending of the function
end
PROGRAM PROCESS :-
- The compiler starts the program.
- Next it checks for the input parameters and values.
- The theta1 and theta2 variables are defined.
- The equations are assigned to the defined variables.
- The calculated values are assigned in array format.
- Finally compiler terminated the program.
FINAL OUTPUT :-
The final output is generated by combining both the programs and whenever the function name is called “ode_func” then the above mentioned function program starts compiling.
GENERATED VIDEO FILE :-
VIDEO FILE LINK :-
The link for video file - https://youtu.be/ZKpCT2KMD_s