# Fourth-order dynamics of the damped harmonic oscillator

John W. Sanders
Department of Mechanical Engineering
California State University, Fullerton
Fullerton, CA 92831, United States

ORCiD: 0000-0003-3059-3815<EMAIL_ADDRESS>
(Received: date / Accepted: date)

## Abstract

It is shown that the classical damped harmonic oscillator belongs to the
family of fourth-order Pais-Uhlenbeck oscillators. It follows that the
solutions to the damped harmonic oscillator equation make the Pais-Uhlenbeck
action stationary. Two systematic approaches are given for deriving the Pais-
Uhlenbeck action from the damped harmonic oscillator equation, and it may be
possible to use these methods to identify stationary action principles for
other dissipative systems which do not conform to Hamilton’s principle. It is
also shown that for every damped harmonic oscillator $x$, there exists a two-
parameter family of dual oscillators $y$ satisfying the Pais-Uhlenbeck
equation. The damped harmonic oscillator and any of its duals can be
interpreted as a system of two coupled oscillators with atypical spring
stiffnesses (not necessarily positive and real-valued). For overdamped
systems, the resulting coupled oscillators should be physically achievable and
may have engineering applications. Finally, a new physical interpretation is
given for the optimal damping ratio $\zeta=1/\sqrt{2}$ in control theory.

##  I Introduction and literature review

Ever since Newton Newton (1687) first laid the foundation for mechanics as a
formal system, the development of analytical mechanics has progressed via the
principle of virtual work and its generalization, d’Alembert’s principle,
d’Alembert (1743) reaching its pinnacle in the principle of stationary action.
Lagrange (1811); Hamilton (1834, 1835) During any time interval
$t\in(t_{1},t_{2})$, the actual motion of a conservative mechanical system is
given by a critical point of the action functional

$\mathcal{S}=\int_{t_{1}}^{t_{2}}Ldt,$ (1)

where the Lagrangian $L=T-V$ is the difference between the system’s total
kinetic energy $T$ and total potential energy $V$. Put simply, the action is
stationary for the actual motion of a conservative system:

$\delta\mathcal{S}=0.$ (2)

While sometimes attributed to Hamilton, Hamilton (1834, 1835) the principle of
stationary action was already well known by Lagrange Lagrange (1811) and
contemporaries. Lanczos (1949); Gray and Taylor (2007) It has been
demonstrated that the action always achieves either a local minimum or a
saddle point, but never a local maximum. Jacobi (1837); Whittaker (1904); Gray
and Taylor (2007)

When dissipative forces are present, Hamilton’s principle (2) may be exchanged
for the d’Alembert-Lagrange principle Lagrange (1811); Kane  _et al._ (2000);
Baddour (2007); Flannery (2011)

$\delta\mathcal{S}+\int_{t_{1}}^{t_{2}}Q_{i}\delta{q}_{i}dt=0,$ (3)

where the $q_{i}$ are the generalized coordinates and the $Q_{i}$ are the
generalized dissipative forces. Like d’Alembert’s principle, d’Alembert (1743)
Hamilton’s principle (2) and the d’Alembert-Lagrange principle (3) can be used
as the basis for numerical solution methods. Gray, Karl, and Novikov (1996);
Kim, Dargush, and Ju (2013) Indeed, much modern work has been devoted to the
development of accurate and efficient variational time integrators based on
these principles. Kane  _et al._ (2000); Marsden and West (2001); Lew  _et
al._ (2003, 2004a, 2004b); Kale and Lew (2007); Ober-Blöbaum, Junge, and
Marsden (2011); Kraus and Maj (2015); Hall and Leok (2015); Lew and Mata
(2016); Schmitt and Leok (2018); Man, Gao, and Zhong (2020); Limebeer, Ober-
Blöbaum, and Farshi (2020); Kraus and Tyranowski (2021) Nevertheless, the fact
remains that the d’Alembert-Lagrange principle is not a true variational
principle, because in general the virtual work terms $Q_{i}\delta{q}_{i}$ are
not the exact variations of a work function. Lanczos (1949) As a result, those
terms must be inserted into the equation manually _after_ the variation
operator $\delta$ has been applied to the action. The _ad hoc_ insertion of
the dissipative terms into an otherwise variational principle is extremely
dissatisfying on physical, mathematical, and aesthetic grounds, and over the
last two centuries many attempts have been made to extend the Lagrangian
framework and the stationary action principle to dissipative systems.

Figure 1: The damped harmonic oscillator.

One of the most ubiquitous models of a dissipative system is given by the
damped harmonic oscillator, as shown in Figure 1. It consists of a mass $m$
attached to a rigid support by a linear spring of stiffness $k$ and a linear
damper of damping coefficient $c$, with $x(t)$ the generalized coordinate
along the single degree of freedom. The governing equation of motion is given
by

$m\ddot{x}+c\dot{x}+kx=0,$ (4)

or equivalently,

$\ddot{x}+2\omega\zeta\dot{x}+\omega^{2}x=0,$ (5)

where we have introduced the natural frequency $\omega=\sqrt{k/m}$ and damping
ratio $\zeta=c/2m\omega$. Henceforth, we shall regard $\omega>0$ and
$\zeta\geq 0$ as real-valued model parameters. Depending on the damping ratio,
the damped harmonic oscillator exhibits different behaviors, including the
absence of damping ($\zeta=0$), underdamping ($0<\zeta<1$), critical damping
($\zeta=1$), and overdamping ($\zeta>1$). While the damped harmonic oscillator
does not conform to Hamilton’s principle (2), there are known variational
principles that reproduce the governing equation (5). As discussed in detail
by Limebeer _et al._ , Limebeer, Ober-Blöbaum, and Farshi (2020) one class of
stationary action principles is based on non-standard Lagrangians. For
example, the dual system method of Bateman, Morse, Feshbach, and Tikochinsky
Bateman (1931); Morse and Feshbach (1953); Feshbach and Tikochinsky (1977)
employs a Lagrangian of the form

