Solving Ordinary Differential Equations for pendulum using ODE45

Wrote a Matlab script to solve the non-linear pendulum equation and simulate it. The equation being solved is:

where b is the damping coefficient (0.05), m is the mass of the pendulum bob (1 kg), L is the length of the pendulum (1 m), theta is the angular displacement and g is the gravitational acceleration (9.82 m/s2). The equation is solved numerically using the ODE45 function with an initial angular velocity of 3 rad/s and angular displacement of 0 radians for 20 sec.

To view the animation please visit:

https://www.youtube.com/watch?v=CByQidAC1V4&feature=youtu.be

clear all
close all
clc

g = 9.81; %Gravitational accelaration (ms^-2)
L = 1;   %Length of the pendulum (m)
m = 1;   %Mass of the pendulum (kg)
b = 0.05; %Damping coefficient 
runtime = [0 20];   %Simulating b/w 0 and 20s
inicondi =[0; 3];    %Initial condition @ 0 angular displacement and 
...angular veliocity of 3 rad/sec

theta = @(t, y) [y(2); (-(g/L)*sin(y(1))-(b/m)*y(2))];  %Fucntion handle for Pendulum Equation
[t, sol] = ode45(theta, runtime, inicondi); 
x_cord = L*sin(sol(:,1));       %Cartesian Coordinates from angular-displacememnt
y_cord = -L*cos(sol(:,1));
x_vel = L.*sol(:,2).*cos(sol(:,1));
y_vel = L.*sol(:,2).*sin(sol(:,1));

figure(1)

%===Plot wall===%
%Horizontal lines
plot([-1.5 1.5],[0 0],\'k\',\'linewidth\',2);
hold on
plot([-1.5 1.5],[0.125 0.125],\'k\',\'linewidth\',1);
plot([-1.5 1.5],[0.25 0.25],\'k\',\'linewidth\',1);
plot([-1.5 1.5],[0.375 0.375],\'k\',\'linewidth\',1);
plot([-1.5 1.5],[0.5 0.5],\'k\',\'linewidth\',2);
%Vertical lines
plot([-1.2 -1.2],[0 0.125],\'k\',\'linewidth\',1);
for lev = 0:3
    for brk = 0:10
                tot = 0;
       if(lev==0 || lev==2)
        if(tot > -1.5 && tot <= 1.5)
            plot([-1.5 -1.5]+(0.3*brk),[0 0.125]+(0.125*lev),\'k\',\'linewidth\',1);
            tot =+ 0.3;            
        end
       else
         if(tot > -1.5 && tot < 1.5)             
            plot([-1.35 -1.35]+(0.3*brk),[0 0.125]+(0.125*lev),\'k\',\'linewidth\',1);
            tot =+ 0.3;           
         end          
       end
       
    end
end

        

%===Plot Pendulum Swing Axis===%
plot (0, 0, \'sk\',\'markersize\',10);
%===Plot Pendulum Vertical Reference Axis===%
plot([0 0], [0.75 -1.25],\'-.k\');
xlim([-1.5 1.5])
ylim([-1.5 1.5])

for i = 1:length(t)
    
figure(1);

%===Plot 1m rod===%
p1 = plot([0 x_cord(i)],[0 y_cord(i)],\'b\',\'LineWidth\',3);


%axis off;


%===Plot Pendulum bob===%
p2 = plot(x_cord(i),y_cord(i),\'o\',\'markersize\',16,\'markerfacecolor\',\'c\',\'MarkerEdgeColor\',\'k\');

%===Velocity vector===%
p3 = quiver(x_cord(i), y_cord(i), 0.1*x_vel(i), 0.1*y_vel(i), \'--r\',\'maxheadsize\',1);
%legend(p3,\'Linear Velocity Vector\');

%===Pendulum Bob Tracer===%
h1 = plot(x_cord(i), y_cord(i),\'.k\',\'markersize\', 0.2);

drawnow;
pause(0.002);

    

%%=====Angular Displacement & Angular Velocity Plots=====%%

%===Ploting using subplot===%

% figure(2);
% subplot(2,1,1)
% plot(t(i),sol(i,1),\'ok\');
% hold on
% xlabel(\'Time\')
% ylabel(\'Angular Displacement\')
% subplot(2,1,2);
% plot(t(i),sol(i,2),\'or\');
% xlabel(\'Time\')
% ylabel(\'Angular Velocity\')
% drawnow;
% hold on


%===Ploting using plot===%

%figure(2);
%plot(t(i), sol(i,2),\'ok\',t(i), sol(i,1),\'or\')
% drawnow;
% hold on


%===Ploting using plotyy===%

%figure(2);
% [hAx, hline1, hline2] = plotyy(t(i), sol(i,1),t(i), sol(i,2));
% set(hAx(1),\'xlim\',[round(min(t(:))),round(max(t(:)))]);
% set(hAx(2),\'xlim\',[round(min(t(:))),round(max(t(:)))]);
% set(hAx(1),\'ylim\',[round(min(sol(:,1))),round(max(sol(:,1)))]);
% set(hAx(2),\'ylim\',[round(min(sol(:,2))),round(max(sol(:,2)))]);
% hline1.Marker = \'*\';
% hline1.MarkerEdgeColor = \'b\';
% hline2.Marker = \'o\';
% hline2.MarkerEdgeColor = \'r\';
% drawnow;
% hold on


F(i) = getframe(gcf);
if i<=(length(t)-1)
delete([p1 p2 p3]);
    end
end

%%===Save the animation in avi format===%%

% video = VideoWriter(\'Pendulum.avi\',\'Uncompressed AVI\');
% video.FrameRate = 15;
% open(video);
% writeVideo(video,F);
% close(video);