A new perspective on the ubiquity of classical harmonic oscillators is presented based on the two-variable Taylor expansion of a perturbed system's total energy E ( q , q ̇ ), where q(t) is the system displacement as a function of time t and q ̇ ( t ) = d q / d t. This generalised approach permits derivation of the lossless oscillator equation from energy arguments only, yielding a universal equation for the oscillation frequency ω = ( 2 E / q 2 ) / ( 2 E / q ̇ 2 ) which may be applied to arbitrary systems without the need to form system-specific linearised models. As illustrated by a range of examples, this perspective gives a unifying explanation for the prevalence of harmonic oscillators in classical physics, can be extended to include damping effects and driving forces, and is a powerful tool for simplifying the analyses of perturbed systems.

The ubiquity of harmonic oscillators is a central theme of classical physics,1–3 with applications ranging from mechanical systems—such as springs,4,5 pendulums,6–8 and rocking bodies9,10—to acoustic resonators,11 electronic circuits,12–14 and laser fluid traps.15,16 Indeed, such is the universality of simple harmonic motion that it can be used to characterise phenomena at almost all scales of interest,17 whether those be modes of molecular vibration,18 the resonant frequencies of bridges,19,20 or perturbations of gravitational orbits.2,3

For undamped mechanical systems involving a fixed mass m subject to conservative forces F, such ubiquity is traditionally understood by considering the time t dependence of small displacements x(t) from a point of stable equilibrium x = 0. In particular, by Taylor expanding the system's potential energy V(x) about x = 0 (where V is assumed to have a minimum), one obtains a local harmonic potential of the form2 
(1)
and thence a restoring force F = d V / d x V ( 0 ) x. It then follows from Newton's second law F = m x ̈, with x ̈ d x ̇ / d t and x ̇ d x / d t, that the perturbed system will undergo simple harmonic motion according to the lossless oscillator equation
(2)
is the natural oscillation frequency.

Since no special assumptions are made about the form of V(x), this perspective on why harmonic oscillations arise in rectilinear mechanical systems is a powerful one. Nevertheless, there are at least three respects in which it fails to account for the ubiquity of harmonic oscillators in classical physics more generally: first (i), equation (2) must be modified on an ad hoc basis to handle multi-directional geometries (such as rotational systems); second (ii), the argument leading to equation (2) requires fixed m, and breaks-down for systems with varying inertia (e.g., multivariable systems subject to constraints); and third (iii), the explanation is a mechanical one, and has no direct application to other contexts (such as oscillating circuits) without recourse to careful argument by analogy.21,22

In this article, we propose a new perspective on the ubiquity of classical harmonic oscillators, one aimed at addressing the problems (i), (ii), and (iii) inherent to the traditional approach. In particular, rather than appealing to a single-variable Taylor expansion of the system's potential energy V, we consider a two-variable expansion of the system's total energy E ( q , q ̇ ), where q(t) is the displacement from equilibrium and q ̇ = d q / d t (Sec. II); in this way, we obtain a universal expression for the lossless oscillator equation
(3)
is a universal equation for the oscillation frequency. Since this equation is derived from energy arguments only, and is expressed in terms of a general coordinate q, it can be used to analyze the stability of arbitrary systems (mechanical, electrical, etc.) without the need to form system-specific linearised models (Secs. III and IV). Moreover, by relaxing the lossless (fixed energy) assumption, Eq. (3) can be further generalised to form a universal driven oscillator equation with damping and forcing effects expressed in terms of power dissipation and delivery (Sec. VII). We illustrate these ideas using a range of examples in Secs. V–VII, including stability results for a novel pendulum system.7,8

Of course, the idea of using energy arguments to treat oscillations in perturbed systems is an old one, and there are close connections between the method introduced here and other energy techniques such as Lagrangian mechanics.1–3 However, our procedure—whereby the lossless oscillator equation (3) is obtained directly from the two-variable Taylor expansion of the energy function E ( q , q ̇ )—does not seem to have been documented elsewhere in the literature, and in this sense may be considered new. Indeed, the universal equation for the natural oscillation frequency ω = ( 2 E / q 2 ) / ( 2 E / q ̇ 2 ) appears to be reported here for the first time. With these points in mind, we consider in Sec. VIII the merits of our approach compared to those of more traditional methods.

The origin of our expression for the lossless oscillator equation (3) may be understood by considering an arbitrary system with total energy E ( q , q ̇ ), where q(t) is a generalised coordinate for the system displacement as a function of time t, and q ̇ = d q / d t. Thus, if the system is in static equilibrium at a point ( q , q ̇ ) = ( q 0 , 0 ), then E has a stationary value there with
(4)
To analyse the behaviour of the system near this equilibrium, we assume that the second-order Taylor expansion of E exists at ( q 0 , 0 ) with a vanishing mixed derivative, that is,
(5)
This assumption applies to a wide range of possible functions E ( q , q ̇ ) but is particularly suitable given the prevalence of classical systems with total energy of the form
(6)
is the kinetic energy, V(q) is the potential energy, and M(q) > 0 is the system inertia, which can, in general, vary with displacement q (see Sec. VI). Energy functions taking this form have a vanishing mixed derivative at ( q 0 , 0 ) because M(q) and V(q) are functions of q only.
Now suppose that the system is perturbed from the equilibrium ( q 0 , 0 ) by supplying (or removing) a small amount of energy E 1, then the total energy becomes
(7)
where the coordinates ( q , q ̇ ) of the perturbed system are
(8)
with q 1 ( t ) and q ̇ 1 ( t ) as perturbations. It then follows from the two-variable Taylor expansion about ( q 0 , 0 ) that23,
(9)
where the derivatives are evaluated at ( q 0 , 0 ) and R3 is the remainder term,23,24 which scales as q 1 3. Here, the first-order derivatives and mixed derivatives vanish according to Eqs. (4) and (5), while for arbitrarily small perturbations ( q 1 0 ), the remainder term R3 may be neglected;24 in this way the expansion simplifies to
(10)
For an undamped, un-driven (constant energy) system, we have d E 1 / d t = 0. Thus, differentiating Eq. (10) with respect to t, for dynamic solutions ( q ̇ 0), we obtain
(11)
This equation is the universal lossless oscillator equation (3) as required, with ω as the natural oscillation frequency, and derivatives evaluated at ( q 0 , 0 ). Notice that if all derivatives of E ( q , q ̇ ) higher than second-order vanish, then R 3 ( q , q ̇ ) = 0, in which case the approximation in Eq. (11) is exact (i.e., an equality) for perturbations of any size.24 

Observe that Eq. (11) is derived using energy arguments only; no assumptions have been made about (i) the system geometry, (ii) whether the system has fixed inertia, or (iii) the system context (e.g., mechanical or electrical). Thus, when compared to the traditional approach outlined in Sec. I, Eq. (11) gives a more widely applicable explanation for the ubiquity of harmonic oscillators in classical physics.25 Indeed, a particular advantage of our approach is that it yields an expression for ω that can be applied to arbitrary systems without the need to form customised linear models (see Secs. V and VI).

The oscillator equation (11) derived in Sec. II may be thought of more generally as a linearised equation of motion for a free system perturbed from equilibrium, and can therefore be used to analyse the linear stability of the stationary point ( q 0 , 0 ). To see this, observe that the general solution to Eq. (11) is
(12)
where q + and q are constants whose values depend on the initial conditions. Thus, if ω 2 > 0, then solutions are stable harmonic oscillations, as discussed above;26 however, if ω 2 < 0, then solutions are unstable, i.e.,
(13)
Notice that because ω 2 = ( 2 E / q 2 ) / ( 2 E / q ̇ 2 ), stability corresponds to a minimum of the energy surface E ( q , q ̇ ) at ( q 0 , 0 ), whereas instability corresponds to a saddle-point.26 These two modes of behaviour can be visualised by sketching the solution trajectories ( q 1 ( t ) , q ̇ 1 ( t ) ) in phase space, yielding what might be termed a universal phase portrait for systems perturbed from equilibrium (see Fig. 1). Indeed, by writing Eq. (10) as
(14)
we see that such trajectories are elliptical (stable) when ω 2 > 0, and hyperbolic (unstable) when ω 2 < 0. Note that the linear analysis leading to Eq. (11) fails if ω = 0, since higher-order terms must then be considered.
Fig. 1.

Universal phase portraits for ( q 1 / q * ) 2 + ( q ̇ 1 / ω q * ) 2 = 1. Orbits about the unit circle describe harmonic oscillations and require ω 2 > 0 and q * 2 > 0. The hyperbolic trajectories correspond to instability in the cases that either: (i) ω 2 < 0 and q * 2 > 0 (branches intersecting the q 1 / | q * | axis); or (ii) ω 2 < 0 and q * 2 < 0 (branches intersecting the q ̇ 1 / | ω q * | axis).

Fig. 1.

Universal phase portraits for ( q 1 / q * ) 2 + ( q ̇ 1 / ω q * ) 2 = 1. Orbits about the unit circle describe harmonic oscillations and require ω 2 > 0 and q * 2 > 0. The hyperbolic trajectories correspond to instability in the cases that either: (i) ω 2 < 0 and q * 2 > 0 (branches intersecting the q 1 / | q * | axis); or (ii) ω 2 < 0 and q * 2 < 0 (branches intersecting the q ̇ 1 / | ω q * | axis).

Close modal

The remainder of this article is focused on illustrating how the ideas developed in Secs. II and III apply to practical systems. To that end, however, it is convenient to summarise the key results in the form of a theorem:

