Motion of Pendulum using ODE45

https://drive.google.com/open?id=1KbbOINw6LcJmgGgy9Iyvzqvq_2on0y3i

Objective

The aim of the project is to write a program in MATLAB to simulate the motion of pendulum with given knowns:

• Gravity,g = 9.81 m⁄s^2
• Length of pendulum,L = 1 m
• Mass of ball, m = 1 kg
• Damping Coefficient, b = 0.05 kg/s
• Total Duration, t = 20 s

Intial Condition at t=0 s,

• Initial Displacement, `θ_o`= 0 rad
• Initial Angular Velocity, θ ̇= 3 `rad⁄s^2`

ODE: A classic example – Pendulum

In engineering, Ordinary differential equations ‘ODE’ is used to describe the transient behaviour of any system. Most simple and classic example seen in textbooks is the pendulum whose ODE has been represented below:

`(d^2 θ)/dt^2 +b/m (dθ)/dt+g/L sin⁡θ=0`

where,

g = gravity in m/s²

L = length of the pendulum in m

m = mass of the ball in kg

b=damping coefficient

The motion of pendulum depends on the Newtons second law. The above shown ODE describes the position of pendulum w.r.t time.

The ODE stated above is the second order differential equation. It can be solved more easily if the order of the equation can be reduced.

Reduction of Order

Any explicit differential equation of order, n can be written as a system of n first-order differential equations by defining a new family of unknown functions, i.e., `y_i = y^(i-1)`

Based on above statement, lets reduce the order of pendulum’s governing ODE.

`(d^2 y)/dt^2 +b/m dy/dt+g/L sin⁡y=0`

Let,

`y=y_1``(d^2 y_1)/dt^2 +b/m (dy_1)/dt+g/L sin⁡y_1 =0`

Rewrite equations as first-order ODE by making substitution: `(dy_1)/dt=y_2`

`(dy_2)/dt+b/m y_2+g/L sin⁡y_1 =0`

`(dy_2)/dt=-b/m y_2-g/L sin⁡y_1`

Resulting system of two first-order ODE,

`y_1^\'=y_2 and y_2^\'=-b/m y_2-g/L sin⁡y_1`

MATLAB Solver: ode45

In order to solve the above reduced ODE, we require the time span and initial conditions.

The reduced ODE can be solved in MATLAB using inbuilt ODE solver, ode45 whose syntax is shown below:

[t,y] = ode45(odefun,tspan,y0)

The above syntax integrates the system of differential equations `y^\'=f(t,y)` from `t_(Initial)` to `t_(Fi nal)`, i.e. tspan with initial conditions, i.e. y0. Thus, giving output as t and y. Here, y represents the angular displacement, θ which is needed to transformed into coordinates in order to be plotted over a plane and visualize the motion

MATLAB Code Structure

MATLAB Code

Main Program:

clear all
close all
clc
b=0.05;     % damping coefficient
m=1;        % mass of ball
g=9.81;     % gravity
l=1;        % length of pendulum
 
t_span=linspace(0,20,100);      % time duration
y_init=[0;3];                   % initial conditions i.e. angular displacement and velocity
 
dydt = @(t,y)pendulum_fun(y,b,m,l,g);
[t, y]=ode45(dydt,t_span,y_init);
 
int_plot(t,y);      % plot of displacement and angular velocity w.r.t. time
 
[y1, x1] = coordinate_fun(y,l);     % coordinate calculator
 
animation_plot(x1,y1,t_span)        % animation generator

Defined Functions:

1. ‘pendulum_fun’ Function: Defined the reduced ODE

function ydot = pendulum_fun(y,b,m,l,g)
    ydot = [y(2);-(b/m)*y(2)-(g/l)*sin(y(1))];
end

2. ‘int_plot’ Function: Plots angular displacement and velocity w.r.t. time

function int_plot(t,y)
    figure (1)
    grid on
    
    subplot(1,2,1)
    plot(t,y(:,1));
    title(\'Displacement vs Time\');
    xlabel(\'Time span in seconds\');
    ylabel(\'Angular displacement in radian\');
    
    subplot(1,2,2)
    plot(t,y(:,2));
    title(\'Angular Velocity vs Time\');
    xlabel(\'Time span in seconds\');
    ylabel(\'Angular velocity in radian per second\');
    
    saveas(figure (1),\'SubPlot\',\'jpg\')
end

3. ‘coordinate_fun’ Function: Conversion from polar to cartesian coordinate system

function [y1, x1]=coordinate_fun(y,l)
    y1=-l*cos(y(:,1));      % y-coordinate for ball\'s position
    x1=l*sin(y(:,1));       % x-coordinate for ball\'s position
end

4. ‘animation_plot’ Function: Generates video for motion of pendulum

function animation_plot(x1,y1,t_span)
    for i = 1:100
        figure (2)
        plot([0 x1(i)],[0 y1(i)],\'linewidth\',2,\'color\',\'b\');
        hold on
        plot(x1(i),y1(i),\'marker\',\'o\',\'linewidth\',12,\'markerfacecolor\',\'r\');
        plot([-1.2 1.2],[0 0],\'linewidth\',2,\'color\',\'g\');
        axis([-1.2 1.2 -1.2 1.2]);
        title([\'Time: \', num2str(t_span(i))]);
        M(i) = getframe(gcf);
        close(2)
    end
 
    movie(M)
    videofile = VideoWriter(\'pendulum.avi\',\'Uncompressed AVI\');
    open(videofile)
    writeVideo(videofile,M)
    close(videofile)
end

Output Results

From the above plots, we can observe the effect of damping coefficient and gravity on the angular displacement and angular velocity w.r.t time. The initial conditions could also be verified from the above shown plots, i.e. angular displacement as 0 rad and angular velocity as 3 rad/s² at 0 seconds.

The below shown video shows the motion of pendulum w.r.t. to time.