$L=\dot{x}\dot{y}+\omega\zeta(x\dot{y}-y\dot{x})-\omega^{2}xy,$ (6)

where $x$ and $y$ are both to be varied. The Euler-Lagrange equations yield
two oscillators: the damped harmonic oscillator (5) and a dual anti-damped
oscillator

$\ddot{y}-2\omega\zeta\dot{y}+\omega^{2}y=0.$ (7)

Another example is given by the Kanai-Caldirola oscillator, Caldirola (1941);
Kanai (1948) which uses a Lagrangian of the form

$L=\frac{1}{2}e^{2\omega\zeta t}(\dot{x}^{2}-\omega^{2}x^{2}),$ (8)

with Euler-Lagrange equation

$e^{2\omega\zeta t}(\ddot{x}+2\omega\zeta\dot{x}+\omega^{2}x)=0.$ (9)

A final approach that will be mentioned here is the Caldeira-Leggett model,
Caldeira and Leggett (1983); Zwanzig (1973); Ford and Kac (1987) which treats
a system as embedded within a heat bath, so that the combination of the system
and its surroundings is conservative. Limebeer, Ober-Blöbaum, and Farshi
(2020) An example is provided by Fukagawa and Fujitani, Fukagawa and Fujitani
(2012) who used Hamilton’s principle coupled with a nonholonomic constraint on
the entropy to incorporate the damping force into the damped harmonic
oscillator equation as well as the viscous terms into the momentum balance
equation for a viscous fluid, one special case being the Navier-Stokes
equations. Stokes (1845)

Quite independent of the damped harmonic oscillator is the Pais-Uhlenbeck
oscillator, Pais and Uhlenbeck (1950) which is used as a representative “toy
model” to study higher-derivative theories for quantum gravity. Mostafazadeh
(2010); Baleanu  _et al._ (2012) In general, the Pais-Uhlenbeck oscillator of
order $2n$ is given by

$\prod_{j=1}^{n}\left(\frac{d^{2}}{dt^{2}}+\omega_{j}^{2}\right)x=0,$ (10)

where the $\omega_{j}$ are real-valued parameters. For $n=2$, this yields the
fourth-order equation

$\ddddot{x}+(\omega_{1}^{2}+\omega_{2}^{2})\ddot{x}+\omega_{1}^{2}\omega_{2}^{2}x=0.$
(11)

In the present work, we will establish that the damped harmonic oscillator
belongs to the family of fourth-order Pais-Uhlenbeck oscillators. To date, it
seems that the connection between the damped harmonic oscillator and the Pais-
Uhlenbeck oscillator has gone unnoticed. This connection leads to two
interesting conclusions: (i) that the damped harmonic oscillator is a critical
point of the Pais-Uhlenbeck action, and (ii) that for every damped harmonic
oscillator $x$, there exists a two-parameter family of coupled dual
oscillators $y$.

The remainder of this paper is organized as follows. In Section II, we show
that every solution to the damped harmonic oscillator equation (5) is also a
solution to a fourth-order equation of the Pais-Uhlenbeck type. In Section
III, we investigate the two-parameter family of coupled-oscillator systems,
leading to new interpretations for the various damping regimes. In Section IV,
we examine the special case of $\zeta=1/\sqrt{2}$, which leads to new insights
into this “optimal control value.” In Section V, we discuss the corresponding
variational formulation based on the Pais-Uhlenbeck action. We also present
two systematic methods for arriving at the Pais-Uhlenbeck action from the
original damped harmonic oscillator equation, which may help identify
variational principles for other dissipative systems. Finally, in Section VI,
we conclude with a brief summary of the present results.

##  II The damped harmonic oscillator as a Pais-Uhlenbeck oscillator

We begin by observing that every solution to (5) is also a solution to the
following fourth-order equation:

$\ddddot{x}+4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\ddot{x}+\omega^{4}x=0,$
(12)

which is a kind of Pais-Uhlenbeck equation in which the coefficient of
$\ddot{x}$ may be positive, negative, or zero, depending on the damping ratio.
This can be demonstrated as follows. Differentiating (5) twice, we obtain

$\ddddot{x}+2\omega\zeta\dddot{x}+\omega^{2}\ddot{x}=0.$ (13)

Substituting

$\dddot{x}=-2\omega\zeta\ddot{x}-\omega^{2}\dot{x}$ (14)

and

$\dot{x}=-\frac{1}{2\omega\zeta}(\ddot{x}+\omega^{2}x)$ (15)

into (13) yields (12). We conclude that the damped harmonic oscillator belongs
to the family of Pais-Uhlenbeck oscillators. Indeed, with the same initial
conditions

$\displaystyle x(0)$ $\displaystyle=x_{0},$ (16) $\displaystyle\dot{x}(0)$
$\displaystyle=v_{0},$ (17) $\displaystyle\ddot{x}(0)$
$\displaystyle=a_{0}\equiv-2\omega\zeta v_{0}-\omega^{2}x_{0},$ (18)
$\displaystyle\dddot{x}(0)$ $\displaystyle=j_{0}\equiv-2\omega\zeta
a_{0}-\omega^{2}v_{0},$ (19)

the uniqueness theorem guarantees that the Pais-Uhlenbeck oscillator (12) will
yield the same solution as the damped harmonic oscillator (5).

Using a trial solution of the form

$x(t)=Ae^{\lambda t},$ (20)

we solve the characteristic polynomial

$\lambda^{4}+4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\lambda^{2}+\omega^{4}=0,$
(21)

yielding four eigenvalues:

$\displaystyle\lambda_{1,2}=$
$\displaystyle\pm\omega\sqrt{-2\left(\frac{1}{2}-\zeta^{2}\right)+\sqrt{4\left(\frac{1}{2}-\zeta^{2}\right)^{2}-1}}=\pm\omega(\zeta+\sqrt{\zeta^{2}-1}),$
(22) $\displaystyle\lambda_{3,4}=$
$\displaystyle\pm\omega\sqrt{-2\left(\frac{1}{2}-\zeta^{2}\right)-\sqrt{4\left(\frac{1}{2}-\zeta^{2}\right)^{2}-1}}=\pm\omega(\zeta-\sqrt{\zeta^{2}-1}).$
(23)