Theorem Let E ( q , q ̇ ) be the energy of an un-damped, un-driven system, where q(t) is the displacement as a function of time t and q ̇ = d q / d t, and let ( q 0 , 0 ) be a stationary point of E ( q , q ̇ ) for which [ 2 E / q q ̇ ] ( q 0 , 0 ) = 0. If the Taylor expansion of E exists at ( q 0 , 0 ) with non-zero second derivatives, then arbitrarily small perturbations q 1 ( t ) = q ( t ) q 0 from q0 obey the lossless oscillator equation
(15)
is the square of the oscillation frequency ω. If all derivatives of E ( q , q ̇ ) higher than second-order vanish, then this approximation is exact for perturbations of any size.

This theorem motivates the following three-step method for analysing free systems perturbed from equilibrium:

  1. Write down an expression for the total energy E ( q , q ̇ ) in terms of q(t) and q ̇ ( t ) = d q / d t.

  2. Determine the stationary points ( q , q ̇ ) = ( q 0 , 0 ) of the energy function E ( q , q ̇ ) by solving
    (16)
  3. Evaluate the second derivatives of E ( q , q ̇ ) at the stationary points ( q 0 , 0 ). If these derivatives are non-zero, with [ 2 E / q q ̇ ] = 0, then arbitrarily small perturbations q 1 ( t ) = q ( t ) q 0 from q0 will satisfy
    (17)

    This equation predicts stable oscillations with frequency ω when ω 2 > 0, and instability with growth-rate γ = | ω | when ω 2 < 0. If all derivatives of E ( q , q ̇ ) higher than second-order vanish, then the approximation in this equation is exact (i.e., an equality) for perturbations of any size.

Further generalisations of this method to account for damping and forcing effects are discussed in Sec. VII.

To demonstrate how our perspective emphasises the ubiquity of harmonic oscillations throughout classical physics, we now apply the method described above to three paradigmatic examples (see Fig. 2): a mass-spring system, an LC-circuit, and a simple pendulum. Each of these systems can, of course, be analysed readily using conventional methods, including more direct energy arguments (see Sec. VIII); however, our purpose in this section is to illustrate the universality of Eq. (17) by using familiar examples. The main practical power of our approach as a way of simplifying problems involving variable inertia is considered in Sec. VI.

Fig. 2.

Three classical oscillators: (a) a mass m on a spring of spring constant k with damping coefficient c; (b) an LCR-circuit with voltage source V S, and current I; and (c) a simple rigid pendulum with rod length l, bob-mass m, and gravity g.

Fig. 2.

Three classical oscillators: (a) a mass m on a spring of spring constant k with damping coefficient c; (b) an LCR-circuit with voltage source V S, and current I; and (c) a simple rigid pendulum with rod length l, bob-mass m, and gravity g.

Close modal

Figure 2(a) depicts a mass-spring system comprising a mass m hanging under gravity g on a spring with spring constant k. We analyse the undamped (c = 0) behaviour of this system using the method described in Sec. IV.

  1. The total energy of the system E ( x , x ̇ ) is given by
    (18)

    where x(t) is the displacement and x ̇ ( t ) = d x / d t.

  2. The stationary points of E ( x , x ̇ ) occur when
    (19)

    that is, at ( x , x ̇ ) = ( x 0 , 0 ), where x 0 = m g / k.

  3. The second-derivatives of E at ( x 0 , 0 ) are
    (20)
    with 2 E / x x ̇ = 0. Hence, by Eq. (17), arbitrarily small perturbations x1 oscillate harmonically about the equilibrium x 0 = m g / k according to
    (21)

    is the oscillation frequency. Since all derivatives of E ( x , x ̇ ) higher than second-order vanish, this approximation is exact for perturbations of any size. We consider the effect of damping ( c 0) in Sec. VII.

Figure 2(b) depicts an LCR-circuit formed by connecting an inductor L, capacitor C, and resistor R in series with a voltage source Vs. We analyse the undamped (R = 0), un-driven ( V s = 0 ) behaviour of this system using the method described in Sec. IV.

  1. The total energy of the system E ( Q , Q ̇ ) is given by
    (22)

    where Q(t) is the charge on the capacitor at time t, and Q ̇ ( t ) = d Q / d t = I is the current through the inductor.

  2. The stationary points of E ( Q , Q ̇ ) occur when
    (23)

    that is, at the point ( Q , Q ̇ ) = ( Q 0 , 0 ) = ( 0 , 0 ).

  3. The second-derivatives of E at ( Q 0 , 0 ) = ( 0 , 0 ) are
    (24)
    with 2 E / Q Q ̇ = 0. Hence, by Eq. (17), arbitrarily small perturbations Q1 oscillate about ( Q 0 , 0 ) = ( 0 , 0 ) according to
    (25)

    is the oscillation frequency. As we found with the mass-spring system, this approximation is exact for perturbations of any size because all derivatives of E ( Q , Q ̇ ) higher than second-order vanish. The damped ( R 0), driven ( V s 0) version of this system is analyzed in Sec. VII.

