Simulating the motion of a pendulum using 2nd order ODE in MATLAB

This project explores the simulation of the motion of a pendulum for the equation of damped vibration of a pendulum.

Below is the equation (taken from "Challenges" page):

An attempt is made to solve the above equation, and then call it using ode45 to get an array of results vs time (position vs time). These positions are coordinates that are plotted over 20 seconds and then animated.

ODE Function:

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


Function here takes time, angle, damping factor, g, length and mass of pendulum as arguments and returns an array of 2 variables that represent angular acceleration and velocity. Here, solving multiple equations concept is used to arrive at variables that are returned and used as arguments for ode45 function in the main code below.


Main Code:

b = 0.05;
g = 9.81;
l = 1;
m = 1;

%initial condition
theta0 = [0; 3];

%time duration
tspan = linspace(0, 20, 500);

[t, results] = ode45(@(t,theta) odefunc(t, theta, b, g, l, m), tspan, theta0);

%subplot(4,4, [9, 10, 13, 13])
plot(tspan, results(:,1), 'color', 'b');
hold on
plot(tspan, results(:,2), 'color', 'r');
hold off

%subplot(4,4, [11, 12, 15, 16])
plot(results(:,1), results(:,2));
xlabel('Position- radians')
ylabel('Angular Velocity')

count = 1;
for i = 1:length(results(:,1))
  x0 = 0;
  y0 = 0;
  x1 = l*sin(results(i,1));
  y1 = -l*cos(results(i,1));
  %subplot(4,4, [1, 2, 3, 4, 5, 6, 7, 8])
  plot([-l l], [0 0], 'linewidth', 3, 'color', 'k')
  hold on
  line([x0 x1], [y0 y1], 'linewidth', 4, 'color', 'b')
  hold on
  plot(x1, y1, '-o', 'markers', 20, 'markerfacecolor', 'r') %plotting bob
  grid on
  axis([-1.5 1.5 -1.5 1.5])
  pause (0.000003)
  hold off
  file_text = sprintf('pendode%d.png',count);
  count = count + 1;

Here, there are 500 instants of time from 0 to 20 seconds for which the motion is simulated.

ode45 uses the function, tspan and initial condition as arguments. Function "odefunc" returns the dtheta1 and dtheta2 values that are solved for to get "results", which is the position of the pendulum with respect to time.

Figure 1:

Plot of results vs time to show the effect of damping.

Plot vs TIme

Figure 2:

Plot of angular velocity vs position.

Ang V. vs Pos

Figure 3:

Here, for loop is used to continually generate plots of the pendulum at different values of results for the length of "results" (which we got from ode45). Line command generates the pendulum string (assumed weightless) at different points through the entire length of results and "markerfacecolor" generates the ball at the end.

Being octave, "sprintf" and "saveas" are used to get the frames, which are compiled using ImageMagick:

magick pendode%d.png[1-500] PendulumOdeAnim.mp4



Pendulum Animation, solved using ODE45 on Octave


Error 1: Attempting lsode in octave

lsode attampt failed, lsode removed

Apparently lsode has been removed in version 5.1+ and ode45 has been enabled. Rest of the code was done using ode45.

Error 2: Syntax error

Syntax error

Second equal sign. Deleted and changed.

Error 3: Subplot experimentation

Subplot error

Was attempting subplot. Was unsure of the syntax until it got clarified. Eventually followed below code to get it right. But eventially subplot was dropped (hence commented out in main code) and normal plots used.

subplot(4,4, [9, 10, 13, 13])

Projects by Arjun Bhat

Curve Fitting using MATLAB
Arjun Bhat · 2019-08-19 20:20:53

Project aims at finding the best fit curve for a given data and finding the error between the original and said curve. The data being used is that of Specific heat (Cp) vs Temeperature. Plot will be made to show variation of Cp with Temp first.  Code: % Data cp Read more

Project involves plotting an Otto Cycle graph in MATLAB The goal is to plot all strokes for a petrol operating by the otto cycle from intake, compression, power to exhaust, including both adiabatic curves. Description: The Otto Cycle is the idealised cycle that descr Read more

Code and Description: The project involves simulating a 2-joint robotic arm. % Lengths l1 = 1; l2 = 0.5; % Angles ang1 = linspace (0, 90, 20); ang2 = linspace (0, 90, 20); count = 1; % loop for i = 1:length(ang1) ANG1 = ang1(i); for j = 1:length(ang2 Read more


The End