Simple Pendulum Motion Using Matlab

OBJECTIVE:

Program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.

What is ODE?

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.

Fig 1 : Pendulum Motion

EQUATION OF MOTION OF SIMPLE PENDULUM:

`d^2theta/dt^2 + b/m* d theta/dt + g/L* sintheta = 0`

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

Considering:

`theta = theta1` equation 1

then

`d theta / dt = d theta1 / dt = theta2` equation 1

`d^2 theta/dt^2 = d^2 theta_1/ dt^2 = d/dt ( d theta_1/dt) = d/dt theta_2` equation 2

putting equation 1 & equation 2 in SIMPLE PENDULUM EQUATION

`d^2 theta_1 / dt^2 + b/m d theta_1/dt + g/L sin theta = 0`

`d theta_2 / dt + b/m theta_2 + g/L sin theta_1 = 0`

`d theta_2 / dt = - b/m theta_2 - g/L sin theta_1`

Creating Matrix Equation:

`d/dt[ [ theta_1 ],[ theta_2]] = [[theta_2] ,[-b/m theta_2 - g/L sin theta]]`

The above is the solution vector and gives the values for theta.

defining the function [dtheta_dt]

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2
dtheta2_dt= -(b/m)*theta2 - (g/l)*sin(theta1)
dtheta_dt = [dtheta1_dt;dtheta2_dt]
end

dtheta_dt is the variable which contains dtheta1_dt & dtheta2_dt

where dtheta1_dt = `d theta_1/dt = theta_2` & dtheta2_dt = `- b/m theta_2 - g/L sin theta_1`

Program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion has been given below:

clear all
close all
clc

b=0.05;
g=9.81;
l=1;
m=0.1;

% initial condition (displacement & velocity)
theta_0=[0;3]

% time points
t_span = linspace(0,10,500)

% solving the ode
[t,x] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span,theta_0)

figure(1)
plot(t,x(:,1),'linewidth',3);
hold on
plot(t,x(:,2),'linewidth',3);
xlabel('time')
ylabel('Amplitude')
hold on ; grid on

%  using for loop
ct=1
ang_vel=x(:,1)
for i=1:length(ang_vel);
    ANG_VEL=ang_vel(i)
    x0=0;
    y0=0;
    x1=-l*sin(ANG_VEL)
    y1=-l*cos(ANG_VEL)
    
    % plot
    figure(2)
    plot([x0 x1],[y0,y1],'linewidth',3)
    hold on
    plot(x1,y1,'marker','o','markerfacecolor','r')
    axis([-1.5 1.5 -1.5 1.5])
    pause(0.03)
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end
movie(M)
videofile=VideoWriter('ode_main.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)



    

• The linspace command will generate 500 values of time between 0 to 10 seconds.
  
  
• The ode45 is an inbuilt function to solve the differential equation. It gives the result at various values of input between time 0 second to 10 seconds. It uses runge kutta 4th order.
  
  
• The inputs of damping coefficient, gravity, length of the pendulum & mass of the ball has been given initially.
• The for loop is a repetition control structure that allows us to efficiently write a loop that needs to execute a specific number of times.

RESULTS:

Fig 2- Displacement vs Velocity

Animation:

Errors made:

Fig 3 - did not use the valid value

Fig 4 - used the correct value