Solving Second order ODE for a Simple Pendulum using MATLAB

Aim:- To study Second-Order ODE of a Simple Pendulum and simulate it in MATLAB

Theory:- Every real-world dynamic system can be expressed mathematically in the form of a differential equation. One such example is the motion of a simple pendulum.

A simple pendulum is an ideal pendulum consisting of a point mass suspended by a weightless inextensible perfectly flexible thread and free to oscillate about its mean position.

Its motion in 2-D space is governed by `2^(nd)` Order Differential Equation, given as;

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

where, `theta` = Angular Displacement from mean position (rad)

b = Damping Coefficient (Ns/m)

m = Mass of point object, i.e. Bob (kg)

g = Acceleration die to Gravity (`m `/` s^2`)

L = Effective length of Pendulum

= length of thread/string + length of hook + radius of spherical bob

In most of the cases mass m, length L, and value of g is constant, but value of Damping Coefficient b varies with the condition of medium, like temperature, moisture, etc.

Fig: Components of Forces Acting on the Pendulum

In order to solve this second-order differential equation, it is broken down into two first-order differential equation, as shown below;

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

Code:- The MATLAB code used for solving these equations for Angular Displacement `theta` and Velocity `omega`.

% Programme to Simulate Motion of Pendulum using 2nd order Differential Equation
clear all
close all
clc

% Pre-Defining Parameters
b = 0.05;
m = 1;
g = 9.81;
l = 1;

% Initial Conditions (position , speed)
theta_0 = [0;3];

% Time Points
t_span = linspace(0,20,100);

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

% Defining Fixed Point
x0 = 1.5;
y0 = 2;

% Column Vector of Angles and Velocity
theta_arr = results(:,1); 
vel_arr = results(:,2);

% Maximum Disp and Velocity
max_disp = max(theta_arr);
max_vel = max(vel_arr);

ct = 1;     % Counter for movie
for (i=1:100)
    
    theta = theta_arr(i);   % Instantaneous Angle
    
    x = x0 + l*sin(theta);  % Corresponding Instantaneous X Co-ordinate
    y = y0 - l*cos(theta);  % Corresponding Instantaneous Y Co-ordinate
    
    figure(ct);
    hold on;
    plot([x0-1,x0+1],[y0,y0],'linewidth',3,'color','g');     % Base of Fixed End
    plot([x0,x],[y0,y],'linewidth',3,'color','b');           % String
    plot(x,y,'o','MarkerFaceColor','red');                   % Marker
    axis([0,3,0,3]);    % Setting Fixed Axis
    grid on;
    
    %Recording Video
    M(ct) = getframe(gcf);
    ct=ct+1;
    %pause(0.0002);
end

close all;  % Closing all the frames

% Displaying Max Displacement and Velocity
fprintf('\nMaximum Angular Displacement = %s deg \n',num2str(max_disp*180/pi));
fprintf('Maximum Angular Velocity     = %s rad/s \n',num2str(max_vel));

% Crteating Plot of Ang Displacement and Ang Velocity against time
figure(ct);
plot(t,results(:,1));    % Position
hold on;
plot(t,results(:,2));    % Velocity
xlabel('Time');
ylabel('Plot');
legend('Angular Position','Angular Velocity');

M(ct) = getframe(gcf);  % Adding last frame of Plot vs time to M
movie(M);
vf = VideoWriter("Challenge_22.avi");
vf.FrameRate = 5;   % Setting Image Rate at 5 Frame Per Second
open(vf);
writeVideo(vf,M);
close(vf);

Explanation of Code:-

