Solving ODE for simple pendulum using MATLAB

Aim: To evaluate the motion of simple pendulum by solving ODE using MATLAB

Theory:

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. 

Consider a pendulum with ball of mass `m` attached to the end string of length `L` with other end being fixed.

where,

`theta` = Angular Displacement

`g` = Gravitational Acceleration

`t` = time

`b` = Damping Coefficient

For solving ODE using MATLAB, higher order differential equations needs to converted into a set of single order differential equations.

A second order differential equation should be convereted into two single order differential equations. This is done by considering,

`theta = theta_1`

`(d theta_1)/dt = theta_2`

and this leads to following differential equations,

`(d theta_1)/dt = theta_2`

`(d theta_2)/dt = -(b/m)*theta_2 - (g/l)*sintheta_1`

Now these equations can be solved using ode45 function of MATLAB.

The parameters choosen for study are as follows:

`L=1 m`

`m=1 kg`

`b=0.05`

`g=9.81 m/s^2`

with initial conditions,

`theta = 0`

`(d theta)/dt = 3`

MATLAB Main Program:

% Program to solve 2nd order ODE for pendulum
close all
clear all
clc

% inputs from challenge
l=1;    % length of pendulum
m=1;    % mass of ball
b=0.05; % Damping Coefficient
g=9.81; % gravitational acceleration

% Initial Condition
theta_0 = [0;3];

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

% ode solution
[t results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

% plot results
figure(1)
plot(t, results(:,1))
hold on
plot(t, results(:,2))
xlabel('time')
ylabel('plot')

% creating animation
ct =1;		% counter variable to keep track of frame number
for i = 1 : length(results(:,1))
    
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    figure(2)
    plot([-1 1], [0 0], 'linewidth', 10, 'color', 'g')		% plot fixed support
    axis([-2 2 -2 2])
    hold on
    line([x0 x1], [y0 y1], 'linewidth', 5, 'color', 'b')	% plot string position
    plot(x1, y1, 'o', 'markersize', 25, 'markerfacecolor', 'r')	% plot ball position
    hold off
    M(ct) = getframe(gcf);	% stores current frame
    ct = ct+1;
end

% Create movie
movie(M)
videofile = VideoWriter('Pendulum_simulation.avi', 'Uncompressed AVI');	% create videofile
open(videofile);	% open videofile
writeVideo(videofile, M)	%  write frames to videofile
close(videofile)		% close video file

ode_func :

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

ode_func converts second order differential equation into two single order differential equations.

Code is explained used comments at each section.

Response Curve:

In plot Blue curve represnts displacemets and Red curve represents Velocity.

Animation: