A new perspective on the ubiquity of classical harmonic oscillators is presented based on the twovariable Taylor expansion of a perturbed system's total energy $ E ( q , q \u0307 )$, where q(t) is the system displacement as a function of time t and $ q \u0307 ( 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 $ \omega = ( \u2202 2 E / \u2202 q 2 ) / ( \u2202 2 E / \u2202 q \u0307 2 )$ which may be applied to arbitrary systems without the need to form systemspecific 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.
I. INTRODUCTION
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 bodies^{9,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}
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 multidirectional geometries (such as rotational systems); second (ii), the argument leading to equation (2) requires fixed m, and breaksdown 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}
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 twovariable Taylor expansion of the energy function $ E ( q , q \u0307 )$—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 $ \omega = ( \u2202 2 E / \u2202 q 2 ) / ( \u2202 2 E / \u2202 q \u0307 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.
II. UNIVERSAL LOSSLESS OSCILLATOR EQUATION
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).
III. APPLICATIONS TO LINEAR STABILITY
IV. METHOD FOR ANALYSING STATIONARY POINTS
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:
This theorem motivates the following threestep method for analysing free systems perturbed from equilibrium:

Write down an expression for the total energy $ E ( q , q \u0307 )$ in terms of q(t) and $ q \u0307 ( t ) = d q / d t$.
 Determine the stationary points $ ( q , q \u0307 ) = ( q 0 , 0 )$ of the energy function $ E ( q , q \u0307 )$ by solving$ [ \u2202 E \u2202 q ] ( q 0 , 0 ) = 0 \u2003 and \u2003 [ \u2202 E \u2202 q \u0307 ] ( q 0 , 0 ) = 0.$
 Evaluate the second derivatives of $ E ( q , q \u0307 )$ at the stationary points $ ( q 0 , 0 )$. If these derivatives are nonzero, with $ [ \u2202 2 E / \u2202 q \u2202 q \u0307 ] = 0$, then arbitrarily small perturbations $ q 1 ( t ) = q ( t ) \u2212 q 0$ from q_{0} will satisfy$ q \u0308 1 + \omega 2 q 1 \u2248 0 , \u2003 where \u2003 \omega 2 = [ \u2202 2 E / \u2202 q 2 \u2202 2 E / \u2202 q \u0307 2 ] .$
This equation predicts stable oscillations with frequency ω when $ \omega 2 > 0$, and instability with growthrate $ \gamma =  \omega $ when $ \omega 2 < 0$. If all derivatives of $ E ( q , q \u0307 )$ higher than secondorder 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.
V. PARADIGMATIC SYSTEMS WITH FIXED INERTIA
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 massspring system, an LCcircuit, 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.
A. Massspring system
Figure 2(a) depicts a massspring 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.
 The total energy of the system $ E ( x , x \u0307 )$ is given by$ E ( x , x \u0307 ) = 1 2 k x 2 \ufe38 elastic \u2009 energy \u2212 mgx \ufe38 gravitational \u2009 energy + 1 2 m x \u0307 2 \ufe38 kinetic \u2009 energy ,$
where x(t) is the displacement and $ x \u0307 ( t ) = d x / d t$.
 The stationary points of $ E ( x , x \u0307 )$ occur when$ \u2202 E \u2202 x = k x \u2212 m g = 0 , \u2003 and \u2003 \u2202 E \u2202 x \u0307 = m x \u0307 = 0 ,$
that is, at $ ( x , x \u0307 ) = ( x 0 , 0 )$, where $ x 0 = m g / k$.
 The secondderivatives of $E$ at $ ( x 0 , 0 )$ are$ \u2202 2 E \u2202 x 2 = k , \u2003 and \u2003 \u2202 2 E \u2202 x \u0307 2 = m ,$with $ \u2202 2 E / \u2202 x \u2202 x \u0307 = 0$. Hence, by Eq. (17), arbitrarily small perturbations x_{1} oscillate harmonically about the equilibrium $ x 0 = m g / k$ according to$ x \u0308 1 + \omega 2 x 1 \u2248 0 , \u2003 where \u2003 \omega = \u2202 2 E / \u2202 x 2 \u2202 2 E / \u2202 x \u0307 2 = k m$
is the oscillation frequency. Since all derivatives of $ E ( x , x \u0307 )$ higher than secondorder vanish, this approximation is exact for perturbations of any size. We consider the effect of damping ( $ c \u2260 0$) in Sec. VII.
B. Oscillating LCcircuit
Figure 2(b) depicts an LCRcircuit formed by connecting an inductor L, capacitor C, and resistor R in series with a voltage source V_{s}. We analyse the undamped (R = 0), undriven $ ( V s = 0 )$ behaviour of this system using the method described in Sec. IV.
 The total energy of the system $ E ( Q , Q \u0307 )$ is given by$ E ( Q , Q \u0307 ) = 1 2 C Q 2 \ufe38 capacitor \u2009 energy + 1 2 L Q \u0307 2 \ufe38 inductor \u2009 energy ,$
where Q(t) is the charge on the capacitor at time t, and $ Q \u0307 ( t ) = d Q / d t = I$ is the current through the inductor.
 The stationary points of $ E ( Q , Q \u0307 )$ occur when$ \u2202 E \u2202 Q = Q C = 0 , \u2003 and \u2003 \u2202 E \u2202 Q \u0307 = L Q \u0307 = 0 ,$
that is, at the point $ ( Q , Q \u0307 ) = ( Q 0 , 0 ) = ( 0 , 0 )$.
 The secondderivatives of $E$ at $ ( Q 0 , 0 ) = ( 0 , 0 )$ are$ \u2202 2 E \u2202 Q 2 = 1 C , \u2003 and \u2003 \u2202 2 E \u2202 Q \u0307 2 = L ,$with $ \u2202 2 E / \u2202 Q \u2202 Q \u0307 = 0$. Hence, by Eq. (17), arbitrarily small perturbations Q_{1} oscillate about $ ( Q 0 , 0 ) = ( 0 , 0 )$ according to$ Q \u0308 1 + \omega 2 Q 1 \u2248 0 , \u2003 where \u2003 \omega = \u2202 2 E / \u2202 Q 2 \u2202 2 E / \u2202 Q \u0307 2 = 1 L C$
is the oscillation frequency. As we found with the massspring system, this approximation is exact for perturbations of any size because all derivatives of $ E ( Q , Q \u0307 )$ higher than secondorder vanish. The damped ( $ R \u2260 0$), driven ( $ V s \u2260 0$) version of this system is analyzed in Sec. VII.
C. Simple rigid pendulum
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.
 The total energy of the system $ E ( \theta , \theta \u0307 )$ is given by$ E ( \theta , \theta \u0307 ) = \u2212 mgl \u2009 cos \u2009 \theta \ufe38 potential \u2009 energy + 1 2 m l 2 \theta \u0307 2 \ufe38 kinetic \u2009 energy ,$
where $ \theta ( t )$ lies in the range $ 0 \u2264 \theta < 2 \pi $.
 The stationary points of $ E ( \theta , \theta \u0307 )$ occur when$ \u2202 E \u2202 \theta = mgl \u2009 sin \u2009 \theta = 0 , \u2003 and \u2003 \u2202 E \u2202 \theta \u0307 = m l 2 \theta \u0307 = 0.$
There are therefore two stationary points $ ( \theta , \theta \u0307 ) = ( \theta 0 , 0 )$, one with $ \theta 0 = 0$ and the other with $ \theta 0 = \pi $.
 At a given stationary point $ ( \theta 0 , 0 )$, the secondderivatives of $E$ are$ \u2202 2 E \u2202 \theta 2 = mgl \u2009 cos \u2009 \theta 0 , \u2003 and \u2003 \u2202 2 E \u2202 \theta \u0307 2 = m l 2 ,$with $ \u2202 2 E / \u2202 \theta \u2202 \theta \u0307 = 0$. Hence, by Eq. (17), arbitrarily small perturbations θ_{1} to the equilibrium θ_{0} satisfy$ \theta \u0308 1 + \omega 2 \theta 1 \u2248 0 , \u2003 with \u2003 \omega 2 = \u2202 2 E / \u2202 \theta 2 \u2202 2 E / \u2202 \theta \u0307 2 = g \u2009 cos \u2009 \theta 0 l .$
