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.

*α*= 1, given the initial conditions $ x ( 0 ) = x 0$ and $ x \u0307 ( 0 ) = v 0$, the underdamped $ ( b / 2 m < \omega 0 )$ solution takes the form

^{1}

*γ*= 0) solution with a time-dependent amplitude. In the remainder of this note, this approximate form is generalized for the case of arbitrary power-law damping functions, $ f \alpha ( x \u0307 ) = \u2212 b x \u0307 | x \u0307 | \alpha \u2212 1$, and the condition which generalizes the linear, “weak-damping” limit $ \gamma \u226a \omega 0$ is discussed.

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 friction^{5–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 $ \alpha > 2$ 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 fluids^{9} sometimes obeys such a nonintegral power law experimentally.^{14} Fractional power-law functions also appear as relaxation functions in certain models of viscoelasticity^{15} 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.

*α*, one may write the exact solution as $ x ( t ) = A ( t ) \u2009 cos \u2009 [ \omega 0 t + \phi ( t ) ]$. We take $ \phi ( 0 ) = 0$, corresponding to $ x \u0307 ( 0 ) = 0$ without loss of generality. For weak damping, it is physically plausible to suppose $ \phi ( t ) = \phi ( 0 )$, but this point will be revisited below. Suppose the phase remains equal to its initial value, the energy contained in the oscillator at some time

*t*is

*α*while also exploring the domain of validity for the approximations used. Indeed, this general method of solution is a particularly intuitive application of multiple scale analysis,

^{11,17}which is essentially an expansion of the solution at long times in powers of the (dimensionless) damping coefficient $ b \omega 0 \alpha / k x 0 1 \u2212 \alpha $. It can be shown that the phase remains equal to its initial value (here zero) to the lowest order in

*b*. Indeed, the case

*α*= 3 was considered as an illustrative example of this method in Ref. 17. Taking $ A ( 0 ) = x 0$, the result is

*α*= 1, corresponding to linear damping, follows from writing $ A 1 \u2212 \alpha = e ( 1 \u2212 \alpha ) \u2009 ln \u2009 A$ and expanding the exponential for small $ \u03f5 = 1 \u2212 \alpha $ so that the limit $ \u03f5 \u2192 0$ may be taken. The remaining work concerns the integral appearing on the right-hand sides of Eq. (10). Due to the periodic, non-negative nature of the integrand, the result should be approximately proportional to time

*t*. Therefore, we may write

*f*. For a periodic function

*f*(

*t*), the average only needs to be computed over a single period. In the case of $ | \u2009 sin \alpha + 1 ( \omega 0 t ) |$, the period is $ \pi / \omega 0$, and this average can be evaluated in terms of the Euler beta function

^{18}

^{,}

*x*

_{0}disappears in the long-time limit, as is known to be the case for quadratic damping.

^{19}As noted previously, the cases for $\alpha =0,\u20091,\u20092$ have been obtained in Ref. 16 by working out the integral in Eq. (12) separately for each case. Note that Eq. (13) contains all of this information as specific cases, and the mathematically brave student can capture all three from a single integration of Eq. (10). With minimal additional effort, one may also obtain

*A*(

*t*) for any real value of damping exponent with $\alpha = 3 2,\u2009\u20093$ depicted as specific examples. It has been claimed

^{13}that the cases $ \alpha > 2$ are not physical, as the amplitude decays too slowly to agree with experimental systems. Indeed, the predicted decay becomes slower with increasing

*α*, since Eq. (13) shows $A t\u2009\u223c t ( \u2212 1 / \alpha \u2212 1 )$ as $ t \u2192 \u221e$. The validity of this claim will be revisited below when multiple damping terms are included.

α . | A(t) . |
---|---|

0 | $ x 0 \u2212 2 b t \pi m \omega 0$ |

1 | $ x 0 e ( \u2212 b t / 2 m )$ |

$ 3 2$ | $ x 0 1 + 3 b t 2 x 0 \omega 0 \Gamma 2 3 4 5 \pi 3 / 2 m \u2212 2$ |

2 | $ x 0 1 + 4 b \omega 0 x 0 t 3 \pi m$ |

3 | $ x 0 1 + 3 b k x 0 2 t 4 m 2$ |

α . | A(t) . |
---|---|

0 | $ x 0 \u2212 2 b t \pi m \omega 0$ |

1 | $ x 0 e ( \u2212 b t / 2 m )$ |

$ 3 2$ | $ x 0 1 + 3 b t 2 x 0 \omega 0 \Gamma 2 3 4 5 \pi 3 / 2 m \u2212 2$ |

2 | $ x 0 1 + 4 b \omega 0 x 0 t 3 \pi m$ |

3 | $ x 0 1 + 3 b k x 0 2 t 4 m 2$ |

*α*= 1 and could be viewed as the weak-damping case for the linear oscillator. Here, we shall see that the restriction of weak damping is generally not equivalent to that of slow dissipation. For $ \alpha \u2260 1$, the approximation will work best at either long or short times. To satisfy $ A \u0307 \u226a \omega 0 A$ for general $ \alpha \u2260 1$, one must have

*α*= 1 reduces directly to $ b / 2 m \u226a \omega 0$, indicating the equivalence of the slow-dissipation limit to the weak-damping limit. For $ \alpha > 1$, Eq. (14) is satisfied in the long-time limit where $ ( \alpha \u2212 1 ) \kappa \omega 0 t \u226b 1$ so that the condition reduces simply to

*b*. An interesting corollary to this statement is that critically damped or overdamped behavior does not occur for any damping strength when $ \alpha > 1$. The absence of critically damped or overdamped dynamics has been demonstrated previously

^{10}for the particular case of

*α*= 2. The conjecture that only underdamped behavior occurs for $ \alpha > 1$ is somewhat speculative, as the asymptotic results are formally correct only to leading order in a hypothetical expansion in powers of the dimensionless parameter $ \kappa \u221d b$. Agreement between the asymptotic results and a numerical integration (see below) turns out to be quite good, even for large values of

*b*with only the phase of the asymptotic solution disagreeing substantially from the numerical solution. It is interesting to note that the fractional Maxwell model for viscoelastic materials also shows an absence of overdamped and critically damped behavior away from the linear regime.

^{15}However, models of viscoelasticity are considerably more complex than those of damped harmonic oscillators, so there is little room for a true analogy.

*α*= 1 represents a sort of crossover case between where slow-dissipation holds at long-times and where it holds at short times. For $ \alpha < 1$, the constant

*κ*scales with the physical constants as $ \kappa \u223c b \omega 0 \alpha / k x 0 1 \u2212 \alpha $, and the slow-dissipation condition yields

*b*, initial amplitude

*x*

_{0}, angular frequency

*ω*

_{0}, and spring stiffness

*k*. Weak damping corresponds to $ b \u226a k x 0 1 \u2212 \alpha / \omega 0 \alpha $, or $ \kappa \u226a 1$. Indeed, for the case of sliding friction (

*α*= 0),

*b*represents the constant sliding friction force, and the weak damping simply means the frictional force is small compared to the scale of the spring's restoring force,

However, even if this weak-damping condition is satisfied, Eq. (16) cannot be satisfied for long times $ \omega 0 t \u226b 1$ since the denominator shrinks with increasing time. Thus for $ \alpha < 1$, 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 $ \mu s > \mu k$ for static (*μ _{s}*) and kinetic (

*μ*) coefficients of friction, so this behavior should be treated more carefully for more physically correct results at long times.

_{k}^{20}Effectively, the model considered treats $ \mu s = \mu k = b / m g$, 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 solutions

^{6}as well as experimental realizations.

^{7,8}For $ \alpha > 1$ and modest values of

*κ*, the agreement between numerical and approximate solutions at long times $ \omega 0 t \u226b 1$ 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}

*α*, $ \alpha \u2032$ and respective amplitudes

*b*, $ b \u2032$. Without loss of generality, we take $ \alpha \u2032 > \alpha $. In the second line, the approximation $ x \u0307 ( t ) \u2243 \omega 0 A \u2009 sin \u2009 ( \omega 0 t )$ has been applied. For $ \alpha \u2260 \alpha \u2032$, the approximate differential equation governing

*A*(

*t*) is no longer separable. However, the heart of the approximation used in the simpler case was to replace oscillatory terms of the form $ | \u2009 sin \alpha + 1 ( \omega 0 t ) |$ by their average values over a single oscillation cycle. Making this approximation and using Eqs. (11) and (12) yields a separable, approximate equation for

*A*(

*t*),

^{22}so let us take the case of

*α*= 1 and $ \alpha \u2032 = 2$ corresponding to a combination of viscous drag and air resistance. In this case the integral in Eq. (20) is straightforwardly evaluated, giving

*α*= 1 and $ \alpha \u2032 = 2$, one has $ c = \pi b k / 2 m$ and $ c \u2032 = 4 b \u2032 \omega 0 k / 3 m$. In the long time limit the exponential in the denominator of Eq. (21) may be neglected, so that the asymptotic form of the amplitude may be written

A few general features can be gleaned from Eq. (22). First, the long-time decay envelope controlled by the exponential factor $ e ( \u2212 b t / 2 m )$ 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 $ \omega 0 x 0 b \u2032$.^{23} In the limit of negligible quadratic damping, $ \omega 0 x 0 b \u2032 \u226a b$, the second term in the denominator may be neglected so one recovers the standard result for linear damping, $A(t)\u2192 x 0 e ( \u2212 b t / 2 m )$. More interesting is the case that $ \omega 0 x 0 b \u2032 \u226b b$ so that the quadratic drag is “stronger.” In this limit $A(t)\u2192 3 b \pi / 8 b \u2032 \omega 0 e ( \u2212 b t / 2 m )$. 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 *x*_{0} 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 $ \alpha > 2$ is nonphysical^{13} due to oscillations which decay too slowly only holds when a single damping term with $ \alpha > 2$ 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 Oscillator^{25} 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

*Classical Dynamics of Particles and Systems*

*Analytical Mechanics*

*Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory*

*Handbook of Mathematical Functions*

*Higher Math for Beginners*

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 $\u27e8v\omega 0x0\u27e9\u223c10\u22125$. 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 $\delta t\u21920$, the motion ceases. A more detailed investigation of this numerical artifact can be found in the supplementary material.^{21}

Since $bx\u0307$ and $b\u2032x\u03072$ both have the same dimensions, *b* and $b\u2032$ differ by a factor of velocity. To compare the quantities, one must introduce an appropriate factor of velocity as $b\u2032\u2192\omega 0x0b\u2032$.