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.

A standard rite of passage for an undergraduate mechanics student is exploring the solution to a harmonic oscillator in the presence of linear damping,

mx¨+kxfα(ẋ)=0,
(1)

for which fα(ẋ)bẋ|ẋ|α1f1(ẋ)=bẋ. In the case α = 1, given the initial conditions x(0)=x0 and ẋ(0)=v0, the underdamped (b/2m<ω0) solution takes the form1 

x(t)=eγt[x0cos[ωt]+(v0+γx0ω)sin[ωt]],
(2)

where γb/2m,ωω02γ2, and ω02k/m. In the weak-damping limit γω0 when v0=0, the solution simplifies substantially to

x(t)x0eγtcos[ω0t](γω0).
(3)

Here, Eq. (3) has the form x(t)A(t)cos(ω0t) with A(t)=x0eγt, corresponding to the undamped (γ = 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α(ẋ)=bẋ|ẋ|α1, and the condition which generalizes the linear, “weak-damping” limit γω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 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 α>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 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.

For an oscillator obeying Eq. (1) with arbitrary α, one may write the exact solution as x(t)=A(t)cos[ω0t+φ(t)]. We take φ(0)=0, corresponding to ẋ(0)=0 without loss of generality. For weak damping, it is physically plausible to suppose φ(t)=φ(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

E(t)=12mẋ2(t)+12kx2(t),
(4)
=12kA2(t)[1+(Ȧ(t)ω0A(t))2cos2(ω0t)(Ȧ(t)ω0A(t))sin(2ω0t)],
(5)

where the second line follows from evaluating ẋ=Ȧ(t)cos(ω0t)A(t)ω0sin(ω0t). We consider the slow dissipation limit for which Ȧω0A, corresponding to a regime in which the timescale over which the amplitude changes appreciably is much longer than the period of oscillation. Neglecting contributions involving Ȧω0A gives

E12kA2(t).
(6)

Differentiating both sides of Eq. (4) with respect to time gives

Ė=ẋ[mx¨+kx]=fα(ẋ)ẋ,
(7)

where the last line follows from using Eq. (1). Taking the time derivative of Eq. (6) gives Ė=kAȦ and equating the two expressions for Ė yields the relationship

kAȦ=fα(ẋ)ẋ.
(8)

Equation (8) becomes quite useful in the limit of slow dissipation (Ȧω0A), where

ẋ(t)Aω0sin[ω0t].
(9)

Using Eq. (8) and taking fα(ẋ)=bẋ|ẋ|α1, one may separate variables and integrate. The basic approach described here was employed previously in Ref. 16 where the exponents α=0,1,2 were treated as separate cases. The aim of the present note is to demonstrate that these cases can be obtained as limits of a more general solution for arbitrary α 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ω0α/kx01α. 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)=x0, the result is

A(t)1αx01α1α=bω01+αk0t|sin1+α(ω0t)|dt(α1).
(10)

The case α = 1, corresponding to linear damping, follows from writing A1α=e(1α)lnA and expanding the exponential for small ϵ=1α so that the limit ϵ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

0t|sinα+1(ω0t)|dt|sin1+α(ωt)|t,
(11)

where f(t) denotes the average value of a function f. For a periodic function f(t), the average only needs to be computed over a single period. In the case of |sinα+1(ω0t)|, the period is π/ω0, and this average can be evaluated in terms of the Euler beta function18 

0t|sinα+1(ωt)|dt2α+1B(α2+1,α2+1)t/π,
(12)

where B(x,y)Γ(x)Γ(y)/Γ(x+y). For α1, the result is

A(t)=x0[1+(α1)b(2ω0)α+1tπkx01αB[α2+1,α2+1]]1/α1(α1).
(13)

The particular expressions which follow for α=0,1,32,2,3 are collected in Table I. For α>1, the dependence of the amplitude on the initial displacement x0 disappears in the long-time limit, as is known to be the case for quadratic damping.19 As noted previously, the cases for α=0,1,2 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 α=32,3 depicted as specific examples. It has been claimed13 that the cases α>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 Att(1/α1) as t. The validity of this claim will be revisited below when multiple damping terms are included.

Table I.

Approximate time-dependence of oscillation amplitude in the low-dissipation limit computed form Eq. (13) for α=0,1,32,2,3.

αA(t)
x02btπmω0 
x0e(bt/2m) 
32 x01+3bt2x0ω0Γ2345π3/2m2 
x01+4bω0x0t3πm 
x01+3bkx02t4m2 
αA(t)
x02btπmω0 
x0e(bt/2m) 
32 x01+3bt2x0ω0Γ2345π3/2m2 
x01+4bω0x0t3πm 
x01+3bkx02t4m2 

Using Eq. (13), we can explore the domain of validity for the slow-dissipation limit with arbitrary damping. The (approximate) ansatz x(t)A(t)cos(ω0t) was inspired by the underdamped solution corresponding to α = 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 α1, the approximation will work best at either long or short times. To satisfy Ȧω0A for general α1, one must have

|Ȧ|A=κω0|1+(α1)κω0t|1ω0,
(14)

where κ2b(2ω0)αB[α/2+1,α/2+1]/(πkx01α) is a dimensionless parameter. The case α = 1 reduces directly to b/2mω0, indicating the equivalence of the slow-dissipation limit to the weak-damping limit. For α>1, Eq. (14) is satisfied in the long-time limit where (α1)κω0t1 so that the condition reduces simply to

|Ȧ|Aω0t11(α1)tω0.
(15)

That is, Eq. (13) becomes valid at times t1/(α1)ω0 regardless of the value of b. An interesting corollary to this statement is that critically damped or overdamped behavior does not occur for any damping strength when α>1. The absence of critically damped or overdamped dynamics has been demonstrated previously10 for the particular case of α = 2. The conjecture that only underdamped behavior occurs for α>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 κ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.

Interestingly, the analytically solvable case α = 1 represents a sort of crossover case between where slow-dissipation holds at long-times and where it holds at short times. For α<1, the constant κ scales with the physical constants as κbω0α/kx01α, and the slow-dissipation condition yields

|Ȧ|Abω0α+1kx01α|1(1α)bω0α+1kx01αt|1ω0.
(16)

Equation (16) is less transparent than the case for α>1, but at short times, it is clear that slow dissipation only occurs when certain conditions are met on damping strength b, initial amplitude x0, angular frequency ω0, and spring stiffness k. Weak damping corresponds to bkx01α/ω0α, or κ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,

bkx0=FfrictionFspring(max)1.
(17)

However, even if this weak-damping condition is satisfied, Eq. (16) cannot be satisfied for long times ω0t1 since the denominator shrinks with increasing time. Thus for α<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 μs>μk 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 μs=μk=b/mg, 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 α>1 and modest values of κ, the agreement between numerical and approximate solutions at long times ω0t1 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 

Fig. 1.

Comparison of the numerical solution of Eq. (1) with α = 0 to the analytic approximation xa(t)=A(t)cos(ω0t) with A(t) given in Table I. Here κ0.042, so that according to Eq. (16), the approximation should work well at short times. Sliding friction allows for the mass to stop away from equilibrium, which is not captured by the approximate solution.

Fig. 1.

Comparison of the numerical solution of Eq. (1) with α = 0 to the analytic approximation xa(t)=A(t)cos(ω0t) with A(t) given in Table I. Here κ0.042, so that according to Eq. (16), the approximation should work well at short times. Sliding friction allows for the mass to stop away from equilibrium, which is not captured by the approximate solution.

Close modal
Fig. 2.

Comparison of the numerical solution of Eq. (1) with α=32 to the analytic approximation xa(t)=A(t)cos(ω0t) with A(t) given in Table I. Here κ0.19. The case shown is fairly generic for α>1 in which the approximation is quite robust, even for a modest κ.

Fig. 2.

Comparison of the numerical solution of Eq. (1) with α=32 to the analytic approximation xa(t)=A(t)cos(ω0t) with A(t) given in Table I. Here κ0.19. The case shown is fairly generic for α>1 in which the approximation is quite robust, even for a modest κ.

Close modal

Finally, we consider briefly the possibility of two distinct types of damping. That is, Eq. (8) is modified as

kAȦ=bẋ2|ẋ|α1bẋ2|ẋ|α1,
(18)
bAα+1ω0α+1|sinα+1(ω0t)|bAα+1ω0α+1|sinα+1(ω0t)|,
(19)

for two power-law functions with exponents α, α and respective amplitudes b, b. Without loss of generality, we take α>α. In the second line, the approximation ẋ(t)ω0Asin(ω0t) has been applied. For αα, 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 |sinα+1(ω0t)| 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),

