Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Exponential solutions to second order linear system

Posted by peeterjoot on October 20, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]


We’re discussing specific forms to systems of coupled linear differential equations, such as a loop of “spring” connected masses (i.e. atoms interacting with harmonic oscillator potentials) as sketched in fig. 1.1.

Fig 1.1: Three springs loop


Instead of assuming a solution, let’s see how far we can get attacking this problem systematically.

Matrix methods

Suppose that we have a set of N masses constrained to a circle interacting with harmonic potentials. The Lagrangian for such a system (using modulo N indexing) is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} \sum_{k = 0}^{2} m_k \dot{u}_k^2 - \frac{1}{{2}} \sum_{k = 0}^2 \kappa_k \left( u_{k+1} - u_k \right)^2.\end{aligned} \hspace{\stretch{1}}(1.1)

The force equations follow directly from the Euler-Lagrange equations

\begin{aligned}0 = \frac{d{{}}}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_{n, \alpha}}}- \frac{\partial {\mathcal{L}}}{\partial {u_{n, \alpha}}}.\end{aligned} \hspace{\stretch{1}}(1.2)

For the simple three particle system depicted above, this is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_0 \dot{u}_0^2 +\frac{1}{{2}} m_1 \dot{u}_1^2 +\frac{1}{{2}} m_2 \dot{u}_2^2 - \frac{1}{{2}} \kappa_0 \left( u_1 - u_0 \right)^2- \frac{1}{{2}} \kappa_1 \left( u_2 - u_1 \right)^2- \frac{1}{{2}} \kappa_2 \left( u_0 - u_2 \right)^2,\end{aligned} \hspace{\stretch{1}}(1.3)

with equations of motion

\begin{aligned}\begin{aligned}0 &= m_0 \dot{d}{u}_0 + \kappa_0 \left( u_0 - u_1 \right) + \kappa_2 \left( u_0 - u_2 \right) \\ 0 &= m_1 \dot{d}{u}_1 + \kappa_1 \left( u_1 - u_2 \right) + \kappa_0 \left( u_1 - u_0 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa_2 \left( u_2 - u_0 \right) + \kappa_1 \left( u_2 - u_1 \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.4)

Let’s partially non-dimensionalize this. First introduce average mass \bar{{m}} and spring constants \bar{{\kappa}}, and rearrange slightly

\begin{aligned}\begin{aligned}\frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_0 &= -\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_1 \right) - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_2 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_1 &= -\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_2 \right) - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_0 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_2 &= -\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_0 \right) - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_1 \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.5)


\begin{aligned}\tau = \sqrt{\frac{\bar{{k}}}{\bar{{m}}}} t = \Omega t\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}\mathbf{u} = \begin{bmatrix}u_0 \\ u_1 \\ u_2\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}B = \begin{bmatrix}-\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \\ \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &-\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \\ \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} & \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} &-\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(

Our system takes the form

\begin{aligned}\frac{d^2 \mathbf{u}}{d\tau^2} = B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(

We can at least theoretically solve this in a simple fashion if we first convert it to a first order system. We can do that by augmenting our vector of displacements with their first derivatives

\begin{aligned}\mathbf{w} =\begin{bmatrix}\mathbf{u} \\ \frac{d \mathbf{u}}{d\tau} \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(

So that

\begin{aligned}\frac{d \mathbf{w}}{d\tau} =\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \mathbf{w}= A \mathbf{w}.\end{aligned} \hspace{\stretch{1}}(1.0.9)

Now the solution is conceptually trivial

\begin{aligned}\mathbf{w} = e^{A \tau} \mathbf{w}_\circ.\end{aligned} \hspace{\stretch{1}}(1.0.9)

We are however, faced with the task of exponentiating the matrix A. All the powers of this matrix A will be required, but they turn out to be easy to calculate

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^2=\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11a)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^3=\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11b)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^4=\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B^2 & 0 \\ 0 & B^2\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

allowing us to write out the matrix exponential

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}B^k & 0 \\ 0 & B^k\end{bmatrix}+\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Case I: No zero eigenvalues

