# Coupled oscillators with imaginary frequency?

A typical problem in Lagrangian mechanics is a system of two coupled bodies under small oscillation. The typical way to do these problems is:

• Write the Lagrangian e.g. $L(\theta_1, \theta_2, \dot \theta_1, \dot \theta_2)$.
• Expand all terms to second order in small quantities.
• Use the Euler-Lagrange equations to find the equation of motion.
• Make the ansatz $\theta_1=Ae^{i\omega t}$ and $\theta_2=Be^{i\omega t}$.
• Solve for $\omega$ by requiring the determinate of the relevant matrix goes to zero.

This procedure usually gives you an expression of the form: $$\omega_\pm^2=\alpha \pm \beta$$ for $\alpha \gt \beta$ (both taken as positive) we then get that $\omega_+$ can be either $\pm |\omega_{+}|$ but these do not lead to linearly independent solutions (since $\sin(|\omega_+|)$ is not linearly independent form $-\sin(|\omega_+|)$), and the solution for $\theta_1$ can then be written as: $$\theta_1=A_+ e^{i|\omega_+|t}+A_- e^{i|\omega_-|t}$$ However when $\alpha \lt \beta$ the solution for $\omega_+$ works just as before, but for $\omega_-$ you have $\omega_-=\pm i |\omega_-|$ which gives you two solutions $\exp(|\omega_-|t)$ and $\exp(-|\omega_-|t)$ which are linearly independent. My question is therefore, in this case how do we write down the solution for $\theta_1$ given that we can only have two constants of integration but have three independent solutions?

• @DanielSank Thanks, I have made the change. It does happen with damping it also happens when working in rotating coordinate systems. But I think my method should still apply either way. – Quantum spaghettification Mar 17 '17 at 16:31
• I think you should discuss the cases under which you get imaginary frequency in the post, because those cases inform the boundary conditions of the problem which in turn tell you which solutions to keep. – DanielSank Mar 17 '17 at 16:36
• @DanielSank Are you sure we need to consider boundary conditions? Boundary conditions are used to determine the integration constants - if we are using them to determine which solution is kept and which is not then we effectively have a third integration constant. – Quantum spaghettification Mar 17 '17 at 16:48
• I suspect that either 1) one of the solutions fails to satisfy the boundary conditions, or 2) The two exponentials are not linearly independent. – DanielSank Mar 17 '17 at 16:50
• In a damped system, the exponentially decreasing solution is allowed, but the growing one is not. – DanielSank Mar 17 '17 at 17:00

## 3 Answers

Fixed points are most elegantly treated in a Hamiltonian formalism so this is what I am going to use here (I will comment briefly how does this translate to Lagrangian formalism in the end).

Consider a Hamiltonian $H(\theta_1,\theta_2,p_1,p_2)$ and a fixed point (equilibrium) at $\theta_{10},\theta_{20},p_{10},p_{20}$, which is defined by $\dot{\theta}_1 = \dot{\theta}_2 = \dot{p}_1 = \dot{p}_2 = 0$. Considering Hamilton's equations this can be also stated as $dH=0$ at the equilibrium point in phase space. If we then take an initial condition with a small perturbation from the equilibrium $\theta_{10} + \delta \theta_1,\theta_{20}+ \delta \theta_2,p_{10}+ \delta p_1,p_{20}+ \delta p_2$ we obtain a set of equations for this perturbation of the form $$\frac{d}{dt} \begin{pmatrix} \delta \theta_1 \\ \delta p_1 \\ \delta \theta_2 \\ \delta p_2 \end{pmatrix} = \begin{pmatrix} \frac{\partial^2 H}{\partial p_1 \partial \theta_1} & \frac{\partial^2 H}{\partial p_1^2 } & \frac{\partial^2 H}{\partial p_1 \partial \theta_2} & \frac{\partial^2 H}{\partial p_1 \partial p_2} \\ -\frac{\partial^2 H}{\partial \theta_1^2} & -\frac{\partial^2 H}{\partial \theta_1 \partial p_1 } & -\frac{\partial^2 H}{\partial \theta_1 \partial \theta_2} & -\frac{\partial^2 H}{\partial \theta_1 \partial p_2} \\ \frac{\partial^2 H}{\partial p_2 \partial \theta_1} & \frac{\partial^2 H}{\partial p_2 \partial p_1 } & \frac{\partial^2 H}{\partial p_2 \partial \theta_2} & \frac{\partial^2 H}{\partial p_2^2} \\ -\frac{\partial^2 H}{\partial \theta_2 \partial \theta_1} & -\frac{\partial^2 H}{\partial \theta_2 \partial p_1 } & -\frac{\partial^2 H}{\partial \theta_2^2} & -\frac{\partial^2 H}{\partial \theta_2 \partial p_2} \end{pmatrix} \begin{pmatrix} \delta \theta_1 \\ \delta p_1 \\ \delta \theta_2 \\ \delta p_2 \end{pmatrix}$$ That looks a little bit overwhelming, but you can substitute your favorite Hamiltonian in there and you will see typically only a few of the components of this perturbation matrix are non-zero. The differential equation above is a linear differential equation of the form $\dot{\mathbf{v}} = \mathbf{A} \mathbf{v}$ whose solution you can find in almost every linear algebra book.

