simple pendulum simulation using MATLAB coding

To write the program to solve the following ODE.This ODE represents the equation of a motion of a damped simple pendulum.

d2θdt2+Cmdθdt+glsinθ=0">d2θdt2+Cmdθdt+glsinθ=0d2θdt2+Cmdθdt+glsinθ=0

θ=θ1,anddθ1dt=θ2∴dθ2dt=d2θ1dt2">θ=θ1,anddθ1dt=θ2∴dθ2dt=d2θ1dt2θ=θ1,anddθ1dt=θ2∴dθ2dt=d2θ1dt2

dθ2dt=-Cmθ2-glsinθ1">dθ2dt=−Cmθ2−glsinθ1dθ2dt=-Cmθ2-glsinθ1

Codes for Calling a function: calling a function means that the function which we have defined with some parameters are going to use in other program, Here this function is going to solve the Second order equation or we can say that its going to integrate the values.

function Dtheta = pendulum_motion_ode(theta,t)
%inputs
b = 0.5;
l = 1;
m = 1;
g = 9.81;
theta1 = theta(1);
theta2 = theta(2);
Dtheta(1) = theta(2);
Dtheta(2) = -(b/m)*theta(2) - (g/l)*sin(theta(1));
Dtheta = [Dtheta(1);Dtheta(2)];
end

Codes for Simulation of Simple Pendulum:

1) we are solving An Ode system and there are two first order ODE's one for velocity and other for displacement. so when we solve the first ode we need to provide initial conditions for the displacement and velocities.

2) Ode 45 is an inbuit solver in matlab and octave which is going to solve the ODE system

3) our Ode depends upon the time, theta and other parameters i.e.

b = 0.5 , damping coefficient.

g = 9.81 m/s ,accceleration due to gravity.

l = 0.9 in m , length.

m = 2 in kg , mass of the bob.

clear all
close all
clc



%initial conditions
theta_0 = [0;3];

%duration
t_span = linspace(0,20,200);
l = 1;

%solving the ode
[results] = lsode('pendulum_motion_ode',theta_0,t_span)



%plotting
plot(t_span,results(:,1))
hold on
plot(t_span,results(:,2))
xlabel('time')
ylabel('amplitude')
hold off
figure(2)
plot(results(:,1),results(:,2))
xlabel('angular displacement')
ylabel('angular velocity')

ct=1;
for i = 1:length(results(:,1))
  position = results(i,1);
  
  x0 = 0;
  y0 = 0;
  
  x1 = l*cos(position-(pi/2));
  y1 = -l*sin(position+(pi/2));
  
  figure(3)
  plot([x0,x1],[y0,y1],'linewidth',3)
  hold on
  plot(x1,y1,'marker','.','markersize',40) 
  plot([-1,1],[0,0],'linewidth',4,'k')
  hold off
  axis([-3,3,-3,3])
  grid on
  pause(0.1)
 

  ct = ct + 1;
  end

Stepwise explaination of the above codes:

1.) In the first code, a function for solving the given ODE is written with variables there and t(time).

2.) Constants of the equation which are coefficient of damping(b), acceleration due to gravity (g), Length of string (L) and mass of the Bob(m) are defined.

3.) Then the second order differential equation is converted into first order differential equation as shown in the code.

4.) In the second part, the initial conditions for angular displacement and angular velocity are input as given in the question. Also the duration for which the simulation has to take place is also specified

(20 seconds).

5.) Using the "lsode" function, the first code which was written to store the variables the differential equations and will be called and will be solved using the initial conditions for the given duration.

The solution will be stored in results vector.

6.) Next step is to plot angular displacement and angular velocity versus time.

7.) Using the for loop, the result of each loop is generated and saved using the sprintf and saveas function. The Bob at the end of the string is created using the marker function.

8.) Finally, all the saved images are stitched together create the required animation using imagemagick.

Plot to show the damping effect:- From the plot we can clearly see that as the damping coefficient increases the damping effect also increases that means there will be less oscilaations.