Provided that B has no zero eigenvalues, we could factor this as

\begin{aligned}\begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}=\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\begin{bmatrix}B^{k + 1/2} & 0 \\ 0 & B^{k + 1/2}\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This initially leads us to believe the following, but we’ll find out that the three springs interaction matrix B does have a zero eigenvalue, and we’ll have to be more careful. If there were any such interaction matrices that did not have such a zero we could simply write

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}\sqrt{B}^{2k} & 0 \\ 0 & \sqrt{B}^{2k}\end{bmatrix}+\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}\sqrt{B}^{2 k + 1} & 0 \\ 0 & \sqrt{B}^{2 k+1} \end{bmatrix}=\cosh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}+ \begin{bmatrix}0 & 1/\sqrt{B} \tau \\ \sqrt{B} \tau & 0\end{bmatrix}\sinh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This is

\begin{aligned}e^{A \tau}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

The solution, written out is

\begin{aligned}\begin{bmatrix}\mathbf{u} \\ \mathbf{u}'\end{bmatrix}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}\begin{bmatrix}\mathbf{u}_\circ \\ \mathbf{u}_\circ'\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

so that

\begin{aligned}\boxed{\mathbf{u} = \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ'.}\end{aligned} \hspace{\stretch{1}}(1.0.11c)

As a check differentiation twice shows that this is in fact the general solution, since we have

\begin{aligned}\mathbf{u}' = \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ + \cosh \sqrt{B} \tau \mathbf{u}_\circ',\end{aligned} \hspace{\stretch{1}}(1.0.11c)


\begin{aligned}\mathbf{u}'' = B \cosh \sqrt{B} \tau \mathbf{u}_\circ + \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ'= B \left( \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ' \right)= B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Observe that this solution is a general solution to second order constant coefficient linear systems of the form we have in eq. 1.5. However, to make it meaningful we do have the additional computational task of performing an eigensystem decomposition of the matrix B. We expect negative eigenvalues that will give us oscillatory solutions (ie: the matrix square roots will have imaginary eigenvalues).

Example: An example diagonalization to try things out


Let’s do that diagonalization for the simplest of the three springs system as an example, with \kappa_j = \bar{{k}} and m_j = \bar{{m}}, so that we have

\begin{aligned}B = \begin{bmatrix}-2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.20)

A orthonormal eigensystem for B is

\begin{aligned}\left\{\mathbf{e}_{-3, 1}, \mathbf{e}_{-3, 2}, \mathbf{e}_{0, 1} \right\}=\left\{\frac{1}{{\sqrt{6}}}\begin{bmatrix} -1 \\ -1 \\ 2 \end{bmatrix},\frac{1}{{\sqrt{2}}}\begin{bmatrix} -1 \\ 1 \\ 0 \end{bmatrix},\frac{1}{{\sqrt{3}}}\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\right\}.\end{aligned} \hspace{\stretch{1}}(1.21)


\begin{aligned}\begin{aligned}U &= \frac{1}{{\sqrt{6}}}\begin{bmatrix} -\sqrt{3} & -1 & \sqrt{2} \\ 0 & 2 & \sqrt{2} \\ \sqrt{3} & -1 & \sqrt{2} \end{bmatrix} \\ D &= \begin{bmatrix}-3 & 0 & 0 \\ 0 & -3 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.22)

We have

\begin{aligned}B = U D U^\text{T},\end{aligned} \hspace{\stretch{1}}(1.23)

We also find that B and its root are intimately related in a surprising way

\begin{aligned}\sqrt{B} = \sqrt{3} iU \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}U^\text{T}=\frac{1}{{\sqrt{3} i}} B.\end{aligned} \hspace{\stretch{1}}(1.24)

We also see, unfortunately that B has a zero eigenvalue, so we can’t compute 1/\sqrt{B}. We’ll have to back and up and start again differently.

Case II: allowing for zero eigenvalues

Now that we realize we have to deal with zero eigenvalues, a different approach suggests itself. Instead of reducing our system using a Hamiltonian transformation to a first order system, let’s utilize that diagonalization directly. Our system is

