# Lagrangian mechanics formulation of a simple free motion of two masses in uniform gravity field

As a part of larger project, I decided to test my Lagrangian formulation of simple system of two rigidly connected point masses as indicated below.

I introduce the generalized coordinates vector $$\bar{q} = [x, y, q]^T$$, where $$x,y$$ are coordinates of position of mass $$M$$ and $$q$$ is the angular deflection from horizontal of the massless rod of length $$l$$, in the middle of which mass $$m$$ is attached. Then, the position coordinates of $$m$$ are given by $$x_m = x + \dfrac{l}{2}\cos(q)$$ $$y_m = y + \dfrac{l}{2}\sin(q)$$ Having this, I defined kinetic ($$T$$) and potential ($$V$$) energies $$T = \dfrac{1}{2}M(\dot{x}^2 + \dot{y}^2) + \dfrac{1}{2}m(\dot{x}_m^2 + \dot{y}_m^2)$$ $$V = g \left( My + m\left(y+\dfrac{l}{2}\sin(q)\right) \right)$$ Then, the Lagrangian $$L$$ and mechanical energy of the system $$E$$ are $$L = T-V$$ $$E = T+V$$

Further, the derivations were done by my by hand (a few times to double check) and using SymPy (symbolic math package of Python) and they match. Upon integrating the EOM, I receive a nicely looking plots for state vector variables $$\bar{Q} = [x,y,q,\dot{x},\dot{y},\dot{q}]^T$$, if my initial conditions involve null $$q$$ nd $$\dot{q}$$. However, if they are non-zero, the motion looks highly non-physical visually (I created a Matlab animation of the motion) and additionally the mechanical energy of the system is not adding up to be constant (see pic below for $$\bar{Q}_0 = [1,1000,1,1,0,-0.1]^T$$).

My question is therefore, are my assumptions and initial formulation correct for the given situation? I am quite new to Lagrangian formalism and afters hours of tackling this seemingly simple problem, started to think there might be some rule of Lagrangian mechanics I violate. All classical mechanics problems I was able to find online are not free motion of multibody systems, therefore I was not able to check with actual solution for problem like this.

To give some more context, my intention is to expand this problem to put the multibody system in orbital flight (central point-mass gravity field) and increase complexity of structure, leading to n-link pendulum attached to bus mass $$M$$. If there is a fundamental error in the proceedings I presented above here, I assume I would not be able to extend the formulation to more complex system I just described.

I will be eternally grateful for any help / advice! I will be also happy to provide more details on my solution.

EDIT To clarify, I post my EL equations: $$M \frac{d^{2}}{d t^{2}} x{\left(t \right)} - 0.5 m \left(l \sin{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} q{\left(t \right)} + l \cos{\left(q{\left(t \right)} \right)} \left(\frac{d}{d t} q{\left(t \right)}\right)^{2} - 2 \frac{d^{2}}{d t^{2}} x{\left(t \right)}\right)= 0$$ $$M \frac{d^{2}}{d t^{2}} y{\left(t \right)} + g \left(M + m\right) + 0.5 m \left(- l \sin{\left(q{\left(t \right)} \right)} \left(\frac{d}{d t} q{\left(t \right)}\right)^{2} + l \cos{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} q{\left(t \right)} + 2 \frac{d^{2}}{d t^{2}} y{\left(t \right)}\right) = 0$$ $$l m \left(0.5 g \cos{\left(q{\left(t \right)} \right)} + 0.25 l \frac{d^{2}}{d t^{2}} q{\left(t \right)} - 0.5 \sin{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} x{\left(t \right)} + 0.5 \cos{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} y{\left(t \right)}\right) = 0$$

$$\begin{bmatrix} M+m & 0 & -\dfrac{l}{2}m\sin(q) \\ 0 & M+m & \dfrac{l}{2}m\cos(q) \\ -\dfrac{l}{2}m\sin(q) & \dfrac{l}{2}m\cos(q) & m\dfrac{l^2}{4} \end{bmatrix} \begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \ddot{q} \end{bmatrix} = \begin{bmatrix} \dfrac{l}{2}m\cos(q)\dot{q}^2 \\ -g(M+m) + \dfrac{l}{2}m\sin(q)\dot{q}^2 \\ -\dfrac{l}{2}gm\cos(q) \end{bmatrix}$$

