Week 3 - Solving second order ODEs

Aim: To simulate pendulum motion of a sphere of mass 1kg, with pendulum length 1m, damping coefficient 0.05 when it’s displaced from equilibrium position at an angular velocity of 3 rad/sec.

Objective: To solve the second order differential equation using MATLAB of a simple oscillating pendulum, also to see how the damping coefficient variation effects the oscillation to slow down.

Procedure:

We have second order ODE to solve for simple oscillating pendulum as follow

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

Where = angle swept by hanging ball from equilibrium position, t = time at which swept angle is,

b = damping coefficient, m = mass of ball, g= 9.81 m/s^2, l = length of rod.

We have l = 1m, m = 1kg and b = 0.05. Also the initial condition is that at = 0 and t = 0, we have angular velocity `dθ/dt` = 3 rad/s

In line 11 of ode_main code, we define initial conditions as theta_0 array which is 2x1 matrix with , in row 1 and angular velocity `dθ/dt` = 3 rad/s in row 2

In line 14 of ode_main code, we define time stamps from 0 to 20s divided into equal 500 intervals to calculate the angular velocity at that point in time as per the governing equation.

In line 1 of ode_function code, we create a 2x1 functional array dtheta_dt to calculate the angular velocity at a given time t and swept angle following the governing ODE. This function arguments are t, theta, b,g,l and m to be provided while calling this function.

In line 3 of ode_function code, variable theta_1 is angle swept by pendulum from normal.

In line 4&5 of ode_function code, variable theta_2 is angular velocity at that particular angle of theta_1 and time t.

In line 6 of ode_function code, differential of theta_2 is equated as per the governing equation as

dtheta_2_dt = -(b/m)*theta_2 -(g/l)*sin(theta_1)

In line 17 of ode_main code, we use ODE solver in Matlab using ode45 command. Here the arguments we call upon are dependent on time t and angle theta, the ode_function which follows the code needed to bring two equations to be solved, the time span and initial boundary condition(theta_0) of angular velocity of 3 rad/s at 0 degrees. This line gives us the results array at time t, which are angular displacement and velocity at that time t.

In line 20&21 of ode_main code, we describe the fixed end of pendulum co-ordinates

In line 24& 25 of ode_main code, we calculate the co-ordinates of end point of pendulum using theta and length ‘l’ of pendulum from geometry concept.

In line 28, we define counter ct as 1 for indexing the plot frames in the movie M

In line 29, we create for loop with index ‘i’ from 1 to 500 for capturing the plots at each time with our time span.

In line 30, we plot the pendulum line using its co-ordinates, round marker shape and increased size.

In line 31, axis command is used to set the figure to be displayed in middle of plot.

In line 32, set command is used to invert the y-axis so that fixed end of pendulum is on top and other end of pendulum indicates the ball position at different timestamps.

In line 33, all the frames of movie M are collected into an array using counter ct with an increment of 1 after each for loop.

In line 39, videofile is created VIdeoWriter command to define video file name and type.

Line 40&42 are helpful in opening and closing videofile. While in 41, writeVideo command is used to create animation on videfile using movie frames M.

Conclusion:

The pendulum animation is created by solving the governing equations for its simple harmonic motion when it’s displaced from neutral position at 3rad/s. This animation clearly shows that the pendulum doesn’t reach its neutral position within 5s. Also when the damping coefficient is increased, the pendulum tends to reach the neutral position quickly.

ode_main.m code

clear all
close all
clc

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

%initial conditions
theta_0 = [0;3];

%time_stamps
t_span = linspace(0,20,500);

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

%Co-ordinates for Fixed end of pendulum
x0 = 0;
y0 = 0;

%Co-ordinates for ball
x1 = l*sin(results(:,1));
y1 = l*cos(results(:,1));

%plotting co-ordinates
ct =1;
for i = 1:500
 plot([x0 x1(i)], [y0 y1(i)],'-o','MarkerSize',10,'linewidth',2)
 axis([-1.5 1.5 -2 1.5])
 set(gca,'YDir','reverse');
 M(ct)=getframe(gcf);
 ct = ct+1;
end

%creating movie animation using frame from each plot
movie(M)
videofile = VideoWriter('Pendulum Motion','Uncompressed AVI');
open(videofile)
writeVideo(videofile, M)
close(videofile)

ode_function.m code :-

function [dtheta_dt] = ode_function(t,theta,b,g,l,m)

theta_1 = theta(1);
theta_2 = theta(2);
dtheta_1_dt = theta_2;
dtheta_2_dt = -(b/m)*theta_2 -(g/l)*sin(theta_1);
dtheta_dt = [dtheta_1_dt;dtheta_2_dt];

end