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

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

In real oscillations, friction, or damping, slows the motion of the system. Due to frictional force, the velocity decreases in proportion to the acting frictional force. These oscillations are technically termed as Damped Oscillations. From Newton’s Second Law, on simplification, a second order differential equation is obtained as a function of displacement and angular velocity.

A simple pendulum in real conditions when made to oscillate, tends to dampen out after some time. The motion of the simple pendulum can be studied by a second order differential equation.

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

Here,

`theta` = angular displacement in radians

'b' = damping factor

`m` = mass of the pendulum

`g` = acceleration due to gravity

`l` = length of the pendulum

The above shown 2nd order differential equation gives the position of the pendulum bob at any instant of time with great accuracy. MATLAB can be used to calculate the oscillation of the pendulum, but it does not have a function to solve a 2nd order differential equation directly. So a 2nd order D.E. can be manipulated to form two 1st order D.E.'s and be solved. Matlab uses the function 'ode45' to solve a 1st order D.E.

Conversion of 2^nd order DE into two 1^st order DE.

Here, `theta` is the angular displacement. So, let `theta` be `theta_1`.

The differential of angular displacement is angular velocity. So,

`(d theta_1)/dt = theta_2`

This is the first differential equation.

The second one is obtained by differentiating `theta_2` with respect to time and we get,

`(d theta_2)/dt = -(b/m)*theta_2 - (g/l)*sintheta_1`

This is the second differential equation.

Both the 1st order differential equation can be represented as a solution matrix.

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

In the matrix form `[[theta_1],[theta_2]]` is the solution matrix, which can be addressed by a single letter, `theta`.

Therefore, `[[theta_1],[theta_2]]` = `theta`, in which,

`theta_1` = the angular position of the pendulum at various instances in time.

`theta_2` = the corresponding angular velocity at the same instances in time.

The function 'ode45' solves the input and returns 2 sets of values. One being the time inputs which has been given as input, and the other being the matrix `theta` as described previously. The second output `theta`, has two columns. The first column gives all the angular displacement `theta_1` in radians and the second column has angular velocity `theta_2`.

Structure of the Matlab Program

Part 1: Writing the main program

• Initialize the parameters of the coefficients like damping factor, acceleration due to gravity, length of the pendulum and mass of the bob. And dependent variables like time interval, initial values of position and initial velocity.
• Call the function 'ode45' and assign the output to two sets of values. In calling the function the main dependent variables are specified. And the function 'ode_function(t,theta,b,g,l,m)', time interval and the initial position and velocity are given as inputs. The function 'ode_function(t,theta,b,g,l,m)' returns a matrix of `theta_2` and `(d theta_2)/dt`.
• The function 'ode45' accepts the values from the inputs and gives two outputs in the form of a matrix calculated by integration. First, the time interval which was given as input. Second, a matrix of angular displacement `theta_1` and angular velocity `theta_2`.
• From the obtained values, a plot of displacement and velocity vs time is made which shows the rate of change of the quantities. The curve can be varied by changing the damping coefficient.
• Now, a counter is assigned to 1. A 'for' loop is opened to plot every position of the pendulum for the specified time interval. A pointer 'i' is assigned to the displacement column of the output matrix. So, whenever the 'for' loop increments the value of the pointer, the successive displacement value of the pendulum is obtained.
• Using trigonometry, the exact position of the pendulum bob is calculated and stored in temporary variables as coordinates.
• A plot of [X,Y] is done and a marker is assigned as per requirements to represent a bob of the pendulum. This plot is stored in an array of the counter, and then the counter is incremented and the loop is closed.
• Outside the loop a movie of the acquired plots is made by writing the following code.

Syntax:

movie(array name)
videofile = VideoWriter('Video name.format','type of file');
open(videofile);
writeVideo(videofile,array name);
close(videofile);

• movie(array name), is a function that takes array name as input and enables the formation of a video.

• VideoWriter('Video name.format', 'type of file'), is a function that takes video name.format and the type of file(i.e., compressed or uncompressed) as input and assigns it to a file.

• open(videofile), it opens the file.

• writeVideo(videofile,array name), writes all the plot frames into the file and saves it.

• close(videofile), closes the file.

Main Program

clear all
close all
clc

b = 0.05; % Damping factor
g = 9.81; % Acceleration due to gravity
l = 1; % Length of the pendulum
m = 1; % Mass of the pendulum

theta_0 = [0;3];% Initial position in radians and velocity

t_span = linspace(0,20,1000);

[t,y] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

figure(1),clf
plot(t,y(:,1));
hold on
plot(t,y(:,2));
xlabel('Time');
ylabel('Parameter');
legend('Angular Displacement','Angular Velocity');

ct = 1;
for i = 1:length(t_span);
    
    X(i) = cos(y(i,1)-pi/2); % position of the pendulum bob
    Y(i) = sin(y(i,1)-pi/2);
    Xo = cos(y(1,1)-pi/2);
    Yo = l-sin(y(1,1)-pi/2);
    
    figure(2)
    clf
    
     plot([Xo X(i)],[ 0 Y(i)],'-');hold on;
     plot(X(i),Y(i),'marker','o','markersize',15,'markeredgecolor','b','markerfacecolor','b');
    axis([-2*l 2*l -2*l l]);
    set(gca,'xtick',[],'ytick',[])
    axis square
    
     M(ct) = getframe(gcf)
     ct = ct+1;
end

movie(M)
videofile = VideoWriter('Damped Harmonic Pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

PART 2: Creating the function 'ode_function(t,theta,b,g,l,m)'

• In a new script, a function is opened with return value being 'dtheta_dt' as shown

function [dtheta_dt] = ode_function(t,theta,b,g,l,m)

end

• The name of the function is 'ode_function', and the input values are t,theta,b,g,l,m.
• In the body of the function, `theta_2` and `(d theta_2)/dt` are calculated and are assigned as elements of the matrix 'dtheta_dt'.
• This matrix is returned to the main program.

Function Body

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

Output

Variation of Angular Displacement and Angular velocity w.r.t. time

The link to the animated output.

https://youtu.be/_Wx3fhoPal4

Errors

Error: In line 14, declaration that the function 'ode45' is dependent on 't' and 'theta' in incorrectly done.

Correction: The declaration is to be done like this, "@(t,theta)" and not only "(t,theta)".

An error is not explicitly shown in this situation but the output is not as per requirement.

Error: The coordinates of the string joining the bob and the rigid surface is incorrectly given in the line35.

Correction: The coordinates are supposed to be given like this,

plot([Xo Xi],[0 Yi],'-');

Thank You.