Solving second order ODE's of Simple Pendulum using MATLAB Program

Objective: To write a program to solve the provided ODE of a simple pendulum and to create an animation File to display the motion of the pendulum

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with the periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained

The second-order ode for simple pendulum is represented as,

`(d ^2 θ)/ (d t^ 2) + b/ m . (d θ)/ (d t )+ g/ l sin θ = 0` d2θdt2+bm.dθdt+glsinθ=0">

• b = damping coefficient 0.05,
• m = mass of ball 1 kg,
• g = acceleration due to gravity 9.81 m/s^2,
• l = length of the string 1 metre.

When these equations are solved, they are converted into two first-order ODEs which are then converted to matrices. These matrices can be solved in MATLAB to obtain the velocity and displacement values with respect to time.

In matrix

`d/ (d t) (( θ_ 1 ),(θ_ 2 )) = [ (θ_ 2),( − b /m θ_ 2 − g /l sin θ _1 ) ]`

by solving these we will get the θ1,θ1 and θ2,θ2 values of the simple pendulum.

MATLAB code for a function to solve the second-order ODEs

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 will create a column vector once the values of dtetha1_dt and dtetha2_dt are defined and assign it to the dtetha_dt.

Input Parameters :

L=1 metre

m=1 kg

b=0.05

g=9.81 m/s^2

Main MATLAB CODE:

clear all ;
close all ; 
clc ; 
b=0.05; 
g=9.8 ;
l=1; 
m =0.1;

theta_0=[0;5];

t_span= linspace(0,20,500);

[t,results]= ode45(@(t,theta) ode_func(t,theta,b,m,g,l),t_span, theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
title('Simple Pendulum')
xlabel('Time in sec')
ylabel('Angular Velocity(rad/s) / Displacement(rad)')
legend('Displacement','Angular velocity')

disp_theta = results(:,1);

ct=1;%counter initiation
for i = 1:length(disp_theta) 
    angle = disp_theta(i); 
    figure(2)
    title('Pendulum Motion')
    
    x0 = 0; %x & y coordinates of the fixed point of pendulum
    y0 = 0;
    
    x1 = l*cos((3*pi/2)+angle); %x & y coordinates of bob
    y1 = l*sin((3*pi/2)+angle);
  
    %plotting 
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g');
    hold on
    plot(x1,y1,'o','markers',15 , 'MarkerFaceColor','g'); %to create the bob 
    hold on
    plot([-3 3],[0 0],'linewidth',3,'color','b'); %fixed point's base
    hold off
    axis([-3 3 -2 2])
    grid on;
    pause(0.03);
    M(ct) = getframe(gcf);
    ct = ct+1;
end
movie(M) 
videofile = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Explanation

First clear all, close all, clc command is used to clear all the data from the workspace and command window and close all this figure which is open its always good to start with these three commands before starting a code

The percent sign is most commonly used to indicate non-executable text within the body of a program. those non-executable texts are for better understanding the code

The input parameter for the equation is defined

To solve the ODE, an inbuilt solver 'ODE45' is used which calls the user-defined function and solves it with respect to the initial boundary conditions and time span

The Important Code to Solve A ODE function

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

[t,results]= Output parameters such as time velocity displacement

ode45 = Function keyword to solve second order ODE

(@(t,theta)= ODE is function of t and theta

ode_func(t,theta,b,m,g,l)= Input parameter for the ode function

t_span = Time interval

theta_0 = initial condition of theta

Using the Plot Command, graphs are plotted

The concept of counter 'ct' is introduced to segregate the frame

For Loop is used to iterate the different positions resulting due to varying values of Angular Displacement

The position of the two ends of the pendulum is defined as coordinates (x0,y0) and (x1,y1) and are plotted

To plot the pendulum bob, the plot command is used along with a marker to define a bob at all points of (x1,y1). The size and color of the marker is also controlled with respective commands

pause command is used to mention the pause time for each iteration results to show

M(ct) = getframe(gcf) is used to capture each frame from the plot

movie(M) command is used to combine those each frame from M(ct) and make a video file

This movie is now written to a video file after creating it and opening it and once written is closed.

The video file is created in the same location as the program.

Graph :

Animation File of the Simple Pendulum