Upgrading to geometric algebra my proof of energy conservation for a rotating rigid body in $D$ dimensions, and solving the Sylvester equation

Question

I wrote a proof from first principles that energy is conserved in a $D$-dimensional rotating rigid body without external forces, and I'd like to ask for some feedback on improving my math with more idiomatic notation, perhaps making use of geometric algebra; my background is in machine learning, meaning I tend to do a lot of linear algebra and statistics.

Say we have a $D$-dimensional rigid body composed of $N$ points of equal mass with coordinates $\mathbf x_i \ \forall i\in\{1\dots N\}$, vectors which together form the matrix $\mathbf X_{D \times N}$. These points are moving and have some velocities and accelerations, similarly forming the matrices $\mathbf V_{D \times N}$ and $\mathbf A_{D \times N}$ respectively.

For the sake of simplicity (and without loss of generality for rotational kinetic energy conservation), let's assume the center of mass of the body is fixed at the origin, and add the assumption that we have no outside forces. Then we can write: $$\mathbf{X1}=\mathbf 0 \\ \mathbf{V1}=\mathbf 0\\ \mathbf{A1}=\mathbf 0$$ where in these equations, $\mathbf 1$ and $\mathbf 0$ are vectors composed of all ones and all zeros respectively.

There are three "first-principles" I will use:

• constant distances between point pairs: $ \frac{\mathrm d}{\mathrm d t} \| \mathbf x_i - \mathbf x_j \| =0 \ \forall i,j $
• all forces are central: $ \mathbf a_i = \sum_j k_{ij} (\mathbf x_i - \mathbf x_j) $
• Newton's (strong?) 3rd Law: $ k_{ij} = k_{ji} $

I'll go over the main results without going into too much detail with the algebra. First, I define the matrices:

$$\mathbf R_{N \times N} = \mathbf{X^T X} \\ \mathbf J_{D \times D} = \mathbf{X X^T} \\ \mathbf L_{D \times D} = \mathbf{X V^T - V X^T} $$

The matrix $\mathbf J$ is homologous to the moment of inertia, but since I'm not using cross-products, I end up with a simpler form. And the matrix $\mathbf L$ is homologous to the angular momentum, which I've recently learned I might be able to write using a wedge product as $$\mathbf L \propto \left<\mathbf x \wedge \mathbf v \right> $$ where I use the stats notation $\left<\mathbf x \right> = \frac{1}{N} \sum_{i=1}^n{\mathbf x_i}$. And I can prove angular momentum is conserved from Newton's 3rd Law and the force centrality.

Here is the first problem I'm encountering: while this pattern with the minus looks like a wedge product in geometric algebra, I also get a lot of a similar pattern but with plus instead of minus; in particular, I use the fact that $$\frac{\mathrm d}{\mathrm d t}\mathbf R = \mathbf{X^T V} + \mathbf{V^T X} = \mathbf 0$$ to prove that there exists a unique skew-symmetric matrix $\mathbf W_{D \times D}$ such that $\mathbf V = \mathbf {WX}$, which I use to represent angular velocity. Plus, I make use of both $[\mathbf{X^T X}]_{ij} = \mathbf{x}_i \cdot \mathbf{x}_j $ as well as $ \mathbf{X X^T} \propto \left<\mathbf{x x^T} \right> $ in my equations, but from what I've read so far, geometric algebra doesn't have both an inner product and an outer product in addition to the wedge product, so I'm unsure how to make all my equation work in the framework of geometric algebra.

Continuing the outline of my proof, I show that $$\mathbf L = -(\mathbf{JW}+\mathbf{WJ}) = \mathbf{JW^T}-\mathbf{WJ^T} $$ I guess this one I can write both as a sum or as a difference, since $\mathbf J$ is symmetric, and $\mathbf W$ is skew-symmetric.

Since angular momentum is conserved, $$\dot{\mathbf L} = -\mathbf{\dot{J}W}-\mathbf{J\dot{W}} - \mathbf{\dot{W}J}-\mathbf{W\dot{J}} = \mathbf 0$$ and I can use this to solve for the angular acceleration $\dot{\mathbf W}$, as I get the Sylvester equation $$\mathbf{J\dot{W}} + \mathbf{\dot{W}J} = -\mathbf{\dot{J}W}-\mathbf{W\dot{J}}$$

Alternatively, I could use the pattern with minuses to write

$$\mathbf{J\dot{W}^T} - \mathbf{\dot{W}J^T} = \mathbf{\dot{J}W^T}-\mathbf{W\dot{J}^T}$$

