I have the Hamiltonian for system of $N$ particles moving in one dimension:
$$H = \sum_{ij} M_{ij} \sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }$$ where  $(q_i,p_i)$ are phase space variables, $i\in\{ 1,\dots, N\}$, and $M_{ij}$ is a constant symmetric matrix. The stationary equilibrium configuration of this system is $q_i=a$ and $p_i=b$ for all $i$ and arbitrary $a$ and $b$. (Static equilibrium requires $b=0$.) How does the system evolve after perturbed infinitesimally from the equilibrium state?

This is the Hamiltonian of vector resonant relaxation of a thin stellar disk. 

**EDIT:**
The equations of motions are
$$\dot{q}_i = \frac{\partial H}{\partial p_i} = \sum_j M_{ij} \frac{p_i - p_j}{\sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }}\,,$$ 
$$\dot{p}_i = -\frac{\partial H}{\partial q_i} = \sum_j M_{ij} \frac{q_j - q_i}{\sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }}\,.$$ 

These equations are remarkably simple in the complex plane. Define $z_i = q_i + i p_i$, then
$$\dot{z}_n = \sum_m i M_{nm} \frac{z_m - z_n}{|z_m - z_n|}\,.$$
Notice that the "effective force" is constant at zero distance it depends on only the absolute value of the relative coordinates, similar to a charge oscillating around an charged insulating plate.

Is there a way to derive the normal mode oscillations of this nonlinear system for "small" perturbations around $z_i=0$?