𝑁 -Dimensional Rigid Body Dynamics
MARC TEN BOSCH, mtb design works, Inc., USA

Fig. 1. A stack of three 4D hypercubes
I present a formulation for Rigid Body Dynamics that is independent of
the dimension of the space. I describe the state and equations of motion of
rigid bodies using geometric algebra. Using collision detection algorithms
extended to 𝑛D I resolve collisions and contact between bodies. My implementation is 4D, but the techniques described here apply to any numberof dimensions. I display these four-dimensional rigid bodies by taking athree-dimensional slice through them. I allow the user to manipulate thesebodies in real-time.

of them. While attention has been given to understanding and visualizing high-dimensional spaces and other abstract mathematicalconcepts, it has most often remained limited to visualizing theseconcepts without any physicality or object-to-object relationships.
Contributions. The contributions of this paper include:(1) extending the geometric algebra-based formulation of classical 3D rigid body dynamics to 𝑛D. By representing geometricalgebra operators as matrices, one can in particular construct,diagonalize and transform the inertia tensor for arbitrary𝑛D simplicial meshes for any 𝑛, in a simple way. This enables formulating the Euler equation in 𝑛D, which allows e.g.
studying the case of the 4D Euler equation under torque-free
conditions.
(2) computing collision and contact processing in 𝑛D, including
static and kinetic friction. I give an 𝑛D formulation for the
Minkowski difference and separating axis theorem collision
detection methods based on geometric algebra.
(3) a method of interacting with 4D objects that is akin to our
3D experience of reality.

CCS Concepts: • Computing methodologies → Physical simulation; Collision detection; Virtual reality.
Additional Key Words and Phrases: Rigid Bodies, 𝑁 -Dimensional, Physics,

Geometric Algebra, Fourth Dimension
ACM Reference Format:
Marc ten Bosch. 2020. 𝑁 -Dimensional Rigid Body Dynamics. ACM Trans.
Graph. 39, 4, Article 1 (July 2020), 6 pages. https://doi.org/10.1145/3386569.
3392483

1

INTRODUCTION

