Week 3 - Solving second order ODEs

Simple Pendulum Simulation using Matlab

Objective: To write a program that solves the second order ODE representing the equation of motion of a simple pendulum.

Equation of Motion representing Pendulum:

Here,

g, gravity =9.81 m/s2,

L, length of the pendulum = 1 m,

m, mass of the ball = 1 kg,

B, damping coefficient 0.05.

Calculation:

To solve the ODE of order 2, it has to be converted into a matrix containing two ODE of order 1.

… (1)

Let

......(2)

Differentiating on both sides,

......(3)

Differentiating eqn. (3)

…..(4)

Sub (2),(3),(4) in (1)

....(5)

In equation 5, term 1 denotes angular displacement and term 2 denotes angular velocity.

Main Program:

%Given variables
L=1;
m=1;
b=0.05;
g=9.81;

%time interval
t= linspace(0,20,100);

%initial conditions at t=0
theta0=[0;3];

% function
[t,output] = ode45(@(t,theta)odefun_second_order_pendulum(theta,L,m,b,g),t,theta0);

% plotting 
figure(1)
plot(t,output(:,1),'color','r','linewidth',3);
hold on;
plot(t,output(:,2),'color','b','linewidth',3);
xlabel('angular displacement, angular velocity');
ylabel('time');
legend('angular displacement','angular velocity');

%movie capture
figure(2)
ct=1;
for i=1:length(t)
  Theta1(i)=output(i,1);
  x0=0;
  y0=0;
  x1=L*sin(Theta1(i));
  y1=-L*cos(Theta1(i));
  plot([x0 x1],[y0 y1],'linewidth',2);
  axis([-2 2  -1.5  0.5]);
  hold on;
  plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','r');
  hold off;
  M(ct)=getframe(gcf);
  ct=ct+1;
end
movie(M)
videofile=VideoWriter('kinematics.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Function Program:

function [dthetabydt]=odefun_second_order_pendulum(theta,L,m,b,g)

dtheta1bydt=theta (2);
dtheta2bydt= (-b*theta(2))/m-(g*sin(theta(1))/L);
dthetabydt= [dtheta1bydt;dtheta2bydt];

end

Result:

Above plot denotes the relation between the X-axis Angular displacement and Angular velocity against Y-axis time.

The above image shows the motion of the pendulum. Complete animation can be viewed at

https://youtu.be/WdbAEEFzHAk