\begin{aligned}\mathbf{u}'' = B \mathbf{u} = U D U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.25)

where D = [ \lambda_i \delta_{ij} ] and

\begin{aligned}\left( U^{-1} \mathbf{u} \right)'' = D \left( U^{-1} \mathbf{u} \right).\end{aligned} \hspace{\stretch{1}}(1.26)


\begin{aligned}\mathbf{z} = U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.27)

so that our system is just

\begin{aligned}\mathbf{z}'' = D \mathbf{z},\end{aligned} \hspace{\stretch{1}}(1.28)


\begin{aligned}z_i'' = \lambda_i z_i.\end{aligned} \hspace{\stretch{1}}(1.29)

This is N equations, each decoupled and solvable by inspection. Suppose we group the eigenvalues into sets \{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0 \}. Our solution is then

\begin{aligned}\mathbf{z} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p+\left( a_z + b_z \tau \right) \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.30)

Transforming back to lattice coordinates using \mathbf{u} = U \mathbf{z}, we have

\begin{aligned}\mathbf{u} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)U \mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)U \mathbf{e}_p+\left( a_z + b_z \tau \right) U \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.31)

We see that the zero eigenvalues integration terms have no contribution to the lattice coordinates, since U \mathbf{e}_z = \lambda_z \mathbf{e}_z = 0, for all \lambda_z = 0.

If U = [ \mathbf{e}_i ] are a set of not necessarily orthonormal eigenvectors for B, then the vectors \mathbf{f}_i, where \mathbf{e}_i \cdot \mathbf{f}_j = \delta_{ij} are the reciprocal frame vectors. These can be extracted from U^{-1} = [ \mathbf{f}_i ]^\text{T} (i.e., the rows of U^{-1}). Taking dot products between \mathbf{f}_i with \mathbf{u}(0) = \mathbf{u}_\circ and \mathbf{u}'(0) = \mathbf{u}'_\circ, provides us with the unknown coefficients a_n, b_n

