Solving a second order ODEs for a Simple Pendulum

Aim

To Simulate a pendulum motion and to solve the equation of the motion of a simple pendulum with damping by using a given parameters on a MATLAB

Simple 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. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth.The oscillations will be reduced due to resistance for the motion by external force that is called Damping.Damping Force will drecreases a angular velocity of the simple pendulum

The motion of a simple pendulum with damping can be studied using second order ODE and the equation used to represent the position of the bob with respect to the time in the real situations

`theta`=Angular displacement in radians

g = gravity in `frac{m}{s^2}`

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To calculate a position of bob using angular displacement

X-axis cordinate=`L cdot cos (theta)`

Y-axis cordinate=`L cdot sin(theta)`

MATLAB program

clc
clear all
close all

%inputs
L=1;     %in Metre
m=1;     %in Kilograms
b=0.05;  %Damping co-efficient
g=9.81;  % in Metre/sec^2
Angular_velocity=3; %in rad/sec
ct=1     %frame counter

%Limits
t_limits=linspace(0,20,600)

%intial condition
theta=[0;3];

%ode solver
[t,results]=ode45(@(t,theta)odesolver(t,theta,b,m,g,L),t_limits,theta)

figure(1)
plot(t,results(:,1),t,results(:,2))
xlabel('time')
legend('Angular Displacement','Angular Velocity')


%To Create animation
figure(2)
 for i=1:1:600
     
    figure(2)
     y=-L*cos(results(:,1))%to get the position using angular displacement
    x=L*sin(results(:,1))
     plot([-3 3],[0 0],[0 x(i,1)],[0 y(i,1)],x(i,1),y(i,1),'o','MarkerSize',20)
     axis([-1.5 1.5 -1 0])
     grid on
     e(ct)=getframe(gcf)
     ct=ct+1
 end
 
%these commands is used to create a video file 
 movie(e)
 videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI')
 open(videofile)
 writeVideo(videofile,e)
 close(videofile)

The Function used for the above program are given below

function [dtheta_dt]=odesolver(t,theta,b,m,g,l)
theta1=theta(1)
theta2=theta(2)
dtheta1_dt= theta2                               %angular velocity = differation of angular displacement with respective time
dtheta2_dt= (-b.*(theta2)./m)-(g.*sin(theta1)./l)% gives a differention of angular velocity with respective time
dtheta_dt=[dtheta1_dt;dtheta2_dt]
end

Code Explanation

• Clc,clear all,close all commands are used to clear a command window,clear all the variables in the workspace,closes all the figure window
• variables are assigned with predescribed values
• mentioned a SI units as a comment
• we have declared a limit for integration that would use by a ODE45 solver
• Initial conditon declared under the name of theta array
• Theta array consist of intial angular displacement on first row and angular velocity on second row

Conversion of 2nd order ODE into two first order ODE

In the matlab ODE45 solver will not able to solve a 2nd order ODE

There are two Initial Conditions

• Angular velocity
• Angular Displacement

Assuming that angular displacemet `theta` would be a `underset(1)(theta)`

`frac{dunderset(1)(theta)}{dt}=underset(2)(theta)`

`underset(2)(theta)` represents a angular velocity

Using the above equation on the Second order ODE

`frac{dunderset(2)(theta)}{dt}=frac{-b cdot underset(2)(theta)}{m}-frac{g cdot sin(underset(1)(theta))}{L}`

Both the 1st order differential equation assinged to the cloumn vector

`frac{d }{dt}[underset(underset(2)(theta)) (underset(1)(theta))]=[underset(frac{-b cdot underset(2)(theta)}{m}-frac{g cdot sin(underset(1)(theta))}{L})(underset(2)(theta))]`

where `theta=[underset(underset(2)(theta))(underset(1)(theta))]`

above described will be used by ODE45 solver as a Expression

• In the Ode45 solver it takes a input as theta and t and expression would next to the input get through a function call and limits and initial conditions are mentioned
• On a function angular displacement and angular velocity are assigned to two different variables
• Two expression is assigned to the two different variables will give the varibles back to the ode45 solver
• A new figure is created.Plot command is used to obtain a graph angular velocity and angular displacement against a time
• For loop is used calculate a x and y cordinates of a bob using an 600 different angular displacement values and to plot a bob position on the new figure window using a plot command
• Each and every time in a for loop will display a one bob postion
• Each and every bob position would capture a data from a figure window by using getframe command and store in a array
• movie command is used to make a movie from a frames
• VideoWriter command is used to create a new video file
• open,writeVideo,close these commands used to open the video file , using a movie object to copy the data to video file and close the video file

Errors

To overcome this error,removed an apostrophe for the number which represents a marker size

To overcome this error,removed a comma in a ode45solver next to the inputs of anonymous function

Plot

Conclusion

Solved an ode of equation of motion of simple pendulum with damping using ODE45 solver and from the plot interpreted is Angular velocity and angular displacement is diminishing with increasing in time because of damper

The animation can be found at the below link

https://youtu.be/speVnHKjmZ0