Pendulum Motion under damping

The above given equation is an ODE of order 2. We can solve it using Octave.

The solution is a 1x2 matrix with y(1) as `theta` and y(2) as `omega`.

Given below is the code in Octave:

clear all
close all
clc
% Solving second order ODE
% Let theta = y and angular velocity = v

t=linspace(0,20,1001);
y0=0;
v0=3;

y=lsode(\"pendulum_fun\",[y0 v0],t);

figure(1)
subplot(3,1,1)
plot(t,y(:,1))
xlabel(\'Time\',\'FontWeight\',\'bold\',\'FontSize\',14)
ylabel(\'Theta\',\'FontWeight\',\'bold\',\'FontSize\',14)

subplot(3,1,2)
plot(t,y(:,2),\'r\')
xlabel(\'Time\',\'FontWeight\',\'bold\',\'FontSize\',14)
ylabel(\"Angular\\nVelocity\",\'FontWeight\',\'bold\',\'FontSize\',14)

velocity=y(:,2)\';
acceleration=diff(velocity)./diff(t);
z=linspace(0,19.98,1000);

subplot(3,1,3)
plot(z,acceleration,\'g\')
xlabel(\'Time\',\'FontWeight\',\'bold\',\'FontSize\',14)
ylabel(\"Angular \\nAcceleration\",\'FontWeight\',\'bold\',\'FontSize\',14)

l=1;
i=1;
while i<1002
theta=y(i,1);
A=l*sin(theta);
B=l*cos(theta);
dx=linspace(0,B,10);
dy=tan(theta)*dx;
f=figure(2);
plot(dy,dx,\'linewidth\',3);
axis([-1 1 0 1.56]);
axis manual;
hold on 
plot(A,B,\'0\',\'MarkerSize\',20)
pause(.01)
i=i+1;
cla;
end

Function:

function [ydot]=pendulum_fun(y,t);

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

ydot=zeros(1,2);
ydot(1)=y(2);
ydot(2)=-(b/m)*y(2)-(g/l)*sin(y(1));
end 

Plots:

Link to animation video:

https://youtu.be/t3FjywBZZ6o