• You have to expand the expression of $T$ to explicitly see $\dot{q}$... Aug 7 '20 at 10:50
• Sure, the kinetic energy is then $$T = \dfrac{1}{2}M(\dot{x}^2 + \dot{y}^2) + \dfrac{1}{2}m\left( \left( \dot{x} - \dfrac{l}{2}\sin(q)\dot{q} \right)^2 + \left( \dot{y} + \dfrac{l}{2}\cos(q)\dot{q} \right)^2 \right)$$ My problem does not come from inability to turn this into a numerically integratable problem. I suppose there is some mistake at the very beginning, in principle of my formulation. Aug 7 '20 at 10:56
• Your theoretical approach is correct. The problems seem to arise when passing to numerical computations. Mechanical energy must be conserved so that the corresponding plot is completely wrong! Aug 7 '20 at 11:39
• Maybe there is a mistake in your Euler-Lagrange equations. Could you post them? By the way, if I am not wrong, these equations are exactly solvable, and if you place a coordinate system at the center of mass of the two mass-points, then the Lagrangian will decouple and it will split into a parabolic motion of the center of mass and a uniform rotation around the center of mass. Aug 7 '20 at 21:27
• Thank you @Futurologist, I've edited the post and added the EL equations. Could you say a bit more about decoupling of equations you mentioned? What steps are necessary to perform it? Aug 8 '20 at 15:53

Remark. This case generalizes to your real goal to have a system of more than two mass-points attached to each other by rods where the rods can rotate relative to each other. Basically, in a uniform force field like the gravitational force field near the surface of the earth the problem becomes equivalent to the geodesic motion of a point onto a torus of dimension equal to the number of rods.

Lagrangian approach. First, you place an inertial coordinate frame (say a frame attached to the ground) in the plane of the system. Denote it by $$O\,\vec{e}_x\,\vec{e}_y$$, where the coordinate vectors are of unit length and perpendicular to each other. Moreover, $$O\,\vec{e}_x\,\vec{e}_y$$ is oriented so that $$\vec{e}_y$$ is the vertical vector along which gravitational acceleration is $$\vec{g} = - \, g \, \vec{e}_y$$. The position vectors of the mass-points, originating from point $$O$$ and ending at the either of the mass points, are $$\vec{r}_M$$ and $$\vec{r}_m$$. The Lagrangian is then $$L\, =\, \frac{M}{2}\,\left|\frac{d\vec{r}_M}{dt}\right|^2 \, + \, \frac{m}{2}\,\left|\frac{d\vec{r}_m}{dt}\right|^2 \, - \, Mg\,\big(\vec{r}_M \cdot \vec{e}_y\big) \, - \, mg\,\big(\vec{r}_m \cdot \vec{e}_y\big)$$ with holonomic constriant $$\big|\vec{r}_M \, - \, \vec{r}_m\big|^2 \, = \, \frac{l^2}{4}$$

