To obtain the modes shapes and resonant frequencies, you start from your equation of motion with no externally applied forces, which is indeed as you've stated.
$$\mathbf M \mathbf{\ddot q} + \mathbf K \mathbf q = \mathbf 0 \qquad (1)$$
For brevity, I've let $\mathbf K = \mathbf K_b + \mathbf K_m$.
Currently, $\mathbf q(t)$ is a function of time. If the system is vibrating in the $i$th mode, everywhere will be vibrating at frequency $\omega_i$ in the same phase, assuming there is no damping. Therefore, we can assume the following:
$$\mathbf q(t) = \mathbf u^{(i)}\exp{(j\omega_i t)}\qquad(2)$$
, where $j$ is the imaginary number defined as $\sqrt{-1}$ and $\mathbf u^{(i)}$ is the $i$th mode shape, which is independent of time. (This assumption may seem bizarre as $\exp{(j\omega_i t)}$ is a complex number, but this does not matter in this scenario, and is used just to show that $\mathbf q(t)$ is sinusoidal in time with frequency $\omega_i$. Equivalently, you could instead use $\mathbf q(t) = \mathbf u^{(i)} (A \sin{\omega_i t} + B \cos{\omega_i t})$.)
Note that:
$$\mathbf {\ddot q}(t) = \frac{d}{dt}\left(\mathbf u^{(i)}\exp{(j\omega_i t)}\right)=-\omega_i^2\mathbf u^{(i)}\exp{(j\omega_i t)}\qquad (3)$$
Substituting equations (2) and (3) into (1) will give you the following, noting that the exponential functions cancel out:
$$\left(-\omega_i^2 \mathbf M + \mathbf K \right) \mathbf u^{(i)}=0 \qquad (4)$$
Equation 4 is a type of problem called the generalised eigenvalue problem, which involves solving equations of the form $\mathbf A \mathbf x = \lambda \mathbf B \mathbf x$ (compare with the standard eigenvalue problem, which involves solving $\mathbf A \mathbf x = \lambda \mathbf x$). In this case, the eigenvalue is $\omega_i^2$ and the eigenvector is $\mathbf u^{(i)}$. To obtain the eigenvalues and eigenvectors involves solving the generalised eigenvalue problem. There are two main approaches: hand calculation or numerical computation. For your problem, the latter will definitely be the recommended option, however, the former should help with smaller problems and lends some insight to how generalised eigenvalue problems are solved.
HAND CALCULATION
Note how equation (4) also has the form $\mathbf A \mathbf u^{(i)} = \mathbf 0$. Obviously, $\mathbf u^{(i)}=\mathbf 0$ would satisfy such an equation, but this is a trivial result that corresponds to no actual vibration. Instead, we are interested in values of $\mathbf u^{(i)}\ne\mathbf 0$ that solve the equation. In order for such a non-zero answer to exist, $\mathbf A$ must have no inverse, i.e. $\mathbf A$ is a singular matrix. (If an inverse did exist, then you could premultiply both sides of the equation with the inverse matrix, getting $\mathbf u^{(i)} = \mathbf A^{-1} \mathbf 0 = \mathbf 0$.) In order for a matrix to be singular, the determinant of the matrix must be zero. i.e.
$$\left| \mathbf A \right| = \left| -\omega_i^2 \mathbf M + \mathbf K \right| = 0$$
This will produce a polynomial in $\omega_i^2$, the roots of which will be each of the resonant frequencies squared. Take care to note that some modes will have the same resonant frequencies, in the form of repeated roots.
For example, say:
$$\mathbf M = \left[\begin{matrix}
m & 0\\
0 & m\\
\end{matrix}\right]
\qquad
\mathbf K = \left[\begin{matrix}
2k & -k\\
-k & 2k\\
\end{matrix}\right]
$$
The determinant of interest would be:
$$
\left|\begin{matrix}
2k-\omega_i^2 m & -k\\
-k & 2k-\omega_i^2 m\\
\end{matrix}\right|
=0
$$
By expanding out the determinant expression, we obtain the polynomial:
$$(2k-\omega_i^2 m)^2-k^2=0$$
The two roots of which are:
$$\omega_i^2 = \frac{k}{m},\frac{3k}{m}$$
So the two resonant frequencies are:
$$\omega_1 = \sqrt{\frac{k}{m}} \qquad \omega_2 = \sqrt{\frac{3k}{m}}$$
Note that you should obtain a number of (not necessarily distinct) eigenvalues equal to the number of rows or columns of your $\mathbf M$ or $\mathbf K$ matrices.
To obtain the eigenvector corresponding to a particular eigenvalue, substitute back into equation (4), to obtain a series of simultaneous equation, the solving of which should determine the eigenvector.
In our example, finding the eigenvector corresponding to the eigenvalue $\omega_1^2$, substitute $\omega_i=\omega_1$ into equation (4)...
$$
\left[\begin{matrix}
2k-k & -k\\
-k & 2k-k\\
\end{matrix}\right]
\left[\begin{matrix}
u^{(1)}_1\\
u^{(1)}_2\\
\end{matrix}\right]
=
\left[\begin{matrix}
k u^{(1)}_1-k u^{(1)}_2\\
-k u^{(1)}_1+k u^{(1)}_2\\
\end{matrix}\right]
=
\left[\begin{matrix}
0\\
0\\
\end{matrix}\right]
$$
From this, we can see that $u^{(1)}_1 = u^{(1)}_2$, and so the eigenvector is therefore:
$$\mathbf u^{(1)} =
\left[\begin{matrix}
1\\
1\\
\end{matrix}\right]
$$
NUMERICAL COMPUTATION
Your problem does not seem practical to do by hand. Thankfully, there are a variety of different eigen-solver algorithms out there that are dedicated to solving this very problem.
A couple of these include:
- Power Iteration method
- Lanczos method
- Krylov-Schur method
- Jacobi-Davidson method
To explain how each of these algorithms work would take more time than I currently have, however I would definitely recommend tracking down these algorithms, reading how they work and using them to solver the generalised eigenvalue problem.
Hope this helps.
