Several coupled oscillators

For the quantum version see quantum harmonic oscillator#Two coupled quantum harmonic oscillators.
Suppose you have $N = 4$ coupled harmonic oscillators with masses $m_{i}$ and spring constants $k_{i j}$ . The equations of motion for the displacements $x_{i} (t)$ can generally be written as:

m_{i} {\ddot{x}}_{i} + \sum_{j} k_{i j} x_{j} = 0,

where the coupling between oscillators is represented by the off-diagonal terms of the stiffness matrix $k_{i j}$ .

The equations of motion can be expressed in matrix form as:

M \ddot{x} + K x = 0,

where $M$ is the mass matrix (diagonal for point masses), $K$ is the stiffness matrix (symmetric for conservative systems), and $x$ is the vector of displacements.

Normal Modes Transformation

To decouple the oscillators, we solve the eigenvalue problem:

K v_{α} = ω_{α}^{2} M v_{α},

where $v_{α}$ are the eigenvectors (normal modes) and $ω_{α}^{2}$ are the eigenvalues (squared angular frequencies of the normal modes). These eigenvectors form a basis in which the motion of the system can be represented.

Using the eigenvectors, we introduce new coordinates $q_{α} (t)$ , called the normal coordinates, defined by:

x (t) = \sum_{α} q_{α} (t) v_{α},

or in matrix form:

x (t) = V q (t),

where $V$ is the matrix of eigenvectors.

Decoupled Equations

Substituting $x (t) = V q (t)$ into the equations of motion and using the orthogonality properties of the eigenvectors, the system transforms into a set of uncoupled equations:

{\ddot{q}}_{α} + ω_{α}^{2} q_{α} = 0.

Each normal coordinate $q_{α}$ behaves like an independent harmonic oscillator with angular frequency $ω_{α}$ . Thus, the normal modes represent a decomposition of the coupled system into $N$ uncoupled harmonic oscillators.

In conclusion, in the transformed coordinates (normal modes), the coupled system indeed behaves as if it consists of 4 independent harmonic oscillators. This equivalence is exact mathematically and illustrates why normal modes are such a powerful concept in analyzing coupled systems.

More detailed approach

Consider the system

\ddot{X} = - M^{- 1} K X

where $M$ and $K$ are symmetric matrices ( $M$ comes from kinetic energy in the Lagrangian and $K$ is the Hessian matrix of a potential function $V$ also in the exact Lagrangian).
To solve it, we construct the associated first order system:

Y = (\begin{matrix} X \\ \dot{X} \end{matrix})

and then

\dot{Y} = (\begin{array}{cc} 0 & I \\ - M^{- 1} K & 0 \end{array}) \cdot Y

\dot{Y} = A \cdot Y

As usual, we use what is said in system of first order ODEs, that is, we find the stationary states of this first-order system. Observe that:

$λ$ is eigenvalue of $A$ if and only if $λ^{2}$ is eigenvalue of $- M^{- 1} K$ . This is so because the block matrix determinant formula

d e t (\begin{array}{cc} H & B \\ C & D \end{array}) = d e t (H - B D^{- 1} C) \cdot d e t (D)

applied to $λ I - A$ .

ξ = (\begin{matrix} ψ \\ η \end{matrix})

is an eigenvector of $A$ associated to $λ$ then $ψ$ is an eigenvector of $- M^{- 1} K$ associated to $λ^{2}$ . This is so because

A ξ = (\begin{array}{cc} 0 & I \\ - M^{- 1} K & 0 \end{array}) \cdot (\begin{matrix} ψ \\ η \end{matrix}) = (\begin{matrix} η \\ - M^{- 1} K ψ \end{matrix}) = (\begin{matrix} λ ψ \\ λ η \end{matrix})

So,

A^{2} ξ = (\begin{array}{cc} - M^{- 1} K & 0 \\ 0 & - M^{- 1} K \end{array}) \cdot (\begin{matrix} ψ \\ η \end{matrix}) = (\begin{matrix} - M^{- 1} K ψ \\ - M^{- 1} K η \end{matrix}) = (\begin{matrix} λ^{2} ψ \\ λ^{2} η \end{matrix})