Our experience of physical space is three dimensional. Consequently,
physically-based simulations (physics engines) have so far been focused on and restricted to the two and three-dimensional cases.
However using the appropriate formulation of the required equations it is possible to extend them to higher dimensions. Geometricalgebra provides a simple dimension-independent formulation. This
allows real time manipulation of 𝑛-dimensional shapes that collidewith each other as if they were real objects, which makes themmuch less abstract, in stark contrast with most people’s experience
Author’s address: Marc ten Bosch, mtb design works, Inc., San Francisco, CA, USA,marc@marctenbosch.com.
Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than theauthor(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or
republish, to post on servers or to redistribute to lists, requires prior specific permission
and/or a fee. Request permissions from permissions@acm.org.
© 2020 Copyright held by the owner/author(s). Publication rights licensed to ACM.
0730-0301/2020/7-ART1 $15.00
https://doi.org/10.1145/3386569.3392483

2

RELATED WORKInteractive Simulation of 3D Rigid Body Dynamics is a broad field. Bender et al. [2014] provide a survey.
[Cameron 1990] has formulated the 3D continuous collision
detection problem as a discrete 4D collision problem by considering
the extrusion of each object over time.
Visualizing 4D objects is an interesting and challenging problem
of its own, with a long history [Abbott 1884; Banchoff 1990; Chu et al.
2009; Hilbert and Cohn-Vossen 1952]. Many ways of manipulating
4D shapes have been proposed [Yan et al. 2012; Zhang and Hanson
2006].
Geometric algebra operates on a space of elements called multivectors, of which vectors are a subspace. While vectors can bethought of as oriented line-segments, other elements can represent
oriented areas (bivectors), volumes (trivectors), and so on. Geometric algebra also defines operations that correspond to translations
and rotations of these elements. The book by Macdonald [2011]
ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

1:2

•

ten Bosch, M.

and the course by Gunn and De Keninck [2019] provide introductions. Dorst et al. [2009] provide an introduction to implementing
geometric algebra in programs. Doran and Lasenby [2003] apply
geometric algebra to the case of 3D rigid bodies.
Cayley [1846] first proposed generalizing the Euler equation
to 𝑛D. The problem has been studied analytically primarily using
matrix analysis, which becomes complex quite quickly (e.g. [Sinclair
and Hurtado 2005]).
I chose to use geometric algebra due to its simplicity, generality,
and coordinate-free nature, which allows the equations of motions
to remain the same as in the 3D case when generalized to 𝑛D.
Projective and conformal geometric algebras allow to combine
the translational and rotational components into a single equation
in a similar way to how homogenous coordinates allow matrices torepresent both these components as single transformations. Gunn
[2011] uses projective geometric algebra to formulate 2D and 3D
rigid body motion in a metric-neutral way that also applies to noneuclidean spaces. In the interest of simplicity I chose not to use thisformulation.

3

BACKGROUNDThis section briefly reviews geometric algebra and its application
to rigid body dynamics in 3D as described by Doran and Lasenby[2003]. However these equations remain the same in 𝑛D.
Using geometric algebra one can write equations for rigid body
evolution in 𝑛 > 1 dimensions as:
𝑥𝑡 = 𝑣

𝑅𝑡 = − 12 𝜔𝑅

𝑣𝑡 = 𝐹 /𝑚

𝐿𝑡 = 𝜏

where 𝑥 and 𝑅 are the position (a vector) and orientation (a rotor), 𝑣
is the velocity (a vector) and 𝜔 is the angular velocity (a bivector),
𝐹 is the net force (a vector), 𝑚 is the mass (a scalar), 𝜏 is the net
torque (a bivector), and 𝐿 is the angular momentum (a bivector). The
t-subscript denotes the time-derivative. The product used for the
angular velocity bivector and the orientation rotor is the geometricproduct, described later.
Angular velocity is represented as a bivector. A bivector is formed
from the exterior product of two vectors:
𝐵 = 𝑎 ∧𝑏𝑛

and has 𝑘 = 2 coordinates, one for each pair of different orthogonalbasis vectors. These components can be stored in memory the same
way as for a vector, in a continuous array of 𝑘 numbers.
If a particle has momentum 𝑝 and position vector 𝑥 from someorigin, the angular momentum of the particle about the origin isdefined as the bivector
𝐿 = 𝑥 ∧ 𝑝.
The usual representation involves a cross product, which is only
defined in 3D. In three dimensions, vectors are “dual” to bivectors:
roughly speaking, in 3D the space orthogonal to a line (vector) isa plane (bivector). This means that in 3D one can use vectors andbivectors somewhat interchangeably. This has led to the wide useof representing angular velocity as vectors. These “axial” vectors
transform differently than regular vectors. Bivectors are a more natural representation because they are straightforwardly transformed
ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

to different coordinate systems using geometric algebra, and makethe equations work in any number of dimensions.
𝑁 -dimensional rotors provide a representation of rotations thatis a replacement for quaternions (3D) and complex numbers (2D)that works in any number of dimensions. The rotor that rotates in
the plane defined by the vectors 𝑎 and 𝑏 by twice the angle between
these two vectors is defined by the geometric product of 𝑎 and 𝑏:
𝑅 = 𝑎𝑏 = 𝑎 · 𝑏 + 𝑎 ∧ 𝑏.
It has a scalar term and a bivector term, similar to the real andimaginary terms of complex numbers and quaternions. This kind of
sum is called a multivector. The geometric product can be extended
to multivectors and hence to rotors themselves. Just like for quaternions, the product of two rotors results in a new rotor that encodes
applying the two rotations one after the other. To rotate a vector
using a rotor, the following formula is used:𝑥 ′ = 𝑅𝑥 𝑅˜
where 𝑅˜ denotes the reverse of 𝑅, which is similar to the conjugate
for quaternions and complex numbers. Products of rotors produce
multivectors that are the sum of 𝑚-vectors where 𝑚 is even and lessthan or equal to 𝑛 (this forms an algebra called the even sub-algebra).
In 3D a general rotor still has only a scalar part and a bivector part(which has 3 components), i.e. the product of two rotors still results
in a rotor that can be represented by a single plane of rotation. But in
4D an object can rotate around two independent planes of rotation
simultaneously: a general rotor has a scalar part, a bivector part(which has 6 components), and a 4-vector part (which has only one
component). I hence store a four-dimensional rotor in memory as acontinuous array of 8 numbers.
The instantaneous velocity of a point 𝑟 on the body, in the referenceframe of the body, is: 𝑣 + 𝑟 · 𝜔. The dot product used (·) is a generalization of the dot product to multivectors, also sometimes called
the left contraction (⌋).
Angular momentum is related to angular velocity by a linear
mapping∫𝑟 ∧ (𝑟 · 𝜔) d𝑉

𝐿(𝜔) =
𝑉

taking bivectors to bivectors. Doran and Lasenby [2003] take thisintegral over the three-dimensional volume of the object.
If a body is rotated by a rotor 𝑅 its angular momentum can be expressed in terms of a time-independent linear mapping I computed
once in the initial local frame of the body:˜𝐿(𝜔) = 𝑅I (𝜔) 𝑅.
(1)
In 3D vectors are dual to bivectors, so the mapping I from bivectors
to bivectors can be turned into a mapping from vectors to vectors
by taking the dual from vectors to bivectors, applying the mapping,then taking the dual again to get back to vectors. This mapping canthen be represented by a matrix: the inertia tensor.
Differentiating Equation (1) with respect to time, keeping in mind
that it is only the rotors that are time dependent gives the familiar
form of Euler’s equation:
𝐿𝑡 (𝜔) = I (𝜔𝑡 ) − 𝜔 × I (𝜔) = 𝜏(2)where × is not the regular cross product but the commutator product𝐴 × 𝐵 = 12 (𝐴𝐵 − 𝐵𝐴), used here for two bivectors 𝐴 and 𝐵. The

𝑁 -Dimensional Rigid Body Dynamics • 1:3

quantities in this equation are expressed in the reference-frame of
the body (which is rotating with the body).

4 𝑛DIMENSIONAL GENERALIZATION
Representing the inertia tensor as a mapping from vectors to vectors
only works in 3D so I do not take the dual.
To represent this bivector mapping, I first define a 𝑘 × 𝑛 matrix[𝑟 ]★ such that𝑟 ∧ 𝑎 = [𝑟 ]★𝑎(3)
for two vectors 𝑟 and 𝑎. It is a generalization of the cross productmatrix [𝑟 ]× such that 𝑟 × 𝑎 = [𝑟 ]×𝑎, but I do not take the dual:
the range of this transformation is a bivector seen as a vector of
dimension 𝑘. Also note that𝑇𝑟 · 𝜔 = [𝑟 ]★
𝜔

where, again, the domain of this matrix transformation is a bivector
seen as a vector of dimension 𝑘.
Using these definitions I define a 𝑘 × 𝑘 matrix 𝐼 that generalizes
the inertia tensor:∫
𝑇
I (𝜔) =
[𝑟 ]★ [𝑟 ]★
d𝑉 𝜔
𝑉
∫
=
Δ𝐼 d𝑉 𝜔𝑉=𝐼𝜔.

The integral is also taken over the 𝑛-dimensional volume instead of
the three-dimensional volume as would normally be done.
I give [𝑟 ]★ and Δ𝐼 in 2, 3, and 4 dimensions in appendix A. Notethat they are given with respect to lexicographically-ordered bivector basis elements (ex: 𝑒𝑥 𝑦 ,𝑒𝑥𝑧 ,𝑒 𝑦𝑧 ), but this is an arbitrary choice.
The inertia tensor’s range and domain are bivectors and transformas such. Therefore I define a 𝑘 ×𝑘 matrix [𝑅]2 such that for a bivector
˜ Then:𝐵 seen as a vector, [𝑅]2 𝐵 = 𝑅𝐵𝑅.
𝐼 ′ = [𝑅]2 𝐼 [𝑅]𝑇2

original rotor using the algorithm by Perwass [2009]. The algorithm
returns normalized vectors, guarantying that the rotor re-formed
from their product will represent a proper rotation. The gyroscopic
term in particular often generates double rotations which cannot
be corrected simply by normalizing.

5

INERTIA TENSOR OF ARBITRARY 𝑛D SIMPLICIALMESH

A 𝑛D simplicial mesh is built from multiple (𝑛 − 1)-simplices, ex: a
3D simplicial mesh is built from multiple triangles, a 4D simplicial
mesh is built from multiple tetrahedra, and so on.
Due to the linearity of the inertia tensor, to compute the inertia
tensor of an arbitrary 𝑛-mesh in any dimension it suffices to sum
the inertia tensors of 𝑛-simplices built from each (𝑛 − 1)-simplex
and the origin of the reference frame.
It is possible to compute the inertia tensor of an arbitrary 𝑛simplex directly from its vertices by evaluating the products ofinertia integrals:∫
𝑃 𝑗𝑘 =𝑗𝑘𝜌 d𝑉
𝑉

where 𝑗 and 𝑘 are coordinates axes. These integrals can be simplified using Gauss’ theorem [Dobrovolskis 1996]. However, generalizing the method of Blow and Binstock [2004] proves to be
simpler. They first compute the covariance matrix of the body, thentransform the covariance matrix into the inertia tensor. Given amatrix 𝑀 that transforms the vertices of a canonical simplex into
those of an arbitrary simplex, the covariance matrix 𝐶 ′ of the arbitrary simplex is related to that of the canonical simplex 𝐶 suchthat 𝐶 ′ = 𝑑𝑒𝑡 (𝑀)𝑀𝐶𝑀𝑇 . The covariance matrix of a canonical 𝑛simplex has two terms of interest 𝐶𝑥𝑥 and 𝐶𝑥 𝑦 . They are the integralof 𝑥 2 and 𝑥𝑦 over the hypervolume of the 𝑛-simplex aligned withthe coordinate axes.

(4)

𝐶𝑥𝑥

=

˜
or, in terms of the mapping itself: 𝑅I (𝑅𝜔𝑅)
𝑅˜ = 𝐼 ′𝜔. This lets me
express equation (2) in a global inertial reference-frame:

𝐶𝑥 𝑦

=

𝐿𝑡 (𝜔) = 𝐼 ′𝜔𝑡 − 𝜔 × 𝐼 ′𝜔 = 𝜏 .

2,(𝑛 + 2)!
1.
(𝑛 + 2)!

(5)

The formula for the time-derivative of a rotor (𝑅𝑡 = − 12 𝜔𝑅) is similar

to and hence can be substituted for the formula for quaternions thatwould be used in 3D systems. My implementation uses a (symplectic
Euler) time-stepping scheme based on the one by Guendelman et al.
[2003] but other time-stepping schemes could be used. I integrate
the gyroscopic term of eq. (5) separately using an implicit Euler
method [Catto 2015].
As the angular velocity is integrated, small errors can accumulate
in the rotor 𝑅 representing the orientation of the body. In three
dimensions, when a quaternion is used instead of a rotor, the effectis corrected by normalizing the quaternion after each time-step.
However this is insufficient in 4D and higher: while simple rotorsthat are formed from the geometric product of two unit vectors
do lie on a (𝑘 + 1)-dimensional sphere within R𝑛 , this is not thecase for rotors which correspond to double rotations (or more). To
correct for small errors I factorize the rotor with respect to thegeometric product to get a set of vectors whose product give the

To transform the covariance matrix into the inertia tensor in 𝑛D,
I take inspiration from Trenkler [2001] and expand [𝑟 ]★ in thecoordinate basis:𝑛−1
Õ[𝑟 ]★ =
𝑟𝑖 [𝑒𝑖 ]★ .
𝑖=0

The inertia tensor can then be derived from the covariance matrix𝐶:𝑇[𝑟 ]★ [𝑟 ]★

=

(

𝑛−1
Õ
𝑖=0

=

𝑛−1
Õ

𝑟𝑖 [𝑒𝑖 ]★)(

𝑛−1
Õ

𝑇
𝑟 𝑗 [𝑒 𝑗 ]★
)

𝑗=0

𝑇
𝑟𝑖 𝑟 𝑗 [𝑒𝑖 ]★ [𝑒 𝑗 ]★

𝑖,𝑗

∫
𝑉

𝑇
[𝑟 ]★ [𝑟 ]★
d𝑉

=

𝑛−1
Õ

𝑇
𝐶𝑖 𝑗 [𝑒𝑖 ]★ [𝑒 𝑗 ]★
.

𝑖,𝑗ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

1:4

•

ten Bosch, M.

It is also needed to compute the mass of the 𝑛-simplex, which can
1 𝑑𝑒𝑡 (𝑀).
be computed from its density and its volume 𝑛!
It is usual and useful to chose a body’s local coordinate-frame
such that its inertia tensor 𝐼 is diagonal. A traditional matrix diagonalization algorithm would produce a rotation matrix, but this(𝑘 × 𝑘) matrix would apply only to bivectors (like in eq. (4)) insteadof vectors, and so would not be useful as is. However, note that if the𝑛 × 𝑛 covariance matrix is diagonal, then the inertia tensor will be
as well. Therefore before computing the inertia tensor I first rotate
the body using the rotation matrix that diagonalizes the covariancematrix.

6

COLLISION RESOLUTION

Using the matrix from (3) I generalize the method of collisions,
including static and kinetic friction, of Guendelman et al. [2003],
but the scheme is general enough to be applied to other methods. I
also use their Contact Graph and Shock Propagation schemes.
After a collision I apply equal and opposite impulses 𝑗 to each
body to obtain 𝑣 0 = 𝑣 ± 𝑗/𝑚 and 𝜔 0 = 𝜔 ± 𝐼 −1 (𝑟 ∧ 𝑗) = 𝜔 ± 𝐼 −1 [𝑟 ]★ 𝑗where 𝑟 points from their respective centers of mass to the collision
location. The new velocities at the point of collision will be 𝑢 0 =𝑇 𝐼 −1 [𝑟 ] with 𝛿 the 𝑛 × 𝑛 identity
𝑢 ± 𝐾 𝑗 where 𝐾 = 𝛿/𝑚 + [𝑟 ]★
★
matrix.
Using this formulation static and kinetic friction can be expressed
the same way as in three dimensions.

7 COLLISION DETECTION
7.1 Hyperspheres with Polytope or Cylinder
To determine if a (𝑛 − 1)-sphere is colliding with an 𝑛-polytope or
𝑛-cylinder I generalize the sphere/convex-object collision detection
algorithm which is based on the Minkowski Difference (for example: [Ericson 2004]). It finds the closest point on the surface of theobject to the center of the hypersphere and checks that it is lessthan one hypersphere radius away. Note that the closest point may
be on a (𝑛 − 1)-cell, or on the cell’s boundary, which it shares with
one or more cells — the dimension of the boundary depends on thenumber of cells that share it (for example in 4D two 3D cells share aface, three cells share an edge, and so on). The projection operator𝑎 ∥ = (𝑎 · 𝐵)/𝐵 of a vector 𝑎 onto a subspace represented by the𝑛-vector 𝐵 is used to find the closest point.

7.2

Two arbitrary convex polytopes

For two arbitrary convex polytopes I generalize the separating axis
theorem collision detection method [Gottschalk 1996] to 𝑛D. A 𝑛polytope 𝑎 contains a certain number of 𝑚-cells, where 𝑚 < 𝑛. Each𝑚-cell 𝑖 is spanned by 𝑚 vectors. The exterior product of these 𝑚vectors is the 𝑚-vector 𝑉𝑎𝑚 (𝑖). Given two 𝑛-polytopes 𝑎 and 𝑏, foreach sets of cells such that 𝑚𝑎 + 𝑚𝑏 = 𝑛 − 1, form all the exterior
products:
𝑉𝑎𝑚𝑎 (𝑖) ∧ 𝑉𝑏𝑚𝑏 ( 𝑗)∀𝑖, 𝑗 .

The potential separating axes are the dual vectors to the resulting(𝑛 −1)-vectors. Note that this contains the case of the empty exterior
product where 𝑚𝑎 = 𝑛 − 1 and 𝑚𝑏 = 0. These are simply the normals
to the (𝑛 − 1)-cells (the faces in 3D) of each polytope.
ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

Fig. 2. A user rolls a duocylinder on top of a layer of hyperspheres. Some ofthem appear to be floating because they are offset in the fourth dimension
from our view.

In 3D this formulation is equivalent to taking the cross product(i.e. the dual of the exterior product of two 3D vectors) of the edge
direction vectors as potential separating axes. This corresponds toedge/edge collisions. In 4D one instead needs to take the dual of theexterior product of all the edge direction vectors of polytope 𝑎 withthe face bivectors of polytope 𝑏 and inversely. This corresponds
to the case of 1D edge / 2D face collisions. The 1D edge / 1D edge
case is subsumed by this case, similarly to how the vertex case issubsumed by the line case in 3D [Ericson 2004].
I implement the specific case of Oriented Boxes by extending
the method of Eberly [2002] to 𝑛D. All scalar triple products are
replaced by exterior products. A 4D hypercube has 8 cells, 24 faces,
32 edges and 16 vertices, but because it is symmetric it has only 4
unique vectors (normals or edge directions) and one can only form 6
unique bivectors. Therefore one needs to check 4 × 2 vectors comingfrom cell normals and 4 × 6 × 2 = 48 vectors coming from exterior
products.

8

INTERACTIVE SIMULATION AND DISPLAY

I display a three-dimensional slice of the interactive four-dimensional
objects. These 4D objects have a three-dimensional surface, represented by a 4D mesh made of simplices (in this case tetrahedra). At
any moment only the intersection between a single 3D plane andthe 4D objects is visible. This is done by slicing the 3D surface ofeach object to get a 2D surface to display, in a manner similar to themethod by Chu et al. [2009], but taking only a single slice. Slicing
each tetrahedron gives a 2D polygon (either a triangle or a quad),and all these polygons together make up the 2D surface to display.
The user can interact with the objects using either a mouse, touch
screen or VR controller. They can tap, drag and rotate objects only
within the confines of their current 3D slice, but obviously the 4D
objects are not confined to that space, and so might disappear outof view. In that case, the user can drag a slider to move the 3D slice
along the fourth axis to find them again. I decided on this display
and interaction method because it lets users play with the 4D objects
in the usual 3D physical environment they are familiar with, whilestill allowing them to find the objects in a very simple way if they
disappear.

𝑁 -Dimensional Rigid Body Dynamics • 1:5

Fig. 3. A users knocks down a stack of hypercubic dominoes. The last dominoes are initially out of view, offset in the fourth dimension.

For the purpose of this simulation, gravity trivially generalizes
to 𝑛D: it points along the normal vector to the ground, towards theground. Air friction also has the same form as in the 3D case.
On touchscreen devices such as phones an accelerometer is usually available and is used to control the direction of gravity. Initiallythe direction of gravity is restricted to the visible 3D slice, but abutton can be pressed to let the gravity vector rotate to face slightly
along the fourth dimension.
Some VR controllers support haptic feedback, and when a 4D
object held by the user hits another 4D object the program slightly
vibrates the controller holding the object.
Any position that the user can see corresponds to a position in
4D, and so the simulation can apply forces or impulses at these
positions to move the objects to follow the motion of the user. To
rotate an object it is only needed to find the 2D plane of rotationin order to form a bivector to apply as torque. In the most generalcase of rotating one 4D frame to another, I first find a rotor that
performs this transformation, then extract its bivector part. To findthis rotor I apply successive rotors to transform corresponding pairsof coordinate axes into each other.
Examples of 4D situations and user interactions can be seen in theaccompanying video and in figures 1, 2, and 3. Commercial softwareis also available [ten Bosch 2017].

9

4D DZHANIBEKOV EFFECT

In this section I briefly look at the case of the 4D Euler equation(5) in torque-free conditions. In 3D a consequence of this equationis the Dzhanibekov effect, first described by Poinsot [1851]. Givena body with three unique moments of inertia, each moment ofinertia corresponds to a specific plane of rotation. Rotations in theplanes with lowest and highest moments of inertia are stable tosmall perturbations in angular momentum, while for rotations in
the intermediate plane a small perturbation gets amplified to a larger,
periodic rotation.
In 4D the 𝜔𝑖 𝑗 component of the equation is:𝐼𝑖 𝑗 𝜔¤ 𝑖 𝑗 = (𝐼 𝑗𝑘 − 𝐼𝑖𝑘 )𝜔𝑖𝑘 𝜔 𝑗𝑘 + (𝐼 𝑗𝑚 − 𝐼𝑖𝑚 )𝜔𝑖𝑚 𝜔 𝑗𝑚
where 𝑘 and 𝑚 are the other two indices besides 𝑖 and 𝑗. Note thepresence of a sum in this equation — unlike in 3D the conditionsfor 4D rotational stability around each plane cannot be definitely
stated solely based on the body having distinct principal moments
of inertia.
The equation for plane (𝑖 𝑗) depends on all the other planes ofrotation (𝑖𝑘,𝑗𝑘,𝑖𝑚,𝑗𝑚) except for the plane perpendicular to it (𝑘𝑚).
This is expected, as 4D objects can rotate in two perpendicularplanes independently.
The equations reduce to the 3D case if the initial angular velocities
lie within a 3D subspace (i.e. if 𝜔𝑖𝑐 = 𝜔𝑐𝑖 = 0 for any specific 𝑐):figure 5c shows the usual unstable rotation when rotating aroundan intermediate plane. The small decrease in angular momentum
over long time scales is due to the implicit integration method.
Double rotations appear to be stable states: figure 5a shows how
a rotation around the (𝑥𝑦) plane and small perturbations aroundthe other planes introduces the dual rotation (𝑧𝑤). The moments ofinertia seem to have little effect on the behavior, a similar graph is
generated for cuboids of size ( 12 , 1, 23 , 2), ( 21 , 12 , 3, 21 ) and others. Even
though rotation around the dual plane is independent (perturbing
only the (𝑧𝑤) plane does not affect the (𝑥𝑦) plane), the equations
are coupled through the other planes, which is what creates theeffect.
Figure 5b shows the effect removing perturbations from the initialstate of figure 5a. Additional oscillations are created but the systemstill ends up in a stable state.
Figure 5d and 4 show the effect of adding an additional perturbation to the initial state of figure 5c. The effect is similar to theprevious case.

1

ωxy

3

0.8

2

0.6

1

0.4

0

0.2

ωzw
5

10

ωxw
ωyw ωxz ωyz

ωzw

−1

0
0

ωxy

ωxz ωxw−20215(a) 𝜔𝑥 𝑦 = 1, 𝜔𝑖 𝑗 = 10−4ωyz ωyw
46810(b) 𝜔𝑥 𝑦 = 3.2, 𝜔𝑥𝑧 = 10−5 ,𝜔𝑥 𝑤 = 10−3

ωyz1

ωxw ωyz

1

ωzw
0

ωxy

0

ωxz

−1

ωxw

ωywωzw

ωyw

ωxy

ωxz

6

8

−1
0

5

10

15

20

25

30

(c) 𝜔𝑥 𝑦 = −0.008, 𝜔 𝑦𝑧 = 1.63

0

2

4

10

(d) 𝜔𝑥 𝑦 = −0.008, 𝜔 𝑦𝑧 = 1.63,
𝜔 𝑦𝑤 = 0.0018

Fig. 5. Angular velocity over time (in seconds) for a 4D cuboid of size
(1.1, 0.3, 2.5, 0.7) in various initial conditions.
ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

1:6

•

ten Bosch, M.

Fig. 4. A hypercuboid rotating under gyroscopic forces.

10

FUTURE WORKThe formulation of rigid body dynamics could be extended to noneuclidean spaces, such as hyperbolic space. Other methods for collision detection, soft bodies, as well as other forces and constraints
could be generalized to higher dimensions. Concerning the interaction method, it would be interesting to allow users to manipulate
4D objects outside of their 3D visible slice.

A

BIVECTOR OPERATOR AND INERTIA TENSOR
MATRICES

2D case:[𝑟 ]★ = −𝑦
Δ𝐼 = 𝑥 2 + 𝑦𝑥2

3D case:−𝑦©
[𝑟 ]★ = ­−𝑧«0

4D case:

0
ª
𝑥®
𝑦¬

𝑥
0
−𝑧

𝑥 2 + 𝑦2©
Δ𝐼 = ­ 𝑦𝑧« −𝑥𝑧

𝑦𝑧
𝑥 2 + 𝑧2
𝑥𝑦

−𝑦
©
­ −𝑧
­
­−𝑤
[𝑟 ]★ = ­
­ 0
­
­ 0
« 0

𝑥
0
0
−𝑧
−𝑤
0

𝑥 2 + 𝑦2
©
­ 𝑦𝑧
­
­ 𝑦𝑤
Δ𝐼 = ­
­ −𝑥𝑧
­
­ −𝑥𝑤
« 0

𝑦𝑧
𝑥 2 + 𝑧2
𝑧𝑤
𝑥𝑦
0
−𝑤𝑥

𝑦𝑤
𝑧𝑤
𝑤2 + 𝑥2
0
𝑥𝑦
𝑥𝑧

−𝑥𝑧
ª
𝑥𝑦 ®
2
2
𝑦 +𝑧 ¬

0
𝑥
0
𝑦
0
−𝑤

0
ª
0®
®
𝑥®
®
0®
®
𝑦®
𝑧¬

−𝑥𝑧
𝑥𝑦
0
𝑦2 + 𝑧2
𝑧𝑤
−𝑦𝑤

−𝑥𝑤
0
𝑥𝑦
𝑧𝑤
𝑤 2 + 𝑦2
𝑦𝑧

0
ª
−𝑥𝑤 ®
®
𝑥𝑧 ®
®
−𝑦𝑤 ®
®
𝑦𝑧 ®
2
2
𝑤 +𝑧 ¬REFERENCES
E.A. Abbott. 1884. Flatland: A Romance of Many Dimensions. Seeley & Co.
T. Banchoff. 1990. Beyond the Third Dimension: Geometry, Computer Graphics, andHigher Dimensions. Scientific American Library. https://books.google.com/books?
id=8n55QgAACAAJJan Bender, Kenny Erleben, and Jeff Trinkle. 2014. Interactive simulation of rigid body
dynamics in computer graphics. Computer Graphics Forum (Print) 33, 1 (2014),
246–270. https://doi.org/10.1111/cgf.12272
Jonathan Blow and Atman J Binstock. 2004. How to find the inertia tensor (or othermass properties) of a 3D solid body represented by a triangle mesh. (2004). http://number-none.com/blow/inertia/
ACM Trans. Graph., Vol. 39, No. 4, Article 1. Publication date: July 2020.

Stephen Cameron. 1990. Collision detection by four-dimensional intersection testing.
Robotics and Automation, IEEE Transactions on 6 (07 1990), 291 – 302. https://doi.
org/10.1109/70.56661
Erin Catto. 2015. Physics for Game Programmers : Numerical Methods. (2015). https://www.gdcvault.com/play/1022197/Physics-for-Game-Programmers-Numerical Presentation at Game Developers Conference 2015.
Arthur Cayley. 1846. Sur quelques propriétés des déterminants gauches. Journal für
die reine und angewandte Mathematik 1846, 32 (1846), 119–123.
Alan Chu, Chi-Wing Fu, Andrew J. Hanson, and Pheng-Ann Heng. 2009. GL4D: A
GPU-based Architecture for Interactive 4D Visualization. IEEE Transactions on
Visualization and Computer Graphics (IEEE Visualization 2009) 15, 6 (Nov.-Dec. 2009),
1587–1594.
Anthony R Dobrovolskis. 1996. Inertia of any polyhedron. Icarus 124, 2 (1996), 698–704.
Chris Doran and Anthony Lasenby. 2003. Geometric algebra for physicists. Cambridge
University Press.
L. Dorst, D. Fontijne, and S. Mann. 2009. Geometric Algebra for Computer Science: An
Object-oriented Approach to Geometry. Elsevier. https://books.google.com/books?
id=yaEqlAEACAAJDavid Eberly. 2002. Dynamic collision detection using oriented bounding boxes. Geometric Tools, Inc (2002).
Christer Ericson. 2004. Real-Time Collision Detection. CRC Press, Inc., USA.
Stefan Gottschalk. 1996. Separating axis theorem, Technical Report TR96-024, Department of Computer Science, UNC Chapel Hill. (1996).
Eran Guendelman, Robert Bridson, and Ronald Fedkiw. 2003. Nonconvex rigid bodies
with stacking. In ACM SIGGRAPH 2003 Papers (San Diego, California) (SIGGRAPH’03). ACM, ACM, New York, NY, USA, 871–878. https://doi.org/10.1145/1201775.
882358Charles Gunn. 2011. Geometry, Kinematics, and Rigid Body Mechanics in Cayley-Klein
Geometries. Doctoral Thesis. Technische Universität Berlin, Fakultät II - Mathematik
und Naturwissenschaften, Berlin. https://doi.org/10.14279/depositonce-3058
Charles G. Gunn and Steven De Keninck. 2019. Geometric Algebra and Computer
Graphics. In ACM SIGGRAPH 2019 Courses (Los Angeles, California) (SIGGRAPH ’19).
Association for Computing Machinery, New York, NY, USA, Article 12, 140 pages.
https://doi.org/10.1145/3305366.3328099
D. Hilbert and S. Cohn-Vossen. 1952. Anschauliche Geometrie. Chelsea Publishing
Company. https://books.google.com/books?id=d6sBd9h1HbMC
A. Macdonald. 2011. Linear and Geometric Algebra. Alan Macdonald. https://books.
google.com/books?id=oxhJYgEACAAJ
Christian Perwass. 2009. Geometric Algebra with Applications in Engineering (1st ed.).
Springer Publishing Company, Incorporated.
Louis Poinsot. 1851. Théorie nouvelle de la rotation des corps. Bachelier.
Andrew J Sinclair and John E Hurtado. 2005. Cayley kinematics and the Cayley formof dynamic equations. Proceedings of the Royal Society A: Mathematical, Physicaland Engineering Sciences 461, 2055 (2005), 761–781.
Marc ten Bosch. 2017. 4D Toys. https://4dtoys.com
Götz Trenkler. 2001. The vector cross product from an algebraic point of view. Discuss.
Math. - General Algebra and Applications 21, 1 (2001), 67–82.
Xiaoqi Yan, Chi-Wing Fu, and Andrew J. Hanson. 2012. Multitouching the Fourth
Dimension. Computer 45, 9 (2012), 80–88. https://doi.org/10.1109/MC.2012.77
Hui Zhang and Andrew J. Hanson. 2006. Physically Interacting with Four Dimensions.
Springer Berlin Heidelberg, Berlin, Heidelberg, 232–242. https://doi.org/10.1007/
11919476_24