The behavior of the solution is determined by the damping coefficient. Just as
in the classical theory of the damped harmonic oscillator, there are four
cases:

  1. 1.

No damping, $\zeta=0$: For undamped oscillators,
$\lambda_{1}=\lambda_{4}=+i\omega$ and $\lambda_{2}=\lambda_{3}=-i\omega$. In
this case, we have

$x(t)=Ae^{+i\omega t}+Bte^{+i\omega t}+Ce^{-i\omega t}+Dte^{-i\omega t}.$ (24)

This reduces to the classical simple harmonic oscillator solution

$x(t)=E\cos{(\omega t)}+F\sin{(\omega t)}$ (25)

upon setting $B=D=0$, $A+C=E$, and $i(A-C)=F$.

  2. 2.

Underdamping, $0<\zeta<1$: In this case, the eigenvalues are complex-valued
and distinct. Writing

$\displaystyle\lambda_{1,2}=$
$\displaystyle\pm\omega(\zeta+i\sqrt{1-\zeta^{2}}),$ (26)
$\displaystyle\lambda_{3,4}=$
$\displaystyle\pm\omega(\zeta-i\sqrt{1-\zeta^{2}}),$ (27)

where $i$ is the imaginary unit, we have

$x(t)=Ae^{\omega\zeta t}e^{i\omega\sqrt{1-\zeta^{2}}t}+Be^{-\omega\zeta
t}e^{-i\omega\sqrt{1-\zeta^{2}}t}+Ce^{\omega\zeta
t}e^{-i\omega\sqrt{1-\zeta^{2}}t}+De^{-\omega\zeta
t}e^{i\omega\sqrt{1-\zeta^{2}}t}.$ (28)

This reduces to the classical damped solution

$x(t)=e^{-\omega\zeta
t}\left[E\cos{\left(\omega\sqrt{1-\zeta^{2}}t\right)}+F\sin{\left(\omega\sqrt{1-\zeta^{2}}t\right)}\right].$
(29)

upon setting $A=C=0$, $B+D=E$, and $-i(B-D)=F$.

  3. 3.

Critical damping, $\zeta=1$: For critically damped systems,
$\lambda_{1}=\lambda_{3}=+\omega$ and $\lambda_{2}=\lambda_{4}=-\omega$. In
this case, we have

$x(t)=Ae^{+\omega t}+Bte^{+\omega t}+Ce^{-\omega t}+Dte^{-\omega t}.$ (30)

This reduces to the classical damped harmonic solution

$x(t)=(C+Dt)e^{-\omega t}$ (31)

upon setting $A=B=0$.

  4. 4.

Overdamping, $\zeta>1$: For overdamped systems, all eigenvalues are real-
valued and distinct, leading to a linear combination of exponential growth and
decay:

$x(t)=Ae^{\lambda_{1}t}+Be^{\lambda_{2}t}+Ce^{\lambda_{3}t}+De^{\lambda_{4}t}.$
(32)

This reduces to the classical damped harmonic solution

$x(t)=Be^{-\omega(\zeta+\sqrt{\zeta^{2}-1})t}+De^{-\omega(\zeta-\sqrt{\zeta^{2}-1})t}$
(33)

upon setting $A=C=0$.

##  III Equivalent coupled-oscillator system

It is well-known that the fourth-order Pais-Uhlenbeck equation (12) is
mathematically equivalent to the following system of two second-order
equations:

$\displaystyle\ddot{x}+\mu_{1}x-\rho_{1}y=$ $\displaystyle 0,$ (34)
$\displaystyle\ddot{y}+\mu_{2}y-\rho_{2}x=$ $\displaystyle 0,$ (35)

Figure 2: Coupled oscillator system equivalent to the fourth-order Pais-
Uhlenbeck oscillator.

where

$\mu_{1}+\mu_{2}=4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right),$ (36)
$\mu_{1}\mu_{2}-\rho_{1}\rho_{2}=\omega^{4}.$ (37)

and both $x(t)$ and $y(t)$ satisfy the Pais-Uhlenbeck equation (12)
individually. This can be thought of as a system of two coupled oscillators,
as shown schematically in Figure 2. The springs in Figure 2 are not meant to
suggest that a damped oscillator is equivalent to two _conservative_
oscillators: as we will see, the stiffnesses of the supporting springs
($k_{1}$, $k_{2}$) and coupling spring ($k^{*}$) are generally not real-valued
and positive, and so the dynamics are not conservative. In general, the masses
and stiffnesses of the coupled system are related to the model parameters as
follows:

$\mu_{1}=\frac{k_{1}+k^{*}}{m_{1}},\quad\rho_{1}=\frac{k^{*}}{m_{1}},\quad\mu_{2}=\frac{k^{*}+k_{2}}{m_{2}},\quad\rho_{2}=\frac{k^{*}}{m_{2}}.$
(38)

whence we obtain

$\frac{m_{2}}{m_{1}}=\frac{\rho_{1}}{\rho_{2}},\quad\frac{k_{1}}{k^{*}}=\frac{\mu_{1}}{\rho_{1}}-1,\quad\frac{k_{2}}{k^{*}}=\frac{\mu_{2}}{\rho_{2}}-1.$
(39)

With initial conditions

$\displaystyle x(0)$ $\displaystyle=x_{0},$ $\displaystyle
y(0)=(a_{0}+\mu_{1}x_{0})/\rho_{1},$ (40) $\displaystyle\dot{x}(0)$
$\displaystyle=v_{0},$
$\displaystyle\dot{y}(0)=(j_{0}+\mu_{1}v_{0})/\rho_{1},$ (41)