And therefore

- M^{- 1} K ψ = λ^{2} ψ

Reciprocally, if $(δ, ψ)$ is a pair such that $- M^{- 1} K ψ = δ ψ$ (of course $δ$ is real and negative) then is easy to check that we have two eigenvalues with two eigenvectors for

(\begin{array}{cc} 0 & I \\ - M^{- 1} K & 0 \end{array}),

the pair

(\sqrt{δ}, (\begin{matrix} ψ \\ \sqrt{δ} ψ \end{matrix}))

and the pair

(- \sqrt{δ}, (\begin{matrix} ψ \\ - \sqrt{δ} ψ \end{matrix}) .

All the eigenvalues of $- M^{- 1} K$ are reals, since is a symmetric matrix and they always can be diagonalized. Moreover, they are all negative numbers, because this matrix reflect a stable equilibrium point, i.e., from it potential energy increases in whatever direction we take. So we conclude:

All the eigenvalues are pure imaginary ones. There is no $e^{α t}$ terms.
For every eigenvalue of $A$ it can be chosen an eigenvector

ψ = (\begin{matrix} ξ \\ η \end{matrix})

such that $ξ$ is real (because is an eigenvector of $- M^{- 1} K$ and their eigenvalues are all real). So

I m (ψ) = (\begin{matrix} 0 \\ I m (η) \end{matrix}) .

So for every eigenvalue-eigenvector of $- M^{- 1} K$ , $(λ_{n}, ξ)$ , we have two solutions:

Ψ^{n} = e^{i ω_{n} t} ξ^{n}

{\tilde{Ψ}}^{n} = e^{- i ω_{n} t} ξ^{n}

where $ω_{n} = \frac{\sqrt{λ_{n}}}{i} \in R$ (remember $λ_{n}$ is negative). It may seem that both solutions have the same initial value, which is a contradiction, but they have different initial velocities!
Moreover, observe that

({\tilde{Ψ}}^{n})^{*} = Ψ^{n},

because we can choose $ξ^{n}$ to be real.

So if we are looking for real solutions we have that

Ψ (t) = \sum_{n} C_{n} Ψ^{n} + {\tilde{C}}_{n} {\tilde{Ψ}}^{n}

(Ψ (t))^{*} = \sum_{n} C_{n}^{*} (Ψ^{n})^{*} + ({\tilde{C}}_{n})^{*} ({\tilde{Ψ}}^{n})^{*}

must be equal, and since the solutions form a vector space we conclude that $C_{n}^{*} = {\tilde{C}}_{n}$

and so

Ψ (t) = \sum_{n} C_{n} Ψ^{n} + C_{n}^{*} (Ψ^{n})^{*} = \sum_{n} 2 R e [C_{n} Ψ^{n}] = \sum_{n} A_{n} c o s (ω_{n} t + ϕ_{n}) ξ^{n}

where $A_{n}$ and $ϕ_{n}$ , determined by initial conditions, are reals and have absorbed the constant to leave $ξ^{n}$ as a real vector.

These fundamental solutions are called normal modes. The constants $ω_{k}$ are called fundamental frequencies.

Continuous limit

Source: ChatGPT.
We consider $N$ oscillators on a 1D lattice with positions $q_{i} (t)$ , where $i = 1, 2, \dots, N$ . The oscillators are coupled to their nearest neighbors, and the Newtonian equation of motion for the $i$ -th oscillator is:

m {\ddot{q}}_{i} = - k (q_{i} - q_{i - 1}) - k (q_{i} - q_{i + 1}) - m ω^{2} q_{i},

where:

$k$ is the spring constant between adjacent oscillators,
$m ω^{2}$ represents a restoring force on each oscillator (like a harmonic potential).

Rewriting, we have:

m {\ddot{q}}_{i} = k (q_{i - 1} - 2 q_{i} + q_{i + 1}) - m ω^{2} q_{i} .

In the continuum limit, the positions $q_{i} (t)$ are replaced by a continuous field $ϕ (x, t)$ , where $x = i a$ and $a$ is the lattice spacing. We make the following approximations:

The displacement difference $q_{i + 1} - q_{i}$ becomes $a \partial ϕ / \partial x$ ,
The discrete second difference $q_{i - 1} - 2 q_{i} + q_{i + 1}$ becomes $a^{2} \partial^{2} ϕ / \partial x^{2}$ .

Substituting into the equation of motion:

m \frac{\partial^{2} ϕ}{\partial t^{2}} = k a^{2} \frac{\partial^{2} ϕ}{\partial x^{2}} - m ω^{2} ϕ .

Divide through by $m$ and define:

c^{2} = \frac{k a^{2}}{m}, m^{2} = \frac{ω^{2}}{c^{2}} .

This gives the equation:

\frac{\partial^{2} ϕ}{\partial t^{2}} = c^{2} \frac{\partial^{2} ϕ}{\partial x^{2}} - m^{2} c^{4} ϕ .

Rewriting the equation in standard form:

\frac{\partial^{2} ϕ}{\partial t^{2}} - c^{2} \frac{\partial^{2} ϕ}{\partial x^{2}} + m^{2} c^{4} ϕ = 0,

which is nothing but Klein--Gordon equation.

Observe that if we had considered a null restoring force we would have obtained the wave equation.

Old stuff, to be incorporated above, maybe.

See the source in anotacioneslatex.tex.

A chain of a huge amount of coupled oscillators: passing to continuous

\begin

\includegraphics[width=15cm]

\end

We are going to study $N$ masses connected to a string with tension. We leave open ends, for the moment.

Let $Ψ_{j} (t)$ be the displacement from the equilibrium of every mass $j$ . Let the mass be $m$ and the distance between the mass be $d$ . Sometimes we will take $x = j d$ and write

Ψ (x, t) = Ψ_{x / d} (t)

Analyzing every mass individually we obtain

{\ddot{Ψ}}_{1} = - \frac{2 T}{m d} Ψ_{1} + \frac{T}{m d} Ψ_{2}

{\ddot{Ψ}}_{2} = \frac{T}{m d} Ψ_{1} - \frac{2 T}{m d} Ψ_{2} + \frac{T}{m d} Ψ_{3}

\dots

In general:

\ddot{Ψ} = \frac{1}{m d} (\begin{array}{cccc} - T & T & 0 & 0 \\ T & - 2 T & T & 0 \\ . . . & . . . & . . . & . . . \\ 0 & 0 & T & - T \end{array}) \cdot Ψ

or in short

\ddot{Ψ} = - M^{- 1} K \cdot Ψ

Observe that the matrix $- M^{- 1} K$ is symmetric and all their eigenvalues are negatives, since it comes from a stable equilibrium point (see section \ref{coupledoscillators}). Moreover, from this section and the previous ones we know that a basis for solutions of this differential equation is

A c o s (ω t + ϕ) ψ

where $- ω^{2}$ is an eigenvalue of $- M^{- 1} K$ ,

$ψ$ is an eigenvector for this eigenvalue, and $A$ and $ϕ$ are real constants determined by the initial conditions: the $A$ 's by the initial positions and the $ϕ$ 's by the initial velocities.

\bigbreak

\textbf{When $N$ is really big} is very difficult to find the eigenvalues by linear algebra, so we can proceed in a different way.

Also, the vector $ψ$ is a very large vector, since the number of masses $N$ will tend to infinity, and we want to study the shape of their components. In fact, when $N$ goes to infinity $ψ$ will be a function of position, giving the displacement of the masses in the initial time.

As usual, we will follow with the complex solution for $Ψ$ and take the real part in the final step. We think in solutions like

Ψ_{j} (t) = e^{i ω t} ξ (j)

where $ξ$ is a complex eigenvector of infinite length (keep an eye: in section \ref{coupledoscillators} we saw that with finite masses we can force $ξ$ to be real, but with infinite masses, a priori, we cannot).

To find infinite length eigenvectors we use the trick of \ref{symmetrytrick}. Since this infinite chain of oscillating masses has translation symmetry, we check that

- M^{- 1} K S = S (- M^{- 1} K)

where

- M^{- 1} K = (\begin{array}{cccccc} ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \\ \dots & - 2 B & C & 0 & 0 & \dots \\ \dots & C & - 2 B & C & 0 & \dots \\ \dots & 0 & C & - 2 B & C & \dots \\ \dots & 0 & 0 & C & - 2 B & \dots \\ ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{array})

and

S = (\begin{array}{cccccc} ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \\ \dots & 0 & 1 & 0 & 0 & \dots \\ \dots & 0 & 0 & 1 & 0 & \dots \\ \dots & 0 & 0 & 0 & 1 & \dots \\ \dots & 0 & 0 & 0 & 0 & \dots \\ ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{array})

What are the eigenvectors for $S$ ? We can take any $β \in \RR$ , and then produce an eigenvector $ξ$ . Observe that

S ξ (j) = ξ (j + 1) = β ξ (j)

so if we take $ξ (0) = 1$ we conclude $ξ (j) = β^{j}$ .

Moreover, since $j \in \ZZ$ we conclude that $| β | = 1$ and therefore must be $α \in [- π, π]$ such that $β = e^{i α}$ . In conclussion, for every $α \in [- π, π]$ we have the eigenpair $(e^{i α}, ξ^{α})$ where $ξ^{α} (j) = e^{i α j}$ , respect to $S$ .

But, what is the eigenvalue respect to $- M^{- 1} K$ ?

(- M^{- 1} K ξ^{α}) (j) = - 2 B e^{i α j} + C e^{i α (j - 1)} + C e^{i α (j + 1)} = (- 2 B + C e^{- i α} + C e^{i α}) e^{i α j}

= (- 2 B + C (e^{- i α} + e^{i α})) e^{i α j} = (2 C \cos (α) - 2 B) e^{i α j}

So $ξ^{α}$ is an eigenvector of $- M^{- 1} K$ associated to the eigenvalue $2 C \cos (α) - 2 B$ , which must be negative (probably $C = B$ , and we would get that).

Now, think that the pair $(2 C \cos (α) - 2 B, ξ^{α})$ produces two solutions for the associated first order system (and then ``truncated''):

Ψ^{α} (t) = e^{i ω_{α} t} ξ^{α}

{\tilde{Ψ}}^{α} (t) = e^{- i ω_{α} t} ξ^{α}

where $ω_{α} = \frac{\sqrt{2 C \cos (α) - 2 B}}{i} \in \RR$ . In the particular case where $B = C = \frac{T}{m d}$ ,

ω_{α} = 2 \sqrt{\frac{T}{m d}} | \sin \frac{α}{2} |

which is known as \textit

\bigbreak

A few remarks:

\begin

\item The general solution is a linear combination of $Ψ^{α}$ 's and ${\tilde{Ψ}}^{α}$ 's. But since $α$ has no restriction, beyond $α \in [- π, π]$ , because we don't have boundary conditions, $α$ varies continuously. So, instead of a linear combination we get

Ψ = \int_{- π}^{π} C (α) Ψ^{α} + \tilde{C} (α) {\tilde{Ψ}}^{α} d α

Observe that this has been developed with a trick: we forgot the initial and final rows of the matrix $- M^{- 1} K$ since we are dealing with a big $N$ . This is the reason why we are finding more than $N$ eigenvalues (infinite, actually): since we are no taking into account the initial and final restriction, it is totally as if we had infinite masses.

\item Let's look for the general \textbf{real} solution. First, observe that:

(Ψ^{α})^{*} = e^{- i ω_{α} t} (ξ^{α})^{*} = e^{- i ω_{α} t} ξ^{- α} = {\tilde{Ψ}}^{- α}

({\tilde{Ψ}}^{α})^{*} = e^{i ω_{α} t} (ξ^{α})^{*} = e^{i ω_{α} t} ξ^{- α} = Ψ^{- α}

Then

Ψ^{*} = \int_{- π}^{π} C^{*} (α) (Ψ^{α})^{*} + {\tilde{C}}^{*} (α) ({\tilde{Ψ}}^{α})^{*} d α = \int_{- π}^{π} C^{*} (α) {\tilde{Ψ}}^{- α} + {\tilde{C}}^{*} (α) Ψ^{- α} d α =

= \int_{- π}^{π} C^{*} (- α) {\tilde{Ψ}}^{α} + {\tilde{C}}^{*} (- α) Ψ^{α} d α

Since it must be $Ψ = Ψ^{*}$ , and so :

C^{*} (- α) = \tilde{C} (α)

(I don't know how to prove this!!, but the idea is that $Ψ^{α}$ and ${\tilde{Ψ}}^{α}$ constitute a basis for the vector space of solutions) and therefore:

Ψ = \int_{- π}^{π} C (α) Ψ^{α} + \tilde{C} (α) {\tilde{Ψ}}^{α} d α = \int_{0}^{π} 2 R e [C (α) Ψ^{α} + \tilde{C} (α) {\tilde{Ψ}}^{α}] d α

\item The parameter $α$ is the seed for the wave number. Let $d$ be the distance between the masses, we can write $x = d \cdot j$ to be the position in the horizontal direction. If we try to write everything int terms of $x$ instead $j$ is useful to choose $κ = \frac{α}{d}$ for the trial solution for the eigenvector. The normal modes would be:

Ψ^{κ} (x, t) = a \cdot c o s (ω_{κ} t + κ x + ϕ)

\item Dispersion relation is a name for the functional relation between wave number $κ$ and frecuency $ω_{κ}$ . In this case is

ω_{κ} = 2 \sqrt{\frac{T}{m d}} s i n (\frac{κ d}{2})

\item The idea for choosing $ξ (j) = A e^{i α j}$ can be seen from other point of view.

The eigenvectors satisfy the relation

- m d ω^{2} ξ (j) = - 2 T ξ (j) + T ξ (j - 1) + T ξ (j + 1)

where we forget the first and the last relation because we assume $N$ very very large

One can observe that if $d$ is very small and taking into account the name change $x = j d$ , equation \ref{eigenvectorrelation} can be interpreted (approximately) as:

- m ω^{2} ξ (x) = T [\frac{ξ (x + d) - ξ (x)}{d} - \frac{ξ (x) - ξ (x - d)}{d}] \approx T [ξ^{'} (x + d / 2)] - ξ^{'} (x - d / 2)]

and so

- m ω^{2} ξ (x) \approx T d ξ^{″} (x)

But since $m = d μ$ , where $μ$ is the density and is supposed to be constant

- \frac{μ ω^{2}}{T} ξ (x) \approx ξ^{″} (x)

We have got an equation quite similar but in the variable $x$ or $j$ if you prefer, so is natural to choice $ξ (j) = A e^{i α j}$ .

\end

\section

We have studied an infinite system. If we want to come back to a finite system we impose boundary conditions, so we reduce the number of fundamental solutions that are allowed.

Particular case: $N = 4$ . We have:

\ddot{Ψ} = \frac{1}{m d} (\begin{array}{cccc} - T & T & 0 & 0 \\ T & - 2 T & T & 0 \\ 0 & T & - 2 T & T \\ 0 & 0 & T & - T \end{array}) \cdot Ψ

We solve using Mathematica and find, evidently, four eigenvalues:

- ω_{1}^{2} = \frac{(- 2 - \sqrt{2}) T}{d m}, - ω_{2}^{2} = - \frac{2 T}{d m}, - ω_{3}^{2} = \frac{(- 2 + \sqrt{2}) T}{d m}, - ω_{4}^{2} = 0

and four eigenvectors

ξ_{1} = (- 1, 1 + \sqrt{2}, - 1 - \sqrt{2}, 1), ξ_{2} = (1, - 1, - 1, 1),

ξ_{3} = (- 1, 1 - \sqrt{2}, - 1 + \sqrt{2}, 1), ξ_{4} = (1, 1, 1, 1)

With pictures ( $T = 2; m = 1; d = 0.5$ ):

\begin{tabular}

\includegraphics[width=60mm]{imagenes/sol1.png} &

\includegraphics[width=60mm]{imagenes/sol2.png} \

\includegraphics[width=60mm]{imagenes/sol3.png} &

\includegraphics[width=60mm]{imagenes/sol4.png} \

\end

Could we obtain the same with the big-N-technique? Let's see. We start with infinite vibrating masses with displacement $Ψ_{j} (t)$ , and so infinite fundamental frequencies

ω_{α} = 2 \sqrt{\frac{T}{m d}} \sin (α / 2)

one for every $α \in [0, π]$ .

If we want to restrict to this particular case, we have to impose some conditions: \textbf{boundary conditions}. For example, if we force $Ψ_{0} = Ψ_{1}$ and $Ψ_{4} = Ψ_{5}$ we are removing the effect of masses 0 and 5 over the masses ranging from 1 to 4 (as desired). Let's study this two conditions: we are going to impose $Ψ_{0} = Ψ_{1}$ and $Ψ_{4} = Ψ_{5}$ and watch what $α$ 's survive.

e^{i ω_{α} t} [C (α) + C (- α)] = e^{i ω_{α} t} [C (α) e^{i α} + C (- α) e^{- i α}]

C (α) + C (- α) = C (α) e^{i α} + C (- α) e^{- i α}

and together with

e^{i ω_{α} t} [C (α) e^{i α 4} + C (- α) e^{- i α 4}] = e^{i ω_{α} t} [C (α) e^{i α 5} + C (- α) e^{- i α 5}]

we arrive to

α = \frac{π n}{4}

Let's study this values one by one:

\begin

\item n=0. Then $α = 0$ , $ω_{α}^{2} = 0$ and $ξ = (1, 1, 1, 1)$ .

Correspond to the case $ω_{4}$ following Mathematica.

\item $n = 1$ . Then $α = \frac{π}{4}$ , and $ω_{α}^{2} = \frac{T}{m d} (2 - \sqrt{2})$ . We are in case 3 from Mathematica. The eigenvector obtained is $ξ = (e^{i \frac{π}{4}}, e^{i \frac{π}{2}}, e^{i \frac{3 π}{4}}, e^{i π})$ , which is complex and do not coincide with the one obtained by Mathematica. But since $- α$ produces the same eigenvalue, we can mix the previous eigenvector with $(e^{- i \frac{π}{4}}, e^{- i \frac{π}{2}}, e^{- i \frac{3 π}{4}}, e^{- i π})$ to make new eigenvectors for this eigenvalue. In fact is a (long) computation to check that the linear combination of this complex vectors that produces the same that the real above eigenvector is

{(x \to - \frac{1}{2} + i (- \frac{1}{2} + \frac{1}{\sqrt{2}}), y \to - \frac{1}{2} + i (\frac{1}{2} - \frac{1}{\sqrt{2}})}

I found it by using Mathematica.

\item $n = 2$ . Then $α = \frac{π}{2}$ , and $ω_{α}^{2} = 2 \frac{T}{m d}$ . This correspond with case $ω_{2}$ above. The eigenvector is not the same, but the phenomenom is the same that for $n = 1$ : we can recover by mixing both eigenvectors corresponding to $α$ and $- α$ .

\item $n = 3$ . Then $α = \frac{3 π}{4}$ , and $ω_{α}^{2} = \frac{T}{m d} (2 + \sqrt{2})$ . Idem.

\item $n = 4$ . Then $α = π$ . Condition $Ψ_{0}^{π} = Ψ_{1}^{π}$ implies $C (π) + C (- π) = 0$ . Moreover, since $e^{i π j} = e^{- i π j}$ , we deduce that $C (π) Ψ^{π} + C (- π) Ψ^{- π} = 0$ , so we can ignore this normal mode.

\item $n = 5$ . Then $α = π + \frac{π}{4}$ . As we reasoned above, this would be the same case as $π + \frac{π}{4} - 2 π = - \frac{3 π}{4}$ which, in fact, is paired with $α = \frac{3 π}{4}$ (case $n = 3$ )

\item All the remaining cases are repeated

\end