Figure 2(c) depicts a simple pendulum formed by a point bob of mass m suspended from a pivot by a light rod of length l in a uniform gravitational field. The pendulum is constrained to move in a vertical plane with θ as the angle between the rod and the direction of gravitational acceleration g. We analyse the behaviour of this system using the method described in Sec. IV.

  1. The total energy of the system E ( θ , θ ̇ ) is given by
    (26)

    where θ ( t ) lies in the range 0 θ < 2 π.

  2. The stationary points of E ( θ , θ ̇ ) occur when
    (27)

    There are therefore two stationary points ( θ , θ ̇ ) = ( θ 0 , 0 ), one with θ 0 = 0 and the other with θ 0 = π.

  3. At a given stationary point ( θ 0 , 0 ), the second-derivatives of E are
    (28)
    with 2 E / θ θ ̇ = 0. Hence, by Eq. (17), arbitrarily small perturbations θ1 to the equilibrium θ0 satisfy
    (29)

There are thus two modes of solution. When the pendulum is hanging downwards from its pivot, i.e., at the stationary point ( θ 0 , 0 ) = ( 0 , 0 ), we have ω 2 = g / l > 0, yielding stable oscillations with frequency ω = g / l. Conversely, when the pendulum is balanced vertically above its pivot, i.e., at the point ( θ 0 , 0 ) = ( π , 0 ), we have ω 2 = g / l < 0, yielding instability with rate γ = g / l.

The fixed-inertia systems discussed in Sec. V are useful for emphasising the universality of Eq. (17). However, the true power of our approach is most apparent when it is used to simplify problems involving variable inertia, for example, multi-dimensional systems subject to constraints. In this section we illustrate such an application, and in so doing obtain new stability results for a constrained pendulum system; we call this system the rigid trapezoidal pendulum.7 

The basic form of the rigid trapezoidal pendulum is depicted in Fig. 3: two simple rigid pendula each of length l > 0 and bob mass m are suspended from fixed points separated by a horizontal distance a, and their bobs joined by a light rod of length b > 0. If b = a, then the equilibrium configuration of the system is rectangular and the constituent pendula oscillate in parallel without exchanging energy. Otherwise, if b a, then the equilibrium configuration resembles a trapezium, and the combined system oscillates with a “rocking” motion.28 

Fig. 3.

Rigid trapezoidal pendulum system. Schematics (a) and (b) show the ϕ 0 = θ 0 equilibrium configuration when either: (a) a > 0 (a to the left of b); or (b) a < 0 (a to the right of b). Schematics (c) and (d) show the corresponding configurations when the system is displaced from equilibrium.

Fig. 3.

Rigid trapezoidal pendulum system. Schematics (a) and (b) show the ϕ 0 = θ 0 equilibrium configuration when either: (a) a > 0 (a to the left of b); or (b) a < 0 (a to the right of b). Schematics (c) and (d) show the corresponding configurations when the system is displaced from equilibrium.

Close modal
Variants of this system have been considered independently by two authors: Ramachandran et al., who treated a system constrained by a light rod of fixed length b = l but allowed l < a < 3 l, and the present author, who treated a system constrained by a light string of variable length b > 0, which necessarily required b < a ( b + 2 l ).7,8 In this article we combine these two systems into a single model whereby a and b both vary subject to the geometric restriction that
(30)
Note that a can be either positive or negative; indeed, labelling the points of suspension as a and b (see Fig. 3), we define a > 0 if a is to the left of b, and a < 0 if a is to the right of b. The a = 0 case corresponds to a situation equivalent to that of a compound pendulum.29 
Observe in Fig. 3 that the two pendula define angles θ and ϕ with the vertical, where θ is the angular displacement of the first pendulum, and ϕ is that of the second. According to Pythagoras' theorem these angles are related according to
(31)
such that ϕ ( θ ) is an implicit function of θ, with
(32)

Notice here that we have assumed the system to be constrained within the equilibrium plane (i.e., the system is not a bifilar one, such as those studied by Cromer27), such that the motion of one pendulum is determined entirely by the motion of the other. In this way, the system has one degree of freedom only, and can be specified completely by the angular displacement θ ( t ) and the velocity θ ̇ ( t ).

