Solving ODE using matlab and simulating the motion of the pendulum

The project involves the simulation of a simple pendulum ODE using MATLAB as the programming language.

AIM: To simulate the transient behavior of a simple pendulum and to create an animation of its motion in MATLAB.

CODE:

clear all 
close all
clc

b = 0.05; %damping co-effecient
g = 9.81; %acceleration due to gravity
l = 1; %length of the pendulum 
m = 0.1; % mass of the pendulum

%initial conditions
theta_0 = [0;3];

%time plots
t_span = linspace(0,20,500); %integrating the system(time array)

%ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
legend('velocity','angular displacement')
xlabel('time')
ylabel('plot')


ct = 1;
for i = 1:length(results)
   THETA = results(i);
   
       x0 = 0;
       y0 = 0;
       x1 = -l* sin(THETA);
       y1 = -l* cos(THETA);
       
       figure(2)
       %plotting
       plot([x0 x1],[y0 y1],'color', 'r', 'linewidth', 4)
       hold on
       plot([-1 1],[0 0], 'color', 'c', 'linewidth', 3)
       plot(x1,y1,'O','markerfacecolor','blue','markersize', 10 );
       hold off
       axis([-1.5 1.5 -1.5 1.5]);
       pause(0.03);
       
       P(ct) = getframe(gcf);
       ct = ct+1;
end
movie(P)
vidfile = VideoWriter('ode_pendulum.avi', 'Uncompressed AVI');
open(vidfile)
writeVideo(vidfile,P)
close(vidfile)

A CONCISE DESCRIPTION OF THE CODE:

The code begins with the standard operation of clear all, close and clc to eliminate and refresh the editor. Followed by defining the inputs for the code such as the damping coefficient 'b', accn due to gravity'g' in m/s2, length of the pendulum'l' in m and the mass of the pendulum 'm' in kg.

The initial conditions of angular displacement '0' and angular velocity of '3' rad/s.

The time array is defined between 0 to 20 sec with intervals of 500.

The ODE is solved using ODE45 and the function defined with the name ode_func along with the input parameters, time array and the initial condition.

The function of the code is as follows:

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

The function defined with regard to the input parameters and breaking down the second-order differential equation to a series of 2 first-order differential equations. As it can be noticed from the program in the form of dtheta_1 and dtheta_2 and having it arranged in the final column matrix which provides the output.

The function is run and the plot of varying velocity and displacement with respect to time can be observed below

The variable ct considers the count of each frame executed for the intervals.

The for loop with 'i' is initiated for every corresponding value of the results with the susbsequent definition of the result of 'i' in THETA.

The initial conditions and the corresponding length of the pendulum in the coordinate system is evaluated.

ERROR encountered:

by considering

x1 = l* sin(THETA);
y1 = l* cos(THETA); in the initial program I happen to come across the error wherein the oscillations occurred in the manner as the image displayed.

Further the coordinate reference was changed to

x1 = -l* sin(THETA);
y1 = -l* cos(THETA); to obtain the satisfactory result.

The plotting was carried out for the oscillations of the mass and the length using the plot command and providing appropriate size and color. Along with the definition of the axis that it is plotted in.

The P(ct) command captures each frame for the loop of values and compiles the frames to generate the movie/animation.

The "vidfile" defines the name and the type of video created by the Videowriter.

The oscillation of the pendulum can be viewed by the video seen below.

Kudos to the mentors at skill-lync.