$x$ will coincide with the damped harmonic oscillator. Now we observe that
there are only two damped harmonic oscillator parameters ($\omega$ and
$\zeta$), but there are four coupled oscillator parameters ($\mu_{1}$,
$\mu_{2}$, $\rho_{1}$, and $\rho_{2}$). It follows that for every damped
harmonic oscillator $x$ there exists a _two-parameter family of dual
oscillators_ $y$ satisfying the Pais-Uhlenbeck equation. For example, provided
$\zeta\neq 1/\sqrt{2}$ (this is a unique case that will be treated in Section
IV), we may set $\mu_{2}=(1/\alpha-1)\mu_{1}$ ($\alpha\neq 0$, $\alpha\neq 1$)
and $\rho_{2}=\beta^{-2}\rho_{1}$ ($\beta\neq 0$), leading to

$\displaystyle\mu_{1}$
$\displaystyle=4\alpha\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right),$
$\displaystyle\mu_{2}$
$\displaystyle=4(1-\alpha)\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right),$ (42)
$\displaystyle\rho_{1}$
$\displaystyle=\beta\omega^{2}\sqrt{16\alpha(1-\alpha)\left(\frac{1}{2}-\zeta^{2}\right)^{2}-1},$
$\displaystyle\rho_{2}$
$\displaystyle=\frac{1}{\beta}\omega^{2}\sqrt{16\alpha(1-\alpha)\left(\frac{1}{2}-\zeta^{2}\right)^{2}-1},$
(43)

where we have chosen the positive radical for $\rho_{1}$ in order to make the
coupling stiffness $k^{*}$ positive when $k^{*}$ is real-valued and $\beta$ is
positive.

###  III.1 Symmetric coupled oscillator

Table 1: Summary of all possible scenarios for a symmetric coupled oscillator ($m_{1}=m_{2}=m\in\mathbb{R}$, $k_{1}=k_{2}=k$). The sign of $\rho$ has been chosen so as to make the coupling stiffness $k^{*}$ positive when it is real-valued. Case | $\mu$ | $\rho$ | $k^{*}=m\rho$ | $k=m\mu-k^{*}$ | $r=k^{*}/k$
---|---|---|---|---|---
$\zeta=0$ | $+\omega^{2}$ | 0 | 0 | $+m\omega^{2}$ | 0
$\zeta\in(0,1/\sqrt{2})$ | real, positive | imaginary | imaginary | complex | complex
$\zeta=1/\sqrt{2}$ | 0 | $+i\omega^{2}$ | $+im\omega^{2}$ | $-im\omega^{2}$ | -1
$\zeta\in(1/\sqrt{2},1)$ | real, negative | imaginary | imaginary | complex | complex
$\zeta=1$ | $-\omega^{2}$ | 0 | 0 | $-m\omega^{2}$ | 0
$\zeta>1$ | real, negative | real, positive | real, positive | real, negative | real, negative
$\zeta\rightarrow+\infty$ | $-\infty$ | $+\infty$ | $+\infty$ | $-\infty$ | $-1/2$

Henceforth, we will restrict attention to the symmetric coupled oscillator, in
which the two masses are identical and real-valued
($m_{1}=m_{2}=m\in\mathbb{R}$) and the two supporting springs have identical
stiffnesses ($k_{1}=k_{2}=k$, not necessarily real or positive), giving
$\rho_{1}=\rho_{2}=\rho$ and $\mu_{1}=\mu_{2}=\mu$. This corresponds to the
particular dual oscillator with $\alpha=1/2$ and $\beta=1$, giving

$\mu=2\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right),$ (44)
$\rho=\omega^{2}\sqrt{4\left(\frac{1}{2}-\zeta^{2}\right)^{2}-1}=2\omega^{2}\zeta\sqrt{\zeta^{2}-1}.$
(45)

The coupling stiffness is given by

$k^{*}=m\rho=2m\omega^{2}\zeta\sqrt{\zeta^{2}-1},$ (46)

and the supporting stiffnesses are given by

$k=m\mu-k^{*}=m(\mu-\rho)=2m\omega^{2}\left(\frac{1}{2}-\zeta^{2}-\zeta\sqrt{\zeta^{2}-1}\right).$
(47)

We define the stiffness ratio as

$r=\frac{k^{*}}{k}=\left(\frac{\mu}{\rho}-1\right)^{-1}=-\left[\frac{\zeta^{2}-\frac{1}{2}}{\zeta\sqrt{\zeta^{2}-1}}+1\right]^{-1}.$
(48)

The real and imaginary parts of the stiffness ratio $r$ are plotted versus
$\zeta$ in Figure 3. Below we consider the various damping cases in turn,
which are summarized in Table 1.

Figure 3: Real and imaginary parts of the stiffness ratio $r$ versus damping
ratio $\zeta$.

  1. 1.

No damping, $\zeta=0$: In the absence of damping, $\mu=+\omega^{2}$, $\rho=0$,
$k^{*}=0$, $k=+m\omega^{2}$, and $r=0$. This corresponds to a case of two
uncoupled simple harmonic oscillators with positive stiffnesses and identical
natural frequencies $\omega$.

  2. 2.

Underdamping, $0<\zeta<1$: For underdamping, the stiffnesses $k$ and $k^{*}$
are non-real. To be more precise, we must make a distinction between damping
ratios less than, equal to, and greater than $\zeta=1/\sqrt{2}$:

    1. (a)

$0<\zeta<1/\sqrt{2}$: We have that $\mu$ is real and positive, $\rho$ is
imaginary, $k^{*}$ is imaginary, $k$ is complex, and the stiffness ratio $r$
is complex.

    2. (b)

$\zeta=1/\sqrt{2}$: We have that $\mu=0$, $\rho=+i\omega^{2}$,
$k^{*}=+im\omega^{2}$, $k=-im\omega^{2}$, and $r=-1$. Again, this case will be
discussed in more detail in Section IV.

    3. (c)

$1/\sqrt{2}<\zeta<1$: We have that $\mu$ is real and negative, $\rho$ is
imaginary, $k^{*}$ is imaginary, $k$ is complex, and the stiffness ratio $r$
is complex.

