Damped Pendulum Simulation using MATLAB   (MATLAB, Skill lync, Easy)

Aim:

To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping.

Solution:

Used parameters

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s^2.

Simulated the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Code:

1) Main code:

clear;
clc;
clf;

t_span = [0 20];
y0 = [0; 3];
[t, y] = ode45('pendulum', t_span, y0);

ang_disp=y(:,1);
ang_vel=y(:,2);
%calculation of co-ordinates
x_coordinate = sin(y(:,1));
y_coordinate =1-(cos(y(:,1)));

%giving a figure handle
figh = figure;
for i=1:length(t)
    clf;
    subplot(2,1,1);
    
    %Values of co-ordinate for particular timestep
    t_k = t(i);
    x_k = x_coordinate(i);
    y_k = y_coordinate(i);
    
    %values of co-ordinates for fixed string of pendulum
    line_x=[0,x_k];
    line_y=[1,y_k];
    
    %plotting pendulum
    plot(x_k,y_k,'ro','LineWidth',3,'MarkerSize',10);
    hold on
    plot(x_coordinate,y_coordinate,'-k','LineWidth',0.5);
    hold on
    
    %plotting string
    plot(line_x,line_y,'-b','LineWidth',2);
    grid on;
    
    %plotting value reference lines
    ref_x1=[0,x_k];
    ref_y1=[y_k,y_k];
    ref_x2=[x_k,x_k];
    ref_y2=[0,y_k];
    hold on
    plot(ref_x1,ref_y1,'k','LineWidth',0.5);
    hold on
    plot(ref_x2,ref_y2,'k','LineWidth',0.5);
    
    %decoration of figure
    xlabel('X co-ordinate');
    ylabel('Y co-ordinate');
    title(['t = ', num2str(t_k),' seconds']);
    axis([-1 1 0 1]);
    axis equal;
    
   
    %Simultaneos Displacement plot
    subplot(2,1,2);
    plot(t,ang_disp,'-k','LineWidth',2);
    hold on
    plot(t(i),ang_disp(i,1),'ob','MarkerSize',10);
    grid on;
    xlabel('Time');
    ylabel('Displacement');
    title('Displacement  vs  Time');
    
    %plotting value reference lines
    refe_x1=[0,t(i)];
    refe_y1=[ang_disp(i,1),ang_disp(i,1)];
    refe_x2=[t(i),t(i)];
    refe_y2=[0,ang_disp(i,1)];
    hold on
    plot(refe_x1,refe_y1,'k','LineWidth',0.5);
    hold on
    plot(refe_x2,refe_y2,'k','LineWidth',0.5);
    
    
    %capturing frames
    movieVector(i) = getframe(figh);
    
end

%making video
myvid = VideoWriter('Pendulum','MPEG-4');
myvid.FrameRate = 10;
open(myvid);
writeVideo(myvid, movieVector);
close(myvid);

2) Code for ODE:

function angle_dot = pendulum(t,y)
    b=0.5;
    m=1;
    l=1;
    g=9.81;
    angle_dot=zeros(2,1);
    angle_dot(1)=y(2);
    angle_dot(2)=-((b/m)*y(2))-((g/l)*sin(y(1)));
end

Corresponding Video link: