Solving second order ODEs to simulate motion of a pendulum using MATLAB

Objective:

To simulate the motion of a pendulum using MATLAB by solving the equations of motion given by second order ODE\'s.

Code:

clear all
close all
clc

b=0.05; % damping coefficient
g=9.81;
l=1; % length of pendulum
m=1 % mass of pendulum

%initial conditions

init_cond= [0;3];

%time duration

t_span = linspace(0,30,1000);

%solve ode

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

%plot results

plot(t,theta(:,1))
hold on
plot(t,theta(:,2))
ylabel(\'Plot\')
xlabel(\'Time\')
legend(\'Displacement\', \'Velocity\')

%% Animation

%Origin

figure(2)
O = [0, 0];
axis([ -1 1 -1.5 0.2])
grid on;

%Loop for animation
ct=1;
for i = 1:length(t)
    %Calculate point of mass
    P=l*[sin(theta(i,1)) -cos(theta(i,1))];
    
    %Define a circle at origin
 
origin_circle=viscircles(O, 0.01);

% define pendulum 

pendulum = line([O(1) P(1)], [O(2) P(2)]);

%represent mass as ball

ball = viscircles(P, 0.05);

pause(0.001) %time interval to update plot

M(ct)=getframe(gcf) %capture frames from the plot window
ct=ct+1; % frame increment

%delete the previous position of the pendulum
if i<length(t)
    delete(pendulum);
    delete(ball);
    delete(origin_circle);
    end

end

%writing and recording movie
movie(M)

videofile = VideoWriter(\'pendulum.avi\',\'Uncompressed AVI\')
open(videofile)
writeVideo(videofile, M)
close(videofile)

Function simplifying the second order equations of motion into a system of first order ode\'s :

%function to calculate velocity and displacement

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

theta1 = theta(1)
theta2 = theta(2)

%breaking down the equation of motion into a set of first order odes

dtheta1_dt = theta2; % change in velocity w.r.t
dtheta2_dt= -(b/m)*theta2 - (g/l)*sin(theta1); % change in displacement w.r.t

dtheta_dt=[dtheta1_dt; dtheta2_dt]; % defining the solution vector
end

Plots and animation:

Figure 1

Animation:

Code explanation:

• The equations of motion of a simple pendulum with damping is given by the following ODE

• The equation consists of a second order differential term which is the change in velocity and a first order differential which is the change in displacement w.r.t.
• In order to write a program to solve this equation, it is necessary to break down the second order ODE into a first order ODE which can be solved by MATLAB.
• A function is defined to hold the set of first order ode to be solved to calculate the displacement and velocity of the pendulum. The change in displacement is defined by dtheta2_dt= -(b/m)*theta2 - (g/l)*sin(theta1) and change in velocity is defined by dtheta1_dt = theta2.
• Both these values are defined in column matrix assigned to a variable called dtheta_dt.
• In order to solve this equation, the necessary inputs are defined at the beginning of the code.
• To solve a differential equation, initial conditions need to be defined. A row vector init_cond is defined holding the values for displacement and angular velocity respectively.
• Since the pendulum needs to be simulated for atleast 20 secs, we define the time period using the linspace command sampling for 500 values between 0-20s to get a smooth representation of the motion.
• The function containing the simplified first ODE\'s is called and defined with an inbuilt matlab syntax for ode45 solver defining the dependant and independant variables along with all the necessary input parameters, time period and initial conditions.
• To visualise the change in displacement and velocity, the results array is plotted as showing in Figure 1.
• To create an animation, we define the origin coordinates where the pendulum is going to be fixed.
• The end point where the mass is going to be suspended is defined P=l*[sin(theta(i,1)) -cos(theta(i,1))] . It is defined within a for loop so that all the values of theta are accounted for the corresponding positions are displayed on the plot.
• To connect both the points together, line command is used to define the initial position of the fixed point and the mass point and the next immediate value for the same.
• viscircles is used in place of a marker to represent a circle around the points.
• M(ct)=getframe(gcf) is defined to record the plot images generated to compile a video later.
• Since we dont need the last recorded values for objects like points and line, an if statement is defined which deletes the points and the lines after each frame is recorded.
• In the final section of the code, a movie file is written and saved as \'pendulum.avi\'.