Week 3 - Solving second order ODEs

Aim: To solve 2^nd order governing ordinary differential equation of a motion of a simple pendulum using MATLAB.

Theory:

A pendulum is nothing, but mass or weight suspended at fixed point (pivot point) by using string of foxed length. If there is no action of external force on pendulum, it remains in equilibrium position. But when it is displaced from its equilibrium position with the application of force, it oscillates at its mean position due to gravity. The time taken to complete one oscillation from left swing to right swing or from right swing to left swing is called as period of pendulum. Period of pendulum is given by,

`T = 2 π √ ( L /g )`

Where L is length of the string of pendulum and g is acceleration due to gravity.

Simple pendulum is a type of pendulum in which a mass (bob) is attached at the end of massless string. String is fixed at pivot which is frictionless. Simple pendulum is shown in figure below.

Figure: Simple pendulum

MATLAB code

1. Function code

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

2. Main code

clear all
close all
clc

% Input of variables 
b=0.05;
g=9.81;
l=1;
m=1;

% Initial condition
theta_0 = [0;3];

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

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

% To defined pivot point of pendulum and plot it.
x0=0;
y0=0;
plot (x0, y0,'marker', 'o')

% To extract the first column of a result and assigned it to 'y'
y =  results(:,1);

% To create images for each time value and stored it in array ‘M’.
ct=1;
for i=1:length(t_span)
    Y1=y(i)
    x1 = l*sin(Y1);
    y1 = -l*cos(Y1);
    plot([x0 x1],[y0 y1],'linewidth',3,'marker','o')
    hold on
    plot(x1, y1, 'marker','o', 'color','g','markersize',30, 'markerfacecolor', 'g')
    hold off 
    title({'Simple harmonic motion'; 'of a pendulum'}, 'color','r')
    axis([-1 1 -2 0])
    pause(0.01)
    M(ct)= getframe(gcf)
    ct=ct+1;
end

%  To create animation of S.H.M. of simple pendulum
movie(M)
videofile = VideoWriter('pendulum.avi', 'Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Errors face while coding

1. pivot point (x0,y0) is not defined

2. Varable 'Y1' is not defined to access indivisual values of y of each time step

3. Improper use of 'hold off' command to erase plot of privious time step

4. Missing of curly bracket in 'title' command

Results:

Youtube link of animation of SHM of simple pendulum: https://youtu.be/LTGZ2dYK5js

Google drive link of animation of SHM of simple pendulum: https://drive.google.com/file/d/11aG_z19XWCimdvw0Y3pYpWic4B8Ut7iX/view?usp=sharing

Conclusion

Governing differential equation of simple pendulum is solved using MATLAB program. Animation is working correctly. Hence it can be concluded that solution obtained by program is correct.