# Problems with modelling a pendulum

I've been trying to model a pendulum using the following diagram:

Where $\theta$ is the angle between the dashed line and $T$. Using this, and assuming $T = mg\cos(\theta)$, I can come up with two equations:

$$a_x = \frac{g}{2}\sin(2\theta)$$ $$a_y = g\sin^2(\theta)$$

Where $a_x$ is the horizontal acceleration, and $a_y$ is the vertical acceleration. However, when I put these equations into my program, the "weight" (circle bit) of the pendulum goes off in a straight line towards the bottom right corner of the screen.

Do the equations work or is it a problem with the program?

I would prefer to do this using Cartesian coordinates, since it is easier for the graphical side of the programming.

• Did you include the fact that the force of tension changes direction? Aug 20 '17 at 14:57
• @garyp Yes, I recalculate $\theta$ each time I move the pendulum Aug 20 '17 at 15:22
• Your assumption about T is incorrect. Aug 20 '17 at 16:27
• @ChesterMiller It is? Why? Aug 20 '17 at 17:01
• Do it in polar coordinates and do the obvious (and fast) coordinate transformation before drawing the thing.
– user107153
Aug 20 '17 at 22:39

In Cartesian coordinates: $$ma_x=m\frac{d^2x}{dt^2}=-T\sin{\theta}=-T\frac{x}{\sqrt{x^2+(h-y)^2}}\tag{1}$$ $$ma_y=m\frac{d^2y}{dt^2}=T\cos{\theta}-mg=T\frac{(h-y)}{\sqrt{x^2+(h-y)^2}}-mg\tag{2}$$where h is the elevation of the ceiling where the rope is attached.

To convert to radial and tangential coordinates, we write: $$x=r\sin{\theta}$$ $$(h-y)=r\cos{\theta}$$So, differentiating to get the velocities, we have: $$\frac{dx}{dt}=r\cos{\theta}\frac{d\theta}{dt}$$ $$\frac{dy}{dt}=r\sin{\theta}\frac{d\theta}{dt}$$ Differentiating again to get the accelerations, we have: $$\frac{d^2x}{dt^2}=-r\sin{\theta}\left(\frac{d\theta}{dt}\right)^2+r\cos{\theta}\frac{d^2\theta}{dt^2}\tag{3}$$ $$\frac{d^2y}{dt^2}=r\cos{\theta}\left(\frac{d\theta}{dt}\right)^2+r\sin{\theta}\frac{d^2\theta}{dt^2}\tag{4}$$ If we substitute Eqns. 3 and 4 into Equns. 1 and 2, we obtain: $$m\left[-r\sin{\theta}\left(\frac{d\theta}{dt}\right)^2+r\cos{\theta}\frac{d^2\theta}{dt^2}\right]=-T\sin{\theta}\tag{5}$$ $$m\left[r\cos{\theta}\left(\frac{d\theta}{dt}\right)^2+r\sin{\theta}\frac{d^2\theta}{dt^2}\right]=T\cos{\theta}-mg\tag{6}$$ If we multiply Eqn. 5 by $\cos{\theta}$ and Eqn. 6 by $\sin{\theta}$, and then add, we obtain: $$mr\frac{d^2\theta}{dt^2}=-mg\sin{\theta}\tag{7}$$ This is the force balance in the tangential direction. If we multiply Eqn. 5 by $-\sin{\theta}$ and Eqn. 6 by $\cos{\theta}$, and add, we obtain: $$mr\left(\frac{d\theta}{dt}\right)^2=T-mg\sin{\theta}\tag{8}$$ This is the force balance in the radial direction.

Instead of using the cathesian coordinates (horizontal and vertical axis), you should use the coordinate system "of the rope". The rope fixes the distance of the mass from the center of the rotation. Hence, decompose the focus into a component parallel and orthogonal to the rope. Only the latter contributes to the acceleration.

Which language do you use. Here the matlab code.

%% Solving the exact pendelum ode

clear all force;

g = 9.81; % gravitational constant [in m/s^2]
L = 1;    % length of the pendelum [in m]

%% define differential equation as a function handle
dPhidt = @(t,phi) [phi(2); -g/L * sin(phi(1))];

tVec    = [0:0.1:20];     % time [in s]
phi0    = [20*pi/180; 0]; % inital conditions
[t,Phi] = ode45(dPhidt, tVec, phi0);

figure(1)
plot(t,Phi(:,1) * 180/pi,'-o',t,Phi(:,2) * 180/pi,'-o')
title('Solution of the exact pendulum equation with ODE45');
xlabel('Time t [in s]');
grid on
legend('\phi [in deg]','\omega [in deg/s]')

%% transform polar coord into carthesian
[x, y] = pol2cart(Phi(:,1)-pi/2, L);

figure(2)
plot(x, y,'-o')
title('Solution of the exact pendulum equation with ODE45');
xlabel('Position A_x [in m]');
ylabel('Position A_y [in m]');
grid on


Note: The pendulum equation $$\ddot \varphi = -\frac{g}{L} \sin{(\varphi)}$$ which its a second order differential equation, is written in a set of two first order differential equations $$\begin{pmatrix} \dot \varphi \\ \dot \omega \end{pmatrix} = \begin{pmatrix} \omega \\ - \frac{g}{L} \sin{[\varphi]} \end{pmatrix} = \begin{pmatrix} phi(2) \\ - \frac{g}{L} \sin{[phi(1)]} \end{pmatrix}$$ where the right hand side displays the matlab code.

• This should be a comment, not an answer. In any event, the OP should respond and tell us if there is some reason he or she wants to do this in Cartesian coordinates. Aug 20 '17 at 14:58
• I use Cartesian coordinates simply because it is much easier for graphics Aug 20 '17 at 15:23
• That's of course true. However, polar coordinates are the "natural" coordinates of this problem, as you get the differential equation $$m \cdot l \cdot \ddot{\varphi}(t) = - m \cdot g \cdot \sin \left(\varphi(t) \right)$$ Transforming the solution back into cartesian coordinates is straight forward. Aug 20 '17 at 15:29

Presumably your weight is attached somewhere (by a rope?).

Consider the situation when the pendulum is not moving. In that case the forces on the bob/weight are (1) gravity directed downward and (2) a tension force directed upward along the strong. Both forces are of equal magnitude and opposite direction, that's why the sum of forces is zero and the pendulum does not move.

In the situation of a swinging pendulum you still have that the tension along the string compensates the component of the gravity along that axis ($mg\cos\theta$) and also the centripetal force. If this was not the case the length of the string would change during swinging.

So in the end, effectively you have only the component $mg\sin\theta$ of gravity acting on the bob, because the other forces cancel each other.

Since the pendulum moves on a circle, its motion is fully described by the angle $\theta$ and it will be much simpler to do calculations in such (polar) coordinates than in Cartesian (x-y) coordinates.

• This is incorrect. There is a centripetal acceleration in the radial direction that needs to be included in the radial force balance. Aug 20 '17 at 22:15
• @ChesterMiller: Thanks for pointing that out. I added the centripetal force in the answer. Aug 21 '17 at 8:50
• See my answer below. Aug 21 '17 at 10:52