We now analyse the rigid trapezoidal pendulum using the method described in Sec. IV.

  1. By combining the energies of the constituent pendula, the total energy E of the system is given by
    (33)

    where for our purposes we assert θ , ϕ [ π / 2 , + π / 2 ].

    Since ϕ ( θ ) varies implicitly with θ, this expression for the energy function may be written as
    (34)

    The energy function is thus of the form in Eq. (6) with inertia M that varies as M ( θ ) = m l 2 ( 1 + ( d ϕ / d θ ) 2 ).

  2. The stationary points of E ( θ , θ ̇ ) occur when both
    (35)
    Hence, for a static equilibrium ( θ , θ ̇ ) = ( θ 0 , 0 ), we require
    (36)
    In principle, Eq. (36) may be solved simultaneously with Eqs. (31) and (32) to yield multiple solutions for θ0; however, this procedure does not look easy in general. To make analytical progress, therefore, we consider the symmetric equilibrium depicted in Fig. 3, i.e.,
    (37)
    In this way Eq. (36) is satisfied straightforwardly with sin θ 0 = sin ϕ 0 and ϕ ( θ 0 ) = 1, where θ0 is found by taking the positive root of Eq. (31), viz.,
    (38)

    [The negative root solution sin θ 0 = ( a + b ) / 2 l is meaningless when ( a + b ) > 2 l and therefore discarded.] We focus our analysis on this symmetric equilibrium.

  3. At the symmetric equilibrium ( θ 0 , 0 ) defined by Eqs. (37) and (38), the second-derivatives of E are
    (39)
    (40)
    and [ 2 E / θ θ ̇ ] ( θ 0 , 0 ) = 0, where we used the fact that ϕ ( θ 0 ) = 1 and ϕ ( θ 0 ) = ( 2 a sin θ 0 / b cos θ 0 ). It therefore follows by Eq. (17) that arbitrarily small perturbations θ1 to the equilibrium θ0 satisfy
    (41)
    where the frequency ω is given by
    (42)

This expression for ω 2 extends the present author's previous analysis to the range θ 0 < 0 (see Fig. 4) and reduces to that of Ramachandran et al. when b = l.7,8 As expected, we recover the result for a simple pendulum if a = b, when θ 0 = 0 and ω 2 = g / l (cf. Sec. V).

Fig. 4.

Plots of ω 2 with sin θ 0 as given by Eq. (42), for five values of b/l. The curves intercept the ω 2 / ( g / l ) = 0 axis at the critical angles θ c ( b , l ) defined by Eq. (43).

Fig. 4.

Plots of ω 2 with sin θ 0 as given by Eq. (42), for five values of b/l. The curves intercept the ω 2 / ( g / l ) = 0 axis at the critical angles θ c ( b , l ) defined by Eq. (43).

Close modal
For configurations satisfying b > 2 l, Eq. (42) predicts ω 2 > 0, in which case the system oscillates harmonically about the equilibrium. However, if b 2 l, then ω vanishes entirely ( ω = 0 ) when
(43)
where θc defines a critical angle. In these cases perturbations to the symmetric equilibrium configuration are stable ( ω 2 > 0) provided that θ 0 > θ c, but are unstable ( ω 2 < 0) if θ 0 < θ c (see Fig. 4). The solutions corresponding to stable oscillations were verified experimentally in our earlier context for configurations with a > b ( sin θ 0 > 0).7,28 Unstable solutions imply transition to an asymmetric equilibrium configuration yet to be fully determined.8 

It is important to emphasise that the method introduced here is substantially simpler than those adopted by previous authors, but nonetheless yields more general results.7,8 For example, the study by Ramachandran et al. is formulated in terms of three coupled Euler–Lagrange equations and two Lagrange multipliers, which are then linearised to yield a 5 × 5 matrix eigenvalue problem.8 In contrast, our method requires mathematics no more advanced than partial differentiation. Indeed, the fundamental power of our approach comes from the linear assumptions built into Eq. (17), which means that the hard work of formulating a customised linear model is avoided entirely.7,8 This property of the method is quite general, and applies to a large class of other problems in classical physics.1–3,5,10,18,20

Since the discussion so far has assumed un-driven, lossless systems, it remains to comment on the extent to which our approach can be generalised to include damping effects and driving forces. To do this, let us consider as before a system having total energy E ( q , q ̇ ) and a stationary point at ( q , q ̇ ) = ( q 0 , 0 ) satisfying [ 2 E / q q ̇ ] = 0. If this stationary point is perturbed by supplying a small amount of energy E 1, then the analysis leading to Eq. (10) goes through as in Sec. II, and arbitrary small perturbations q 1 ( t ) to q0 satisfy
(44)
In Sec. II, we assumed constant energy ( E 1 ̇ = 0); however, if the system is instead subject to damping effects and driving forces, then the rate-of-change of E 1 is the difference between the power supplied due to forcing E ̇ F and the power dissipated due to damping E ̇ D, i.e.,
(45)
Hence, differentiating Eq. (44) with respect to time, we have
(46)
This expression may be rearranged in the form of what we refer to as the universal driven oscillator equation
(47)
where the damping rate α and natural frequency ω are
(48)
respectively, and [ δ 2 q / δ t 2 ] F denotes the acceleration due to forcing, and is defined as
(49)
In the case of damped-driven systems, therefore, step 3 of our method in Sec. IV must be modified to:
  1. Evaluate the second derivatives of E ( q , q ̇ ) at the stationary point ( q 0 , 0 ). If these derivatives are non-zero, with [ 2 E / q q ̇ ] = 0, then arbitrarily small perturbations q1 from q0 will satisfy
    (50)

    where the damping rate α, frequency ω, and forcing [ δ 2 q / δ t 2 ] F are defined in Eqs. (48) and (49). As in Sec. IV, if all derivatives of E ( q , q ̇ ) higher than second-order vanish, then this approximation is exact (an equality) for perturbations of any size.

