Normal modes in physics - solution of set of 2nd order linear differential equations In physics we study normal modes for vibrating system. Mathematically speaking we have a set of linear differential equations. For example we have something like this:
$$
\begin{align*}
m \frac{dx_1}{dt^2} &= k(x_2-x_1) \\ 
M \frac{dx_2}{dt^2} &= k(x_3-x_2) - k(x_2-x_1) \\
m \frac{dx_3}{dt^2} &= -k(x_3-x_2)
\end{align*}
$$
We seek solutions for $x_i(t)$. A general solution is to assume to look for these functions to be sines, cosines or exponents for all $x_i(t)$ having the same frequency $\omega$. For above example our solution would be something like:
$$
\begin{align*}
x_1 (t) &= A e^{\omega t + \varphi} \\ 
x_2 (t) &= B e^{\omega t + \varphi} \\
x_3 (t) &= C e^{\omega t + \varphi}
\end{align*}
$$
I could not find a clear explanation why we assume that we seek solution with the same $\omega$ (normal modes), and later why the general solution is the superposition of all such solutions (we obviously find . Many thanks for comments and explanation.
 A: The system modeled as
$$
\begin{align*}
m \frac{d^2x_1}{dt^2} &= k(x_2-x_1) \\ 
M \frac{d^2x_2}{dt^2} &= k(x_3-x_2) - k(x_2-x_1) \\
m \frac{d^2x_3}{dt^2} &= -k(x_3-x_2)
\end{align*}
$$
can be represented shortly as
$$
\ddot x = A x
$$
Normally in this kind of problems $A$ can be represented as $T^{-1}\Lambda T$ with $T^{-1}T = I$ where $\Lambda$ is a diagonal matrix with the eigenvalues from $A$. So we can follow with
$$
\ddot y = \Lambda y,\ \ \ y = T x
$$
This system doesn't contain attenuation terms so it represents an oscillatory system. This system will remain forever oscillating along their natural frequencies  associated to the eigenvalues $\lambda_k$. The most plausible system with the properties you cited is an actuated dissipative system with the structure
$$
\ddot x + 2B \dot x + A x = u
$$
with $u = f(\omega t)$ in which $f$ represents a sinusoidal adequate function. The solution for this system can be written as
$$
y_k(t) = e^{-b_k t}\left(c_k e^{-\sqrt{b_k^2-\lambda_k}t}+\bar c_k e^{\sqrt{b_x^2-\lambda_k}t}\right)+ y_p(\omega t)
$$
For passive systems, the homogeneous component goes asymptotically to zero remaining the particular $y_p(\omega t)$