There are thus two modes of solution. When the pendulum is hanging downwards from its pivot, i.e., at the stationary point $ ( \theta 0 , 0 ) = ( 0 , 0 )$, we have $ \omega 2 = g / l > 0$, yielding stable oscillations with frequency $ \omega = g / l$. Conversely, when the pendulum is balanced vertically above its pivot, i.e., at the point $ ( \theta 0 , 0 ) = ( \pi , 0 )$, we have $ \omega 2 = \u2212 g / l < 0$, yielding instability with rate $ \gamma = g / l$.
VI. NOVEL SYSTEM WITH VARIABLE INERTIA
The fixedinertia 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, multidimensional 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}
A. Model for the rigid trapezoidal pendulum
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 \u2260 a$, then the equilibrium configuration resembles a trapezium, and the combined system oscillates with a “rocking” motion.^{28}
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 Cromer^{27}), 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 $ \theta ( t )$ and the velocity $ \theta \u0307 ( t )$.
B. Stability analysis
We now analyse the rigid trapezoidal pendulum using the method described in Sec. IV.
 By combining the energies of the constituent pendula, the total energy $E$ of the system is given by$ E = \u2212 mgl ( cos \u2009 \theta + cos \u2009 \varphi ) \ufe38 potential \u2009 energy + 1 2 m l 2 ( \theta \u0307 2 + \varphi \u0307 2 ) \ufe38 kinetic \u2009 energy ,$
where for our purposes we assert $ \theta , \varphi \u2208 [ \u2212 \pi / 2 , + \pi / 2 ]$.
Since $ \varphi ( \theta )$ varies implicitly with θ, this expression for the energy function may be written as$ E ( \theta , \theta \u0307 ) = \u2212 mgl ( cos \u2009 \theta + cos \u2009 \varphi ( \theta ) ) + 1 2 m l 2 ( 1 + ( \varphi \u2032 ) 2 ) \theta \u0307 2 .$The energy function is thus of the form in Eq. (6) with inertia M that varies as $ M ( \theta ) = m l 2 ( 1 + ( d \varphi / d \theta ) 2 )$.
 The stationary points of $ E ( \theta , \theta \u0307 )$ occur when both$ \u2202 E \u2202 \theta = mgl ( sin \u2009 \theta + \varphi \u2032 \u2009 sin \u2009 \varphi ) + m l \varphi \u2032 \varphi \u2033 \theta \u0307 2 = 0 , and \u2003 \u2202 E \u2202 \theta \u0307 = m l 2 ( 1 + ( \varphi \u2032 ) 2 ) \theta \u0307 = 0.$Hence, for a static equilibrium $ ( \theta , \theta \u0307 ) = ( \theta 0 , 0 )$, we require$ [ \u2202 E \u2202 \theta ] ( \theta 0 , 0 ) = mgl ( sin \u2009 \theta 0 + \varphi \u2032 ( \theta 0 ) \u2009 sin \u2009 \varphi ( \theta 0 ) ) = 0.$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.,$ \theta = \theta 0 = \varphi 0 = \varphi ( \theta 0 ) .$In this way Eq. (36) is satisfied straightforwardly with $ sin \u2009 \theta 0 = sin \u2009 \varphi 0$ and $ \varphi \u2032 ( \theta 0 ) = \u2212 1$, where θ_{0} is found by taking the positive root of Eq. (31), viz.,$ \u2212 1 \u2264 sin \u2009 \theta 0 = ( a \u2212 b ) 2 l \u2264 1.$
[The negative root solution $ sin \u2009 \theta 0 = ( a + b ) / 2 l$ is meaningless when $ ( a + b ) > 2 l$ and therefore discarded.] We focus our analysis on this symmetric equilibrium.
 At the symmetric equilibrium $ ( \theta 0 , 0 )$ defined by Eqs. (37) and (38), the secondderivatives of $E$ are$ [ \u2202 2 E \u2202 \theta 2 ] ( \theta 0 , 0 ) = mgl ( \varphi \u2033 ( \theta 0 ) \u2009 sin \u2009 \varphi + ( 1 + \varphi \u2032 ( \theta 0 ) 2 ) \u2009 cos \u2009 \theta 0 ) = mgl ( 2 b + 4 l \u2009 sin \u2009 3 \theta 0 ) / b \u2009 cos \u2009 \theta 0 ,$$ [ \u2202 2 E \u2202 \theta \u0307 2 ] ( \theta 0 , 0 ) = m l 2 ( 1 + \varphi \u2032 ( \theta 0 ) 2 ) = 2 m l 2 ,$and $ [ \u2202 2 E / \u2202 \theta \u2202 \theta \u0307 ] ( \theta 0 , 0 ) = 0$, where we used the fact that $ \varphi \u2032 ( \theta 0 ) = \u2212 1$ and $ \varphi \u2033 ( \theta 0 ) = ( 2 a \u2009 sin \u2009 \theta 0 / b \u2009 cos \u2009 \theta 0 )$. It therefore follows by Eq. (17) that arbitrarily small perturbations θ_{1} to the equilibrium θ_{0} satisfy$ \theta \u0308 1 + \omega 2 \theta 1 \u2248 0 ,$where the frequency ω is given by$ \omega 2 = [ \u2202 2 E / \u2202 \theta 2 \u2202 2 E / \u2202 \theta \u0307 2 ] ( \theta 0 , 0 ) = g l ( b + 2 l \u2009 sin 3 \theta 0 b \u2009 cos \u2009 \theta 0 ) .$
This expression for $ \omega 2$ extends the present author's previous analysis to the range $ \theta 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 $ \theta 0 = 0$ and $ \omega 2 = g / l$ (cf. Sec. V).
C. Remarks on the analysis
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}
VII. DAMPING EFFECTS AND DRIVING FORCES
 Evaluate the second derivatives of $ E ( q , q \u0307 )$ at the stationary point $ ( q 0 , 0 )$. If these derivatives are nonzero, with $ [ \u2202 2 E / \u2202 q \u2202 q \u0307 ] = 0$, then arbitrarily small perturbations q_{1} from q_{0} will satisfy$ q \u0308 1 + 2 \alpha q \u0307 1 + \omega 2 q 1 \u2248 [ \delta 2 q \delta t 2 ] F ,$
where the damping rate α, frequency ω, and forcing $ [ \delta 2 q / \delta t 2 ] F$ are defined in Eqs. (48) and (49). As in Sec. IV, if all derivatives of $ E ( q , q \u0307 )$ higher than secondorder vanish, then this approximation is exact (an equality) for perturbations of any size.
A. Damped massspring system
B. Driven LCRcircuit
VIII. CONCLUSION
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 \u0307 )$ 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 \u0307 = 0$. In Sec. VII we relaxed this assumption and considered the possibility of time dependence in $ E 1$ by setting $ E 1 \u0307 = E \u0307 F \u2212 E \u0307 D$, where $ E \u0307 F$ is the power delivery due to system forcing and $ E \u0307 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 \u0307 D$. Preliminary investigation suggests that this approach may be helpful for analysing dampeddriven systems when power inputs and power losses are wellcharacterised.
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 multidirectional 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 \u0307 )$.^{3} For example, Lagrangian mechanics is typically studied at an advanced undergraduate level, whereas our method is based on firstyear 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 massspring system and the LCcircuit 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 \u0307 )$ will usually only yield an oscillator equation when the inertia is fixed, and the potential and kinetic energy terms are proportional to q^{2} and $ q \u0307 2$ respectively. In contrast, our method is universal because the twovariable 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 \u0307 )$ is locally paraboloidal.
In summary, we have used the twovariable Taylor expansion of a system's total energy $ E ( q , q \u0307 )$ 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 $ \omega = ( \u2202 2 E / \u2202 q 2 ) / ( \u2202 2 E / \u2202 q \u0307 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.