在Mathematica中创建一个符号正交矩阵

Question

我需要在Mathematica中创建一个3x3的实正交符号矩阵。如何做到呢？

- marcellus

---

什么类型的矩阵？Mathematica内置了旋转矩阵和反射矩阵函数，两者都是正交的。 - Niki

我想构建一个符号矩阵，使得该矩阵在后续计算中始终被视为正交矩阵。Mathematica中的RotationMatrix命令无法实现此功能。 - marcellus

在科学计算堆栈交换网站上有一个相关的问题：http://scicomp.stackexchange.com/questions/74/symbolic-software-packages-for-matrix-expressions - MRocklin

Answer 1

虽然我不建议这样做，但是...

m = Array[a, {3, 3}];
{q, r} = QRDecomposition[m];
q2 = Simplify[q /. Conjugate -> Identity]

假设我们在实数域上工作，那么q2是一个符号正交矩阵。

- Daniel Lichtblau

---

Daniel，感谢您的回答。您是否看到有任何直接实现正交性的六个条件和行列式为1的条件的方法，而不涉及QR分解？谢谢。Marcello - marcellus

Answer 2

我觉得您想在Mathematica中使用SO(3)群参数化。由于向量之间的相互正交和模长相等，您只有3个独立符号（变量），因此有6个约束条件。一种方法是构造绕三个轴线的独立旋转，并将这些矩阵相乘。以下是实现该方法可能过于复杂的代码:

makeOrthogonalMatrix[p_Symbol, q_Symbol, t_Symbol] :=
  Module[{permute, matrixGeneratingFunctions},
    permute =  Function[perm, Permute[Transpose[Permute[#, perm]], perm] &];
    matrixGeneratingFunctions = 
       Function /@ FoldList[
            permute[#2][#1] &,
            {{Cos[#], 0, Sin[#]}, {0, 1, 0}, {-Sin[#], 0, Cos[#]}},
            {{2, 1, 3}, {3, 2, 1}}];
    #1.#2.#3 & @@  MapThread[Compose, {matrixGeneratingFunctions, {p, q, t}}]];

这是它的工作原理：

In[62]:= makeOrthogonalMatrix[x,y,z]
Out[62]= 
{{Cos[x] Cos[z]+Sin[x] Sin[y] Sin[z],Cos[z] Sin[x] Sin[y]-Cos[x] Sin[z],Cos[y] Sin[x]},
 {Cos[y] Sin[z],Cos[y] Cos[z],-Sin[y]},
 {-Cos[z] Sin[x]+Cos[x] Sin[y] Sin[z],Cos[x] Cos[z] Sin[y]+Sin[x] Sin[z],Cos[x] Cos[y]}}

您可以通过对各列（或行）点积进行Simplify，检查矩阵是否正交。

- Leonid Shifrin

---

Leonid，感谢您的回答。实际上，我想定义一个通用的符号SO（3）矩阵，而不是像您一样从欧拉角开始。基本上，我想设置一个通用矩阵（例如使用mat = Table [Subscript [m，i，j]，{i，3}，{j，3}]），并强制要求该矩阵的元素始终满足正交性条件和行列式= 1条件，而无需在以后指定这一点。Marcello - marcellus

@user1087909，那么Daniel的答案应该是正确的方法。 - Leonid Shifrin

出奇地，在回复我的评论/问题时，我要说Leonid的方法是正确的。 - Daniel Lichtblau

@Daniel 确实很有趣。我重新查看了我的代码，意识到我不应该使用 FoldList 等复杂的东西，而是应该明确列出所有三个矩阵 - 我的代码没有任何好处，非常令人困惑。很快会改成这个版本。 - Leonid Shifrin

Answer 3

我找到了一种“直接”的方法来实现特殊正交性。见下文。

(*DEFINITION OF ORTHOGONALITY AND SELF ADJUNCTNESS CONDITIONS:*) 
MinorMatrix[m_List?MatrixQ] := Map[Reverse, Minors[m], {0, 1}] 
CofactorMatrix[m_List?MatrixQ] := MapIndexed[#1 (-1)^(Plus @@ #2) &, MinorMatrix[m], {2}] 
UpperTriangle[ m_List?MatrixQ] := {m[[1, 1 ;; 3]], {0, m[[2,   2]], m[[2, 3]]}, {0, 0, m[[3, 3]]}}; 
FlatUpperTriangle[m_List?MatrixQ] := Flatten[{m[[1, 1 ;; 3]], m[[2, 2 ;; 3]], m[[3, 3]]}];
Orthogonalityconditions[m_List?MatrixQ] := Thread[FlatUpperTriangle[m.Transpose[m]] == FlatUpperTriangle[IdentityMatrix[3]]]; 
Selfadjunctconditions[m_List?MatrixQ] := Thread[FlatUpperTriangle[CofactorMatrix[m]] == FlatUpperTriangle[Transpose[m]]]; 
SO3conditions[m_List?MatrixQ] := Flatten[{Selfadjunctconditions[m], Orthogonalityconditions[m]}]; 

(*Building of an SO(3) matrix*) 
mat = Table[Subscript[m, i, j], {i, 3}, {j, 3}]; 
$Assumptions = SO3conditions[mat]

那么

Simplify[Det[mat]] 

提供 1；...以及

MatrixForm[Simplify[mat.Transpose[mat]]

提供了单位矩阵； ...最后

MatrixForm[Simplify[CofactorMatrix[mat] - Transpose[mat]]]

返回一个零矩阵。

========================================================================

这正是我在提问时寻找的内容！然而，请告诉我您对这种方法的看法。

Marcellus

- marcellus

Answer 4

Marcellus，您需要使用SO(3)的一些参数化方法，因为您的一般矩阵必须反映该组的RP³拓扑结构。没有单一的参数化方法可以涵盖整个群，而不会出现多值性或奇异点。维基百科有一个关于SO(3)上各种图表的好页面。

也许最概念上最简单的之一是来自李代数so(3)的指数映射。定义一个反对称的实数A（它跨越so(3)）。

A = {{0, a, -c},
     {-a, 0, b},
     {c, -b, 0}};

那么MatrixExp[A]是SO(3)的一个元素。 我们可以使用以下方法来验证：

Transpose[MatrixExp[A]].MatrixExp[A] == IdentityMatrix[3] // Simplify

如果我们写成t^2 = a^2 + b^2 + c^2，我们可以简化矩阵指数。

{{   b^2 + (a^2 + c^2) Cos[t]  , b c (1 - Cos[t]) + a t Sin[t], a b (1 - Cos[t]) - c t Sin[t]}, 
 {b c (1 - Cos[t]) - a t Sin[t],    c^2 + (a^2 + b^2) Cos[t]  , a c (1 - Cos[t]) + b t Sin[t]}, 
 {a b (1 - Cos[t]) + c t Sin[t], a c (1 - Cos[t]) - b t Sin[t],    a^2 + (b^2 + c^2) Cos[t]}} / t^2

请注意，这基本上是与RotationMatrix相同的参数化。 与输出进行比较。

RotationMatrix[s, {b, c, a}] // ComplexExpand // Simplify[#, Trig -> False] &;
% /. a^2 + b^2 + c^2 -> 1

- Simon

---

非常好的回答。然而，我更喜欢直接定义一个满足SO（3）条件的矩阵...请参见我的上面的答案，并告诉我您的想法。Marcellus - marcellus

Answer 5

尽管我非常喜欢Marcellus对自己问题的回答，但它并不完全正确。 不幸的是，他得出的条件也会导致

Simplify[Transpose[mat] - mat]

评估为零矩阵！这显然是不正确的。这里有一种既正确又更直接的方法：

OrthogonalityConditions[m_List?MatrixQ] := Thread[Flatten[m.Transpose[m]] == Flatten[IdentityMatrix[3]]];
SO3Conditions[m_List?MatrixQ] := Flatten[{OrthogonalityConditions[m], Det[m] == 1}];

即将一个旋转矩阵乘以它的转置矩阵会得到单位矩阵，而旋转矩阵的行列式为1。

- twan