x0A(t)dAcAα+cAαtπk,
(20)

where c=b(2ω0)α+1B(α/2+1,α/2+1) and c=b(2ω0)α+1B(α/2+1,α/2+1). The general case does not lend itself to transparent results,22 so let us take the case of α = 1 and α=2 corresponding to a combination of viscous drag and air resistance. In this case the integral in Eq. (20) is straightforwardly evaluated, giving

A(t)cλe(ct/πk)1λce(ct/πk),
(21)

where λx0/c+cx0. For α = 1 and α=2, one has c=πbk/2m and c=4bω0k/3m. 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(t)ω0t1x01+8bω03bπx0e(bt/2m).
(22)

A few general features can be gleaned from Eq. (22). First, the long-time decay envelope controlled by the exponential factor e(bt/2m) 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 ω0x0b.23 In the limit of negligible quadratic damping, ω0x0bb, the second term in the denominator may be neglected so one recovers the standard result for linear damping, A(t)x0e(bt/2m). More interesting is the case that ω0x0bb so that the quadratic drag is “stronger.” In this limit A(t)3bπ/8bω0e(bt/2m). 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 α>2 is nonphysical13 due to oscillations which decay too slowly only holds when a single damping term with α>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.

Fig. 3.

Comparison of numerical solution of Eq. (1) with the combined total damping force bẋ|ẋ|α1bẋ|ẋ|α1 for α = 1, α′ = 2, and b=ω0x0b. The analytic expression A(t)cos[ω0t] with A(t) given by Eq. (21) provides an excellent approximation in the slow-dissipation regime. At long times, the simpler long-time limit (LTL) in Eq. (22) becomes an accurate representation.

