Why do systems of $n$ coupled oscillators have $n$ normal modes? Consider a linear system of $n$ differential equations with constant coefficients corresponding to a physical scenario where I have $n$ coupled oscillators (like $n$ masses attached by springs in series or something like that, also no friction). I was told that if I have $n$ oscillators in my coupled system, then I will always have $n$ distinct normal modes corresponding to this system.
For that to be the case, though, my coefficient matrix would need to have $n$ distinct eigenvalues which correspond to each normal mode. Going further, if my system had $n$ distinct normal modes it would mean that all of my normal modes could be decoupled since my coefficient matrix could be diagonalized by a suitable change of basis, thus suggesting that all linear coupled systems with constant coefficients have normal modes which are independent from one another.
My question pertains to how we could say in general that our matrices will have $n$ distinct eigenvalues as from what I am aware, generally speaking that may not be the case. Any help would be appreciated in understanding this generalization.
 A: The easiest way to solve the problem is to take advantage of the Lagrangian formalism. Denoting the $n$ generalized coordinates by $x_1, \ldots x_n$, the most general form of a Lagrangian of the system under consideration is given by $$ L(x_1, \ldots x_n, \dot{x}_1, \ldots \dot{x}_n)= \frac{1}{2}\sum\limits_{i, j} \dot{x}_i M_{ij} \dot{x}_j - \frac{1}{2} \sum_{i,j} x_i V_{i j} x_j = \frac{1}{2} \dot{x}^TM\dot{x}-\frac{1}{2}x^T V x$$ with a strictly positive symmetric $n \times n$ mass matrix $M^T=M \gt 0$ and a nonnegative symmetric matrix $V^T= V \ge 0$.
In a first step, the mass matrix $M$ can be diagonalzed by a real orthogonal transformation $R$, such that  $M = R^T \hat{M} R$ with the diagonal matrix $\hat{M} = {\rm diag}(m_1, \ldots, m_n)$ (where $m_i \gt 0$) and $R^T = R^{-1}$. Introducing the new generalized coordinates $y= \hat{M}^{1/2} R \, x \Leftrightarrow x= R^T \hat{M}^{-1/2}  y$, the Lagrangian  takes the form $$L(y, \dot{y}) = \frac{1}{2} \dot{y}^T \dot{y}- \frac{1}{2}y^T K y$$ with the symmetric non-negative coupling matrix $K= \hat{M}^{-1/2}R V R^T \hat{M}^{-1/2} $.
In the final step, the matrix $K$ is diagonalized, $K =S^T \hat{K} S$, with $\hat{K} = {\rm diag} (\omega_1^2, \ldots, \omega_n^2)$ (with $\omega_\alpha^2 \ge 0$ as $K$ is non-negative) and a real orthogonal matrix $S$. Introducing the normal coordinates $Q= S y = S \hat{M}^{1/2} R \, x $, the Lagrangian takes the final form $$L(Q, \dot{Q})= \frac{1}{2} \dot{Q}^T \dot{Q} -\frac{1}{2} Q^T \hat{K} Q= \frac{1}{2} \sum\limits_{\alpha = 1}^n (\dot{Q}_\alpha^2 - \omega_\alpha^2 \, Q_\alpha^2).$$ Depending on the physical context, it is possible that some of the eigenvalues $\omega_\alpha$ might vanish (so-called zero modes), which occurs e.g. in the situation where two mass points are connected by a spring, which, however, does not affect the linear motion of the center-of-mass coordinates. Likewise, some of the eigenvalues may be degenerate, which allows you to arbitrarily rotate the normal coordinates in the corresponding eigenspace.