Figure 4 shows representative coupled oscillator dynamics for an underdamped
case with $\zeta=0.2$, $\omega=1$, $x_{0}=1$, and $v_{0}=0$. In this case, the
damped harmonic oscillator $x(t)$ is real-valued, while the dual oscillator
$y(t)$ is imaginary.

Figure 4: Representative coupled oscillator dynamics for an underdamped system
with $\zeta=0.2$, $\omega=1$, $x_{0}=1$, and $v_{0}=0$. The damped harmonic
oscillator $x(t)$ is real-valued, while the dual oscillator $y(t)$ is
imaginary.

  3. 3.

Critical damping, $\zeta=1$: For critical damping, $\mu=-\omega^{2}$,
$\rho=0$, $k^{*}=0$, $k=-m\omega^{2}$, and $r=0$. As in the absence of
damping, this corresponds to another case in which the coupling spring
stiffness is zero. The difference is that here the supporting spring
stiffnesses are _negative_.

  4. 4.

Overdamping, $\zeta>1$: All quantities are real-valued. We have that $\mu$ is
negative, $\rho$ is positive, $k^{*}$ is positive, $k$ is negative, and $r$ is
negative. In the limit as $\zeta\rightarrow+\infty$, $k^{*}\rightarrow+\infty$
and $k\rightarrow-\infty$, and $r\rightarrow-1/2$. This leads to an
interesting interpretation of overdamping. As the system becomes more and more
overdamped, the coupling spring becomes ever more positively stiff, and the
supporting springs become ever more negatively stiff, with their ratio
approaching $-1/2$ asymptotically. Figure 5 shows representative coupled
oscillator dynamics for an overdamped case with $\zeta=1.5$, $\omega=1$,
$x_{0}=1$, and $v_{0}=0$. In this case, the damped harmonic oscillator $x(t)$
and the dual oscillator $y(t)$ are both real-valued.

Figure 5: Representative coupled oscillator dynamics for an overdamped system
with $\zeta=1.5$, $\omega=1$, $x_{0}=1$, and $v_{0}=0$. The damped harmonic
oscillator $x(t)$ and the dual oscillator $y(t)$ are both real-valued.

It should be noted that spring-like elements with effectively negative
stiffnesses are physically achievable. Lee, Goverdovskiy, and Temnikov (2007)
Thus, it should be possible to test the results presented here for overdamped
oscillators experimentally. Coupled oscillators mimicking the overdamped
harmonic oscillator may even find applications to vibration suppression
systems.

##  IV Optimal control: $\zeta=1/\sqrt{2}$

Figure 6: Plot of normalized response amplitude $X(\nu)$ versus frequency
ratio $\nu$ for various damping ratios $\zeta$. The local maximum does not
exist for $\zeta>1/\sqrt{2}\approx 0.7071$.

Consider a forced damped harmonic oscillator of the form

$\ddot{x}+2\omega\zeta\dot{x}+\omega^{2}x=f_{0}\sin{(\Upomega t)},$ (49)

where $f_{0}$ is the forcing magnitude and $\Upomega$ is the forcing
frequency. At steady state, the solution has normalized magnitude

$X(\nu)=\frac{1}{\sqrt{(1-\nu^{2})^{2}+(2\zeta\nu)^{2}}},$ (50)

where $\nu=\Upomega/\omega$ is the frequency ratio. A plot of $X(\nu)$ versus
$\nu$ for various values of $\zeta$ is shown in Figure 6. It is
straightforward to show that $X(\nu)$ has a local maximum at
$\nu=\sqrt{1-2\zeta^{2}}$, provided $\zeta<1/\sqrt{2}$. For
$\zeta>1/\sqrt{2}$, there is no longer a local maximum. The critical value
$\zeta=1/\sqrt{2}$ at which the local maximum disappears happens to be the
“optimal control value” in terms of the tradeoff between rise time and percent
overshoot in response to a step input. Jacobs (1965)

Looking at the second-order damped harmonic oscillator equation (5),
$\zeta=1/\sqrt{2}$ does not appear to have a special role; it is simply one
case of underdamping. However, it is _immediately_ clear from the Pais-
Uhlenbeck equation (12) that $\zeta=1/\sqrt{2}$ is significant. The inertial
term,

$4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\ddot{x},$ (51)

is the only term in the Pais-Uhlenbeck equation (12) that involves the damping
ratio. We might call this the “inertial damping” term. Interestingly, the
inertial damping term vanishes for $\zeta=1/\sqrt{2}$, yielding

$\ddddot{x}+\omega^{4}x=0.$ (52)

This leads to a physical (rather than a purely control-theoretical)
interpretation for the optimal control value. When $\zeta=1/\sqrt{2}$, the
oscillator’s “jounce” (time rate-of-change of jerk) is not affected by its
current acceleration, only by its current displacement from equilibrium. In
other words, there is no inertial damping for the optimal control value.

Additionally, $\zeta=1/\sqrt{2}$ leads to a degenerate case in which the
coupled oscillator parameters $\mu_{1}=-\mu_{2}$, according to (36). Setting
$\mu_{1}=\gamma$ and $\rho_{2}=\beta^{-2}\rho_{1}$ ($\beta\neq 0$), we have

$\displaystyle\mu_{1}$ $\displaystyle=\gamma,$ $\displaystyle\mu_{2}$
$\displaystyle=-\gamma,$ (53) $\displaystyle\rho_{1}$
$\displaystyle=\beta\sqrt{-\gamma^{2}-\omega^{4}},$ $\displaystyle\rho_{2}$
$\displaystyle=\frac{1}{\beta}\sqrt{-\gamma^{2}-\omega^{4}}.$ (54)

In Section III.1, when we considered the case of $\zeta=1/\sqrt{2}$, we were
assuming identical masses and identical supporting springs, which amounts to
setting $\gamma=0$ and $\beta=1$ here. That was the only underdamped value for
which the stiffness ratio (48) was real-valued ($r=-1$), giving

$k^{*}=+im\omega^{2},\quad k=-im\omega^{2}.$ (55)

Thus, optimal control is the only case in which the symmetric coupled
oscillator with real masses has purely imaginary spring stiffnesses. Figure 7
shows representative coupled oscillator dynamics for an optimal control case
with $\zeta=1/\sqrt{2}$, $\omega=1$, $x_{0}=1$, and $v_{0}=0$. In this case,
the damped harmonic oscillator $x(t)$ is real-valued, while the dual
oscillator $y(t)$ is imaginary.

Figure 7: Representative coupled oscillator dynamics for an optimal control
case with $\zeta=1/\sqrt{2}$, $\omega=1$, $x_{0}=1$, and $v_{0}=0$. The damped
harmonic oscillator $x(t)$ is real-valued, while the dual oscillator $y(t)$ is
imaginary.

##  V Variational formulation

The fact that every solution to the damped harmonic oscillator equation (5)
also satisfies the fourth-order Pais-Uhlenbeck equation (12) also leads to a
new variational formulation for the damped harmonic oscillator. The Pais-
Uhlenbeck action is given by

$\mathcal{S}=\int_{t_{1}}^{t_{2}}\left[\ddot{x}^{2}-4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\dot{x}^{2}+\omega^{4}x^{2}\right]dt,$
(56)

and the Pais-Uhlenbeck equation (12) can be recovered by taking
$\delta\mathcal{S}=0$. It follows that every solution to the damped harmonic
oscillator equation (5) constitutes a critical point of the Pais-Uhlenbeck
action. It is instructive to compare this result to previous variational
formulations of the damped harmonic oscillator. Like the dual system approach
of Bateman, Morse, Feshbach, and Tikochinsky, Bateman (1931); Morse and
Feshbach (1953); Feshbach and Tikochinsky (1977) the Pais-Uhlenbeck action
introduces dual oscillators. The difference here is that there is actually a
two-parameter family of dual oscillators $y$; moreover, the dual oscillators
do not appear explicitly in the action and are not varied independently of
$x$. Unlike the Kanai-Caldirola Lagrangian, Caldirola (1941); Kanai (1948) the
Pais-Uhlenbeck Lagrangian does not depend explicitly on time. And unlike the
Caldeira-Leggett approach Caldeira and Leggett (1983); Zwanzig (1973); Ford
and Kac (1987) taken by Fukagawa and Fujitani, Fukagawa and Fujitani (2012)
the Pais-Uhlenbeck action does not require the use of non-holonomic
constraints to include the dissipative terms.

It may be that variational formulations can be found for other dissipative
systems by considering higher-order equations. The problem becomes how to
identify those equations. One possible method is to manipulate the second-
order equation directly, as we did in Section II. Having arrived at (12), we
may multiply by $\delta x$,

$\ddddot{x}\delta
x+4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\ddot{x}\delta
x+\omega^{4}x\delta x=0,$ (57)

integrate by parts judiciously,

$\int_{t_{1}}^{t_{2}}\left[\ddot{x}\delta\ddot{x}-4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\dot{x}\delta\dot{x}+\omega^{4}x\delta
x\right]dt=0,$ (58)

thus completing the variation and arriving at the Pais-Uhlenbeck action
principle

$\delta\int_{t_{1}}^{t_{2}}\left[\ddot{x}^{2}-4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\dot{x}^{2}+\omega^{4}x^{2}\right]dt=0.$
(59)

Another method is to start from the d’Alembert-Lagrange principle (3) and look
for a particular variation that puts it in variational form. For the damped
harmonic oscillator, the d’Alembert-Lagrange principle gives

$-\int_{t_{1}}^{t_{2}}(m\ddot{x}+c\dot{x}+kx)\delta{x}dt=0.$ (60)

Both the inertial term $m\ddot{x}\delta x$ and the spring term $kx\delta x$
can be put into variational form, because they each contain an even number of
time derivatives. However, the damping term $c\dot{x}\delta x$ has an odd
number of derivatives and therefore cannot be expressed as the exact variation
of a work function. If we had an additional term proportional to
$x\delta\dot{x}$, it could cancel out the damping term. Thus, we might look
for a variational principle of the form

$\int_{t_{1}}^{t_{2}}(m\ddot{x}+c\dot{x}+kx)(\delta{x}+\alpha\delta\dot{x})dt=0.$
(61)

However, this will give a term proportional to $\ddot{x}\delta\dot{x}$, which
also has an odd number of derivatives. Thus, in general we must look for a
variational principle of the form

$\int_{t_{1}}^{t_{2}}(m\ddot{x}+c\dot{x}+kx)(\delta{x}+\alpha\delta\dot{x}+\beta\delta\ddot{x})dt=0.$
(62)

Expanding, we have

$\int_{t_{1}}^{t_{2}}(m\ddot{x}\delta{x}+\alpha m\ddot{x}\delta\dot{x}+\beta
m\ddot{x}\delta\ddot{x}+c\dot{x}\delta{x}+\alpha c\dot{x}\delta\dot{x}+\beta
c\dot{x}\delta\ddot{x}+kx\delta{x}+\alpha kx\delta\dot{x}+\beta
kx\delta\ddot{x})dt=0.$ (63)

The terms with odd numbers of derivatives are $\alpha m\ddot{x}\delta\dot{x}$,
$c\dot{x}\delta{x}$, $\beta c\dot{x}\delta\ddot{x}$, and $\alpha
kx\delta\dot{x}$. Upon integration by parts, we want $\alpha
m\ddot{x}\delta\dot{x}$ to cancel $\beta c\dot{x}\delta\ddot{x}$, and for
$c\dot{x}\delta{x}$ to cancel $\alpha kx\delta\dot{x}$. This requires

$\alpha m=\beta c\quad\text{and}\quad c=\alpha k,$ (64)

which gives $\alpha=c/k$ and $\beta=m/k$. Hence, we are looking for a
variational principle of the form

$\int_{t_{1}}^{t_{2}}(m\ddot{x}+c\dot{x}+kx)\left(\delta{x}+\frac{c}{k}\delta\dot{x}+\frac{m}{k}\delta\ddot{x}\right)dt=0.$
(65)

