Solving 2nd order ODE of a simple pendulum and animating its transient behaviour

This program solves the second order ODE which represents the transient behaviour of simple pendulum by converting it into two first order ODE\'s. It also plots the angular displacement and velocity vs time graph and creates an animation of pendulum with damped oscillations.

First the values of length, mass, acceleration due to gravity and damping cofficient are given and then initial values of position and velocity are provided.

After that time span for the plot is provided.

Then a function is created to define the two first order ODE\'s which is then used to solve the equations using ode45 command. This gives us \'results\' which have the position (1st column) and velocity (2nd column) at different time values.

These values are used to plot angular displacement and velocity vs time curve and then an animation is created to simulate the motion of damped pendulum.

Youtube link for the animation is https://youtu.be/BQP-lpSZ7To

% solving 2nd order ODE for a simple pendulum

clear all
close all
clc

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

% initial conditions
theta0 = [0 3];   % initial displacement = 0, initial velocty = 3

% time span
t_span = linspace(0,20,500);    % plot is made for 20 sec


% ODE function
function [dtheta_dt] = ode_f(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];
  endfunction

% solving ODE

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

% plotting displacement and velocity vs time

figure(1)
plot(t,results(:,1))    % plots angular displacement vs time
hold on
plot(t,results(:,2))    % plots angular velocity vs time
xlabel(\'time\')
ylabel(\'Angular displacement and velocity\')
legend(\'Angular displacement\',\'Angular velocity\')

% making animation of pendulum

figure(2)

ct=1
for i = 1:length(t_span)
  THETA = results(i,1);
  
  x_pos = l*sin(THETA);
  y_pos = l*cos(THETA);
  
  plot([-1 1],[0.15 0.15],\'linewidth\',50,\'color\',\'b\')
  hold on
  plot([0 x_pos], [0 -y_pos],\'linewidth\',4,\'color\',\'g\')
  plot(x_pos,-y_pos,\'o\',\'markersize\',40,\'markerfacecolor\',\'magenta\',\'markeredgecolor\',\'black\')
  plot(0,0,\'o\',\'markersize\',10,\'markeredgecolor\',\'black\')
  hold off
  
  axis([-1.25 1.25 -1.25 0.5])
  pause(0.001)
  
  text = sprintf(\'zzpendulum%d.png\',ct);
  saveas(gca,text);
  ct=ct+1;

end