MATLAB Scripting to simulate motion of simple pendulum using ODE45

Aim:- To write a program that solves the ODE which represents the equation of motion of a simple pendulum.

Governing equations:-

The ODE which represents motion of simple pendulum

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

Where,

`g=`acceleration due to gravity`(m/s^2)`

`L=`length of the pendulum`(m)`

`m=`mass of the ball`(kg)`

`b=`damping coefficient

Software used:- MATLAB R2019a

Main code:-

%solving 2nd order ODE(simple pendulum)
close all
clear all
clc

%inputs
b=0.05;   %damping coefficient
g=9.81;   %acceleration due to gravity
l=1;      %length of string
m=1;      %mass of bob

%intial conditions
theta_0=[0;3];      %[displacement ; velocity]

%time points
t_span = linspace(0,20,200);

%solving ode by using inbuilt function
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

%Plot of displacement and angular velocity w.r.t time 
figure(1)
hold on
plot(t,results(:,1),'linewidth',1,'color','b');
plot(t,results(:,2),'linewidth',1,'color','r');
legend('Angular displacement','Angular velocity')
xlabel('Time(s)')
ylabel('Angular displacement & velocity')
title('Angular displacement & velocity vs Time')
grid on;

%simulation of pendulum
figure(2)
ct=1;
for i=1:length(results(:,1))
    
    x0=0;
    y0=0;
    
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    
    plot([-1 1],[0 0],'linewidth',5,'color','k');
    hold on
    line([x0 x1],[y0 y1],'linewidth',3,'color','b');
    plot(x1,y1,'o','markers',20,'markerfacecolor','r');
    axis([-1.5 1.5 -1.5 .5]);
    title('Simulation of pendulum')
    grid on
    pause(.1);
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%creating movie
movie(M)
videofile=VideoWriter('pendulum_.avi','Motion JPEG AVI');
videofile.FrameRate=10
open(videofile)
writeVideo(videofile,M)
close(videofile)

• The inputs are given i.e. b,g,l,m .Then the initial displacement and angular velocity are '0' & '3' resp. The time is given as 20 seconds with 200 intervals using linspace command.
• The ODE is solved using inbuilt 'ode45' function. As it solves only 1st order differential eqns. the give ode is converted into 1st order by using another function called 'ode_func' which is explained below.
• This gives two results i.e. displacement , angular velocity. with help of plot command these results are plotted in figure(1).
• By using 'for' loop the frames are generated to create a movie. Using the marker command created a sphere and one end of the pendulum is fixed the other end moves according to the solution of the ODE solved.
• Used 'markersize' command to set the bob size, 'o' for circular shape of bob and 'markerfacecolor' for the color of the bob.
• Using movie command we stitch the frames and make the whole animation through it.

Syntax of 'ode45':-

[t,y] = ode45(odefun,tspan,initial condtions)

where, odefun is the differential eqn to be solved and tspan =[ti tf] integrates the system of differential equations y′=f(t,y) from ti to tf with initial conditions.

Function code:-

% function for converting given ODE into 1st order ODE
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

In this function the given ode is converted into 1st order and the output of angular displacement & velocity is given.

`(d theta_1)/dt=theta_2`

`(d theta_2)/dt=-b/mtheta_2-g/L sintheta_1`

Output:-

We can see that the displacement and velocity decreases with time as it is damped motion.

Animation video link :- https://drive.google.com/open?id=1abCuy5hJXN8w7XZaYKjCD5_TXLCgHvCK