The full set of solutions is found by an Ansatz of the form $\mathbf{v}(t) = \mathbf{v}_0 e^{\lambda t}$. From that you get that the solutions have to have $\mathbf{v}_0$ an eigenvector of the matrix $\mathbf{A}$ with a possibly complex eigenvalue $\lambda$. There will be generally four eigenvalues and four corresponding linearly independent eigenvectors of $\mathbf{A}$.

Now, the mathematical theory of classical mechanics tells us that due to the "symplectic structure" of classical mechanics, the eigenvalues will be grouped in something called a loxodromic set. A loxodromic set of numbers in the complex plane is a set of numbers which is symmetric under reflections with respect to both the real and imaginary axes. For instance, this is a loxodromic set of four numbers $$\{\alpha + i\beta, -\alpha + i\beta, \alpha - i \beta, - \alpha - i \beta \}$$ where $\alpha,\beta$ are real numbers. A fixed point with this kind of set of eigenvalues corresponds to a system with a runaway oscillation $e^{(\alpha+i\beta) t}$. This case can happen only in two degrees of freedom and higher, you will not get intuition for it from one degree of freedom! One particular case where this can happen are rotating systems, where the particle is forced to corotate and oscillate by some attractive potential but also spirals out due to the centrifugal force.

One of the other possibilities for the eigenvalues is a "degenerate" loxodromic set of the form $$\{\alpha,-\alpha, i\beta,-i\beta\}$$ with $\alpha,\beta$ again real. These correspond to a system where a perturbation in one direction causes a stable oscillation whereas in the other direction it is unstable and the trajectory starts drifting away. The loxodromic symmetry of the eigenvalues corresponds to the reversibility of the classical-mechanical system. This means that for every motion in one general direction, you can also find the inverse in some sense; for every $e^{\alpha t}$ runaway, there is an $e^{-\alpha t}$ slow approach to the equilibrium (this corresponds to a particle which has just the right energy to reach a potential maximum and is slowly crawling up towards it).

Now for your question: your Ansatz $\theta_1 = A e^{i\omega t}, \theta_2 = B e^{i\omega t}$ corresponds to the more general Ansatz $\mathbf{v}(t) = \mathbf{v}_0 e^{\lambda t}$ in the Hamiltonian formalism and should also give you a loxodromic set of frequencies.

But you should now see that if we have an initial condition at a general position and velocity close to the equilibrium, it will be a linear combination of the basis of all the four linearly independent eigendirections - independently of whether the loxodromic set is degenerate or non-degenerate. (Remember that you have to give four numbers $\theta_1,\theta_2,\dot{\theta}_1, \dot{\theta}_2$ as your initial condition!!)

If your initial condition (and dynamical system) is real, the loxodromy of the set ensures that your initial conditions will actually combine in a real basis of evolutions written for a non-degenerate set as follows $$\frac{1}{2i} (e^{(\alpha + i \beta)t} - e^{(\alpha - i \beta})t), \frac{1}{2} (e^{(\alpha + i \beta)t} + e^{(\alpha - i \beta})t), \frac{1}{2i} (e^{(-\alpha + i \beta)t} - e^{(-\alpha - i \beta})t), \frac{1}{2} (e^{(-\alpha + i \beta)t} + e^{(-\alpha - i \beta)t})$$ and similarly for the degenerate sets. So even though the frequencies can get mixed-complex (note purely imaginary or real) for a system with two or more degrees of freedom, the symmetries of the mechanical system save the day and you still get proper non-complexified evolution.

## Free oscillations

In the case of free small oscillation around a stable equilibrium, the characteristic frequencies of the system are always positive, $\omega^2=\alpha\pm\beta>0$.

The lagrangian of a system with $n$ degrees of freedom and with a stable equilibrium at $q_0=(q_{0,1},\ldots,q_{0,n})$ can be written as $$L=\frac 12 \dot\eta^TM(q_0)\dot\eta+\frac 12 \eta^TV(q_0)\eta,$$ where $M_{ij}=\sum_{a=1}^N m_a\frac{\partial \vec r_a}{\partial q_i}\cdot\frac{\partial \vec r_a}{\partial q_j}$ and $V_{ij}=\frac{\partial^2V}{\partial q_i\partial q_j}$ are $n\times n$ matrices and $\eta,\dot\eta$ are $n$-component column matrices denoting the small displacements from equilibrium. The first term in the lagrangian corresponds to the kinetic energy so it is always positive for any $\dot\eta\neq 0$, hence $M$ is positive definite. The condition for a stable equilibrium implies that the matrix $V(q_0)$ is also positive definite, so the second term in the lagrangian is positive for $\eta\neq 0$.

