Week 3 - Solving second order ODEs

A simple pendulum consists of mass which is attached by the string of length L from a pivot point P. when displaced by initial force the mass gets moved back and forth performing the periodic motion. By applying Newton's second law motion we can obtain the equation of motion for a simple pendulum which is given as.

PROGRAM:

% A Program for Solving second order ODEs
clear all
close all
clc

% Inputs
% Demping coefficient
b = 0.05;
% Acceleratin due to gravity
g = 9.81;
% Length of Pendulum in meters
l = 1;
% Mass
m = 1;

% Time duration for simulation
t_span = linspace(0, 20, 100);
% Providing Initial Condition as Angular displacement & Velocity
y_initial = [0;3];
% Calculation using ODE45

a = @(t,y)ode_second_order_pendulum(y,b,m,1,g);
[t,y] = ode45(a,t_span,y_initial);


% Plot of Angular displacement & Velocity with respect to Time function
figure(1);
grid on
subplot(1,2,1)
plot(t,y(:,1));

title('Displacement vs Time');
xlabel('Time in seconds');
ylabel('Angular displacement in rad');
subplot(1,2,2)
plot(t,y(:,2));

title('Angular velocity Vs Time');
xlabel('Time in seconds');
ylabel('Angular velocity in rad/sec');
grid off
saveas(figure(1), 'subplot', 'jpg')

% Polar Co-ordinate To Cartesian Co-ordinate
y1 = -1*cos(y(:,1)); 
x1 = 1*sin(y(:,1));

% Video for Pendulum motion
for i = 1:100
    figure(2)
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','r');
    hold on
    plot(x1(i),y(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0 0],'linewidth',3,'color','b');
    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('singlependulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

% The function is given as
function ydot = ode_second_order_pendulum(y,b,m,l,g)
ydot = [y(2);-(b/m)*y(2)-(g/1)*sin(y(1))];
end

EXPLANATION OF CODE:

• For clearing command window & workspace, closing all plots.
• The code starts with the inputs which are provided in the problems. then we make a time domain for which the whole code is going to run.
• It was done by the 'linspace' command from 0 to 20 seconds with 100 values in between them.
• After that, we provide the initial condition for the angular displacement and velocity in column vector form, because we need to calculate both values from the equation of motion for plotting purposes.
• The code comprises of function which solves the second-order differential equation which is taken inside the main code.
• The main code has another function that calculates the array of data for time and displacement.
• After calculation plotting is done in which subplot command is used, in the first plot first column is plotted and in the second subplot, the second column is plotted concerning time.
• Then we are changing the polar co-ordinate to cartesian co-ordinate because we have used the 'sin' command instead of 'sind' command in function.
• Then we have created a for loop for which we are creating the frame motion for the simple pendulum and saving it in the movie.
• saving the movie in .avi format.

• I have found a few errors and solutions were found during execution.
• Use of the semi-colon(;) in coding results in no output in the command window.

As we can see from the plots the displacement is getting reduced as we move further into time because of the damping factor. From the plot, we can clearly observe displacement is starting from 0 as an initial condition. In the angular velocity vs time plot, it has started from 3 rad/sec and then further gets decreased.

• https://youtu.be/vFO2RQyrJZ8