We observe that $\delta{x}+(c/k)\delta\dot{x}+(m/k)\delta\ddot{x}$ is actually
a variation of the governing equation, up to a factor of $k$. Expanding, we
have

$\int_{t_{1}}^{t_{2}}\left(m\ddot{x}\delta{x}+\frac{mc}{k}\ddot{x}\delta\dot{x}+\frac{m^{2}}{k}\ddot{x}\delta\ddot{x}+c\dot{x}\delta{x}+\frac{c^{2}}{k}\dot{x}\delta\dot{x}+\frac{mc}{k}\dot{x}\delta\ddot{x}+kx\delta{x}+cx\delta\dot{x}+mx\delta\ddot{x}\right)dt=0.$
(66)

Judicious integration by parts yields

$\int_{t_{1}}^{t_{2}}\left(-m\dot{x}\delta\dot{x}+\frac{mc}{k}\ddot{x}\delta\dot{x}+\frac{m^{2}}{k}\ddot{x}\delta\ddot{x}+c\dot{x}\delta{x}+\frac{c^{2}}{k}\dot{x}\delta\dot{x}-\frac{mc}{k}\ddot{x}\delta\dot{x}+kx\delta{x}-c\dot{x}\delta{x}-m\dot{x}\delta\dot{x}\right)dt=0,$
(67)

which simplifies to

$\int_{t_{1}}^{t_{2}}\left[\frac{m^{2}}{k}\ddot{x}\delta\ddot{x}+\left(\frac{c^{2}}{k}-2m\right)\dot{x}\delta\dot{x}+kx\delta{x}\right]dt=0.$
(68)

Each term can now be put into variational form, yielding

$\delta\int_{t_{1}}^{t_{2}}\left[\frac{m^{2}}{2k}\ddot{x}^{2}+\frac{1}{2}\left(\frac{c^{2}}{k}-2m\right)\dot{x}^{2}+\frac{1}{2}kx^{2}\right]dt=0,$
(69)

or equivalently,

$\delta\int_{t_{1}}^{t_{2}}\left[\ddot{x}^{2}-4\omega^{2}\left(\frac{1}{2}-\zeta^{2}\right)\dot{x}^{2}+\omega^{4}x^{2}\right]dt=0.$
(70)

Future work may consider whether either of the above approaches generalizes to
other kinds of dissipative systems.

##  VI Summary and conclusion

The main result of this paper is that the classical damped harmonic oscillator
belongs to the family of fourth-order Pais-Uhlenbeck oscillators. Three
immediate corollaries follow:

  1. (I)

For every damped harmonic oscillator $x$, there exists a two-parameter family
of dual oscillators $y$ satisfying the Pais-Uhlenbeck equation. The damped
harmonic oscillator and any of its duals constitute a system of two coupled
oscillators with non-standard spring stiffnesses. The various cases are
summarized in Table 1 for coupled oscillators with identical masses and
identical supporting springs. Of particular interest is the case of
overdamping $\zeta>1$, for which the stiffness of the coupling spring is
positive while the stiffnesses of the supporting springs are negative. This
scenario should be physically achievable Lee, Goverdovskiy, and Temnikov
(2007) and may even find applications to vibration suppression systems.

  2. (II)

The optimal control value $\zeta=1/\sqrt{2}$, which gives the best tradeoff
between rise time and percent overshoot in response to a step input, Jacobs
(1965) appears prominently in the fourth-order formulation and has a clear
physical interpretation. For this value of the damping ratio, there is no
inertial damping.

  3. (III)

The solutions of the damped harmonic oscillator make the Pais-Uhlenbeck action
stationary, yielding a new variational formulation for the damped harmonic
oscillator. Two systematic methods have been given for arriving at the Pais-
Uhlenbeck action starting from the damped harmonic oscillator equation. It may
be possible to generalize these approaches and identify stationary actions for
other kinds of dissipative systems which do not conform to Hamilton’s
principle.

## Declarations

### Conflicts of interest

The author declares that he has no conflicts of interest or competing
interests.

