The Swinging Pendulum

Introduction

---

Mathematical Model

---

Method

---

Numerical Results and Visualization

Case I Without Damping Coefficient:

At a step size (h) = 0.05s and a number of points = 201 point:

• As can be seen from the previous two figures, Euler's method is numerically unstable (ill-conditioned) as its step size is does not account for the sudden change at the slope. Thus, it will yield inaccurate results. On the other hand, Runge-Kutta of order four provided us with stable and precise approximations, even with a large step size.

At a step size (h) = 0.005s and a number of points = 2001:

• As can be seen from the above two figures, even with a ten times smaller step size Euler's method provide us with a precise values but not stable since the amplitude of the displacement and velocity is increasing. While Runge-Kutta yields the same precise results from the previous step size (h= 0.05), proving that Runge-Kutta is superior to Euler's method.

Displacement Relative Error:

• In this section, we have compared the results of the displacement from both methods for each logarithmic spaced step size from 1 to 1E-5. We notice that the results would become more similar when the step size is smaller. However, the precision would become better at the expense of the computational time, as the last tested step size (h = 1E-5) took approximately 30 minutes to compute the results for both methods. Therefore, we picked a step size around 1E-2 for the optimal value.

Finite Difference:

At step size = 0.005s and T = 2.0061 s

• As can be seen in the above Figure, Finite difference's method can also be used to approximate the values of the displacement of the Pendulum with a good precision similar to Runge-Kutta's method which implies stability. Apart from this, Finite difference approach is limited to the fact that it most be attached to a couple of boundary conditions which are most probable not available.

Cord Tension:

• As seen in the above figure which illustrates the variation in the cord tension from zero to nearly 150 and it goes back and forth. This happens because we exclude the air resistance in this case which will cause continuous, unvarying oscillations. The tension will be minimum when theta equals pi/2 and omega equals zero which will make the tension to be zero. Whereas, it will reach the maximum when theta equals zero and the omega is maximum. Also, we can see the black line which represent the weight of the object.

---

Case II: With Damping Coefficient:

At a step size (h) = 0.005s, a number of points = 2001, and Cd = 0.5:

• Given that the pendulum will experience air friction that will result in damping, we included a damping coefficient in the differential equation to simulate the results of the induced friction. We observe above that the angular displacement, velocity, and acceleration will decrease with time till it reaches zero.

Cord Tension:

• Cord tension in the second case where we considered the air resistance which will add extra term (damping coefficient). This will cause a change in omega and that will affect theta. Therefore, the cord tension will be changed in every oscillation until it becomes as the weight of the object because the theta will be zero and the omega will also be zero. So, the maximum value will be varying, and the cord tension will be zero just in the first moment when theta equals pi/2 and omega equals zero.

---

Pendulum Animations

Without damping coefficient:

With damping coefficient:

When length of cord is doubled:

When gravity is doubled:

When mass of object is doubled:

With initial angular velocity equals -4.5:

Discussion

• We observe that the mass has no affect on the period of the oscillation. While the length of the cord is proportional to the period of the pendulum and the gravity is inversely proportional to the time it takes to complete one cycle.

• If the pendulum experiences drag, it will slow down until it reaches equilibrium.

• When the pendulum has enough kinetic energy, it can complete a full revolution with the condition of no drag.

Conservation of Energy

No drag:

With drag:

Discussion

• When the pendulum experiences no drag, we observe that the potential energy will be maximized when it reaches the highest point, while the kinetic energy will be maximized when it reaches the highest velocity (when it is completely vertical). The energy will be conserved between the two energy forms.

• If the pendulum experiences drag, the loss will increase at the cost of the kinetic energy and will stagnate if only the potential energy is available.

---

Summary

Finally, let us take a quick look to what we did in this project. As we have seen in this project that the swinging pendulum problem requires a good effort to be solved due to the complexity in the analytical solution. So, we choose the three methods to solve it numerically (RK4, Euler, and finite difference methods). We construct our mathematical model to illustrate the problem and makes it easier for us to solve and we list the needed equations, assumptions, and methods. By that, we have taken a clear illustration of the problem. Then, we started our work to get the numerical results. We have put two cases. The first case is without damping coefficient and the second is with. In our work, we have seen the inaccuracy in the Euler method in the first case but RK4 method have worked very well. Also, we calculated the cord tension to give a more understanding for the problem. In the second case, we have added a damping coefficient that would simulate factors like air friction. After we got the numerical results, we illustrated our work in animation to make it clearer and we have shown many different examples. We concluded our work by showing the conservation of energy in both cases.

References

- Pendulum Numerical Animation Tutorial

- Steven C. Chapra, Raymond P. Canale - Numerical Methods for Engineers, 6th Edition -

McGraw-Hill Higher Education (2009) (pp. 824-837)

- Steven C. Chapra Dr. - Applied Numerical Methods with MATLAB for Engineers and Scientists-McGraw-Hill Education (2017)

---

Authors:

Haroun Alattas

Mohammed Alabdullah

Mohammad Alnegheimish

Appendix :

Code:

Differential equations function:

function dy = dydtsys(t,y,g,L,Cd)
dy = [y(2); (-g/L)*sin(y(1))-Cd*y(2)]; % dy(1) = omega    dy(2) = alpha
end

Main script:

clc
clearvars
close all

% assumptions 
m = 5; %5Kg
L = 1; %1 meter
g = 9.81; % m/s^2             
% g2 on mars or any other planet
t = [0 10]; % seconds
h = 0.005; % step size
iter = t/h; % number of iteration
Cd = 0.5; % damping coefficient 
% Initial conditions   
theta0 = pi/2;
omega0 = 0;
y0 = [theta0 omega0];
                  