• The programme starts by defining the four parameters, namely Damping Coefficient (b), Mass (m), Acceleration due to gravity (g) and Effective length of the Pendulum (l), as specified in the challenge.
• After this, a variable named theta_0 has been declared, which is basically a two-element row-vector. Its first element sets the initial Angular Displacement from mean and the later specifies the Angular Velocity at this position. Here in this code, these parameters are set to 0 rad and 3 rad/s respectively.
• Further linspace() function is used to assign a time-point array for 20 sec with 100 set-points to the variable
• An in-built function ode45() is used to solve the required differential equation. Its parameter is mentioned as below:
  • First parameter @(t,theta), specifies the independent variable time (t) and the dependent variable (theta).
  • This parameter is followed by the actual differential equation (ode_func) mentioned as a separate function. It will be discussed in detail at the later stage.
  • The third parameter is the range of the independent variable within which the equation has to be solved. In this case, it is the t_span.
  • The fourth and last parameter is the set of initial conditions, being specified by theta_0.
• The fixed ends (`x_0,y_0`) of the pendulum is set to (1.5,2).
• The results array is sliced as
  • theta_arr is the column vector of Angular Displacement obtained by array slicing get the first column of results as results(:,1).
  • Similarly, results(:,2) returns the second column of results to variable vel_arr, a column vector containing Angular Velocity corresponding to theta
• max() function returns the maximum displacement and velocity of theta_arr and vel_arr to max_disp and max_vel
• ct variable is initiated with value as ‘1’ to facilitate making M, movie frame array.
• To produce animation’s static frames, for loop is used.
• theta_arr(i) command assigns the instantaneous value of displacement to variable
• This value of theta is then used to calculate corresponding X and Y coordinates (x,y) with respect to the fixed point (x0,y0) using trigonometry relation as
• `x = x_0 + L sin(theta) ; y = x_0 + L cos(theta)`
• Explanation of formula
  • Since pendulum oscillates on both side of fixed x – coordinate (x0) and sin() is an odd function, therefore ‘+’ sign with respect to x0 gives required x-coordinate
  • Since pendulum oscillates only on one side of fixed y – coordinate (y0) and cos() is an even function, therefore ‘+’ sign with respect to y0 gives required y-coordinate.
• The first line after (hold on;) command creates the Horizontal base of the Pendulum.
• The next command (%String) creates the instantaneous image string, by joining the fixed point (x0,y0) and the calculated instantaneous position of pendulum (x,y).
• (%Marker) command makes a red marker at the point (x,y), and axis is set so as to x and y range from 0 to 3.
• These individual frames are added to the movie vector using ct counter
• After the loop ends, (close all; ) command closes all the individual frames. After this, fprintf() function is used to print maximum Displacement and Velocity onto the console.
• Next 7 lines of code plots the Angular Displacement and Angular Velocity along Y-Axis against the time (t) along X-Axis, for the specified duration. This figure is then added to the movie vector M as the last frame.
• vf handle open a Video object, which is used to create an animation/movie from M with an image rate set to 5 fps(frames per second).

The Function Used (ode_func):-

function [dtheta_dt] = ode_func(t,theta,b,m,g,l)
    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

This function along with four basic parameters also takes the t and theta, where theta is the column vector mentioned on LHS of the equation.

• The first two line specifies the elements of theta as `theta_1` and `theta_2` as mentioned in the equation.
• Next two lines represent the two first-order differential equation as discud=ssed in the theory.
• These two equations are then clubbed to return dtheta_dt as a column vector back to the main code.

Input Parameters:- Parameters taken for the Simulation are as following

Damping Factor	b	0.05 Ns/m
Mass of Bob	m	1 kg
Acceleration due to Gravity	g	9.81 m/s2
Effective Length	l	1 m
Initial Conditions
time	t	0 sec
Angular Position	θ	0 rad
Angular Velocity	ω	3 rad/s
Time Range		20 sec

Output Console:- The following console shows the Maximum Angular Displacement in degree and Angular Velocity in rad/s

Following figures shows the three frames of Left Displacement, Mean position and Right Position respectively of the Pendulum.

Video:- Youtube link of an animation of oscillating Pendulum for given specification for first 20 sec (as specified in the challenge)

https://youtu.be/-a-W57suVO4

Error:- No such error was encountered while writing this code. A slight help form documentation helped to control the frame rate of the movie.