Solving Second Order ODE for the Damped Oscillations of a Simple Pendulum

Aim: To write a program using MATLAB to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

Software used: Octave

Introduction:

In mechanics, simple harmonic motion (SHM), is a special type of periodic motion where the restoring force is directly proportional to the displacement and acts in the direction opposite to that of the displacement and also no other forces are involved.

A simple pendulum, in real-time conditions, when made to oscillate tends to dampen out after some time.

Newton's second law determines the movement of a pendulum. When this law is written down, we get a second order Ordinary Differential Equation which describes the position of the 'ball' w.r.t time.

Here, the position(theta) of the bob w.r.t time is determined by using the second-order ODE which is;

where, b= damping coefficient

g= gravity in m/s^2

m= mass of the ball in Kg

L= length of the pendulum in m

θ= angular displacemnet in radians.

Conversion of Second order ODE into First order:

Here, θ is the angular displacement. So, let θ be θ₁. Now, the differential of angular displacement will be angular velocity.

So, `(d θ₁)/dt= θ₂`

This, is the first differential equation.

The second one is obtained by differentiating θ₂ with respect to time and thereafter we get,

`(d θ₂)/dt=-(b/m).ϴ₂-(g/l).sin θ₁`

This is the second differential equation.

Both the first order differential equations can be now represented as a solution matrix.

In the matrix form `[■(θ₁@ θ₂)]` is the solution matrix, which can be addressed by a single letter ϴ.

Therefore,`[■(θ₁@ θ₂)]=ϴ` , in which

`θ₁=`angular position of the pendulum at various instances of time

`θ₂=`corresponding angular velocity at the same instances of time.

The function ‘ode45’, solves the input and returns two sets of values. One being the time inputs which has been given as input, and the other being the matrix `θ` as described previously.

The second output `θ` has two columns. The first column gives all the angular displacement `θ₁` in radians and the second column has angular velocity `θ₂`.

CODES:

1) Program to create a function that solves ODEs

function[dtheta_dt]=ode_func(theta,t,b,g,l,m)
b=0.05;
g=9.81;
l=1;
m=1;

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

2) Octave Code for generating an animation of a simple pendulum by solving second order ODE

close all
clear all
clc

%Defining the inputs or constants
b=0.05;
g=9.81;
l=1;
m=1;

%Initial conditions
theta_0=[0 3];  %Initial condition of [position;velocity]
t_span=linspace(0,20, 150); %Defining the time domain

%calling the inbuilt ode solver
[Results]=lsode('ode_func',theta_0,t_span);

%Angular Dispacment Plot
plot(t_span,Results(:,1),'color','b')
hold on

%Angular Velocity
plot(t_span,Results(:,2),'color','r')
hold off 

xlabel('Time')
ylabel('Plot')
pause()

theta=Results(:,1);

ct=1;
for i=1:length(theta)
    Theta=theta(i);
    
    x0=0;
    y0=0;
    x1=l*sin(Theta);
    y1=-l*cos(Theta);
    
    figure(ct)
    plot([-0.5 0.5],[0 0],'linewidth',5,'color','b')
    hold on
    plot(x0,y0,'.','color','b') 
    hold on
    axis([-1.5 1.5 -1.5 1.5])  %defining the axis
    
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g')
    hold on
    
    plot(x1,y1,'o','markers',20,'markerfacecolor','r')
    hold off
    
    pause(0.003)
    txt=sprintf('ode_function%d.png',ct);
    
    saveas(gca,txt);
    ct=ct+1;
    
end   

EXPLANATION OF THE CODE:

Firstly, the given details are assigned as variables in the code.Then time span of the simple pendulum is set by using linspace command then theta is set by using linspace command then theta value is given as a column matrix. After this we use ode45 function and the corresponding "ode_func" function call in main program. Ode45 is a solver used which will yield a much accurate value.

Then we use the plot command to plot the graphs and the legend command is used for identifying each value in the graph.

We make use of the looping statements, like for loop to generate an accurate simulation corresponding to the second-order ODE.

Now, we fix the coordinates for the pendulum.

Marker command is used for making bob of the pendulum. The axis is fixed using the axis command.

OUTPUT:

ANIMATION:

CONCLUSION:

Successfully plotted a graph of Time vs Plot for a simple pendulum and also generated an animation depicting the motion of a simple pendulum.