Solving simple pendulum

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the mtion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec2 at time t=0.

Objective: To understand and develop the code for the pendulum oscillation. Animate the pendulum oscillation. Plot the displacement and angular velocity graph.

Results:

The Animation video link is attached along with the code. The output graph image is attached.

Successfully animated the pendulum with a damping coefficient of 0.05 and angular velocity=3rad/s2 for the time range of 0-20s.

Video link:

https://www.youtube.com/watch?v=e1IZqraeKUc

Function code:

%damped pendulum function program(d^2theta/dt^2)+(b/m)*(dtheta/dt)+(g/l)*(sin(theta))
function [dtheta]=ode_func(t,dtheta,L,m,b,g)
theta1=dtheta(1);
theta2=dtheta(2);
theta_1=theta2;
theta_2=-(b/m)*theta2-(g/L)*sin(theta1);
dtheta=[theta_1;theta_2];
end

clear all
close all
clc
%time dependent damped pendulum simulation (d^2theta/dt^2)+(b/m)*(dtheta/dt)+(g/l)*(sin(theta))
L=1; %Length of pendulum in m
m=1; %Mass of ball in kg
b=0.05; %Damping coefficient
g=9.81; %Gravity in m/s2
theta_0=[0;3]; %intial condition %displacement=0 & angular velocity=3 rad/s
n=500; %step size
tf=20; %final time in s
dt=tf/n; %time of each step
t_span=linspace(0,tf,n); %time array
%solve ODE
[t,r]=ode45(@(t,theta) ode_func(t,theta,L,m,b,g),t_span,theta_0);
%plotting graphs
subplot(6,5,[3,4,5,8,9,10,13,14,15])
plot(t,r(:,1),'b');
hold on;
plot(t,r(:,2),'r');
xlabel('time');
ylabel('plot');
legend('Displacement','Velocity');
hold off;
subplot(6,5,[24,25,29,30]);
plot(r(:,1),r(:,2));
xlabel('position-radians');
ylabel('angular velocity');
%pendulum animation
%loop for animation
for i=1:length(r(:,1))
    x0=0;
    y0=0;
    x1=L.*sin(r(i,1)); 
    y1=-(L.*cos(r(i,1)));
    subplot(6,5,[1,2,6,7,11,12,16,17,21,22,26,27]);
    plot([-0.2,0.2],[0,0],'linewidth',2,'color','k');
    hold on;
    plot([x0 x1],[y0 y1],'linewidth',4,'color','b');
    plot(x1,y1,'-or','markersize',15,'MarkerFaceColor','r');
    %to show the title with time stamp
    titlestring = ['Damped pendulum motion at t =',num2str(i*dt),'seconds'];
    title(titlestring ,'fontsize',14);
    axis([-1 1 -2 1]);
    grid on;
    M(i)= getframe(gcf);
    hold off;
end
%save the video
movie(M)
videofile=VideoWriter('SPendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)