/// 注：本文虽然大篇幅讨论的虚构四维世界中的虚构的物理规律，但也对纯数学的单位四维2-向量的均匀随机分布规律进行了讨论。
内容概要
- 四维贾尼别科夫效应
- 有能量交换的多体旋转系统
- 角动量的随机分布律
- 四维刚体的的主轴定理
贾尼别科夫效应
一个前苏联的宇航员在太空中把弄一个T字形零件，他偶然发现让零件漂浮着快速旋转起来，零件的旋转方向却突然时不时在空中翻转，后来就以他的名字(Джанибеков)命名了这个贾尼别科夫效应(Dzhanibekov Effect)。
解释这个现象的方法很多，它的根源来源于角动量守恒与那个把手的转动惯量的各向异性。我们先从最简单的牛顿第一定律开始，它说物体不受外力作用一定会保持匀速直线运动或静止不动。然而这没有考虑物体的旋转。当物体不受外力时是不是应该保持匀速角速度旋转呢？并不是。其实，物体不受外力时并不是速度与角速度不变，而是动量与角动量不变，由于动量等于质量乘速度，一般物体的质量不会发生变化，所以能推出物体的速度也恰好不变，但物体的角动量与角速度之间的关系就复杂了。就拿一个网球拍来说，由于它各个方向的尺寸不同，绕不同的轴旋转的难易程度也不同，即每个方向的转动惯量不同。角动量等于转动惯量乘以角速度，当物体旋转起来后，这些轴也在跟着旋转，每个方向的转动惯量也在变，要想保持输出的总角动量不变，角速度就必须跟着转轴变才能抵消掉转动惯量的变化，这就是零件不受外力旋转时自己会翻转的原因。如果旋转物体是一个均匀对称的图形，比如正方体球体这种，它各个方向的转动惯量都相同，就只会老老实实地匀速旋转，不再有贾尼别科夫效应。
其实这个效应的出现条件还挺苛刻的，它还要求物体的初始旋转轴是介于长轴和短轴之间的中间轴，这是什么意思呢？由三维刚体力学的分析可知，任何刚体的转动惯量都等价于一个三轴椭球或长方体的惯量，只要描述了三个主轴的转动惯量就能确定整个物体的转动行为。因此我们可以把网球拍等效为一块长方体的板子，上图分别标出了网球拍绕长轴$a$，短轴$c$与中间轴$b$的转动，只有绕中间轴（上图中标红的轴）的转动不稳定才产生贾尼别科夫效应，这个关于旋转稳定性的定理亦称网球拍定理。关于为什么只有中间轴旋转不稳定的，有基于力学的解释，也有基于能量的解释，还有基于欧拉旋转微分方程稳定性分析的解释等等，这不是本文的重点，就不再赘述了。
四维情况
虽然我们无法实际在想象的四维空间中做实验来目睹四维效应，但可以通过Tesserxel中的物理引擎进行模拟。具体引擎中使用的物理公式我放在了附录中。首先要想旋转具有不稳定性，就必须要找四个“轴”都不等的物体，最简单就是超长方体了。不同于三维，四维空间的旋转发生在平面上，主轴之间的两两组合有六个方向，初始旋转位置也有六种可能性。
设超长方体有4条不等长的边$a$、$b$、$c$、$d$，且$a < b < c < d$。我通过数值模拟计算发现了以下规律：
- 初始旋转平面选在$ab$、$bc$或$cd$上都是稳定的。$ab$类似三维的最长轴，它的转动惯量最小；$cd$类似三维的最短轴，它的转动惯量最大；虽然旋转平面$bc$是一种“中间轴”，但却是稳定的。点击这里在Tesserxel中在线模拟$bc$平面的旋转。
- 初始旋转平面选在$ac$或$bd$时会出现类似三维的中间轴不稳定的情况，且产生几乎跟三维情况相同的贾尼别科夫效应。点击这里在Tesserxel中在线模拟$ac$平面的旋转。
- 初始旋转平面选在$ad$上时旋转也不稳定，且会有两种周期的反向旋转交替进行，这是四维贾尼别科夫效应独有的一种现象。由于每个周期都会让旋转方向翻转，两个不同的周期共同作用导致刚体局部坐标系中的角动量变化曲线看起来像是一种脉宽调制(PWM)波形，点击这里在Tesserxel中在线模拟$ad$平面的旋转。
以上的在线模拟中可以按键盘L在惯性坐标系与刚体的局部坐标系之间切换，刚体的轴从小到大对应的坐标分别是$x$、$y$、$z$、$w$。可以看出:
1.在惯性坐标系中：任何情况下，物体的角动量（字母J表示）的各个分量都不变，是一条直线，代表着角动量绝对守恒，而在有贾尼别科夫效应的场景中，可以看到除了初始角速度对应的分量不怎么变以外，在旋转方向发生颠倒时会临时出现一些较小的其它角速度分量（字母W表示）。
2. 将坐标系切换到局部坐标系后可以看到在每次物体反转时初始角速度对应的分量的方向也会跟着翻转。
这两种坐标系下角速度的行为为什么会截然不同呢？据说当时宇航员贾尼别科夫在太空偶然发现了这个效应后，他非常担心地球自转着某天也会突然翻转过来。当然后面我们知道地球的自转轴是三轴椭球最短的方向，旋转是稳定的，不用担心发生这种事。不妨假设地球的自转轴就是中间轴，本来地球是自西向东转的，太阳东升西落，当地球南北极颠倒后，由于总角动量的方向不会变，从宇宙的惯性系中看地球还是在自西向东转动，然而在地球上的人的感觉就不同了：原来的东西南北方已经由于地球的颠倒变成了实际的西东北南方，地面的人们会觉得太阳西升东落，但太空中的人们还是认为太阳东升西落，想必现在读者应该能理解惯性系和局部系之间角速度一个几乎不变、另一个周期性反向了吧。
稳定性规律总结
观察上面的三种情况不难发现稳定性的判断有以下规律：若平面是由按大小顺序排列相邻的边张成的则是稳定旋转，若中间相隔一个边则产生单周期的贾尼别科夫效应，若中间相隔两个边则产生双周期的贾尼别科夫效应，此外我还发现了一些关于等角旋转的规律：
- 若初始的旋转是等角双旋转则旋转是稳定的，比如同时按相同角速度旋转$ac$与$bd$平面：点击这里在Tesserxel中在线模拟。
- 若是非等角的双旋转，则看两个绝对垂直的旋转的角速度哪个大哪个的效应就更明显，且初始旋转越接近完美的等角双旋转，翻转周期就越长。下面这个在线模拟中，双旋转的两绝对垂直的旋转速度比例是$15/16$，已经相当接近等角了，因此方向翻转的周期特别长，需要等半分多钟才能看到角动量的翻转：点击这里在Tesserxel中在线模拟。
小插曲
我第一次听说贾尼别科夫效应其实是Marc ten Bosch在他的论文《𝑁-Dimensional Rigid Body Dynamics》中提到的四维空间中的贾尼别科夫效应，然后我再去搜三维情形才知道苏联宇航员的那个故事，然后我才发现我之前些的四维引擎根本没考虑不均匀的转动惯量引起的刚体进动项，我赶忙补上了引擎中的这个大漏洞。
Marc ten Bosch的那篇论文里的结论跟我这篇文章的结论并不一样，他的模拟结果表明，四个轴都不相等的超长方体最终都会趋于双旋转这一种稳定状态，然而我认为该结果有误：试想一个比较接近超立方体的超长方体，它的各方向转动惯量虽有不同但差异不大。若初始旋转是一个简单旋转，对应角速度是简单2-向量，虽然总角动量不一定是简单2-向量，但应该不会偏离太多；如果演化到了等角双旋转，则角速度就是自对偶或反自对偶2-向量，跟简单2-向量的差异太大，不可能满足总角动量守恒。注意到Marc自己也说他模拟中出现了隐式求解器数值误差积累导致角速度明显减小（Tesserxel使用的是精度很高的四阶龙格库塔法），由于我不知道他是怎样算的，只能猜测很可能是同样原因导致了错误收敛到等角旋转的结果。我在数学软件Mathematica中和自己用typescript写的tesserxel物理引擎中求解都表明四维自由旋转的刚体不会自主趋于稳定的等角双旋转状态。
有能量交换的多体旋转系统
刚才我们看到的是孤零零的刚体的自由旋转。如果刚体内部有损耗，比如一个盛着液体的杯子，维基百科上说它将在保持总角动量不变的情况下逐渐趋于能量最低的稳定状态。那么一个有趣的问题来了: 一般四维星球的自转到底最可能是怎样的呢？是接近简单旋转可能性大，还是接近等角双旋转或一般双旋转的可能性大呢？Higher space论坛中有人认为，等角双旋转的星球上各点速度大小一样，是能量最均分的稳定状态，简单旋转中的赤道速度最高，极道速度为零，并不是一种平衡态，因此最终所有星球都会演化至等角双旋转。然而跟Marc ten Bosch的疏忽一样，对于孤立的星球在总角动量的约束下是没办法从简单旋转直接演化到等角旋转的状态的，唯一的机会只能在星球形成的星系的前期演化过程中。由于四维天体力学存在很多设定上的硬伤Bug，星系的轨道系统在三次方反比衰减的引力场中不稳定，因此这里我们暂不讨论具体星球如何获得初始角动量，而是来看一个热力学“玩具模型”：想象有一个绝对光滑的超球形的腔体，里面充满了一堆初速度随机分布的小球（其实是超球），它们在里面像空气分子做布朗运动那样自由碰撞，小球表面有摩擦力使得它们可以交换动量与角动量，但没法从超球壁获得额外角动量，因此可以保持总角动量守恒。经过一段时间达到平衡态时这些小球的角速度的分布规律是怎样的呢？我在Tesserxel中放置了40个小球，并标出了每个小球的角速度的等角程度——自对偶为1，简单为0，反自对偶为-1，具体计算公式见这里。点击这里在Tesserxel中在线模拟
由于它们在超球壳内碰撞且超球壳无法提供角动量，因此总角动量会保持守恒，但不同于单个孤立的刚体，整个系统的角动量由两部分组成，一个是轨道角动量，一个是自转角动量。这两部分各自是不守恒的，只有加起来的总角动量才守恒。给个简单的例子来加深理解: 想象真空中有两个质量几乎相同的宇航员，他们之间距离很近，保持相对静止。这时他们伸出手，以下图绿色箭头标出的方式击掌，相互作用力让他们开始旋转且互相远离。不难发现，他们都朝着同一个方向旋转，看似违背了总角动量守恒，其实除了自转有角动量，他们质心的运动也有轨道角动量，且方向是顺时针与自转的逆时针方向相反，所以总角动量还是守恒。
通过Tesserxel的模拟发现，这些小球确实有时会共同趋于等角双旋转，若开始自转角动量自对偶部分大些，则最终全都会趋于自对偶的角动量，若初始旋转反自对偶分量大些，则最终全都趋于反自对偶角动量。
若初始角动量是简单的，则最终的稳态自对偶与反自对偶都有可能，有时则是完全均匀随机分布，看不出往哪边收敛，不排除模拟时间短，后续还是会收敛到对偶状态。我还偶然发现了一次同时向着两端收敛的情况，然而后面几乎没有再复现该情况。
虽然我不敢保证上面的规律不是计算误差导致的假像，但这些现象至少跟总角动量守恒并不矛盾，因为自转多出来的等角分量可以被轨道角动量多出方向相反的分量来抵消掉，因此从这个意义上说星球演化趋于等角双旋转是有道理的。然而四维天体轨道的稳定性不允许我们等足够的时间让星球充分碰撞交换角动量，因此实际星球很难说都是等角双旋转的。既然无法想象合理的星系演化过程来推导星球的角动量，如果直接假设星球有随机方向的角动量，简单与复合2-向量怎样分布才算真正的随机呢？下面我们假设角动量的大小为1，探究如何均匀地随机生成任意方向的角速度2-向量。
角动量的随机分布律
如何随机生成一个任意的单位2-向量呢？下面给出我想到的几种方法:
- 这里已经介绍了基于左右等角分解的生成随机的简单单位2-向量的算法，可它不能生成复合的2-向量。有理由相信，将这些足够多的随机的2-向量叠加起来的方向就能逼近得到任意的均匀分布的随机2-向量。其实这相当于一个三维空间中的随机行走模型：我们每次叠加一个随机的简单2-向量相当于在三维的自对偶和反自对偶空间中都各进行了一步相同距离的随机行走，复合得到的新2-向量的对偶性取决于两个空间中随机行走的距离之比。
- 这里已经介绍了基于对偶分解的可生成四维空间中的随机旋转的算法，对随机旋转取对数得到生成它的类似角动量的生成元2-向量，将它们单位化即可得到方向随机分布的2-向量。
- 还有种直接粗暴的做法：2-向量不过是有六个自由度的向量，我们要生成单位2-向量其实就是在六维空间中的五维超球面上生成均匀分布的点。
以上三种方法生成的2-向量的分布规律是不是一样的呢？显然这些方法在四维空间中都是各向同性的，因此差异就体现在2-向量的对偶程度的分布上，即自对偶分量与反自对偶分量的分布比例。我专门定义了一种对偶程度来定量描述它：自对偶的对偶程度为1，反自对偶为-1，简单则为0。
选读：对偶程度函数的定义
由于自对偶与反自对偶2-向量的另一种分量为0，直接定义比例会带来无穷大的问题，因此我将比值通过反正切函数转换成了角度，并进行归一化。设$A$是一个单位2-向量，则它的对偶程度定义为
$$\mathrm{duality}(A)={4 \over\pi}(\arctan(|A+A^*|,|A-A^*|)-{\pi\over 4})$$
它的几何意义是把单位圆周上的点的$x$、$y$坐标之比对应自对偶分量与反自对偶分量的比例，用弧长参数来衡量离哪个状态更接近。
下面是我模拟的结果：三种生成算法分别生成了两亿个2-向量，绘制出了三种算法得到的2-向量的对偶程度分布的直方图，其中红色对应生成100个简单2-向量叠加起来的“随机行走”算法，绿色对应生成随机旋转取对数的算法，蓝色对应五维超球面上取随机点的算法。可以看到随机旋转取对数生成的分布要向代表着简单2-向量的中部更集中些，而另外两种算法得到的分布几乎是一样的，我们可以确定“随机行走”算法与五维超球算法得到的分布才是真正的均匀分布。
为什么随机旋转取对数得到的结果更偏向于简单2-向量呢？这是因为简单2-向量只能生成简单旋转，而复合2-向量有时也能生成简单旋转造成的。考虑2-向量${11\pi/ 6}e_{xy}+2\pi e_{zw}$生成的旋转：由于$2\pi e_{zw}$旋转了整整一周，相当于没旋转，因此整个旋转等效为沿$xy$方向旋转-30°的简单旋转，而这个2-向量其实非常接近于自对偶，取对数时只会给出最简单直接的旋转角度小于180度的生成元，不会恢复出原先的复合旋转，这样就导致了简单2-向量的比重被放大了。
四维刚体的的主轴定理
不知你有没有这样的疑惑，既然四维旋转的自由度是6，这意味着转动惯量有六个参数来描述（转动惯量是对称矩阵，一定可对角化因此需要6个参数），但刚才讨论的超长方体的尺寸只有四个，少了两个自由度，我们会不会遗漏了一些旋转稳定性分析的情况呢？比如会不会存在着一些比超长方体都还不对称的图形，它们有着更复杂的旋转规律呢？答案是没有。转动惯量的六个参数是不独立的，独立的自由度就是4，它们均可等效为超长方体的转动惯量。如果读者不关心具体细节，本文就到这里结束了。下面附上四维刚体转动惯量的推导以及主轴定理的证明。
选读1：四维转动惯量张量的推导
刚体动量公式
我们先来看刚体动量公式的推导作为热身。质点的动量定义为$p=mv$，其中$m$为质点的质量，$v$为质点的速度。刚体是许多质点的聚合体，于是刚体的总动量是各质点动量之和：
$$p=\sum_i m_iv_i$$定义总质量$m$为$m=\sum_i m_i$，则刚体的总动量相当于按各质点的质量比例对速度求加权平均的速度再乘总质量，即$$p=m\sum_i {m_i\over m} v_i$$令刚体的平均速度$v$为$$v=\sum_i {m_i\over m} v_i$$于是我们得到刚体的动量公式为$p=mv$。简直跟质点的动量公式一模一样，更进一步，根据刚体重心的性质我们还可以知道这里的平均速度$v$正好就是重心的运动速度。这个性质在任何维度下都成立。
刚体角动量公式
下面来看角动量。三维空间中质点的角动量向量定义为$$J=r\times p=mr\times v$$其中$r$是质点的位置，$p=mv$是质点的动量。但是我们知道，四维空间没有叉乘运算，角动量的方向不能再说是哪根轴，而应该说在哪个平面中，因此需要把定义中的叉乘换成楔积：$$J=r\wedge p = m r \wedge v$$刚体的总角动量$J$是各质点角动量之和：$$J=\sum_i J_i = \sum_i m_i r_i \wedge v_i$$
刚体有个特点，那就是为了保持刚体不形变，它各部分的速度之间是有强制约束关系的，具体体现在可以只用质心$r$处的速度$v$和绕质心转动的角速度$\omega$两个全局物理量来描述所有质点的速度分布，即：$$v_i=v+(r_i-r).\omega$$注意这里位置乘上角速度得到线速度中的乘法运算不再是叉乘而是一种新的“点乘”，点乘法法则如下：
$$e_{y}.e_{xy}=-e_{x}$$$$e_{x}.e_{xy}=e_{y}$$$$e_{z}.e_{xy}=0$$其中$x$、$y$、$z$可以是任意坐标轴对应的字母。它其实是几何代数中推出的一种运算，读者只要清楚角速度$e_{xy}$方向代表着$xy$平面上发生的从$x$正半轴向$y$正半轴的旋转方向就可以很容易验证这与原先定义的角速度公式是一样的。
我们将新的速度关系式代入角动量表达式中，得到
$$J=\sum_i m_i r_i \wedge (v+(r_i-r).\omega)$$
我们将两项拆开，第一项就是轨道角动量$L$，第二项就是自转角动量$S$：
$$J=L+S=(\sum_i m_i r_i \wedge v)+(\sum_i m_i r_i \wedge( (r_i-r).\omega))$$
注意质心$r$的定义恰好是质点位置对质量的加权平均，因此做变形：$$L=m\sum_i {m_i\over m} r_i \wedge v$$容易发现最终可化简为$L=m r \wedge v$，即轨道角动量计算公式与单个质点的公式相同。
第二项自转角动量$S$要麻烦一些：$$S = \sum_i m_i r_i \wedge ((r_i-r).\omega)$$我们给前面一项配一个项，把所有位置坐标都转换到相对质心来:$$S = \sum_i m_i (r_i-r) \wedge ((r_i-r).\omega) + \sum_i m_i r \wedge ((r_i-r).\omega)$$后面补的那项其实为零，因为有：$$\sum_i m_i r \wedge ((r_i-r).\omega)=mr\wedge((r-r).\omega)=0$$我们可以看出，质心的位置只影响轨道角动量而不影响自转角动量。因此不妨可假设质心位于原点，前面这个表达式可写成$$S=\sum_i m_i r_i \wedge (r_i.\omega)$$到这一步无法再继续化简了。我们设单个质点的坐标为$r_i=(r_{xi},r_{yi},r_{zi},r_{wi})$，把楔积和点积都按坐标展开通过冗长的计算可得上式等价于下面的一个矩阵表达式：
$$S = \mathrm{I}\omega$$其中$$\mathrm{I}=\sum_i m_i\begin{pmatrix}r_{xi}^2+r_{yi}^2&r_{yi}r_{zi}&r_{yi}r_{wi}&-r_{xi}r_{zi}&-r_{xi}r_{wi}&0\\r_{yi}r_{zi}&r_{xi}^2+r_{zi}^2&r_{zi}r_{wi}&r_{xi}r_{yi}&0&-r_{xi}r_{wi}\\r_{yi}r_{wi}&r_{zi}r_{wi}&r_{xi}^2+r_{wi}^2&0&r_{xi}r_{yi}&r_{xi}r_{zi}\\-r_{xi}r_{zi}&r_{xi}r_{yi}&0&r_{yi}^2+r_{zi}^2&r_{zi}r_{wi}&-r_{yi}r_{wi}\\-r_{xi}r_{wi}&0&r_{xi}r_{yi}&r_{zi}r_{wi}&r_{yi}^2+r_{wi}^2&r_{yi}r_{zi}\\0&-r_{xi}r_{wi}&r_{xi}r_{zi}&-r_{yi}r_{wi}&r_{yi}r_{zi}&r_{zi}^2+r_{wi}^2\end{pmatrix}$$我们把2-向量$S$与$\omega$看成是以{$e_{xy},e_{xz},e_{xw},e_{yz},e_{yw},e_{zw}$}为基的六维向量。那一大坨矩阵$\mathrm{I}$就是所谓的四维转动惯量矩阵。转动惯量是种张量：它本质是把角速度映射到角动量的线性映射，张量与矩阵的关系类似向量与坐标的关系：同样的张量选择了不同的坐标系后会有不同的矩阵形式。
最后我们整理一下刚体角动量的公式，这个公式在任何维度下都成立，只不过惯量矩阵越来越庞大，其边长随维度平方级增长，整体元素随维度四次方增长，我们把2-向量看作二阶张量， 其实是个四阶张量。
$$J=L+S=(mr\wedge v)+\mathrm{I}\omega$$
选读2：主轴定理的证明
主轴与矩阵对角化
三维世界的长方体的转动惯量是一个对角矩阵，而一般图形的惯量矩阵由于都是实对称的，根据线性代数一定能够找到一个旋转变换$R$将该矩阵对角化，即$\mathrm{I}’=\mathrm{RIR}^{-1}$，且旋转矩阵$\mathrm{R}$由$\mathrm{I}$的特征向量组成。
同理容易知道四维的超长方体的转动惯量也是一个对角矩阵，我们可以用同样的方法将其对角化：$$\mathrm{I}’=\mathrm{RIR}^{-1}=\mathrm{diag}(I_{xy},I_{xz},I_{xw},I_{yz},I_{yw},I_{zw})$$其中$\mathrm{R}$为$\mathrm{I}$的特征向量组成的六维正交矩阵，$I_{ij}$为对应的特征值。
但注意现在矩阵是六维的，就算得到了变换矩阵$R$也是在六维的2-向量空间中的正交变换，但四维旋转的自由度是6，六维旋转的自由度却是15，因此求出了六维空间中的正交变换不一定能对应到合法的四维空间旋转，问题似乎在这里卡住了。
等角坐标下的惯量矩阵
从转动惯量的矩阵表达式中我们能够发现一些规律，比如矩阵是对称的，且副对角线元素始终为0（由于绝对垂直的旋转不会相互影响造成的）。我们还知道左右等角旋转是互不干扰的，这启发我们如果转换到“等角坐标系”中转动惯量矩阵形式说不定能够解耦成更小的分块矩阵处理，从而得到转动惯量的更多性质，解耦还有个好处就是可以降低找主轴对角化矩阵的维数，消除冗余自由度，或许能避开刚才的困境。
我们将原来的2-向量的六个坐标面组成的基$$e_{xy},e_{xz},e_{xw},e_{yz},e_{yw},e_{zw}$$变换成三个自对偶与三个反自对偶单位2-向量组成的基：$${\sqrt 2\over 2}(e_{xy}+e_{zw}),{\sqrt 2\over 2}(e_{xz}-e_{yw}),{\sqrt 2\over 2}(e_{xw}+e_{yz}),{\sqrt 2\over 2}(e_{xy}-e_{zw}),{\sqrt 2\over 2}(e_{xz}+e_{yw}),{\sqrt 2\over 2}(e_{xw}-e_{yz})$$则在新的坐标系下，转动惯量矩阵$\mathrm{I}$变成了
$$\mathrm{I}’=\sum_i m_i\begin{pmatrix}
\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)&0&0&\frac{1}{2}(r_{xi}^2+r_{yi}^2-r_{zi}^2-r_{wi}^2)&r_{yi}r_{zi}-r_{xi}r_{wi}&r_{xi}r_{zi}+r_{yi}r_{wi}\\0&\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)&0&r_{xi}r_{wi}+r_{yi}r_{zi}&\frac{1}{2}(r_{xi}^2-r_{yi}^2+r_{zi}^2-r_{wi}^2)&r_{zi}r_{wi}-r_{xi}r_{yi}\\0&0&\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)&r_{yi}r_{wi}-r_{xi}r_{zi}&r_{xi}r_{yi}+r_{zi}r_{wi}&\frac{1}{2}(r_{xi}^2-r_{yi}^2-r_{zi}^2+r_{wi}^2)\\\frac{1}{2}(r_{xi}^2+r_{yi}^2-r_{zi}^2-r_{wi}^2)&r_{xi}r_{wi}+r_{yi}r_{zi}&r_{yi}r_{wi}-r_{xi}r_{zi}&\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)&0&0\\r_{yi}r_{zi}-r_{xi}r_{wi}&\frac{1}{2}(r_{xi}^2-r_{yi}^2+r_{zi}^2-r_{wi}^2)&r_{xi}r_{yi}+r_{zi}r_{wi}&0&\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)&0\\r_{xi}r_{zi}+r_{yi}r_{wi}&r_{zi}r_{wi}-r_{xi}r_{yi}&\frac{1}{2}(r_{xi}^2-r_{yi}^2-r_{zi}^2+r_{wi}^2)&0&0&\frac{1}{2}(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)\end{pmatrix}
$$这个矩阵有点庞大，但它确实能写成分块矩阵的形式：
$$\mathrm{I}’=\begin{pmatrix}
\lambda \text{I}_{3\times 3} & \text{P} \\ \text{P}^T & \lambda \text{I}_{3\times 3}
\end{pmatrix}$$其中$\lambda=\frac{1}{2}\sum_i m_i(r_{xi}^2+r_{yi}^2+r_{zi}^2+r_{wi}^2)$，而$\text{P}$通常为一个非对称的三阶矩阵。如果原来的矩阵是一个对角阵，变换到等角坐标后会怎样呢？设$$\mathrm{I}_1=\mathrm{diag}(I_{xy},I_{xz},I_{xw},I_{yz},I_{yw},I_{zw})$$将其变换至等角坐标，得到$$\mathrm{I}’_1=\frac{1}{2}\begin{pmatrix}I_{xy}+I_{zw}&0&0&I_{xy}-I_{zw}&0&0\\0&I_{xz}+I_{yw}&0&0&I_{xz}-I_{yw}&0\\0&0&I_{xw}+I_{yz}&0&0&I_{xw}-I_{yz}\\I_{xy}-I_{zw}&0&0&I_{xy}+I_{zw}&0&0\\0&I_{xz}-I_{yw}&0&0&I_{xz}+I_{yw}&0\\0&0&I_{xw}-I_{yz}&0&0&I_{xw}+I_{yz}\end{pmatrix}$$跟前面的$\mathrm{I}’$矩阵表达式对比可知，等角坐标下的主对角线上元素必须相等，于是有$$I_{xw}+I_{yz}=I_{xy}+I_{zw}=I_{xw}+I_{yz}$$从这里我们能够看出，虽然六阶矩阵对角化后有六个自由度，但转动惯量的几何性质让这个自由度只剩4个，即仅剩超长方体的四个尺寸的自由度。要证主轴定理只差最后一步：证明任何矩阵$\mathrm{I}’$都能通过四维空间的旋转变换变成$\mathrm{I}’_1$的形式。对比$\mathrm{I}’$与$\mathrm{I}’_1$可知，原来的对角矩阵转换至等角坐标下却不再是对角阵，而是$\mathrm{I}’$中的子阵$\mathrm P$也是对角阵的分块矩阵，我们只要找到合适的旋转变换能将$\text{P}$对角化即可找到旋转至主轴坐标系的变换。
求解四维旋转变换
四维空间中的旋转变换作用在2-向量上也相当于一个六维的矩阵，它在等角坐标下的形式非常简单：由于左右等角旋转的对易性，任意的四维2-向量旋转矩阵转化到等角坐标系下将会变成含两个三阶正交矩阵（$\text{A}$、$\text{B}$）的分块矩阵：
$$\text{R}=\begin{pmatrix}
\text{A} & 0 \\ 0 & \text{B}
\end{pmatrix}$$我们正是要靠这个旋转矩阵来把矩阵$\text{P}$化成某个对角阵$\Lambda$：$$\begin{pmatrix}
\text{A} & 0 \\ 0 & \text{B}
\end{pmatrix} \begin{pmatrix}
\lambda \text{I}_{3\times 3} & \text{P} \\ \text{P}^T & \lambda \text{I}_{3\times 3}
\end{pmatrix} \begin{pmatrix}
\text{A}^T & 0 \\ 0 & \text{B}^T
\end{pmatrix} =\begin{pmatrix}
\lambda \text{I}_{3\times 3} & \Lambda \\ \Lambda & \lambda \text{I}_{3\times 3}
\end{pmatrix}$$
计算左边矩阵乘法，对比等式两边可得出$\Lambda=\text{AP} \text{B}^T$。由于$\text{A}$、$\text{B}$都是正交阵，我们可以通过奇异值分解这个经典算法来求解$\text{A}$、$\text{B}$，求解完后则得到了等角坐标系下的完整的四维2-向量旋转矩阵$\text{R}$。
注意矩阵$\text{R}$是六阶的，所以还需要求出原来的四维旋转矩阵或旋转的旋量/四元数表示。最优雅的方式是直接在等角坐标系下完成这个转换工作：正交阵$\text{A}$、$\text{B}$代表的其实就是该旋转的左右手等角旋转部分，把这些三阶正交矩阵表示的旋转直接转换成四元数表示的旋转就能得到整个旋转的四元数表示，最后要注意的是由于四元数右乘的顺序问题，矩阵B需要取一次转置再映射到四元数，正好将奇异值分解得到的$\text A$与$\text B^T$直接分别转成左乘与右乘的四元数即得到最终的四元数版旋量表示的四维旋转。