The inverse problem of rotation of stress matrix

Question

I want to find the degree of rotation between a stress state and the initial stress state.

em[α_, β_, γ_] := 
 EulerMatrix[{α, β, γ}]\[Transpose].{{10, 0, 
    0}, {0, 20, 0}, {0, 0, 
    25}}.EulerMatrix[{α, β, γ}]
res = NMinimize[
  Norm[EulerMatrix[{Pi/2, 0, 0}]\[Transpose].{{10, 0, 0}, {0, 20, 
       0}, {0, 0, 25}}.EulerMatrix[{Pi/2, 0, 0}] - 
    Transpose[em[a, b, c]], "Frobenius"], {a, b, c}, WorkingPrecision->50]

I use the above code to get the solution {a - > 1.16727157, B - > 1.41653756 * 10 ^ - 12, C - > 0.403524756}. But the error between C which is about 0.4 * 180 and the real rotation angle 0 is large. I want to get {Pi/2, 0, 0} or a numerical solution with a small error with {Pi/2, 0, 0}. What shall I do to get a more accurate answer?

Response to comments:

Even if I limit the three variables from - Pi to Pi, the error of the result (c->-0.221782576) is still large. I feel that the norm selection is not accurate, but I don't know how to further improve this code to get more accurate results ( I want to get {Pi/2, 0, 0} or a numerical solution with a small error with {Pi/2, 0, 0}).

em[α_, β_, γ_] := 
 EulerMatrix[{α, β, γ}]\[Transpose].{{10, 0, 
    0}, {0, 20, 0}, {0, 0, 
    25}}.EulerMatrix[{α, β, γ}]
res = NMinimize[{Norm[
    EulerMatrix[{Pi/2, 0, 0}]\[Transpose].{{10, 0, 0}, {0, 20, 0}, {0,
         0, 25}}.EulerMatrix[{Pi/2, 0, 0}] - Transpose[em[a, b, c]], 
    "Frobenius"], (0 <= a <= Pi) && (-Pi <= b <= Pi) && (-Pi <= c <= 
      Pi)}, {a, b, c}]

In other words, I want to use Mathematica to solve the following matrix equation precisely:

A\[Transpose].{{10, 0, 0}, {0, 20, 0}, {0, 0, 25}}.A == {{35/2, (
   5 Sqrt[3])/2, 0}, {(5 Sqrt[3])/2, 25/2, 0}, {0, 0, 25}}

The referenced answer of matrix A is EulerMatrix[{Pi/3, 0, 0}].

If you don't constrain your angles, don't expect a unique solution. — J. M.'s missing motivation, Mar 25 '20 at 05:31
@J.M.'stechnicaldifficulties I have limited the range of variables, but the error of this result is still large. How can I further modify it? — A little mouse on the pampas, Mar 25 '20 at 07:53
@J.M.'stechnicaldifficulties I see you answer a similar question here . Thank you for your answer. — A little mouse on the pampas, Apr 05 '20 at 02:13

Henrik Schumacher · Accepted Answer · 2020-03-30T05:39:15.167

So the problem is as follows: Given two symmetric matrices A and B, find a rotation R that minimizes Norm[A - Transpose[R].B.R]. It is well-known that such a rotation has to map the principal axes = eigenspaces of A onto those of B. More precisely, R maps the eigenspace of the smallest eigenvalue of A onto that one of the smallest eigenvalue of B and so on.

Hence, this problem can be solved by using Eigensystem.

A = #\[Transpose].# &@RandomReal[{-1, 1}, {3, 3}];
R = RandomVariate[CircularRealMatrixDistribution[3]];
B = R\[Transpose].A.R;

{λ, U} = Eigensystem[A];
U = Normalize /@ U; (* only necessary for exact and symbolic A*)
{μ, V} = Eigensystem[B];
V = Normalize /@ V; (* only necessary for exact and symbolic B*)

Now the rows of U and V correspond to unit vectors in the eigenspaces. However, each of these unit vectors has two possible directions and it cannot be told a priorily which one is chosen by Eigensystem. If the eigenvalues of A are pairwise distinct (and thos of B, too) R must be among one of the rotations in the list canditates:

signs = Select[Tuples[{1, -1}, 3], Times @@ # == Det[U] &];
canditates = Table[U\[Transpose].(s V), {s, signs}]

I point out htat each matrix S in the list canditates is a rotation that satisfies B == S\[Transpose].A.S, so this problem definitely has more than one solution.

Finally, if really required, you can obtain the Euler angles with

EulerAngles /@ canditates

Thank you very much for your answer. Your code is accurate and fast. — A little mouse on the pampas, Mar 29 '20 at 09:21
I have another question, how to use your algorithm to solve this matrix equation: A\[Transpose].{{10, 0, 0}, {0, 20, 0}, {0, 0, 25}}.A == {{35/2, ( 5 Sqrt[3])/2, 0}, {(5 Sqrt[3])/2, 25/2, 0}, {0, 0, 25}}. A's reference answer is EulerMatrix[{Pi/3, 0, 0}]. — A little mouse on the pampas, Mar 30 '20 at 05:27
Ah good point. Here the input matrices are exact and Eigensystem does not normalize the eigenvectors in this case. I added a fix for that. — Henrik Schumacher, Mar 30 '20 at 05:40

score 1 · Answer 2 · answered Mar 30 '20 at 07:25

It can be solved in this way, but the warning information will be output:

Solve[Thread[
   EulerMatrix[{a, b, c}]\[Transpose].{{10, 0, 0}, {0, 20, 0}, {0, 0, 
       25}}.EulerMatrix[{a, b, c}] == {{35/2, (5 Sqrt[3])/2, 0}, {(
      5 Sqrt[3])/2, 25/2, 0}, {0, 0, 25}}], {a, b, c}] // FullSimplify

score 0 · Answer 3 · answered Apr 05 '20 at 01:59

r = {{10, 0, 0}, {0, 20, 0}, {0, 0, 25}};
rt = RotationMatrix[{{1, 0, 0}, {0, 1, 0}}]\[Transpose].{{10, 0, 
      0}, {0, 20, 0}, {0, 0, 
      25}}.RotationMatrix[{{1, 0, 0}, {0, 1, 0}}] // FullSimplify;

fg = FindGeometricTransform[r, rt, TransformationClass -> "Rigid"]

EulerAngles[Drop[TransformationMatrix[Last[fg]], -1, -1]]

The inverse problem of rotation of stress matrix

3 Answers3

Linked