NSolve/FindRoot matrix equation and conditions

Question

I have an eigenvalue (vector) problem (please see the theory, maybe I missed something in the post or did not understand Eqs (3.5.24) ) how to obtain two unknown matrices (36 unknowns) to satisfy all given conditions. Unknown matrices are

 XX = Table[X[10 i + j], {i, 6}, {j, 3}];
 YY = Table[Y[10 i + j], {i, 6}, {j, 3}];

which should be joined and present the 36 unknowns

 unknown = ArrayFlatten[({{XX, YY}})];

Let's suppose that we have two matrices in the following form:

JJ = {{0.`, 0.`, 0.`, 1.`, 0.`, 0.`}, {0.`, 0.`, 0.`, 0.`, 1.`, 
0.`}, {0.`, 0.`, 0.`, 0.`, 0.`, 1.`}, {-1.`, 0.`, 0.`, 0.`, 0.`, 
0.`}, {0.`, -1.`, 0.`, 0.`, 0.`, 0.`}, {0.`, 0.`, -1.`, 0.`, 0.`, 
0.`}};

 HH={{5.33449*10^7, -3.14159*10^7, 3.99684*10^6, 0., 199.842, 
 0.}, {-3.14159*10^7, 2.82358*10^7, -200000., -200., 
 0., -199.842}, {3.99684*10^6, -200000., 4.23581*10^6, 0., 199.842, 
 0.}, {0., -200., 0., 0.01, 0., 0.}, {199.842, 0., 199.842, 0., 
 0.00999211, 0.}, {0., -199.842, 0., 0., 0., 0.00999211}};

We obtained eigenvalues s1, s2 and s3, (s4=-s1, s5=-s2, s6=-s3) and take theirs imaginary parts

q0 = Eigenvalues[N[{HH, JJ}, 5]];

s1 = N[Im[q0]][[5]];
s2 = N[Im[q0]][[3]];
s3 = N[Im[q0]][[1]]; 

to construct the diagonal matrix s in the form

s = {{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}};

Now, we need to satisfy (find unknowns) with NSolve or FindRoot all conditions where

 zero = Table[0, {i, 3}, {j, 3}];
 onem = IdentityMatrix[3];


 one = Transpose[XX].JJ.XX;
 two = Transpose[YY].JJ.YY;
 three = Transpose[YY].HH.XX;
 four = Transpose[XX].HH.YY;

 five = Transpose[XX].HH.XX;
 six = Transpose[YY].HH.YY;
 seven = Transpose[XX].JJ.YY;
 eight = -Transpose[YY].JJ.XX;


 FindRoot[{one == zero, two == zero, three == zero, four == zero, 
 five == s, six == s, seven == onem, 
eight == onem}, {{X[11], 1}, {X[12], 1}, {X[13], 1}, {X[21], 
1}, {X[22], 1}, {X[23], 1}, {X[31], 1}, {X[32], 1}, {X[33], 
1}, {X[41], 1}, {X[42], 1}, {X[43], 1}, {X[51], 1}, {X[52], 
1}, {X[53], 1}, {X[61], 1}, {X[62], 1}, {X[63], 1}, {Y[11], 
1}, {Y[12], 1}, {Y[13], 1}, {Y[21], 1}, {Y[22], 1}, {Y[23], 
1}, {Y[31], 1}, {Y[32], 1}, {Y[33], 1}, {Y[41], 1}, {Y[42], 
1}, {Y[43], 1}, {Y[51], 1}, {Y[52], 1}, {Y[53], 1}, {Y[61], 
1}, {Y[62], 1}, {Y[63], 1}}]

but we have more equations then unknowns. If I fill use only first four equations, I can get the solution, but how to satisfy the rest number of equations?

Also, I choose roots around 1. How to get one set of solutions to satisfy all 8 matrix equations?

Answer 1

The book's author's primary goal seems to be his self-entertainment.

He wants to solve the first order ODE

$$\dot{\mathbf{u}}(t) = \mathbf{J} \, \mathbf{H} \, \mathbf{u}(t)$$

and applies an exponential ansatz

$$ \mathbf{u} = \sum_{r=1}^{2 n} \hat{\mathbf{u}}_r \operatorname{e}^{\lambda_r \,t + c_r}$$

with integration constants $c_r$ and solutions $(\hat{\mathbf{u}}_r,\lambda_r)$ of the eigenvalue problem

$$\mathbf{J} \, \mathbf{H} \, \mathbf{u}_r = \lambda_r \, \mathbf{u}_r.$$

So. Let's. Just. Solve. That.

n = 3;
JJ = ArrayFlatten[{{ConstantArray[0., {n, n}], 
     N@IdentityMatrix[n]}, {-N@IdentityMatrix[n], 
     ConstantArray[0., {n, n}]}}];
HH = N[{{4.9348022005446793094172455`5.*^7, -3.1415926535897937`*^7, 
     0, 0, 0, 0}, {-3.1415926535897937`*^7, 
     2.0238963636413604`*^7, -200000.`5., 0, 0, 0}, {0, -200000.`5., 
     238963.6364136009748945761`5., 0, 0, 0}, {0, 0, 0, 0.01`5., 0, 
     0}, {0, 0, 0, 0, 0.0099921105457425255`5., 0}, {0, 0, 0, 0, 0, 
     0.0099921105457425255`5.}}];

{λ, vecs} = Eigensystem[JJ.HH];

Bamm. That's it.

Let's synthesize the general solution $\mathbf{u}(t)$ of the ODE and test it:

u = t \[Function] Evaluate[Simplify[Total[Exp[t λ + Array[C, Length[λ]]] vecs]]];
u'[t] - JJ.HH.u[t] // Simplify // Chop

{0, 0, 0, 0, 0, 0}

Done.

Afterwards, the author means to say (but he did't):

$$\mathbf{x}_r = \Re(\hat{\mathbf{u}}_e) \quad \text{and} \quad \mathbf{y}_r = \Im(\hat{\mathbf{u}}_r).$$

In Mathematica terms, this is

x = Re[vecs];
y = Im[vecs];

Afterwards, he collects a bunch of equations that are automatically satisfied. Here are some of them (3.5.12 and 3.5.13):

ppos = Flatten[Position[λ/I // Chop, _?Positive]];
npos = Flatten[Position[λ/I // Chop, _?Negative]];
ω = 0.5 (λ[[ppos]]/I - λ[[npos]]/I) // Chop;

Table[HH.x[[ppos[[i]]]] - (ω[[i]] JJ.y[[ppos[[i]]]]), {i, 1, n}] // Chop
Table[HH.x[[npos[[i]]]] - (-ω[[i]] JJ.y[[npos[[i]]]]), {i, 1, n}] // Chop

{{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}}

{{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}}

The matrices $\mathbf{X}$ and $\mathbf{Y}$ can be obtained from that by

XX = Transpose[x[[1 ;; n]]];
YY = Transpose[y[[1 ;; n]]];

I leave checking all the other conditions to you.