Perform the linear chainge of coordinates: \begin{align} &\vec{r}_G \, =\, \frac{M}{M+m} \, \vec{r}_M \, + \, \frac{m}{M+m} \, \vec{r}_m\\ &\vec{r}_l \, = \, \vec{r}_M \, - \, \vec{r}_m \end{align} with inverse \begin{align} &\vec{r}_M \, =\, \vec{r}_G \, + \, \frac{m}{M+m} \, \vec{r}_l\\ &\vec{r}_m \, = \, \vec{r}_G \, - \, \frac{M}{M+m} \, \vec{r}_l \end{align} In the new coordinates, the Lagrangian becomes \begin{align} L \, =& \, \frac{M}{2}\,\left|\frac{d\vec{r}_G}{dt} \, + \, \frac{m}{M+m} \, \frac{d\vec{r}_l}{dt}\right|^2 \, + \, \frac{m}{2}\,\left|\frac{d\vec{r}_G}{dt} \, - \, \frac{M}{M+m} \, \frac{d\vec{r}_l}{dt}\right|^2 \\ &\, - \,g\,\Big(\big(M\,\vec{r}_M \, + \, m \, \vec{r}_m \big)\cdot \vec{e}_y\Big) \\ \, =&\, \frac{M}{2} \left|\frac{d\vec{r}_G}{dt}\right|^2 \, + \, \frac{M\,m}{(M+m)}\left(\frac{d\vec{r}_G}{dt}\cdot\frac{d\vec{r}_l}{dt}\right) \, + \, \frac{M\,m^2}{2(M+m)^2}\, \left|\frac{d\vec{r}_l}{dt}\right|^2\\ &\, + \, \frac{m}{2} \left|\frac{d\vec{r}_G}{dt}\right|^2 \, - \, \frac{M\,m}{(M+m)}\left(\frac{d\vec{r}_G}{dt}\cdot\frac{d\vec{r}_l}{dt}\right) \, + \, \frac{M^2\,m}{2(M+m)^2} \left|\frac{d\vec{r}_l}{dt}\right|^2\\ &\, - \,(M+m)\,g\,\Big(\vec{r}_G \cdot \vec{e}_y\Big) \end{align} and after some regrouping and simplification the Lagrangian and the holonomic constraint become: $$L \, = \, \frac{M\,m}{2(M+m)}\, \left|\frac{d\vec{r}_l}{dt}\right|^2 \, + \, \frac{(M+m)}{2}\, \left|\frac{d\vec{r}_G}{dt}\right|^2 \, - \, (M+m)\,g\,\Big(\vec{r}_G \cdot \vec{e}_y\Big)$$ $$\big|\vec{r}_l\big|^2 = \frac{l^2}{4}$$ Introduce the generalized coordinates: $$\vec{r}_l \, = \, \frac{l}{2}\, \cos(q) \, \vec{e}_x \, + \, \frac{l}{2}\, \sin(q) \, \vec{e}_y$$ $$\vec{r}_G \, = \, x_G \, \vec{e}_x \, + \, y_G\, \vec{e}_y$$ which parametrixe the holonimc constrints (they tirivialize it). Then, the Lagrnagian simplifies to $$L \, = \, \frac{M\,m\ l^2}{8(M+m)}\, \left(\frac{dq}{dt}\right)^2 \, + \, \frac{(M+m)}{2}\, \left(\frac{dx_G}{dt}\right)^2 \, + \, \frac{(M+m)}{2}\, \left(\frac{dy_G}{dt}\right)^2 \, - \, (M+m)\,g\,y_G$$ Consequently, the Euler-Lagrange equations are \begin{align} &\frac{d^2{x}_G}{dt^2} \, = \, 0 \\ &\frac{d^2{y}_G}{dt^2} \, = \, - \, g\,\\ & \frac{d^2q}{dt^2} \, = \, 0 \end{align} You can solve them right away \begin{align} &x_G \, = \, x_G(0) \, + \, v_x\, t \\ &y_G \, = \, y_G(0) \, + \, v_y\, t \, - \, \frac{g}{2} \, t^2\,\\ &q \, = \, q(0) \, + \, \omega_0 \, t \end{align} where $$x_G(0), \, y_G(0), \, q(0)$$ are constants describing the initial configuration of the rod and $$v_x, \, v_y, \, \omega_0$$ are constants describing the initial velocities of the rod.

The rigid body interpretation. The special case you are modelling in your OP can be thought of as a solid body, basically a rigid rod, which moves in a plane and whose mass is concentrated at its ends. First, you place an inertial coordinate frame (say a frame attached to the ground). Denote it by $$O\,\vec{e}_x\,\vec{e}_y\,\vec{e}_z$$, where the coordinate vectors are of unit length and perpendicular to each other. Moreover, $$O\,\vec{e}_x\,\vec{e}_y\,\vec{e}_z$$ is oriented so that $$\vec{e}_y$$ is the vertical vector along which gravitational acceleration is $$\vec{g} = - \, g \, \vec{e}_y$$, the system moves in the coordinate plane $$O\,\vec{e}_x\,\vec{e}_y$$, just like on your picture, and $$\vec{e}_z$$ is pointing perpendicularly to the picture.