### Data availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * Newton (1687) I. Newton,  _Philosophiae Naturalis Principia Mathematica_ (Royal Society of London, 1687).
  * d’Alembert (1743) J. L. R. d’Alembert,  _Traité de dynamique_ (J. B. Coignard, 1743).
  * Lagrange (1811) J. L. Lagrange,  _Mécanique Analytique_ (Gauthier-Villars, 1811).
  * Hamilton (1834) W. R. Hamilton, “On a general method in dynamics,” Philosophical Transactions of the Royal Society 124, 247–308 (1834).
  * Hamilton (1835) W. R. Hamilton, “Second essay on a general method in dynamics,” Philosophical Transactions of the Royal Society 125, 95–144 (1835).
  * Lanczos (1949) C. Lanczos,  _The Variational Principles of Mechanics_ (Dover Publications, 1949).
  * Gray and Taylor (2007) C. G. Gray and E. F. Taylor, “When action is not least,” American Journal of Physics 75, 434–458 (2007).
  * Jacobi (1837) C. G. J. Jacobi, “Zür Theorie der Variations-Rechnung und der Differential-Gleichungen,” Journal für die reine und angewandte Mathematik 1837, 68–82 (1837).
  * Whittaker (1904) E. T. Whittaker,  _A Treatise on the Analytical Dynamics of Particles and Rigid Bodies_ (Cambridge University Press, 1904).
  * Kane  _et al._ (2000) C. Kane, J. E. Marsden, M. Ortiz, and M. West, “Variational integrators and the newmark algorithm for consevative and dissipative mechanical systems,” International Journal of Numerical Methods in Engineering 49, 1295–1325 (2000).
  * Baddour (2007) N. Baddour, “Leading-edge applied mathematical modeling research,” (Nova Science Publisher, 2007) Chap. Hamilton’s Principle for the Derivation of Equations of Motion, pp. 155–182.
  * Flannery (2011) M. R. Flannery, “The elusive d’Alembert-Lagrange dynamics of nonholonomic systems,” American Journal of Physics 79, 932–944 (2011).
  * Gray, Karl, and Novikov (1996) C. G. Gray, G. Karl, and V. A. Novikov, “Direct use of variational principles as an approximation technique in classical mechanics,” American Journal of Physics 64, 1177–1184 (1996).
  * Kim, Dargush, and Ju (2013) J. Kim, G. F. Dargush, and Y.-K. Ju, “Extended framework of Hamilton’s principle for continuum dynamics,” International Journal of Solids and Structures 50, 3418–3429 (2013).
  * Marsden and West (2001) J. Marsden and M. West, “Discrete mechanics and variational integrators,” Acta Numerica 10, 357–514 (2001).
  * Lew  _et al._ (2003) A. Lew, J. E. Marsden, M. Ortiz, and M. West, “Asynchronous variational integrators,” Archive for Rational Mechanics and Analysis 167, 85–146 (2003).
  * Lew  _et al._ (2004a) A. Lew, J. E. Marsden, M. Ortiz, and M. West, “Finite element methods: 1970’s and beyond,” (CIMNE, 2004) Chap. An overview of variational integrators, pp. 98–115.
  * Lew  _et al._ (2004b) A. Lew, J. E. Marsden, M. Ortiz, and M. West, “Variational time integrators,” International Journal of Numerical Methods in Engineering 601, 153–212 (2004b).
  * Kale and Lew (2007) K. G. Kale and A. J. Lew, “Parallel asynchronous variational integrators,” International Journal for Numerical Methods in Engineering 70, 291–321 (2007).
  * Ober-Blöbaum, Junge, and Marsden (2011) S. Ober-Blöbaum, O. Junge, and J. E. Marsden, “Discrete mechanics and optimal control: An analysis,” Control, Optimisation, and Calculus of Variations 17, 322–352 (2011).
  * Kraus and Maj (2015) M. Kraus and O. Maj, “Variational integrators for nonvariational partial differential equations,” Physica D: Nonlinear Phenomena 310, 37–71 (2015).
  * Hall and Leok (2015) J. Hall and M. Leok, “Spectral variational integrators,” Numerische Mathematik 30, 681–740 (2015).
  * Lew and Mata (2016) A. J. Lew and P. Mata, “A brief introduction to variational integrators,” in  _Nonlinear Structural Dynamics and Flexible Multibody Dynamics_ (2016) pp. 201–291.
  * Schmitt and Leok (2018) J. M. Schmitt and M. Leok, “Properties of hamiltonian variational integrators,” IMA Journal of Numerical Analysis 38, 377–398 (2018).
  * Man, Gao, and Zhong (2020) S. Man, Q. Gao, and W. Zhong, “Variational integrators in holonomic mechanics,” Mathematics 8 (2020).
  * Limebeer, Ober-Blöbaum, and Farshi (2020) D. J. N. Limebeer, S. Ober-Blöbaum, and F. H. Farshi, “Variational integrators for dissipative systems,” IEEE Transactions on Automatic Control 65, 1381–1396 (2020).
  * Kraus and Tyranowski (2021) M. Kraus and T. M. Tyranowski, “Variational integrators for stochastic dissipative hamiltonian systems,” IMA Journal of Numerical Analysis 41, 1318–1367 (2021).
  * Bateman (1931) H. Bateman, “On dissipative systems and related variational principles,” Physical Review 38, 815–819 (1931).
  * Morse and Feshbach (1953) P. M. Morse and H. Feshbach,  _Methods of Theoretical Physics_ (McGraw-Hill, 1953).
  * Feshbach and Tikochinsky (1977) H. Feshbach and Y. Tikochinsky, “Quantization of the damped harmonic oscillator,” Transactions of the New York Academy of Sciences 38, 44–53 (1977).
  * Caldirola (1941) P. Caldirola, “Forze non conservative nella meccanica quantistica,” Il Nuovo Cimento 18, 393–400 (1941).
  * Kanai (1948) E. Kanai, “On the quantization of the dissipative system,” Progress in Theoretical Physics 3, 440–442 (1948).
  * Caldeira and Leggett (1983) A. O. Caldeira and A. J. Leggett, “Path integral approach to Brownian motion,” Physica A 121, 587–616 (1983).
  * Zwanzig (1973) R. Zwanzig, “Nonlinear generalized Langevin equations,” Journal of Statistical Physics 9, 215–220 (1973).
  * Ford and Kac (1987) G. W. Ford and M. Kac, “On the quantum langevin equation,” Journal of Statistical Physics 46, 803–810 (1987).
  * Fukagawa and Fujitani (2012) H. Fukagawa and Y. Fujitani, “A variational principle for dissipative fluid dynamics,” Progress of Theoretical Physics 127, 921–935 (2012).
  * Stokes (1845) G. Stokes, “On the theories of the internal friction of fluids in motion, and of the equilibrium and motion of elastic solids,” Transactions of the Cambridge Philosophical Society 8, 287–319 (1845).
  * Pais and Uhlenbeck (1950) A. Pais and G. E. Uhlenbeck, “On field theories with non-localized action,” Physical Review 79, 145–165 (1950).
  * Mostafazadeh (2010) A. Mostafazadeh, “A Hamilton formulation of the Pais-Uhlenbeck oscillator that yields a stable and unitary quantum system,” Physics Letters A 375, 93–98 (2010).
  * Baleanu  _et al._ (2012) D. Baleanu, I. Petras, J. H. Asad, and M. P. Velasco, “Fractional Pais-Uhlenbeck oscillator,” International Journal of Theoretical Physics 51 (2012).
  * Lee, Goverdovskiy, and Temnikov (2007) C.-M. Lee, V. N. Goverdovskiy, and A. I. Temnikov, “Design of springs with “negative” stiffness to improve vehicle driver vibration isolation,” Journal of Sound and Vibration 302, 865–874 (2007).
  * Jacobs (1965) O. L. R. Jacobs, “The damping ratio of an optimal control system,” IEEE Transactions on Automatic Control 10, 473–476 (1965).