but this is not a Sylvester equation anymore, and I don't know how to solve it. Is there an idiomatic way to write this equation using wedge products? And how would I solve it once it's in terms of wedge products?

Finally, I can complete my proof of energy conservation by showing that I can compute the total energy of the system using the matrix trace in the equation $$\mathrm{Tr}[\mathbf{LW}]=2\mathrm{Tr}[\mathbf{V^T V}]$$ then taking the derivative of the left side to show that $$\frac{\mathrm d}{\mathrm d t}\mathrm{Tr}[\mathbf{LW}]=\mathrm{Tr}[-\mathbf{(WW^T)L}]=0$$ since the trace of the product of a symmetric matrix with a skew-symmetric matrix is zero, but I'm not using the wedge-like looking pattern of something minus its transpose anymore. And I can find the acceleration matrix for the individual particles with $$\mathbf{A} = \mathbf{WV} + \dot{\mathbf{W}}\mathbf{X}$$

My two questions are:

• Can I rewrite all of this using geometric algebra, despite making use of both inner and outer products?
• What happens to my Sylvester equation once I use geometric algebra? How do I solve it?

Answer 1

Technically, Lagrange/Hamiltonian formulations do not require energy conservation. However, if the equations of motion are independent with time, then they do require it by Noether's theorem. From this perspective, the key underlying axiom when defining your system is that it is autonomous, which Lagrangian/Hamiltonian mechanics emphasise more transparently. Granted, the mere assumption that there exists a Lagrangian is quite restrictive for the equations of motion.

Furthermore, you are assuming Newton's laws of motion, which may be a bit restrictive. The Lagrangian perspective explains where they come from using general symmetry arguments and allow more general laws of motion.

Yes, Clifford algebras can certainly be used in this case, but the true structure that you are interested in is the Lie algebra of $SO(D)$ which is a special case. You depending on your applications, it might be computationally expensive to use the former.

The formalism can be sketched as the following. Your configuration space can be more compactly described with $SO(D)$, the set of rigid rotations in $D$ dimensions, rather than $\mathbb R^{ND}$. Note that the dimension is now independent of $N$, and it is an improvement (being smaller) when $D\leq 2N+1$. Of course, the $N$ dependence enters again when calculating the moment of inertia, but it can be either seen as an amortised cost since you need do it once for a simulation. A configuration will therefore be represented by a rotation $R$ with $RR^T = 1$ (which can be in turns be represented as a matrix, Euler angles, unit quaternions when $D=3$ etc.).

For the Lagrangian formulation, you assume the existence of the Lagrangian of the form $L(R,\dot R,t)$. From the assumption of autonomy, you can remove the explicit $t$ dependence. Note that this leaves more alternatives for that free rotation is. Typically, you would need to impose invariance by left multiplication of the Lagrangian which will give conservation of angular momentum. Now, it is of the form $L(\Omega)$ with: $$ \Omega = R^T\dot R $$ which is an element of the Lie algebra of $SO(D)$. Physically, it is angular velocity in the frame of the moving rigid body. Newton's laws only considers a quadratic Lagrangian in $\Omega$, the "prefactor" being the moment of inertia. However, it could be of higher order, and Newton's laws can be seen as a harmonic approximation of a more general theory.

Hope this helps.

Answer to comment

For your first comment, I will detail some computations. Actually your $W$ related to my $\Omega$, which physically represents angular velocity. You'll need to be a bit careful though. There are two different angular velocities: the one in the reference frame of the moving solid or the one in the inertial reference frame. Mathematically, from the group defining relation: $$ RR^T = 1 $$ You get by derivation: $$ R\dot R^T+\dot R R^T = 0 $$ So you can naturally define two skew symmetric matrices (guaranteeing existence and unicity): $$ \Omega = R^T\dot R \quad \omega = \dot R R^T $$ the former being the intrinsic angular momentum and the latter being the extrinsic angular momentum (which ressembles more your $W$). This interpretation is justified by the fact that they are related by a conjugation by $R$. Mathematically, they represent the push forward of the Lie group's tangent spaces to the Lie algebra in the first case by using left multiplication and in the second case using right multiplication.

