Solving a Second Order ODE for a pendulum motion and Generating its Animation

Aim:

Write a program in MATLAB that can solve a second order ode and plot the animation of a pendulum.

Objective:

Solve a second order ode.

Plot position and velocity of pendulum vs time.

Generate an animation of the pendulum’s motion.

The workflow of the program is as follows:

• Create a function which can easily input an ode in form of matrix.
• Define parametric definitions like mass, length, damping coefficient and acc. Due to gravity.
• Input the initial conditions
• Input the time span over which the ode is integrated.
• Now solve the ode using “ode45” command and within this command use the previously written function to input the ode in form of matrix.
• Now, plot time vs the position and velocity of the pendulum.
• Now, plot the pendulum for each value of position in a “for” loop.
• Collect the each frame in an array.
• Now generate the .avi file of pendulum’s motion by combining all the frames together.

Program:

Here is a program for a function that inputs the ode’s in form of matrix:

%% Description 
%
% ode_define is a function which defines the ode in form of matrix
%
%% Author
%
% Umang Tyagi
% Date of Creation : 27 July 2020
%
%% Program
function [dtheta_dt] = ode_define(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

Now, here is the actual program that solves the ode and generates the required plot and the animation:

%% Description
%
% This is a program to solve an ODE for a pendulum's motion and generate an animation of its motion. 
%
%% Output
%
% 1) Plot of bob's position and its angular velocity w.r.t. time.
% 2) Animation file of motion of pendulum for 20 seconds.
%
%% Author
%
% Umang Tyagi
% Date of Creation : 27 July 2020
%
%% Program
clear all
close all
clc

%% 1) Parameter Definition
l = 1;                % length of pendulum (m)
m = 1;                % mass of bob (kg)
b = 0.05;             % damping coefficient
g = 9.81;             % acceleration due to gravity (m/s^2)

%% 1.1) Initial Conditions
theta_0 = [0;3];                       % first column contains position and second colum contains velocity

%% 1.2) Time Span
t_span = linspace(0,60,500);           % it is span over which the ode is integrated

%% 2) Solving ODE
% in order to solve the ode, ode45 function is used. 
% ode function takes the input of ide in matrix form, the t_span and initial conditions.
% ode45 returns two solutions in case of 2nd order ode.
[t, result] = ode45(@(t,theta) ode_define(t,theta,b,g,l,m), t_span, theta_0);

% output 't' is an array of values of time.
% output 'result' is a 2 column array. 
% the first column conatins the values of theta w.r.t time.
% the second column contains the valuesd of velocity w.r.t time.

%% 3) Plotting Position and Velocity vs Time
figure(1);                                    % generate an empty figure 

%% 3.1) Plot time(t) vs position(theta)
plot(t,result(:,1),'linewidth',1);            % plot t vs theta
hold on                                       % hold on command is used to hold the current plot on the same figure.
%% 3.2) Plot time(t) vs velocity
plot(t,result(:,2),'linewidth',1 );           % plot t vs velocity 
xlabel('time')                                % label the x_axis
ylabel('position and Velocity');              % label the y-axis
legend('Theta','Velocity')                    % insert the legends on the plot
hold off                                      % turn the hold command off 
pause(2);

%% 4) Plotting Pendulum 
figure(2);                                    % generate another empty figure
%% 4.1) Looping 
for i = 1:length(result)                      % looping all values of theta(first column of result)
    plot([-0.5 0.5],[0 0],'color','k','linewidth',0.8);      % plot the base line
    hold on                                   % hold all further plots on same figure
   
    x = l*sin(result(:,1));                   % define x coordinate of the bob
    y = -l*cos(result(:,1));                  % define y coordinate of the bob
    plot([0 x(i)],[0 y(i)],'linewidth',0.8,'color','k');     % plot the pendulum's string


    plot(x(i),y(i),'.','MarkerSize',50,'color','r');         % plot the bob
    axis([-2 2 -2 2]);                        % define the axis
   M(i) = getframe(gcf);                      % store each plot as a frame in array M
   hold off                                   % hold off the plot, to get a new plot on next iteration
end

%% 5) Generating Animation
movie(M);                                                           % combine all frames   
pendulum_animation = VideoWriter('pendulum_animation.avi');         % generate video in avi format
open(pendulum_animation);                                           % open video
writeVideo(pendulum_animation,M);                                   % write video using M array
close(pendulum_animation);                                          % close video

Output:

Here is the plot of time vs the position and velocity of theta:

Here is the link to the animation file of the pendulum’s motion:

https://drive.google.com/file/d/15XnU_QWHUra8bMm-qDQ-ePAIYz-iulHH/view?usp=sharing

Conclusion:

• The required plots of time vs position and velocity of pendulum is generated correctly.
• The animation of pendulum’s motion is also generated correctly.