Solving ODE of simple pendulum motion using Octave / MATLAB

Our objective is to write a program that solves the following ODE.

The following ODE represents the equation of motion of a simple pendulum with damping.

`(d^2theta)/(dt^2)+b/m (d theta)/dt + g/L sintheta = 0`

where, 

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

Given data:-> 

L = 1 meter,

m = 1 kg,

b = 0.05.

g = 9.81 m/s2.

At initial condition, time = 0, angular displacement = 0,  angular velocity = 3 rad/s^2.

Time span = 0 to 20 seconds.

MATLAB/ Octave has standard solver for solving ordinary differential equations (ODEs) the function is "ODE45". This function calculate solutions of ODEs by using Numerical method i.e. Runge-Kutta method for a time domains at variable time steps. Function syntax is as follows.

Output = ode45("Function", "time", "Initial condition") 

For coding we are using ODE45 solver. In following code an anonymous function f(t,theta) is used.

here, consider angular displacement = theta = theta(1)

for angular velocity i.e.`(dθ)/dt = theta_(2)= (d theta_1)/(dt)` 

For angular acceleration, `(d^2θ)/(dt^2)=  -(b/m(theta_(2)))-(g/L sintheta_(1))`

And initial condition is given by `theta_0`

Octave code:-

clear all
close all
clc

l=1; %Length of pendulum string in meter
m=1; %Mass of pendulum in kg
b=0.05; %Damping coefficient
g=9.81; %Gravitational acceleration in (N/m^2)

t=linspace(0,20,100); %time span 0 to 20 seconds with 100 time steps.
theta0=[0 3]; %initial conditions [angular displacement Angular velocity]

%Ordinary differential equation function 
[t, theta]= ode45(@(t,theta) [theta(2) ; -(b*theta(2)/m)-(g*sin(theta(1))/l)], t, theta0);
%Here theta is angular displacement which is shown by theta(1)
%Rate of change of theta is angular velocity which is shown by theta(2)
%Rate of change of theta(2) is angular acceleration

% Creation of graph resulting magnitude of angular velocity & displacement w.r.t time
figure(1);
plot(t,theta(:,1),'color','b'); %plotting time v/s angular displacement
hold on;
plot(t,theta(:,2),'color','r'); %plotting time v/s angular velocity
legend('Angular displacement','Angular velocity');
xlabel('Time (Seconds)');
ylabel('Magnitude');
title('Pendulum plot for angular displacement & angular velocity')

for i= 1:length(theta(:,1))
  figure(2)
  clf;
  %origin of pendulum (x0,y0)(Constrain point)
  x0=0;
  y0=0;
  % Other end of pendulum (x1,y1) (Centre of bob)
  x1=l*sin(theta(i,1));
  y1=l*cos(theta(i,1));
  plot([-1 1], [0 0]);
  hold on;
  line([x0 x1], [-y0 -y1]);
  plot(x1,-y1,"0",'markersize',40,'markerfacecolor','g'); % Plotting of bob
  grid on;
  axis([-1.5 1.5 -1.5 0.5]); %axis dimensions
  title('Pendulum animation');
  printf('Creating GIF- Progrss... %d/%dn', i, length(t));
  img = print('-RGBImage');
  imwrite(img, 'animation.gif', 'DelayTime', 0.1,'Compression','bzip','WriteMode','append');
endfor

After running above we get output plot of time v/s magnitude of angular velocity & angular displacement

And we get animation of pendulum for time scale 0 to 20 seconds.