From the Lagrange formulation, you obtain another skew symmetric matrix (or abstractly another element of the Lie algebra): $$M = \frac{\partial L}{\partial \Omega}$$ which you identify as angular momentum which depends on $\Omega$ (in your case, you are assuming a linear relation, but it can be more complicated). The Euler Lagrange equation gives you (Euler's equation): $$ \dot \Omega = [\Omega,M] $$ with $[\Omega,M]$ the Lie bracket. If you are only interested in the case where they are represented by skew symmetric matrices, $[\Omega,M] = \Omega M-M\Omega$. You can deduce from this the conservation of the Hamiltonian: $$ H = \text{Tr}(M\Omega^T)-L $$ as well as conservation of extrinsic angular momentum (which corresponds to your $L$, and can be obtained abstractly from Noether's theorem using invariance by left multiplication): $$ m = RMR^T $$ If you want more details on the abstract Lie algebra approach, check out Arnold's Mathematical Methods of Classical Mechanics (appendix 2, there is also an application to fluid dynamics).

For your second comment, my discussion is purely theoretical and is axiomatic as well. I just chose a different set of axioms. I take for granted the existence of a Lagrangian, the same way you take for granted Newton's laws. We both make assumptions, and, thankfully for this problem, we arrive at the same conclusion. It is therefore enlightening to have these two perspectives as they emphasise different aspects of the problem.

For your final comment, it is true that the Lagrangian formalism has a harder time dealing with non-conservative forces. As I said before, the existence of a Lagrangian is quite restrictive, and in some sense, the "game is rigged" in favour for energy conservation. However, once you've established that you can describe your system with a Lagrangian, it is a powerful method of resolution, especially when going to Hamiltonian mechanics. Furthermore, it is the default choice in physics when trying to infer new laws of motion.

Answer 2

In the Geometric Algebra picture, the primary operator of interest between two vectors is the geometric product that is defined as the sum of the inner product and the wedge (exterior) product of the two vectors: $$ab=a\cdot b+a\wedge b$$ One can use the fact that reversing the operators yields, $$ba=b\cdot a+b\wedge a$$ then since $a\cdot b\equiv b\cdot a$ and $a\wedge b=-b\wedge a$, then we can define the inner and wedge products in terms of the geometric product: \begin{align} a\cdot b&=\frac{1}{2}\left(ab+ba\right) \\ a\wedge b&=\frac{1}{2}\left(ab-ba\right) \end{align} You can also take a projection of the blade (i.e., the grade) to yield these terms as well: \begin{align} a\cdot b&=\langle ab\rangle_0 \\ a\wedge b&=\langle ab\rangle_2\\ \end{align} where the $0$-subscript indicates the scalar component and $2$ indicates the bivector component of the multivector. It is also somewhat common that having no subscript indicates the grade-0 projection (i.e., $\langle ab\rangle\equiv\langle ab\rangle_0$).

So despite your insistence, the inner and wedge (exterior) products are not only well-defined in the geometric algebra picture, but are requisite.

---

In terms of the kinematics of a rotating rigid body, this is well covered in Section 3.4 of Doran & Lasenby's Geometric Algebra for Physicists. Essentially, you build up the algebra to discuss the rotor, such that a rotation of a vector $a$ about the bivector (plane) $\hat{B}\equiv\hat{\mathrm{e}}_i\hat{\mathrm{e}}_j$ is, $$a'=Ra\tilde{R}$$ with $R=\exp(-\hat{B}/2)$.

If you have a frame $\hat{e}_k$ in a fixed plane, you can then apply the time-dependent rotor $R(t)$ to obtain the rotated plane, $$\hat{f}_k=R(t)\hat{e}_k\tilde{R}(t)$$ which has a time-derivative, $$\dot{\hat{f}_k}=\dot{R}\hat{e}_k\tilde{R}+R\hat{e}_k\dot{\tilde{R}}=\dot{R}\tilde{R}\hat{f}_k+\hat{f}_kR\dot{\tilde{R}}$$ which, since $\dot{R}\tilde{R}=-R\dot{\tilde{R}}=-\left(\dot{R}\tilde{R}\right)^\sim$ means that $2\dot{R}\tilde{R}=-\Omega$ determines your angular momentum. Note also that multiplying on the right by $R$ yields, $$\dot{R}=-\frac{1}{2}\Omega R$$ so that if $\Omega=\text{const}$, then $R(t)=R_0\exp(-\Omega t/2)$.

An earlier section of the text covers the motions of a collection $N$ particles in the geometric algebra framework, utilizing a center-of-mass approach, if I recall correctly, so I suspect that the text would be worth your time investment in reading.

---

For your starting point, you probably want to completely erase the notion of matrices from your mathematical background and stick with vectors and multivectors, which constitute the basis of the algebra. As such, there would be no such thing as a Sylvester equation to solve, since we have no matrices.