Project

Solving second order ODEs- SIMPLE PENDULUM

TToom Nithin Chandra10 Mar 2021

MATLAB

AIM

To write a code in MATLAB which solve the simple pendulum motion equation considering damping coefficient

THEORY

simple harmonic motion equation of pendulum is derived from newtons second law it can be given as

`(d^2(theta))/(dt^2)+ (g/l)sin(theta) =0` - here the damping effect is not considered

By considering the damping coefficient equation becomes `(d^2(theta))/(dt^2)+ (b/m)(d(theta))/dt +(g/l)sin(theta) =0`

where, `theta` = angular displacement , b = damping coefficient, g= gravity, l =length of pendulum, m= mass of ball

The above equation is a 2nd order ODE to solve this equation we need to covert this 2nd order equation to a number of 1st order ODE,

let `theta` = `theta1`, and `(d(theta))/dt = (d(theta1))/dt =theta2` - let it be equation 1 , `(d^2(theta))/dt^2 = (d^2(theta1))/dt^2 = (d(theta2))/dt`

then 2nd order equation can be written as `(d(theta2))/(dt) = -(b/m)theta2 -(g/l)sin(theta)` - let it be equation 2

from both equation 1 and 2 we can write vector form `d/dt [(theta1),(theta2)] = [(theta2),(-(b/m)theta2 - (g/l)sin(theta))]`

solution of this equation gives us velocity and displacement

PROCEDURE

first, inputs like damping coefficient(b),gravity(g),length of pendulum(l),mass(m) values are given
and the initial conditions i.e, angular displacement as '0 rad' and velocity as '3 rad/s' are given
time-space is defined here in this we generated motion for 0-20 seconds
now use Matlab inbuilt ODE solver 'ode45' to solve our equation, which takes inputs function which should be solved, time interval and initial conditions and dependent and independent variables t, theta. so, wrote a function where the two equations are placed and named it as 'ode_func' then provided the time interval and initial conditions
then the solution obtained is stored in the results array where it gives two values wrt to time angular displacement and velocity thus w.r.t time angular displacement and velocity graphs are plotted
to create an animation, for loop is created to plot the position of the pendulum for each time in a time interval, to plot the pendulum position fixed coordinates are defined as (0,0) and other end coordinates vary as `lsin(theta),lcos(theta)`
these points are joined by the straight line using the plot command and at the end, a marker symbol is given and its size is defined
all images are collected using getframe command and by using the movie command a video file is created and all it is written with the frames that collected one by one thus each image is displayed and created as a video

CODE

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

MAIN CODE

clear all
close all
clc

%inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%intial conditions
theta_0 = [0 3];

%time points
t_span= linspace(0,20,200);
%ode solving
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,1));
hold on
plot(t,results(:,2));
ylabel('angular velocity and displacement');
xlabel('time');
legend('angular displacement','velocity');
%creating animation
ct=1;
for i = 1:length(t)
    x0 =0;
    y0 =0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    figure(2)
    plot([x0,x1],[y0,y1],'color','r')
    hold on
    
    plot([-1.5 1.5 1.5 -1.5],[0 0 0.1 0.1],'b')
    
    plot(x1,y1,'.','MarkerSize',50,'color','g')
    axis([-1.5 1.5 -1.5 0.5]);
    hold off
    M(ct)=getframe(gcf);
    ct = ct+1;
end
movie(M)
videofile = VideoWriter('pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

OUTPUT