A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action
$$S_L[x_i] = \int_{t_0}^{t_1}\mathbb{d} t \left\{ \sum_{i=0}^{N-1} \frac{1}{2}m \dot x_i^2 - \frac{1}{2}m\omega^2\left(x_i - x_{i+1}\right)^2 \right\}\,,$$
where $x_N \coloneqq x_0$.

This system can be solved in terms of travelling waves, obtained by discrete Fourier transform. Initial condition in #35 is desired, see below.

Unfortunately discrete Fourier transform is a transformation in the complex field. We first write the system in matrix form
$$S_L[x_i] = \int_{t_0}^{t_1}\mathbb{d} t \left\{ \frac{1}{2}m \dot x_i \delta_{ij} \dot x_j - \frac{1}{2}m x_i A_{ij} x_j\right\}\,,$$
where $A_{ij} / \omega^2$ is equal to $(-2)$ if $i=j$, $1$ if $|i-j|=1$ or $|i-j|=N$, and $0$ otherwise.

Then we can extend it with non-dynamic phases $\phi_i \in \mathbb{R}$, and introduce $X_i \coloneqq x_i \mathbb{e}^{-\phi_i}$, $X^\ast_i \coloneqq x_i \mathbb{e}^{+\phi_i}$ so that
$$S_L[x_i] = S_L[x_i, \phi_j] \equiv S_L[X^\ast_i, X_j] = \int_{t_0}^{t_1}\mathbb{d} t \left\{ \frac{1}{2}m \dot X^\ast_i \delta_{ij} \dot X_j - \frac{1}{2}m X^\ast_i A_{ij} X_j\right\}\,.$$

$A_{ij}$ can be diagonalised by the inverse discrete Fourier transform (motivation / proof ?)
$$X_i = (F^{-1})_{ik} Y_k = \frac{1}{\sqrt{N}}\sum_k \mathbb{e}^{i \frac{2\mathbb{\pi}}{N} k\mathbb{i}} Y_k\,.$$

Note $F^{-1}$ is unitary, but not special unitary. Can we make it special unitary? Maybe use the fact that $\sqrt{N} F^{-1}$ makes a Vandermonde determinant?

Calculating gives (motivation / proof ?)
$$S_L[X^\ast_i, X_j] = S_L[Y^\ast_i, Y_j] = \sum_{k=0}^{N-1} \int_{t_0}^{t_1}\mathbb{d} t \left\{ \frac{1}{2}m \dot Y^\ast_k \dot Y_k - \frac{1}{2}m \omega^2\cdot4\sin^2\frac{2\mathbb{\pi}k}{N} Y^\ast_k Y_k\right\}\,.$$
Using the same transformation to separate the non-dynamic phases, we can arrive at a real action
$$S_L[y] = \sum_{k=0}^{N-1} \int_{t_0}^{t_1}\mathbb{d} t \left\{ \frac{1}{2}m \dot y_k^2 - \frac{1}{2}m \omega^2\cdot4\sin^2\frac{2\mathbb{\pi}k}{N} y_k^2\right\}\,.$$

We can solve $N$ independent oscillators
$$\dot y_k^2 + 4\omega^2\sin^2\frac{2\mathbb{\pi}k}{N} y_k^2 \equiv 4\omega^2\sin^2\frac{2\mathbb{\pi}k}{N} y_{k0}^2$$
with initial conditions $(y_{k0}, t_{k0})$ as in in #35, find suitable constant phases so that $x_i$'s are real, and generate the data.