Fig. 3.

Comparison of numerical solution of Eq. (1) with the combined total damping force bẋ|ẋ|α1bẋ|ẋ|α1 for α = 1, α′ = 2, and b=ω0x0b. The analytic expression A(t)cos[ω0t] with A(t) given by Eq. (21) provides an excellent approximation in the slow-dissipation regime. At long times, the simpler long-time limit (LTL) in Eq. (22) becomes an accurate representation.

Close modal

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.

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.

1.
F. J.
McCormack
, “
Smooth transitions between the three damping cases for the harmonic oscillator
,”
Am. J. Phys.
63
(
12
),
1151
(
1995
).
2.
J. R.
Taylor
,
Classical Mechanics
(
University Science Books
,
Saulsalito, CA
,
2005
).
3.
S. T.
Thornton
and
J. B.
Marion
,
Classical Dynamics of Particles and Systems
, 5th ed. (
Thomson-Brooks/Cole
,
Belmont, CA
,
2004
).
4.
G. R.
Fowles
and
G. L.
Cassiday
,
Analytical Mechanics
, 7th ed. (
Thomson-Brooks/Cole
,
Belmont, CA
,
2005
).
5.
C.
Barratt
and
G. L.
Strobel
, “
Sliding friction and the harmonic oscillator
,”
Am. J. Phys.
49
(
5
),
500
501
(
1981
).
6.
A.
Marchewka
,
D. S.
Abbott
, and
R. J.
Beichner
, “
Oscillator damped by constant-magnitude friction force
,”
Am. J. Phys.
72
(
4
),
477
483
(
2004
).
7.
P.
Onorato
,
D.
Mascoli
, and
A.
De Ambrosis
, “
Damped oscillations and equilibrium in a mass-spring system subject to sliding friction forces: Integrating experimental and theoretical analyses
,”
Am. J. Phys.
78
(
11
),
1120
1127
(
2010
).
8.
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
).
9.
F. M.
White
,
Fluid Mechanics
, 7th ed. (
McGraw-Hill
,
New York
,
2009
).
10.
B. R.
Smith
, Jr.
, “
The quadratically damped oscillator: A case study of a nonlinear equation of motion
,”
Am. J. Phys.
80
(
9
),
816
824
(
2012
).
11.
P. B.
Kahn
and
Y.
Zarmi
, “
Weakly nonlinear oscillations: A perturbative approach
,”
Am. J. Phys.
72
(
4
),
538
552
(
2004
).
12.
S. J.
Elliott
,
M.
Chandchi Tehrani
, and
R. S.
Langley
, “
Nonlinear damping and quasi-linear modelling
,”
Philos. Trans. R. Soc. A
373
,
20140402
(
2015
).
13.
G. D.
Quiroga
and
P. A.
Ospina-Henao
, “
Dynamics of damped oscillations: Physical pendulum
,”
Eur. J. Phys.
38
(
6
),
065005
(
2017
).
14.
G.
Terenzi
, “
Dynamics of SDOF systems with nonlinear viscous damping
,”
J. Eng. Mech.
125
(
8
),
956
963
(
1999
).
15.
R. H.
Pritchard
and
E. M.
Terentjev
, “
Oscillations and damping in the fractional Maxwell materials
,”
J. Rheol.
61
(
2
),
187
203
(
2017
).
16.
X. J.
Wang
,
C.
Schmitt
, and
M.
Payne
, “
Oscillations with three damping effects
,”
Eur. J. Phys.
23
(
2
),
155
164
(
2002
).
17.
C. M.
Bender
and
S. A.
Orszag
,
Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory
(
Springer-Verlag
,
New York
,
1999
).
18.
Handbook of Mathematical Functions
, edited by
M.
Abramowitz
and
I.
Stegun
(
Dover
,
New York
,
1965
).
19.
Y. B.
Zeldovich
and
I. M.
Yaglom
,
Higher Math for Beginners
, translated by Eugene Yankovsky (
Mir Publishers
,
Moscow
,
1987
).
20.

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 vω0x0105. 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 δt0, the motion ceases. A more detailed investigation of this numerical artifact can be found in the supplementary material.21 

21.
See supplementary material at http://dx.doi.org/10.1119/10.0001794 for the Python-based Jupyter notebooks with the basic computational scheme.
22.
The general case of this integral can be evaluated formally in terms of the hypergeometric function 2F1(a,b;c;x) as
dxcxα+cxα=x1α2F1(1,α1αα;1+α2ααα;ccxαα)c(1α),
but one can argue that this expression is not much in the way of progress toward a physically transparent expression.
23.

Since bẋ and bẋ2 both have the same dimensions, b and b differ by a factor of velocity. To compare the quantities, one must introduce an appropriate factor of velocity as bω0x0b.

24.
R.
Hauko
and
R.
Repnik
, “
Damped harmonic oscillation: Linear or quadratic drag force?
,”
Am. J. Phys.
87
(
11
),
910
914
(
2019
).
25.
Further information about the torsional oscillator is available at the TeachSpin website
, <http://www.teachspin.com/torsional-oscillator.html>.
26.
P. F.
Hinrichsen
and
C. I.
Lardner
, “
Combined viscous and dry friction damping of oscillatory motion
,”
Am. J. Phys.
86
(
8
),
577
584
(
2018
).

Supplementary Material