\begin{aligned}\mathbf{u}(\tau)= \sum_{ \lambda_n < 0, \lambda_p > 0 }\left( (\mathbf{u}_\circ \cdot \mathbf{f}_n) \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_n}{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( (\mathbf{u}_\circ \cdot \mathbf{f}_p) \cosh \sqrt{\lambda_p} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_p}{\sqrt{-\lambda_p} } \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p.\end{aligned} \hspace{\stretch{1}}(1.32)

Supposing that we constrain ourself to looking at just the oscillatory solutions (i.e. the lattice does not shake itself to pieces), then we have

\begin{aligned}\boxed{\mathbf{u}(\tau)= \sum_{ \lambda_n < 0 }\left( \sum_j \mathbf{e}_{n,j} \mathbf{f}_{n,j}^\text{T} \right)\left( \mathbf{u}_\circ \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' }{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right).}\end{aligned} \hspace{\stretch{1}}(1.33)

Eigenvectors for eigenvalues that were degenerate have been explicitly enumerated here, something previously implied. Observe that the dot products of the form (\mathbf{a} \cdot \mathbf{f}_i) \mathbf{e}_i have been put into projector operator form to group terms more nicely. The solution can be thought of as a weighted projector operator working as a time evolution operator from the initial state.

Example: Our example interaction revisited


Recall that we had an orthonormal basis for the \lambda = -3 eigensubspace for the interaction example of eq. 1.20 again, so \mathbf{e}_{-3,i} = \mathbf{f}_{-3, i}. We can sum \mathbf{e}_{-3,1} \mathbf{e}_{-3, 1}^\text{T} + \mathbf{e}_{-3,2} \mathbf{e}_{-3, 2}^\text{T} to find

\begin{aligned}\mathbf{u}(\tau)= \frac{1}{{3}}\begin{bmatrix}2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2\end{bmatrix}\left( \mathbf{u}_\circ \cos \sqrt{3} \tau + \frac{ \mathbf{u}_\circ' }{\sqrt{3} } \sin \sqrt{3} \tau \right).\end{aligned} \hspace{\stretch{1}}(1.34)

The leading matrix is an orthonormal projector of the initial conditions onto the eigen subspace for \lambda_n = -3. Observe that this is proportional to B itself, scaled by the square of the non-zero eigenvalue of \sqrt{B}. From this we can confirm by inspection that this is a solution to \mathbf{u}'' = B \mathbf{u}, as desired.

Fourier transform methods

Let’s now try another item from our usual toolbox on these sorts of second order systems, the Fourier transform. For a one variable function of time let’s write the transform pair as

\begin{aligned}x(t) = \int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}\tilde{x}(\omega) = \frac{1}{{2\pi}} \int_{-\infty}^{\infty} x(t) e^{i \omega t} dt\end{aligned} \hspace{\stretch{1}}(

One mass harmonic oscillator

The simplest second order system is that of the harmonic oscillator

\begin{aligned}0 = \dot{d}{x}(t) + \omega_\circ^2 x.\end{aligned} \hspace{\stretch{1}}(

Application of the transform gives

\begin{aligned}0 = \left( \frac{d^2}{dt^2} + \omega_\circ^2 \right)\int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega= \int_{-\infty}^{\infty} \left( -\omega^2 + \omega_\circ^2 \right)\tilde{x}(\omega) e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(

We clearly have a constraint that is a function of frequency, but one that has to hold for all time. Let’s transform this constraint to the frequency domain to consider that constraint independent of time.

\begin{aligned}0 =\frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i \omega t}\int_{-\infty}^{\infty} \left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') e^{-i \omega' t} d\omega'=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i (\omega -\omega') t}=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \delta( \omega - \omega' )=\left( -{\omega}^2 + \omega_\circ^2 \right)\tilde{x}(\omega).\end{aligned} \hspace{\stretch{1}}(

How do we make sense of this?

Since \omega is an integration variable, we can’t just mandate that it equals the constant driving frequency \pm \omega_\circ. It’s clear that we require a constraint on the transform \tilde{x}(\omega) as well. As a trial solution, imagine that

\begin{aligned}\tilde{x}(\omega) = \left\{\begin{array}{l l}\tilde{x}_\circ & \quad \mbox{if latex \left\lvert {\omega – \pm \omega_\circ} \right\rvert < \omega_{\text{cutoff}}$} \\ 0 & \quad \mbox{otherwise}\end{array}\right.\end{aligned} \hspace{\stretch{1}}($

This gives us

\begin{aligned}0 = \tilde{x}_\circ\int_{\pm \omega_\circ -\omega_{\text{cutoff}}}^{\pm \omega_\circ +\omega_{\text{cutoff}}}\left( \omega^2 - \omega_\circ^2 \right)e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(

Now it is clear that we can satisfy our constraint only if the interval [\pm \omega_\circ -\omega_{\text{cutoff}}, \pm \omega_\circ + \omega_{\text{cutoff}}] is made infinitesimal. Specifically, we require both a \omega^2 = \omega_\circ^2 constraint and that the transform \tilde{x}(\omega) have a delta function nature. That is

\begin{aligned}\tilde{x}(\omega) = A \delta(\omega - \omega_\circ)+ B \delta(\omega + \omega_\circ).\end{aligned} \hspace{\stretch{1}}(

Substitution back into our transform gives

\begin{aligned}x(t) = A e^{-i \omega_\circ t}+ B e^{i \omega_\circ t}.\end{aligned} \hspace{\stretch{1}}(

We can verify quickly that this satisfies our harmonic equation \dot{d}{x} = -\omega_\circ x.

Two mass harmonic oscillator

Having applied the transform technique to the very simplest second order system, we can now consider the next more complex system, that of two harmonically interacting masses (i.e. two masses on a frictionless spring).


Our system is described by

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_1 \dot{u}_1^2+ \frac{1}{{2}} m_2 \dot{u}_2^2-\frac{1}{{2}} \kappa \left( u_2 - u_1 \right)^2,\end{aligned} \hspace{\stretch{1}}(

and the pair of Euler-Lagrange equations

\begin{aligned}0 = \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_i}}- \frac{\partial {\mathcal{L}}}{\partial {u_i}}.\end{aligned} \hspace{\stretch{1}}(

The equations of motion are

\begin{aligned}\begin{aligned}0 &= m_1 \dot{d}{u}_1 + \kappa \left( u_1 - u_2 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa \left( u_2 - u_1 \right) \end{aligned}\end{aligned} \hspace{\stretch{1}}(


\begin{aligned}\begin{aligned}u_1(t) &= \int_{-\infty}^\infty \tilde{u}_1(\omega) e^{-i \omega t} d\omega \\ u_2(t) &= \int_{-\infty}^\infty \tilde{u}_2(\omega) e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(

Insertion of these transform pairs into our equations of motion produces a pair of simultaneous integral equations to solve

\begin{aligned}\begin{aligned}0 &= \int_{-\infty}^\infty \left( \left( -m_1 \omega^2 + \kappa \right) \tilde{u}_1(\omega) - \kappa \tilde{u}_2(\omega) \right)e^{-i \omega t} d\omega \\ 0 &= \int_{-\infty}^\infty \left( \left( -m_2 \omega^2 + \kappa \right) \tilde{u}_2(\omega) - \kappa \tilde{u}_1(\omega) \right)e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(

As with the single spring case, we can decouple these equations with an inverse transformation operation \int e^{i \omega' t}/2\pi, which gives us (after dropping primes)

\begin{aligned}0 = \begin{bmatrix}\left( -m_1 \omega^2 + \kappa \right) &- \kappa \\ - \kappa &\left( -m_2 \omega^2 + \kappa \right) \end{bmatrix}\begin{bmatrix}\tilde{u}_1(\omega) \\ \tilde{u}_2(\omega)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(

Taking determinants gives us the constraint on the frequency

\begin{aligned}0 = \left( -m_1 \omega^2 + \kappa \right) \left( -m_2 \omega^2 + \kappa \right) - \kappa^2=m_1 m_2 \omega^4-(m_1 + m_2) \omega^2 =\omega^2 \left( m_1 m_2 \omega^2 -\kappa (m_1 + m_2) \right).\end{aligned} \hspace{\stretch{1}}(

Introducing a reduced mass

\begin{aligned}\frac{1}{{\mu}} = \frac{1}{m_1}+\frac{1}{m_2},\end{aligned} \hspace{\stretch{1}}(

the pair of solutions are

\begin{aligned}\begin{aligned}\omega^2 &= 0 \\ \omega^2 &=\frac{\kappa}{\mu}\equiv \omega_\circ^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(

As with the single mass oscillator, we require the functions \tilde{u}{\omega} to also be expressed as delta functions. The frequency constraint and that delta function requirement together can be expressed, for j \in \{0, 1\} as

\begin{aligned}\tilde{u}_j(\omega) = A_{j+} \delta( \omega - \omega_\circ )+ A_{j0} \delta( \omega )+ A_{j-} \delta( \omega + \omega_\circ ).\end{aligned} \hspace{\stretch{1}}(

With a transformation back to time domain, we have functions of the form

\begin{aligned}\begin{aligned}u_1(t) &= A_{j+} e^{ -i \omega_\circ t }+ A_{j0} + A_{j-} e^{ i \omega_\circ t } \\ u_2(t) &= B_{j+} e^{ -i \omega_\circ t }+ B_{j0} + B_{j-} e^{ i \omega_\circ t }.\end{aligned}\end{aligned} \hspace{\stretch{1}}(

Back insertion of these into the equations of motion, we have

\begin{aligned}\begin{aligned}0 &=-m_1 \omega_\circ^2\left( A_{j+} e^{ -i \omega_\circ t } + A_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( A_{j+} - B_{j+} \right) e^{ -i \omega_\circ t } + \left( A_{j-} - B_{j-} \right) e^{ i \omega_\circ t } + A_{j0} - B_{j0} \right) \\ 0 &=-m_2 \omega_\circ^2\left( B_{j+} e^{ -i \omega_\circ t } + B_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( B_{j+} - A_{j+} \right) e^{ -i \omega_\circ t } + \left( B_{j-} - A_{j-} \right) e^{ i \omega_\circ t } + B_{j0} - A_{j0} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(

Equality requires identity for all powers of e^{i \omega_\circ t}, or

\begin{aligned}\begin{aligned}0 &= B_{j0} - A_{j0} \\ 0 &= -m_1 \omega_\circ^2 A_{j+} + \kappa \left( A_{j+} - B_{j+} \right) \\ 0 &= -m_1 \omega_\circ^2 A_{j-} + \kappa \left( A_{j-} - B_{j-} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j+} + \kappa \left( B_{j+} - A_{j+} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j-} + \kappa \left( B_{j-} - A_{j-} \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(

or B_{j0} = A_{j0} and

\begin{aligned}0 =\begin{bmatrix}\kappa -m_1 \omega_\circ^2 & 0 & - \kappa & 0 \\ 0 & \kappa -m_1 \omega_\circ^2 & 0 & - \kappa \\ -\kappa & 0 & \kappa -m_2 \omega_\circ^2 & 0 \\ 0 & -\kappa & 0 & \kappa -m_2 \omega_\circ^2 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(

Observe that

\begin{aligned}\kappa - m_1 \omega_\circ^2=\kappa - m_1 \kappa \left( \frac{1}{{m_1}} + \frac{1}{{m_2}} \right)=\kappa \left( 1 - 1 - \frac{m_1}{m_2} \right)= -\kappa \frac{m_1}{m_2},\end{aligned} \hspace{\stretch{1}}(

(with a similar alternate result). We can rewrite eq. as

\begin{aligned}0 =-\kappa\begin{bmatrix}m_1/m_2 & 0 & 1 & 0 \\ 0 & m_1/m_2 & 0 & 1 \\ 1 & 0 & m_2/m_1 & 0 \\ 0 & 1 & 0 & m_2/m_1 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(

It’s clear that there’s two pairs of linear dependencies here, so the determinant is zero as expected. We can read off the remaining relations. Our undetermined coefficients are given by

\begin{aligned}\begin{aligned}B_{j0} &= A_{j0} \\ m_1 A_{j+} &= -m_2 B_{j+} \\ m_1 A_{j-} &= -m_2 B_{j-} \end{aligned}\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}\boxed{\begin{aligned}u_1(t) &= a+ m_2 b e^{ -i \omega_\circ t }+ m_2 c e^{ i \omega_\circ t } \\ u_2(t) &= a- m_1 b e^{ -i \omega_\circ t }- m_1 c e^{ i \omega_\circ t }.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(

Observe that the constant term is not really of interest, since it represents a constant displacement of both atoms (just a change of coordinates).


\begin{aligned}u_1(t) - u_2(t) = + (m_1 + m_2)b e^{ -i \omega_\circ t }+ (m_1 + m_2)c e^{ i \omega_\circ t },\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}m_1 \dot{d}{u}_1(t) + \kappa (u_1(t) - u_2(t) )=-m_1 m_2 \omega_\circ^2 \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)+(m_1 + m_2) \kappa \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \omega_\circ^2 + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \kappa \frac{m_1 + m_2}{m_1 m_2} + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)= 0.\end{aligned} \hspace{\stretch{1}}(


We’ve seen that we can solve any of these constant coefficient systems exactly using matrix methods, however, these will not be practical for large systems unless we have methods to solve for all the non-zero eigenvalues and their corresponding eigenvectors. With the Fourier transform methods we find that our solutions in the frequency domain is of the form

\begin{aligned}\tilde{u}_j(\omega) = \sum a_{kj} \delta( \omega - \omega_{kj} ),\end{aligned} \hspace{\stretch{1}}(1.63)

or in the time domain

\begin{aligned}u(t) = \sum a_{kj} e^{ - i \omega_{kj} t }.\end{aligned} \hspace{\stretch{1}}(1.64)

We assumed exactly this form of solution in class. The trial solution that we used in class factored out a phase shift from a_{kj} of the form of e^{ i q x_n }, but that doesn’t change the underling form of that assumed solution. We have however found a good justification for the trial solution we utilized.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: