An approximate solution is presented for simple harmonic motion in the presence of damping by a force which is a general power-law function of the velocity. The approximation is shown to be quite robust, allowing for a simple way to investigate amplitude decay in the presence of a general type of weak, nonlinear damping.
Students and instructors might wonder why linear damping is the only case investigated in detail in most texts on advanced, undergraduate-level, classical mechanics.2–4 The exclusive focus of most authors on linear, viscous damping is likely due to the fact that Eq. (1) is a linear equation—and therefore solvable in terms of elementary functions—only for the case α = 1. It is worth noting that other values of α also lead to physically plausible dissipative forces, since α = 0 represents a constant force of “dry,” sliding friction5–8 while α = 2 corresponds to turbulent air resistance.9,10 Additionally, a cubic drag term, which has attracted significant attention as a nontrivial mathematical problem,11 has appeared as an important feature in studies of energy-harvesting devices and cochlear implants.12 Large values of have also been explored as drag terms for the nonlinear pendulum.13
As no fundamental law forbids noninteger values of α, it is perhaps unsurprising that noninteger values of α also find applications. The effective drag force on objects moving through fluids9 sometimes obeys such a nonintegral power law experimentally.14 Fractional power-law functions also appear as relaxation functions in certain models of viscoelasticity15 relevant to the study of biological tissues. Such systems are substantially more complicated than harmonic oscillators. Still, aside from the purely mathematical curiosity of nonlinear damping, intuition gleaned from investigating an oscillator with nonlinear damping could be useful to students who wish to explore viscoelasticity and its role in biomechanics. The aim of this note is to summarize a unified, analytic treatment of the asymptotic behavior predicted by Eq. (1) for all possible values of α. Additionally, we consider briefly the situation in which more than one damping term is present. The general approach employed throughout this note serves as a useful example of how one may extract meaningful predictions from analytically intractable systems by using asymptotic methods. Few interesting problems in physics are exactly solvable, and students should benefit from early exposure to such methods in a familiar setting such as the harmonic oscillator.
However, even if this weak-damping condition is satisfied, Eq. (16) cannot be satisfied for long times since the denominator shrinks with increasing time. Thus for , the asymptotic results for A(t) are valid for the weak-damping regime at short times. Figure 1 depicts the gradual breakdown for the case α = 0. For certain parameter choices, it is possible for the sliding friction to overwhelm the restoring force of the spring when the amplitude has decayed to a sufficiently small value. Generally, one has for static (μs) and kinetic (μk) coefficients of friction, so this behavior should be treated more carefully for more physically correct results at long times.20 Effectively, the model considered treats , and our main interest is in demonstrating the validity of the asymptotic result at short times. It should be noted that examples in which the mass comes to rest away from equilibrium have been published for exact solutions6 as well as experimental realizations.7,8 For and modest values of κ, the agreement between numerical and approximate solutions at long times is generally quite good, as shown in Fig. 2. A Python implementation of the numerical approach used to generate these figures is included in the supplementary material.21
A few general features can be gleaned from Eq. (22). First, the long-time decay envelope controlled by the exponential factor which also describes the case of only linear damping. Rather than some crossover/intermediate behavior, we find the long-time decay is dominated by the term with the smallest value of α. So far we have made no restrictions whatsoever on the relative sizes of b and .23 In the limit of negligible quadratic damping, , the second term in the denominator may be neglected so one recovers the standard result for linear damping, . More interesting is the case that so that the quadratic drag is “stronger.” In this limit . While the decay at long times is still dominated by the exponential envelope characteristic of linear drag, the dependence of the instantaneous amplitude on its initial value x0 disappears in the long-time limit, as in the case of purely quadratic drag.19 Comparison of the asymptotics to a numerical integration of the equation of motion is shown in Fig. 3. A Python implementation of the numerical approach is available in the supplementary material.21 Due to the exponential decay at long times, it should be emphasized that the argument that terms with is nonphysical13 due to oscillations which decay too slowly only holds when a single damping term with is present. Additionally, it is known that distinguishing between (for example) linear and quadratic functions as the dominant drag term from simulated experimental data is not always possible.24 Experimentally, one might expect the difficulty in determining the drag force to be compounded by the potential existence of more than one power-law term, perhaps with one much “weaker” than the other. It might prove worthwhile to explore the feasibility of using Eq. (21) to fit the relative amplitudes rather than distinguishing between two pure, power laws.
The methods described in this note allow one to investigate the amplitude decay in oscillating systems with any type of power-law damping. This approach provides a concise theoretical complement to the experimental investigation of nonlinear damping in oscillating systems. In particular, one could create a fairly rich lab experiment with strong theoretical and computational components by using the TeachSpin Torsional Oscillator25 which offers ways of producing linear (α = 1), quadratic (α = 2), and sliding (α = 0) drag in a mechanical system with a nearly linear restoring force. Additionally, the TeachSpin apparatus allows one to combine these types of drag, and the results of this note could form the basis of a theoretical model for such a setup which has no exact solutions. Finally, the availability of affordable and accurate accelerometers allows one to collect detailed data in oscillating systems which can be analyzed to reconstruct the effective drag force.26 The results outlined in this note could provide a simplified theoretical component for interpreting such experimental data.
ACKNOWLEDGMENTS
The author thanks Robert Hauko and the anonymous reviewers for the constructive criticisms and valuable insights. The author also thanks Joseph M. Starobin for the illuminating discussions and for introducing him to Ref. 19.
REFERENCES
The issue of the mass stopping away from its equilibrium location due to friction is somewhat subtle. Once the mass reaches zero velocity, it will stop if the available static friction is greater in magnitude than the restoring force at that location. Within a numerical integration scheme, the velocity is never strictly zero, so a careful implementation of the correct switchover between kinetic and static friction would be somewhat tedious to simulate. In Fig. 1, it is observed that the mass appears visually to stop when the speed approaches an average value . This small, nonzero average results in motion which is not visible in the figure and represents a slight error in the numerical integration. As the step size , the motion ceases. A more detailed investigation of this numerical artifact can be found in the supplementary material.21
Since and both have the same dimensions, b and differ by a factor of velocity. To compare the quantities, one must introduce an appropriate factor of velocity as .