[t, y1] = eulersys(@dydtsys,t,y0,h,g,L,Cd);
[t, y2] = rk4sys(@dydtsys,t,y0,h,g,L,Cd);



alpha_RK4 = (-g/L).*sin(y2(:,1))-Cd.*y2(:,2);
alpha_Euler = (-g/L).*sin(y1(:,1))-Cd.*y1(:,2);

% RK4
figure(1)
%Displacement
subplot(3,1,1);
plot(t,y2(:,1),'red');
legend('Displacement (\theta)')
title('RK4 Approximation');
xlabel('Time (s)');
ylabel('Radian');
%Velocity
subplot(3,1,2);
plot(t,y2(:,2),'blue');
legend('Velocity (\theta/t)')
% title('RK4 Approximation');
xlabel('Time (s)');
ylabel('Radian/s');
% Acceleration
subplot(3,1,3);
plot(t,alpha_RK4,'green');
legend('Acceleration (\theta/t^2)')
% title('RK4 Approximation');
xlabel('Time (s)');
ylabel('Radian/s^2');
%------------------------------------------------
% Euler
figure(2)
%Displacement
subplot(3,1,1);
plot(t,y1(:,1),'red');
legend('Displacement (\theta)')
title('Eulers Approximation');
xlabel('Time (s)');
ylabel('Radian');
%Velocity
subplot(3,1,2);
plot(t,y1(:,2),'blue');
legend('Velocity (\theta/t)')
% title('Eulers Approximation');
xlabel('Time (s)');
ylabel('Radian/s');
% Acceleration
subplot(3,1,3);
plot(t,alpha_Euler,'green');
legend('Acceleration(\theta/t^2)')
% title('Eulers Approximation');
xlabel('Time (s)');
ylabel('Radian/s^2');
%------------------------------------------------
% Finite difference approach
%Without Drag
if Cd == 0
    T = 2*pi*sqrt(L/g);
    p = @(x) 0;
    q = @(x) (-g/L);
    r = @(x) 0;
    a = 0; c = T;
    x0 = [a c];
    dx = h;
    alpha = pi/2; beta = pi/2;
    bc =[alpha beta];
    [y,x] = bvpD(p,q,r,dx,x0,bc);
    figure(3)
    plot(x,y);
    xlim([min(x),max(x)])
    legend('Displacement (\theta)')
    title('Finite Difference approximation');
    xlabel('Time (s)');
    ylabel('Radian');
end
%------------------------------------------------


% tension in cable = m*g*cos(theta)+m*L*omega^2
T_Euler = (m*g).*cos(y1(:,1))+(m*L).*(y1(:,2).^2);
T_RK4 = (m*g).*cos(y2(:,1))+(m*L).*(y2(:,2).^2);

%plot tension in cable (Euler)
figure; hold on
plot(t,T_Euler','c.');
plot(t,m*g, 'k.');
xlabel('iteration');
ylabel('Tension (N)');
title('Tension in a pendulum cable as a function of time (Euler)');
hold off

%plot tension in cable (RK4)
figure; hold on
plot(t,T_RK4','c.');
plot(t,m*g, 'k.');
xlabel('iteration');
ylabel('Tension (N)');
title('Tension in a pendulum cable as a function of time (RK4)');
hold off
%------------------------------------------------
% Height difference
height = L-L.*cos(y2(:,1));
% Kinetic Energy 
ke = 0.5*m*(L.*y2(:,2)).^2;
pe = m*g*height;
loss = max(pe+ke) - (ke+pe);



% % Animation of RK4
fig4 = figure('Position', [100, 100, 768, 300]);

x =  L*sin(y2(:,1));
y = -L*cos(y2(:,1));


subplot(1,2,1)

hh1 = plot([0, x(1)],[0, y(1)], 'black.-', 'MarkerSize', 40, 'LineWidth', 2);
axis equal
axis([-1.5*L 1.5*L -1.5*L 0.3*L]);
ht = title(sprintf('Time: %0.2f sec', t(1)));

names = categorical({'PE','KE','Loss'});
names = reordercats(names,{'PE','KE','Loss'}); 
enrgy = [pe,ke,loss];
subplot(1,2,2)
hh2 = bar(names,enrgy(1,:),'FaceColor','flat');
hh2.CData(1,:) = [1 0 0];
hh2.CData(3,:) = [0 0.5 0];
ylim([0,max(enrgy(:,1))]);
title('Conservation of Energy')

% Get figure size
pos = get(gcf, 'Position');
width = pos(3);
height = pos(4);
framenum = 1:3:length(t);
% Preallocate
mov = zeros(height, width, 1, length(framenum), 'uint8');
frame = 1;
tic
for i =framenum
    set(hh1(1), 'XData', [0, x(i)]  , 'YData', [0, y(i)]);
    set(hh2, 'YData', enrgy(i,:));
    set(ht, 'String', sprintf('Time: %0.2f sec', t(i)));
%     drawnow;

   % Get frame as an image
   f = getframe(gcf);

   % Create a colormap for the first frame. For the rest of the frames, use
   % the same colormap
   if i == 1
      [mov(:,:,1,i), map] = rgb2ind(f.cdata, 256, 'nodither');
   else
      frame = frame+1;
      mov(:,:,1,frame) = rgb2ind(f.cdata, map, 'nodither');
   end
end
toc
imwrite(mov, map, 'animationENGDRAG.gif', 'DelayTime', 0, 'LoopCount', inf)