Notice that Eq. (50) is not in fact a linear differential equation unless E ̇ D q ̇ 1 2 and E ̇ F q ̇ 1, a property which reveals useful insights into the meaning of α and [ δ 2 q / δ t 2 ] F. Indeed, if the energy function E ( q , q ̇ ) takes the form of Eq. (6), then
(51)
where M ( q 0 ) is the inertia. Thus: (i) the damping rate α represents the ratio between the dissipated power E ̇ D and the kinetic energy; and (ii) the acceleration due to forcing [ δ 2 q / δ t 2 ] F represents the ratio between the supplied forcing power E ̇ F and the system momentum. With these ideas in mind, let us briefly illustrate how Eq. (50) works when applied to some canonical examples.
Figure 2(a) shows a mass m hanging under gravity g on a spring with spring constant k, coupled to a viscous damping pot with damping coefficient c. In Sec. V we found that for c = 0 (i.e., no damping), perturbations x1 to the system equilibrium x 0 = ( m g / k ) satisfy
(52)
where 2 E / x 2 = k and 2 E / x ̇ 2 = m are determined by the expression for the energy E ( x , x ̇ ) in Eq. (18). Recall that the approximation in this case is exact, since the higher-order derivatives of E ( x , x ̇ ) vanish (see Sec. IV).
If we now suppose that the system is damped ( c 0 ), then the viscous damping force is F v = c x ̇, such that the dissipative power E ̇ D = ( force × velocity ) is
(53)
where we note that the perturbed system has x ̇ = x ̇ 1. According to Eq. (50), therefore, in the presence of system damping Eq. (52) is modified to
(54)
where the damping rate α is given by Eq. (48), i.e.,
(55)
and where we have assumed no driving forces (i.e., E ̇ F = 0). As with Eq. (52), because all of the second derivatives of E ( x , x ̇ ) vanish, the approximation in Eq. (54) can be replaced with an equality (see Sec. IV). This discussion shows that routine application of Eq. (50) yields the standard equation of motion for a mass-spring system with viscous damping coefficient c.1,3
In Sec. V we considered an oscillating circuit formed by connecting an inductor L in series with a capacitor C, and showed that charge perturbations Q1 satisfy
(56)
where 2 E / Q 2 = 1 / C and 2 E / Q ̇ 2 = L are given by the expression for the energy E ( Q , Q ̇ ) in Eq. (22). As above, the approximation here is exact since the higher-order derivatives of E ( Q , Q ̇ ) vanish (see Sec. IV).
We now suppose that the system is modified by adding a resistor R and voltage source V S ( t ) in series to form an LCR-circuit, as depicted in Fig. 2(b). In this way the energy function E ( Q , Q ̇ ) is unchanged, but energy can be dissipated by the resistor R with damping power
(57)
and supplied by the voltage source V S ( t ) with power
(58)
where I = d Q / d t = Q ̇ 1 is the current, and we used the fact that the perturbed system has Q ̇ = Q ̇ 1.
According to Eq. (50), therefore, charge perturbations Q 1 ( t ) in a driven LCR-circuit satsify
(59)
where the damping rate α is
(60)
and we have written the forcing term [ δ 2 Q / δ t 2 ] F as
(61)
This is the standard equation for driven oscillations in an LCR-circuit.1,13 Because all of the second derivatives of E ( Q , Q ̇ ) vanish, this equation is exact (see Sec. IV) and can be used to describe large amplitude perturbations, such as those arising due to forcing with resonant power.14 
We have proposed a new perspective on the ubiquity of classical harmonic oscillators based on the two-variable Taylor expansion of a perturbed system's total energy E ( q , q ̇ ), where q(t) is the system displacement as a function of time t and q ̇ ( t ) = d q / d t (Sec. II). In particular, we have shown that if E ( q , q ̇ ) has a stationary point at some static equilibrium ( q 0 , 0 ), then perturbations q 1 ( t ) to this equilibrium satisfy the lossless oscillator equation
(62)
is a universal equation for the oscillation frequency.

Unlike the traditional explanation for the prevalence of simple harmonic oscillators—which is based on conservative forces and is limited to rectilinear mechanical systems with fixed inertia (see Eq. (2))—our perspective uses energy arguments only, making it more universal in at least three key respects. First, Eq. (62) is framed in terms of generalised coordinates ( q , q ̇ ) and applies to arbitrary geometries (e.g., rotational systems) without the need for ad hoc modifications. Second, Eq. (62) provides a formal framework for handling problems involving variable inertia (see, e.g., Sec. VI). And third, Eq. (62) is derived without reference to Newtonian forces, meaning that it encompasses a wider range of contexts, including electromagnetic systems (see Secs. V and VI).

Our analysis of lossless systems assumes a fixed energy perturbation E 1 satisfying E 1 ̇ = 0. In Sec. VII we relaxed this assumption and considered the possibility of time dependence in E 1 by setting E 1 ̇ = E ̇ F E ̇ D, where E ̇ F is the power delivery due to system forcing and E ̇ D is the power dissipation due to damping. In this way we extended our analysis to obtain a universal driven oscillator equation from energy arguments only, including an expression for the damping rate α in terms of E ̇ D. Preliminary investigation suggests that this approach may be helpful for analysing damped-driven systems when power inputs and power losses are well-characterised.

Equation (62) is essentially a linearised equation of motion for a system perturbed from equilibrium (Sec. III), and for this reason may also be used to study instability (Sec. IV). In particular, the linear assumptions built into Eq. (62) make it possible to determine stability conditions (and natural frequencies) without the need to form customised linear models. This feature is especially useful for simplifying problems involving variable inertia, such as multi-directional systems subject to constraints; for example, in Section VI we showed how equation (62) can be used to obtain new stability results for a constrained pendulum system, using mathematics no more advanced than partial differentiation.7,8

It should be emphasised that our method offers several advantages in comparison to more traditional energy techniques, such as a solution by the Euler–Lagrange equations or direct differentiation of the energy function E ( q , q ̇ ).3 For example, Lagrangian mechanics is typically studied at an advanced undergraduate level, whereas our method is based on first-year calculus, and can be used to tackle “complicated” problems at an earlier stage (see, e.g., Sec. VI C). Likewise, although direct differentiation of the energy function is effective for very simple systems (such as the mass-spring system and the LC-circuit described in Sec. V), in other contexts the resulting equation of motion is often difficult or impossible to solve.5,7,20 Indeed, without an additional linearisation process, differentiation of E ( q , q ̇ ) will usually only yield an oscillator equation when the inertia is fixed, and the potential and kinetic energy terms are proportional to q2 and q ̇ 2 respectively. In contrast, our method is universal because the two-variable Taylor expansion leading to Eq. (62) imposes parabolic approximations on both the potential energy and the kinetic energy from the outset (see Eq. (10)). These approximations are useful for the reasons discussed above (since they allow us to treat arbitrary potentials, and variable inertia problems without forming customised linear models); however, they also reveal a conceptual subtlety that implies the conventional explanation for the ubiquity of classical harmonic oscillators is incomplete. In particular, our analysis suggests that harmonic oscillations arise in classical systems near equilibrium not because the potential energy function V(q) is locally parabolic, but because the total energy function E ( q , q ̇ ) is locally paraboloidal.

In summary, we have used the two-variable Taylor expansion of a system's total energy E ( q , q ̇ ) to present a fresh perspective on the ubiquity of harmonic oscillators in classical physics, and in so doing we have derived a universal equation for the natural oscillation frequency ω = ( 2 E / q 2 ) / ( 2 E / q ̇ 2 ). Crucially, this perspective yields a routine method for characterising the stability of perturbed systems, and can render complicated problems amenable to analysis by elementary methods. Indeed, we look forward to exploring further applications of our approach to a range of new contexts in future publications.

