Linear ODE and Fourier Series

Question

Let $m,k_0,k$ be positive real numbers and $x_1$, $x_2$ be real-valued functions of time.

Suppose we have following system of two coupled ODEs ( motivated by a coupled oscillator with two masses and three springs ) :

$m(x_1)'' = -k_0x_1 - k(x_1-x_2) $
$m(x_2)'' = -k_0x_2 + k( x_1 - x_2)$

By doing some algebraic manipulations i ended up finding out that the solution to both $x_1(t)$ and $x_2(t)$ results in a linear combination of two cosine terms, of two different frequencies.
More precisely :
$$x_1(t) = c_1\cos(w_0t + \phi_1) + c_2\cos(w_1t + \phi_2)$$
with $$w_0=\sqrt{\frac{k}{m}}$$
and $$w_1=\sqrt{\frac{k_0 +2k}{m}} $$

I was intrigued with the fact that possible highly chaotic solutions $x_1(t), x_2(t)$ ( the motions of each mass after arbitrary initial conditions on the coupled oscillator ) could be described so nicely as a linear combination of two simple solutions.
I wanted to understand why there exists such a nice simplification.
After asking in physics.stackexchange, a suggestion was to check Fourier series. Since the solutions $x_1(t),x_2(t)$ are periodic continuous functions ( no damping in our oscillator) , it's possible to decompose it into a fourier series, i was told.

For the past day, i've been reading and understanding Fourier Series but i still lack to see the exact connection in this case.

For instance, decomposing a periodic function into a fourier series might get us infinitely many cosine terms ( infinitely many coefficients ) ...
Why do we have only two in that case ?

Secondly, the frequency of the cosine terms in a fourier series are multiples of the lowest frequency ( fundamental frequency ).
Why the frequency of the fastest cosine term is not a multiple of the frequency of the slowest cosine term, but rather a function of $k_0$.

In summary, i want to understand how such chaotic motions of each mass after an arbitrary set of initial conditions could be so nicely described as a linear combination of two simple harmonic motions.
I can see how this might have to do with Fourier series, but i still lack the accurate understanding of such connection.

Thanks a lot in advance.

Answer 1

The key is to look at the form of the ODE, which comes from Newton's second law:

$$m_ia_i = m_i\ddot{x}_i = F_i(x_1, \ldots, x_n).$$

There are two observations to make here. First, if you have several particles and each of them obeys Newton's second law, you can group the equations all together into a configurational equation of motion $$M\ddot{x} = F(x)$$ where $M$ is a $n\times n$ positive-definite mass matrix (typically diagonal) and $x$ and $F$ are $n$-vectors of particle positions and forces. (Note that this extends trivially to the case where the particles are in 3D or other ambient dimension).

Next, for conservative forces (which most forces are, including the springs in your example) you can write the force as coming from a scalar potential function $V:\mathbb{R}^{n}\to\mathbb{R}$, with $F= -DV$ where $D$ is the Jacobian. Thus $$M\ddot{x} = -DV(x).$$ Now suppose $x_0$ is a minimizer of the potential energy. Then $DV(x_0)=0$ and $x_0$ is an equilibrium point: a point at which all net forces on all particles are zero. Moreover, since $x_0$ is a minimum, the second derivative $HV(x_0)$ is positive-definite (bowl-shaped).

Finally, let's expand the right-hand side of Newton's law using Taylor's theorem, for small perturbations away from the equilibrium point $x(t) = x_0 + \delta x(t)$: $$M\ddot{\delta x} = - HV(x_0)\delta x + O(\|\delta x\|^2) \qquad (*)$$ where we have used that $DV(x_0)=0$. For $\delta x$ small, the correction term is negligible and $$\ddot{\delta x} = -M^{-1} HV(x_0)\delta x.$$ Here is where the Fourier analysis comes in, which we will use here in another form. Let's expand $\delta x(t)$ in the eigenbasis of $M^{-1} HV(x_0)$: $\delta x(t) = \sum \alpha_i(t) v_i$. Then $$\sum \ddot{\alpha_i}v_i = -\sum \lambda_i \alpha_i v_i,$$ and since the eigenvectors are orthogonal, $$\ddot{\alpha_i} = -\lambda_i \alpha_i.$$ Notice that since $HV$ is positive-definite, the $\lambda_i$ are positive. This is a one-dimensional ODE that even I can solve: $$\alpha_j(t) = Ae^{i \sqrt{\lambda_j} t} + Be^{-i \sqrt{\lambda_j} t}$$ or, if we assume $\delta x(0) = 0$, $$\alpha_j(t) = A_j\sin\left(\sqrt{\lambda_j} t\right).$$

Whew! Finally we put the pieces together to get a full solution for $x(t)$:

$$x(t) = x_0 + \sum A_i \sin\left(\sqrt{\lambda_i} t\right)v_i.$$

Here's the punchline: near their equilibrium points, all conservative physical systems behave as a superposition of harmonic modes. This is why different shaped glasses ring at different frequencies, for instance; when you strike a wine glass, it vibrates at a superposition of the eigenmodes of the elastic potential energy of the glass, and the tones you hear are the eigenfrequencies $\sqrt{\lambda_i}$.

Finally, when the potential $V$ is quadratic (as is the case for springs), something special happens: the equation $(*)$ is valid not only for small $\delta x$, but for all $\delta x$. This is why the motion of coupled harmonic oscillators isn't chaotic at all: for all initial condition it's a superposition of harmonic modes. In your case, since you have two degrees of freedom $HV$ is $2\times 2$ and has two eigenvectors, which is why your solution has (only) two terms.

---

Postscript for the attentive reader: the above moved very quickly over a lot of ground to try to keep the answer concise; one detail that might reasonably be objectionable is why the non-symmetric $M^{-1}HV$ should have an orthogonal eigenbasis with positive eigenvalues. The right way to think about the eigenspace of $HV$ is in terms of the generalized eigenvalue problem $$HVv_i = \lambda_i Mv_i$$ which in turn can be interpreted as an eigendecomposition using the nonstandard inner product $\langle\rangle_M$. It can be shown that since $M$ and $HV$ are both positive-definite, a full set of generalized eigenvectors (caution: this term can also have an entirely different sense) exists that are orthonormal under the inner product $\langle v_i, v_j\rangle_M = \delta_{ij}$, and have positive eigenvalues. Incidentally, thinking about the configuration of a set of particles as a high-dimensional Riemannian manifold with metric given by the mass matrix is an extremely fruitful idea throughout classical mechanics.