Denote by $$\vec{r}_M$$ the position vector, pointing from the origin $$O$$ of the coordinate frame to the mass-point $$M$$ and by $$\vec{r}_m$$ the position vector, pointing from the origin $$O$$ of the coordinate frame to the mass-point $$m$$. The condition that the two mass points are connected by a mass-less rod of length $$\frac{l}{2}$$, can be written as the quadratic equation. $$\big|\vec{r}_M - \vec{r}_m\big|^2 = \frac{l^2}{4}$$ Let $$G$$ be the center of gravity of the rod, i.e. the center of gravity, which is $$\vec{r}_G = \frac{M}{M+m} \, \vec{r}_M \, + \, \frac{m}{M+m} \, \vec{r}_m$$ The force that acts on the system of the two mass points on the rod is gravity and is applied to the center of mass $$\vec{r}_G$$ of the rod. The general equations of motion of a rigid body can be written as \begin{align} &\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\ & \frac{d}{dt}\,I \,\omega\, \vec{e}_z \, = \, \vec{T} \end{align} where $$\omega \, \vec{e}_z$$ is the angular velocity of the rod and $$I = \frac{M\,m\,l^2}{4\,(M+m)}$$ is its moment of inertia along the $$O\, \vec{e}_z$$ axis with respect to the center of mass $$G$$ of the rod. The other moments of inertia are zero because this is a one dimensional rod. The vector $$\vec{T}$$ is the sum of the torques of all forces acting on the rod, calculated with respect to the center of mass $$G$$. However, there is only one force, the gravity force $$-(M+m)\, g \, \vec{e}_y$$, acting on the rod and it is applied to the center of gravity $$G$$. Since $$\vec{T} = \vec{GG} \times \big(-(M+m)\, g \, \vec{e}_y\big) = \vec{0}\times \big(-(M+m)\, g \, \vec{e}_y\big) = \vec{0}$$ i.e. the torque is zero, the equations of motion become \begin{align} &\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\ & \frac{d}{dt}\,I \,\omega\, \vec{e}_z \, = \, \vec{0} \end{align} As you can see, these equations are decoupled and can be written as \begin{align} &\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\ & \,I \, \frac{d\omega}{dt} \, = \, {0} \end{align} But, the angular velocity, in your notations, is simply $$\omega = \frac{dq}{dt}$$ and also
$$\vec{r_G} = x_G \, \vec{e}_x + y_G\, \vec{e}_y$$ so \begin{align} &\frac{d^2{x}_G}{dt^2} \, = \, 0 \\ &\frac{d^2{y}_G}{dt^2} \, = \, - \, g\,\\ & \,I \, \frac{d^2q}{dt^2} \, = \, {0} \end{align} You can solve them right away \begin{align} &x_G \, = \, x_G(0) \, + \, v_x\, t \\ &y_G \, = \, y_G(0) \, + \, v_y\, t \, - \, \frac{g}{2} \, t^2\,\\ &q \, = \, q(0) \, + \, \omega_0 \, t \end{align} where $$x_G(0), \, y_G(0), \, q(0)$$ are constants describing the initial configuration of the rod and $$v_x, \, v_y, \, \omega_0$$ are constants describing the initial velocities of the rod.

If I have time, I can get to these equations from Lagrangian point of view too. Basically you just have to change coordinates of your Lagrangian and the holonomic constrant $$\big|\vec{r}_M - \vec{r}_m\big|^2 = \frac{l^2}{4}$$ like this: \begin{align} &\vec{r}_G \, =\, \frac{M}{M+m} \, \vec{r}_M \, + \, \frac{m}{M+m} \, \vec{r}_m\\ &\vec{r}_l \, = \, \vec{r}_M \, - \, \vec{r}_m \end{align}

Remark. This case does not generalize to your real goal to have a system of more than two mass-points attached to each other by rods where the rods can rotate relative to each other.

• I do understand that the OP asked you to elaborate on the change of variables in a comment, and that it would be impossible to do so by writing another comment, but does this in any way (numeric, Lagrangian, possible mistakes) answer the actual question of the OP? Aug 8 '20 at 18:08
• @GiorgioComitini pointed out correctly - the answer of Futurologist is very valuable and insightful (thanks a lot for that!), however I am mostly interested in Lagrangian approach to this problem and extending the problem to (what you correctly understood) multiple point masses attached to massless rods. Aug 8 '20 at 21:22
• @AleksanderFiuk I added a Lagrangian approach that simplifies the system and can generazie to your goal. Aug 9 '20 at 19:01