1.
H. J. J.
Braddick
,
Vibrations, Waves, and Diffraction
(
McGraw-Hill
,
New York
,
1965
).
2.
T. W. B.
Kibble
,
Classical Mechanics
(
McGraw-Hill
,
New York
,
1966
).
3.
M.
Lunn
,
A First Course in Mechanics
(
Oxford U. P
.,
Oxford
,
2004
).
4.
R.
Hauko
,
D.
Andreevski
,
D.
Paul
,
M.
Šterk
, and
R.
Repnik
, “
Teaching of the harmonic oscillator damped by a constant force: The use of analogy and experiments
,”
Am. J. Phys.
86
(
9
),
657
662
(
2018
).
5.
M. J.
O'Shea
, “
Harmonic and anharmonic behaviour of a simple oscillator
,”
Eur. J. Phys.
30
,
549
558
(
2009
).
6.
G. L.
Baker
and
J. A.
Blackburn
,
The Pendulum: A Case Study in Physics
(
Oxford U. P
.,
Oxford
,
2005
).
7.
J. J.
Bissell
and
S. K.
Bhamidimarri
, “
Introducing the trapezoidal pendulum: Dynamics of coupled pendula suitable for distance teaching
,”
Phys. Educ.
55
,
065008
(
2020
).
8.
P.
Ramachandran
,
S. G.
Krishna
, and
Y. M.
Ram
, “
Instability of a constrained pendulum system
,”
Am. J. Phys.
79
(
4
),
395
400
(
2011
).
9.
D. T.
Gillespie
, “
Simple harmonic motion of a round body rolling on a concave curve
,”
Am. J. Phys.
52
(
2
),
180
182
(
1984
).
10.
M. J.
Mazzoleni
,
M. B.
Krone
, and
B. P.
Mann
, “
Dynamics of rocking semicircular, parabolic, and semi-elliptical disks: Equilibria, stability, and natural frequencies
,”
J. Vibration Acoust.
137
(
4
),
041017
(
2015
).
11.
J. T.
Wilkinson
,
C. B.
Whitehouse
,
R. F.
Oulton
, and
S. D.
Gennaro
, “
An undergraduate experiment demonstrating the physics of metamaterials with acoustic waves and soda cans
,”
Am. J. Phys.
84
,
14
20
(
2016
).
12.
S. S.
Ivković
,
M. Z.
Marković
,
D. Ž.
Ivković
, and
N.
Cvetanovic̀
, “
LCR circuit: New simple methods for measuring the equivalent series resistance of a capacitor and inductance of a coil
,”
Eur. J. Phys.
38
,
055705
(
2017
).
13.
M. C.
Faleski
, “
Transient behavior of the driven RLC circuit
,”
Am. J. Phys.
74
,
429
437
(
2006
).
14.
F. B.
Kneubil
, “
Driven series RLC circuit and resonance: A graphic approach to energy
,”
Phys. Teach.
58
,
256–259
(
2020
).
15.
O.
Isaksson
,
M.
Karlsteen
,
M.
Rostedt
, and
D.
Hanstorp
, “
An optical levitation system for a physics teaching laboratory
,”
Am. J. Phys.
86
(
2
),
135
142
(
2018
).
16.
J. T.
Marmolejo
,
O.
Isaksson
,
R.
Cabrera-Trujillo
,
N. C.
Giesselmann
, and
D.
Hanstorp
, “
A fully manipulable damped driven harmonic oscillator using optical levitation
,”
Am. J. Phys.
88
(
6
),
490
498
(
2020
).
17.
R.
Marshall
, “
Seek it here, find it there, shm is everywhere
,”
Phys. Educ.
49
(
2
),
255
258
(
2014
).
18.
S.
Kryhin
and
O.
Orlyansky
, “
Different approaches and methods for solving oscillator problems
,”
Eur. J. Phys.
40
,
025006
(
2019
).
19.
S. H.
Strogatz
,
D. M.
Abrams
,
A.
McRobie
,
B.
Eckhardt
, and
E.
Ott
, “
Crowd synchrony on the Millennium Bridge
,”
Nature
438
,
43
44
(
2005
).
20.
T.
Troy
,
M.
Reiner
,
A. J.
Haugen
, and
N. T.
Moore
, “
Small oscillations by conservation of energy
,”
Phys. Educ.
52
,
065009
(
2017
).
21.
R. B.
Abbott
and
C. G.
Fry
, “
Acoustical, mechanical and electrical analogies
,”
Am. J. Phys.
5
,
166
167
(
1937
).
22.
A.
Bloch
, “
Electromechanical analogies and their use for the analysis of mechanical and electromechanical systems
,”
J. Inst. Electr. Eng. - Part I: General
92
(
52
),
157
169
(
1945
).
23.
K. F.
Riley
,
M. P.
Hobson
, and
S. J.
Bence
,
Mathematical Methods for Physics and Engineering
, 3rd ed. (
Cambridge U. P
.,
Cambridge
,
2006
).
24.
The second-order Taylor expansion of E ( q , q ̇ ) has remainder term R 3 ( q , q ̇ ) given by (Ref. 23) R 3 ( q , q ̇ ) = 1 3 ! ( q 1 3 3 E q 3 + 3 q 1 2 q ̇ 1 3 E q 2 q ̇ + 3 q 1 q ̇ 1 2 3 E q q ̇ 2 + q ̇ 1 3 3 E q ̇ 3 ) ,with derivatives evaluated at ( q , q ̇ ) = ( η , ζ ), where η [ q 0 , q 0 + q 1 ] and ζ [ 0 , q ̇ 1 ]. Provided q ̇ 1 O ( q 1 / τ ), where τ is a finite timescale (typically of order 1 / ω), the terms in this expression scale as q 1 3; in this way R3 may be neglected from Eq. (9) as q 1 0 when compared to terms in q 1 2. Notice that if all derivatives of E ( q , q ̇ ) higher than second-order vanish, then R 3 ( q , q ̇ ) = 0, and Eq. (10) is exact for q1 of any size.
25.
As expected, Eq. (11) reduces to Eq. (2) for rectilinear, fixed-mass energy functions of the form E ( x , x ̇ ) = V ( x ) + 1 2 m x ̇ 2.
26.
The condition ω 2 > 0 can, in principle, correspond to a maximum of E ( q , q ̇ ) at ( q 0 , 0 ), but this would require negative inertia.
27.
A.
Cromer
, “
Many oscillations of a rigid rod
,”
Am. J. Phys.
63
,
112
121
(
1995
).
28.
J. J.
Bissell
and
S. K.
Bhamidimarri
, YouTube video of the trapezoidal pendulum available at <https://youtu.be/GZAkr-OyTbM>.
29.
T. H.
Richardson
and
S. A.
Brittle
, “
Physical pendulum experiments to enhance the understanding of moments of inertia and simple harmonic motion
,”
Phys. Educ.
47
,
537–544
(
2012
).