Solve Second Order ODE of Simple Pendulum and Animation of its motion

AIM: Write a program in MATLAB that solves the second order ODE of pendulum and animate the motion of pendulum.

Introduction

In Engineering, ODE (ordinary differential equation) is used to describe the transient behavior of a system. A Pendulum is simplest example of the ODE.

The way the pendulum moves depends on the Newtons second law. When this law is written down foe pendulum, we get a second order ODE that describes the position of the "ball" w.r.t time which as follow:

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/L sin⁡θ=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.

MATLAB have ODE45 command to solve the ODE. But for the second order ODE first we have to convert into first order ODE.

So, convert our second order ODE in first order

First, we plot the transient behavior of pendulum than we create the animation of the motion of the simple pendulum.

Step-wise explanation of code for plot the graph of transient behavior of pendulum:

1. First, we have to assign the required input which is damping coefficient(b), Length of the pendulum(L), mass of the ball(m) and gravity(g).
2. We have to assign value of theta in matrix form because it contains 2 values which is angular displacement and angular velocity.
3. We have to give the time span and convert into this span into equal interval with the use of linspace command.
4. We have to use ode45 command to solve first order ODE with the help of sub function ode_func with give the angular displacement and angular velocity at any instant of time.
5. In the subfunction we solve our 2-differential equation with given input.
6. Now the plot the graph of angular displacement of ball with respect to time and angular velocity with respect to time.
7. We plot these two graphs in figure so we use legend command to easy understand of angular displacement and angular velocity curve.
8. We use title command give the title of the graph.

For simulation of motion of pendulum,we have to add some steps

1. First for the show motion of pendulum in different figure. we assign the figure.
2. We use the for loop for the show change in displacement with respect to time.
3. We assign the coordinate of string at fixed end (x₀, y₀) = (0,0)
4. We need to calculate the coordinate of free end of string (x₁, y₁) according to the value of angles which depends in time.
5. We have to plot line for pivot the string in it.
6. Then we have to plot another line for showing equilibrium position of ball.
7. Now plot the string of the graph.
8. We have to show the ball in animation so we use the marker command to plot the sphere.
9. We have to give the proper the value of X and Y axis.
10. For save every position on the ball we use pause command.
11. We have to make the new variable M to store the frame of the current figure.
12. We have to give increment to counter command.
13. Then after the for loop, we have to create the movie of variable which is use to store the frame. So, the final animation is created.

The errors faced while programming & steps taken to overcome those errors:

1. I make mistake in subfunction where I put the theta1 in place of theta2 in the line no 6. So, I can’t get the right graph of transient behavior.

Updates: I correct the mistake and get right graph

2. First, I am not using hold off command and also the write hold on command at right place.so I get the ugly animation.

Updates: I correct the position of hold on command and put the hold off command at proper position.

Screenshots showing errors thrown in command window

1.

2.

Program for the plot transient behavior and animation of simple pendulum

clear all 
close all 
clc

%inputs
b=0.05; %damping co-efficient
L=1;    %Length of the Pendulum in m
m=1;    %mass of the ball in Kg
g=9.81;  %gravity in m/s^2

%initial vaule
theta_0 =[0;3]; %initial anglar displacement is zero and initial angular velocity is 3 rad/s

%time points
t_span = linspace(0,20,500); %total time for motion is 20 sec

%solve ODE
[t results]= ode45(@(t,theta) ode_func(t,theta,b,g,L,m),t_span,theta_0); %ode45 to solve ODE with using ode_func

hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('Time')
ylabel('Anglar displacement & Angular velocity')
legend('Angular displacement','Angular velocity')
title('Transient behaviour of Simple Pendulum')
hold off

figure(2)
ct = 1
for i = 1:length(t_span)
    
    x0 = 0;
    y0 = 0;
    
    x1 = 1*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    
   
    plot([-1 1],[0 0],'linewidth',10,'color','k')
    hold on
    plot([0 0],[0 -1.2],'linestyle','--','color','k')
    
    plot([x0 x1],[y0 y1],'linewidth',5,'color','r')
    plot(x1,y1,'marker','o','markerFacecolor','b','markersize',30)
    title('Simple Pendulum motion')
    
    axis([-1.5 1.5 -1.5 0.2])
    pause(0.0001)
   
    
    M(ct)= getframe(gcf);
    hold off
    ct=ct+1;
    
end

movie(M)
Videofile= VideoWriter('Pendulam.avi','uncompressed AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)

Program of sub function

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

Graph of transient behavior of simple pendulum

Animation of the simple pendulum