This is my solution for your problem:

The Kinetic Energy is:

$$T=\frac{1}{2}\,M(\dot x^2+\dot y^2)+\frac 12\,m\,\left[\left(\dot x -\frac 12 l \sin(\varphi)\,\dot \varphi\right)^2 +\left(\dot y +\frac 12 \,l\cos(\varphi)\dot \varphi\right)^2\right]$$

and the potential energy

$$U=M\,g\,y+m\,g\,\left(y +\frac 12 l\sin(\varphi)\right)$$

the generalized coordinates are:

$$\vec q=\begin{bmatrix} x \\ y \\ \varphi \\ \end{bmatrix}$$

The EOM's:

$$\ddot x-\frac 12 \frac{\cos(\varphi)\,m\,l}{M+m}\,\dot \varphi^2=0$$

$$\ddot y-\frac 12 \frac{\sin(\varphi)\,m\,l}{M+m}\,\dot \varphi^2+g=0$$ $$\ddot \varphi=0$$

the numerical simulation give you for $$E=T+U$$

the initial conditions are all zero except for $$~\varphi(0)=0.3$$

thus: if you have this equations, the total energy $$E=T+U~$$ is constant as it should be.

edit

in case that $$\varphi(0)=\varphi_0$$ and all other initial conditions are zero, the solution of the equations of motion are:

$$x(\tau)=0~,y(\tau)=-\frac 12\,g\,\tau^2~,\varphi(\tau)=\varphi_0$$ and $$E=\frac 12 m\,g\,l\,\sin(\varphi_0)=~\text{constant}$$

The equations of motion:

from Euler Lagrange you obtain this three equations

$$\left( M+m \right) { \ddot x}-1/2\,ml\sin \left( \varphi \right) \ddot \varphi -1/2\,ml\cos \left( \varphi \right) {\dot \varphi }^{2} =0$$

$$\left( M+m \right) {\ddot y}+1/2\,ml\cos \left( \varphi \right) \ddot\varphi -1/2\,ml\sin \left( \varphi \right) {\dot \varphi }^{2}+M g+mg=0$$

and

$$-1/4\,ml \left( 2\,\sin \left( \varphi \right) { \ddot x}-2\,\cos \left( \varphi \right) {\ddot y}-\ddot\varphi \,l-2\,g\cos \left( \varphi \right) \right) =0$$

solve those equations for $$\ddot x~,\ddot y~,\ddot \varphi$$ you obtain my equations above

• Could you clarify how did you arrive to $\ddot{\varphi} = 0$? Did you follow EL equation $\dfrac{\partial L}{\partial q_i} - \dfrac{\partial}{\partial t}\dfrac{\partial L}{\partial \dot{q}_i} =0$ ? I agree with the energy plot, mine should look like that as well in the end. Actually, if my initial conditions are like yours, I also obtain constant mechanical energy, but maybe could you try with non-zero $q(0) \neq 0$and $\dot{q}(0) \neq 0$? Aug 7 '20 at 13:56
• I follow EL equation and then solve for $\ddot{x}~,\ddot{y}~,\ddot{\varphi}$. with arbitrary IC what so ever, I obtain constant total energy as it should be!, this is also a good test, if not something is wrong with your EOM's
– Eli
Aug 7 '20 at 14:33
• I have tried to recreate your solution but I get lost at some point. I have actually managed to get your solution by mistake, when my $q$ (your $\varphi$) was not a differentiable variable, but a constant. Indeed, the energy is conserved in that case. Could you please verify if you have not made a similar mistake? If not, could you elaborate a bit more on how you arrived to the final equations? Aug 8 '20 at 15:55
• @AleksanderFiuk I wrote you the equations
– Eli
Aug 8 '20 at 21:18
• Thank you, @Eli. I arrive to the same equations, which basically solves the issue of this post. There must be a silly mistake somewhere on the way to numerical implementation. Aug 8 '20 at 21:34