Equations of motion reads $$M(q_0)\ddot\eta+V(q_0)\eta=0.$$ We look for solutions of the type $\eta=\rho\cos(\omega t+\phi)$ and then we obtain that the characteristic vectors $\rho$ satisfy $$\left( V(q_0)-\omega^2 M(q_0)\right)\rho=0.$$ Multiplying on the left by $\rho^T$ and on the right by $\rho$ we obtain $$\omega^2=\frac{\rho^TV(q_0)\rho}{\rho^TM(q_0)\rho}.$$ Since both $V(q_0)$ and $M(q_0)$ are positive definite matrices, then $\omega^2>0$ for any non-trivial solution $\rho$.

If the equilibrium is neutral, i.e., there are restoring generalized forces in all direction of the configuration space except for some directions where the force is zero, then the matrix $V(q_0)$ is only positive semidefinite. In that case, $\omega^2\geq 0$.

## Damped oscillations

Lagrangian formalism is not suited for dissipative systems since it assumes that the system is described given its kinetic and potential energies. However, for some systems one can suitably change the formalism in order to include dissipative forces. In particular, dissipative forces proportional to the particles velocities can be derived from the so called dissipation function $\mathcal F$. Euler-Lagrange equations then read $$\frac{d}{dt}\frac{\partial L}{\partial \dot \eta_i}-\frac{\partial L}{\partial \eta_i}+\frac{\partial\mathcal F}{\partial \eta_i}=0.$$

The dissipation function is quadratic in the velocities and is related to the rate of dissipated energy which means it is non-negative. For small oscillations around equilibrium can be written as $$\mathcal F=\frac 12 \dot\eta^T F(q_0)\eta,$$ and the equations of motions are $$M(q_0)\ddot\eta+F(q_0)\dot\eta+V(q_0)\eta=0.$$ Unlikely the free case, this system is not easy to solve because in general the matrices $M$, $F$ and $V$ cannot be simultaneously diagonalized. An exception happens when the damping also depends on the particle masses, the matrices can be simultaneously diagonalized and the equations of motion decouple in the normal coordinates, $$\ddot\zeta_i+f_i\dot\zeta_i+\omega_i^2\zeta_i=0,$$ where $f_i\geq 0$, $i=1,\ldots,n$ are the eigenvalues of $F(q_0)$. These equations can be easily solved using complex functions $z=z_0\exp(-i\Omega_i t)$ and at the end taking only its real part. Plugging $z$ into the equation of motion gives $$\Omega_i=\pm\sqrt{\omega_i^2-\frac{f_i^2}{4}}-i\frac{f_i}{2}.$$ If the damping is small, we can neglect the quadratic term in $f_i$ so that the solution has the nice form $$\zeta_i=C_ie^{-if_it/2}\cos(\omega_it+\varphi_i).$$ If $f_i^2/4\geq\omega_i^2$ then the system does not oscillate at all and these regime are called critically damped and overdamped.

• See the comments above, we get negative $\omega^2$ with the presence of damping or when working in a rotating coordinate system. I am not sure the former is easy to deal with using Lagrangian mechanics but the case of rotating coordinate systems is. – Quantum spaghettification May 30 '17 at 15:34
• @Quantumspaghettification All that procedure you are describing is valid for conservative systems. In that case the characteristic frequencies are always positive for the reason I showed. When there is damping things are not so easy since Euler-Lagrange equations is assumed for non-dissipative systems. The lagrangian itself assumes monogenic system (system described by potentials). In case of rotating coordinates, you can use any generalized coordinates you wish but you must write the lagrangian in terms of an inertial system. – Diracology May 30 '17 at 16:40
• @Diracology The point is obviously that the equilibrium is not stable, stability is not discussed anywhere in the OP and, thus, it is not necessary that $\omega^2 > 0$. – Void Jun 4 '17 at 12:08
• @Diracology You can write your Lagrangian any way you like and treat it and its corresponding equations purely mathematically, monogenicity is only required for equivalence with d'Alembert's principle. The transformation into a corotating system allows you to reduce the time-dependence of e.g. a rotating potential in the Lagrangian into a time-independent description, you just have to find the right generalized potential which captures all the rotational effects (it does exist). – Void Jun 4 '17 at 12:09
• @Void Wouldn't stability be a necessary condition in order to study oscillations? At least neutral stability (vanishing force in one direction of the configuration space while the others have restoring force) would be required. In the case of neutral stability, $\omega^2\geq 0$. – Diracology Jun 4 '17 at 14:32

If there are only two independent solutions (the other two solutions are a growing function and have no physical meaning) with the initial conditions respective two independent decisions and the coefficient before the growing solution is zero. If the eigenvalue has a negative real part, then the solution is damped. If the real part of the eigenvalue is zero, then oscillating.