Suppose we have a rigid body extended in R3 and want to characterize its dynamic properties. Where xj are the 3 coördinates for 0 ≦ j < 3:
m = ∫ dm is the total mass,
mj = ∫ xjdm are the 3 first moments of the mass distribution.
mjk = ∫ xjxkdm are the 6 independent second moments of the mass distribution.

See note on summation convention. The coördinates of the center of mass (CM) are at mj/m.
These 10 numbers characterize the body if it is rigid. In the absence of external forces, these 10 together with initial position (3 numbers), velocity (3 numbers) and either angular velocity or angular momentum (3 numbers each), are enough to predict the future position, orientation, velocity and angular velocities. The velocity of the CM and angular momentum are constant, of course. The orientation can be given by an orthogonal matrix but really needs only 3 independent numbers. The development below assumes that the point of the body at the origin is constrained not to move. This matches the tumbling problem above if mj = 0, a situation that may be brought about by translating the coördinate system.

In 3 dimensions there are 3 independent values in an antisymetric tensor and there are 3 independent values in a vector. Axial vectors are thus able to code the information that requires an antisymetric tensor in higher dimensions. The tensor of inertia, Iij = δijmkk – mij, derived from the 2nd moments, fits well with these axial vectors for angular velocity, ω, and momentum, L, to make the convenient and simple formula “L = Iω”. The tensor of inertia looses this magic beyond 3D however but the 2nd moments remain in control.

This code corroborates that L = Iω is correct in 3D when I is viewed as a symmetric matrix and L and ω are both (axial) vectors. Here is a generalized formula in n dimensions.

This program corroborates the conventional 3D relation L = Iω.

This differs from this only in style and notation.

The equations so far look very simple but the tensor of inertia must presumably be expressed in the laboratory frame and thus varies with time. If Aij is the orthogonal matrix specifying the current orientation of the rigid body then Jij = AikAjlIkl is the inertia tensor in lab coordinates which varies with time.

The time rate of change of Aij depends on the angular velocity ωik thus: d/dt Aij = Aikωjk

A calculation cycle can proceed thus: We know L and I which are constant.

I should write the code!