Solution (source code)

= Solution

The variables $y_j=I_je^{i\Omega_j}$ are <complex inclinations>, despite the question's description as eccentricities. Diagonalize the <Laplace-Lagrange inclination matrix>: let real <eigenvectors> $\mathbf v_k$ satisfy $B\mathbf v_k=\lambda_k\mathbf v_k$, and form $V=(\mathbf v_1,\mathbf v_2,\mathbf v_3)$. The physical matrix is diagonalizable by <positive diagonal symmetrization of a matrix>, as shown below.

The <matrix exponential> solution is
$$
\mathbf y(t)=e^{iBt}\mathbf y(0)
=\sum_{k=1}^3 c_k\mathbf v_ke^{i\lambda_kt},\qquad
(c_1,c_2,c_3)^T=V^{-1}\mathbf y(0).
$$
Write $c_k=I_ke^{i\gamma_k}$ with $I_k\ge0$, and define $I_{jk}=I_kv_{jk}$. This gives
$$
\boxed{y_j(t)=\sum_{k=1}^3I_{jk}e^{i(\lambda_kt+\gamma_k)}.}
$$
The <eigenvalues> determine nodal precession frequencies, the <eigenvectors> determine relative inclinations within each <secular eigenmode>, and the initial <complex inclinations> determine their amplitudes and phases. The real $I_{jk}$ can have either sign, encoding opposite nodal phases. The zero-frequency mode is a common fixed tilt.