# Linear Fractionally Damped Oscillator

Authors: Mark Naber
Journal: International Journal of Differential Equations
 (2010)

---

## Abstract

The linearly damped oscillator equation is considered with the damping term generalized to a Caputo fractional derivative. The order of the derivative being considered is0≤v≤1. At the lower end (v=0)
the equation represents an undamped oscillator and at the upper end (v=1)
the ordinary linearly damped oscillator equation is recovered. A solution is found analytically, and a comparison with the ordinary linearly damped oscillator is made. It is found that there are nine distinct cases as opposed to the usual three for the ordinary equation (damped, over-damped, and critically damped). For three of these cases it is shown that the frequency of oscillation actually increases with increasing damping order before eventually falling to the limiting value given by the ordinary damped oscillator equation. For the other six cases the behavior is as expected, the frequency of oscillation decreases with increasing order of the derivative (damping term).


---

## Body

## 1. Introduction

In this paper the linearly damped oscillator equation is considered with the damping term replaced by a fractional derivative [1] whose order, ν, will be restricted to, 0≤ν≤1,(1.1)Dt2x+λ0Dtνx+ω2x=0.
Burov and Barkai [2, 3] examined such an equation in connection with critical behavior.  They were able to determine a solution in terms of generalized Mittag-Leffler functions.  Nonlinear fractional oscillators have been studied numerically by Zaslavsky et al. [4].  He was primarily interested in chaotic behavior.  It is hoped that a careful study of the analytic solution to the linear fractionally damped equation will help shed light on properties of the nonlinear equation and be of use for direct applications of fractionally damped oscillations (see, e.g., [5, 6]).In this paper the Caputo formulation of the fractional derivative will be used.  The Caputo derivative is preferred over the Riemann-Liouville derivative for physical reasons.  Consider the Laplace transform of the two formulations of the fractional derivative for0<ν<1(1.2)ℒ(Dtν0RLf(t))=sνF(s)-Dtν-10RLf(t)|0,ℒ(Dtν0Cf(t))=sνF(s)-f(t)|0.
The constant term arising from the Laplace transform of the Caputo derivative is merely the initial value of the function.  For the Riemann-Liouville derivative this is not the case.  The constant term arising from the Laplace transform currently has no simple physical interpretation.  Hence the Caputo fractional derivative seems to be more useful for modeling physical systems.If the order of the fractional damping term is allowed to become3/2 (outside the range of values considered in this paper) the equation is usually referred to as the Bagley-Torvik equation (see, e.g., [1, 7, 8]).  The solution of this equation exhibits damped oscillatory behavior similar to what we expect to find for the equation studied in this paper.  The Bagley-Torvik equation was originally derived to study the motion of a rigid plate in a Newtonian fluid [7].The analytic solution to the fractionally damped equation is found by means of Laplace transform.  For the sake of clarity, and for pointing out some unique difficulties with the factional equation, a comparison with the Laplace transform method as applied to the nonfractional case is made.  It is found that there are nine distinct cases for the fractionally damped equation as opposed to the usual three cases for the Nonfractional equation.  In six of the nine cases the results are as expected; increasing the order of the factional derivative increases the effects of the damping (i.e., the frequency of the damping slows as the order of the derivative increases).  However, in three cases, the frequency of the damping actually increases as the order of the fractional derivative increases until a peak value is reached after which the frequency falls to its Nonfractional limit.  The physical reason for this increase in the oscillation frequency is not yet clear.


## 2. The Nonfractional Case

Before a solution to the linear fractionally damped oscillator equation is constructed it will be useful to review the Laplace transform method of solution for the linearly damped oscillator equation(2.1)Dt2x+λDtx+ω2x=0.
The constants λ and ω are taken to be real and positive. λ is the damping force per unit mass. ω2 is the restoring force per unit mass.  In both cases (fractional and Nonfractional) the following initial conditions will be used:(2.2)x(0)=x0,Dtx(0)=x1.
Transforming (2.1) together with the initial conditions (2.2) gives the following:(2.3)s2X(s)-sx0-x1+λ(sX(s)-x0)+ω2X(s)=0,
or(2.4)X(s)=sx0s2+λs+ω2+x1+λx0s2+λs+ω2.
Equation (2.4) can be inverted using tables, however, to shed light on a problem that will happen later, (2.4) will be inverted via the complex inversion integral.  The exponents on the s variable in both terms are whole numbers, hence there will not be a branch cut in the contour integral and the Bromwich contour can be used(2.5)x(t)=Residue-12πi∫estX(s)ds.
Recall that the Bromwich contour begins at γ-i∞ goes vertically up to γ+i∞  (where γ is chosen so that all poles will lie to the left of the vertical contour line and thus all poles will be captured within the contour) and then travels in a half circle (to the left, counter clockwise) back to γ-i∞.  For this problem there is no contribution from the contour integral.  The only contribution comes from the residue.  The residue is generated from the roots of the following quadratic equation(2.6)s2+λs+ω2=0.
There are three different cases.(1)λ>2ω; 2 unequal real roots that are negative(2.7)s1,2=-λ±λ2-4ω22.(2)λ=2ω; 2 repeated real roots that are negative(2.8)s3=-λ2=-ω.(3)λ<2ω; 2 complex roots whose real parts are negative(2.9)s4,5=-λ±i4ω2-λ22.
See Figure 1 for a graphical representation of the location of the roots.Figure 1Note that in cases one and three the poles will be of order one and in case two the pole will be of order two.  Note also that if there were no damping the poles would be on the imaginary axis at±iω.Computing the residue for case one gives(2.10)Residue=lims→s1,2(s-s1,2)est(sx0s2+λs+ω2+x1+λx0s2+λs+ω2),
or(2.11)x(t)=es1t2s1+λ(s1x0+x1+λx0)+es2t2s2+λ(s2x0+x1+λx0).
As s1 and s2 are both negative this solution will decay exponentially.  This is usually referred to as the over-damped case.Case three is computed the same way as case one.  Now the poles are complex so the exponential function can be expressed using sine and cosine with an over-all exponential damping factor(2.12)x(t)=e-αt(x0cos(ρt)+(2x1+λx0)2ρsin(ρt)),
where  ρ=ω2-λ2/4and  α=λ/2Notice that the presence of damping causes the effective angular frequency, ρ, to be smaller than the un-damped angular frequency; that is, the oscillations go slower, as one might expect if there were damping to impede the motion.  This is usually referred to as the under-damped case.  By comparison, cases one and two could be viewed as having a zero frequency or an infinite period.For case two the pole is of order 2, and the residue is given by(2.13)lims→-ωdds((s+ω)2est(sx0+x1+λx0s2+λs+ω2)).
Recall that λ=2ω, this allows the denominator to be factored.  The limit then becomes(2.14)lims→-ωdds(est(sx0+x1+2ωx0)),
or(2.15)x(t)=e-ωt(t(ωx0+x1)+x0).
This is usually called the critically damped case.  Graphs of sample solutions to these three cases can be found in any introductory book on differential equations.


## 3. The Fractional Case

Now consider (1.1).  In this case, λ has units of time raised to the power ν-2.  Hence the over all units of the second term remain the same as in (2.1).  There are two cases to consider, 0<ν<1 and 1<ν<2.  We will consider the former in this paper.  The Laplace transform of (1.1) is(3.1)s2X-sx0-x1+λ(sνX-sν-1x0)+ω2X=0,
or(3.2)X=sx0+x1+λsν-1x0s2+λsν+ω2.
If (3.2) is inverted using a contour integral a branch cut is needed on the negative real axis due to the fractional exponents on the complex variable s.  Hence, a Hankel contour will be used.  This contour starts at γ-i∞ and goes vertically up to γ+i∞ (where γ is again chosen so that all poles will lie to the left of the vertical contour line) then travels in a quarter circle arc (to the left) to just above the negative real axis (i.e., -∞).  The contour then has a cut that goes into the origin (following the negative real axis), around the origin in a clockwise sense (to just below the negative real axis) and then back out to -∞.  The contour then has another quarter circle arc to γ-i∞.Now the question is, where are the poles?  This is a somewhat more involved question than in the standard linearly damped model.  To find the poles the following equation needs to be solved:(3.3)s2+λsν+ω2=0,
Which, for an arbitrary ν, is not a trivial problem.  To determine if there are solutions, and if so how many, let s=reiθ then (3.3) breaks into 2 equations, a real and an imaginary part(3.4)r2cos(2θ)+λrνcos(νθ)+ω2=0,r2sin(2θ)+λrνsin(νθ)=0.
Could there be a solution on the positive real axis?  No, in this case θ=0 and the first equation of (3.4) would be the sum of three positive nonzero terms, which would never be zero.  Could there be a solution on the negative real axis?  No, in this case θ=π and the second term of the second equation of (3.4) would never be zero.  Using similar arguments we can show that there are no solutions on the positive or negative imaginary axes, recall 0<ν<1.  It can also be shown that no solutions are in the right half plane (both terms of the second equation would always be positive).  If there are solutions, they should be in pairs, complex conjugates, with π/2<θ<π and -π/2>θ>-π.  To attempt to find a solution first solve the second equation of (3.4) for r and substitute this into the first equation (only look for positive θ values first)(3.5)(-λsin(νθ)sin(2θ))2/(2-ν)cos(2θ)+λ(-λsin(νθ)sin(2θ))ν/(2-ν)cos(νθ)+ω2=0.
The reader may be worried about the negative sign and the fractional exponent in (3.5), however, for the restricted angular range being considered, π/2<θ<π, sin(2θ) is always negative.  So, the argument of the root will always be positive.Given values forν,λ,andω it would appear to be impossible to solve (3.5) for θ.  Equation (3.5) can be simplified to a more aesthetically pleasing form(3.6)((sin(νθ))ν(sin(2θ))2)1/(2-ν)sin((2-ν)θ)=(ωλ1/(2-ν))2.Now it needs to be seen if there is aθ value that will satisfy (3.6).  For this equation to be true sin((2-ν)θ) needs to be positive.  This will only happen on the restricted domain π/2<θ<π/(2-ν).  Now the question becomes, on this restricted domain can we pick a θ value that will make the left-hand side of (3.6) as large or as small as we wish?  Thus ensuring that no matter what the values of ω, λ, and ν we are given we can always find a θ value that will satisfy (3.6).  Consider the two limits(3.7)limθ→π/2+((sin(νθ))ν(sin(2θ))2)1/(2-ν)sin((2-ν)θ)=∞,limθ→π/(2-ν)-((sin(νθ))ν(sin(2θ))2)1/(2-ν)sin((2-ν)θ)=0.
Since the left-hand side of (3.6) is continuous in θ and we have the two limits above, (3.7), it is guaranteed that there will be at least one solution to (3.6) and hence there will be at least two poles for the residue calculation.  If we can show that the left hand side of (3.6) decreases monotonically in θ over the restricted domain then we know that there will be only one solution to (3.6), and thus only two poles in the residue calculation.  To show that the left hand side of (3.6) decreases monotonically in θ we need to show that the derivative of the left hand side of (3.6) with respect to θ is always negative, that is,(3.8)∂∂θ{((sin(νθ))ν(sin(2θ))2)1/(2-ν)sin((2-ν)θ)}<0.
Computing the derivative, doing some algebra, and throwing away over all factors that are always positive we have(3.9)ν2sin2(2θ)-4νsin(2θ)sin(νθ)cos((2-ν)θ)+4sin2(νθ)>0.
On the restricted domain sin(νθ)>0, sin(2θ)<0, and cos((2-ν)θ)≤1.  This reduces (3.9) to(3.10)ν2sin2(2θ)-4νsin(2θ)sin(νθ)+4sin2(νθ)>0.
Equation (3.10) can now be factored into a perfect square and prove the assertion made in(3.8)(3.11)(νsin(2θ)-2sin(νθ))2>0.
Hence the left hand side of (3.6) will decrease monotonically on the restricted domain with the upper bound being ∞ and the lower bound being 0.  To summarize, it has just been shown that there is always one solution, with a positive angle, to (3.6) and this solution must be such that π/2<θ<π/(2-ν).  Consequently there will be two poles for the residue calculation and they will be complex conjugates of each other.  Notice that for the fractionally damped equation repeated roots are not possible.  Repeated roots can only happen when the order of the derivative becomes one.  See Figure 2 for a graphical representation of the location of the roots.Figure 2Now the question of the poles that has been settled the solution to (1.1) can be generated.  Denote the two poles as(3.12)s6,7=β±iσ=re±iθ,
where β and σ are determined from the r and θ values that satisfy (3.6) in the usual way, r=β2+σ2 and tan(θ)=σ/β.  Note that β is negative, the two solutions are in the second and third quadrants and s7 is the complex conjugate of s6.  Note also that β plays the role of -λ/2from the Nonfractional case.  The poles are of order one and the residue is given by(3.13)Residue=lims→s6(s-s6)est(sx0+x1+λsν-1x0s2+λsν+ω2)+lims→s7(s-s7)est(sx0+x1+λsν-1x0s2+λsν+ω2)=es6t(s6x0+x1+x0λs6ν-12s6+νλs6ν-1)+es̅6t(s̅6x0+x1+x0λs̅6ν-12s̅6+νλs̅6ν-1),
where s7 has been replaced by s̅6.  After some algebra this can be reduced to(3.14)2eβtcos(σt)x0(2r2+νλ2r2ν-2+λrν(ν+2)cos(θ(ν-2)))+𝔄4r2+4νλνrνcos((2-ν)θ)+ν2λ2νr2ν-2+2eβtsin(σt)x0(λrν(ν-2)sin((ν-2)θ))+x1(2rsin(θ)+νλrν-1sin(θ(ν-1)))4r2+4νλrνcos((2-ν)θ)+ν2λ2r2ν-2,
where 𝔄 denote to x1(2rcos(θ)+νλrν-1cos(θ(ν-1))).For the contour integral the only contributions come from the paths along the negative real axis(3.15)λπ∫0∞(Rx0-x1)sin(νπ)+(x0/R)(R2+ω2)sin(π(ν-1))(R2+ω2)2+2λRν(R2+ω2)cos(νπ)+(λRν)2e-RtRνdR.
The solution to (1.1) is then (3.15) subtracted from (3.14).  This may look overly complicated but the solution does have the general form of(3.16)x(t)=Aeβtcos(σt)+Beβtsin(σt)-Decayfunction.
Notice that the decay function, (3.15), goes to zero if ν goes to zero or one; that is, if (1.1) goes to its Nonfractional limits the decay function goes away, as expected.  The damping factor eβt is similar to the damping factor for the Nonfractional case, e-λt/2.  Notice that since the poles, for the residue calculation, have nonzero imaginary and nonzero real parts we will not have the same three distinct cases as we did for the Nonfractional case (critically damped, over-damped, and under-damped).


## 4. The Oscillation Frequency

Consider the frequency of the oscillation component of the solution,σ=Im(s6).  One question we might ask is: how does the frequency change as we change the order of the fractional damping?  When ν is set to zero we have an un-damped oscillator with frequency(4.1)σ=λ+ω2.
When ν is set to one we have the three cases given in Section 2: over-damped, critically-damped, and under-damped. So, the frequency may be zero or nonzero, that is,(4.2)σ=0λ≥2ω,σ=4ω2-λ22λ<2ω.
For 0<ν<1 there will always be a nonzero frequency.  Note that 0≤4ω2-λ2/2<λ+ω2.In the Nonfractional case increasingλ causes the frequency of oscillation to become smaller, monotonically, until the critical cases are reached and the oscillation period becomes infinite (these are the critical and over-damped cases).  In the fractional case the frequency of oscillation, σ=Im(s6), now depends on the order of the derivative, ν, as well as λ and ω.  To try to determine how σ depends on these three parameters consider s to be a function of ν, on 0≤ν≤1, implicitly defined by(4.3)s2+λsν+ω2=0,
for fixed values of λ and ω (both being positive).  Let us restrict our attention to the upper half plane for s.  As such s will be one to one on 0≤ν<1.  Due to the fractional exponent causing a branch cut on the negative real axis s will not be one to one at ν=1.Now the question arises, doesσ fall monotonically with respect to ν?  To get the answer to this consider the derivative of (4.3) with respect to ν and isolate ds/dν (remember, λ and ω are being held fixed)(4.4)dsdν=-λsνln(s)s2s2+λsνν.
The imaginary part of this equation is(4.5)dσdν=Im(dsdν).
Specifically, consider this equation at ν=0(4.6)dσdν|ν=0=Im((s2+ω2)ln(s)s2s2-ν(s2+ω2))|ν=0=λ(ln(λ+ω2))4λ+ω2.
This gives three initial slopes for the rate of change of σ with respect to ν.λ+ω2>1⇒ The frequency initially increases with increasing damping order.λ+ω2=1⇒The frequency initially is not changing with increasing damping order.λ+ω2<1⇒ The frequency initially decreases with increasing damping order.This is not entirely what might have been expected.  In the first case the oscillation frequency actually increases before falling.  Hence there will be some values ofν for which the fractional damping will actually cause the oscillations to go faster than the un-damped oscillator (the damping will still cause the amplitude to decrease).  Each of the above three cases can become any of the three Nonfractionally damped cases by letting ν→1 ((2.7), (2.8), and (2.9)).  Hence, there are nine cases for the linear fractionally damped oscillator.There are some graphs of solutions to the imaginary part of (4.3) (the oscillation frequency) for various values of ν, λ, and ω.  In all three graphs the oscillation frequency is on the vertical axis and the order of the derivative is on the horizontal axis.  The three graphs for each case correspond to what would be under-damped, critically damped, and over-damped for a damped oscillator with whole order derivatives.Figure3 is a representative graph of case one, λ+ω2>1.  The green graph is for λ=ω=1, the black graph is for λ=2andω=1, and the red graph is for λ=3andω=1.Figure 3Figure4 is a representative graph from case two, λ+ω2=1, a flat start. The red graph is for λ=2(2-1)andω=λ/2, the green graph is for λ=1/2andω=1/2, and the black graph is for λ=15/16andω=1/4.Figure 4Figure5 is a representative graph for case three, a decreasing start.  The green graph is for λ=1/2andω=1/8, the black graph is for λ=1/2andω=1/4, and the red graph is for λ=ω=1/2.Figure 5


## 5. Conclusion

In this paper the linear fractionally damped oscillator equation was solved analytically.  It was found that the solution is very similar to the Nonfractional case (decayed oscillations but with the inclusion of an additional decay function).  It was found that there are nine distinct cases, as opposed to the usual three for the ordinary damped oscillator.  An unexpected result was that for three of the cases the oscillation frequency actually increases with increasing order of derivative of the damping term (till a peak value is reached, then the frequency decreases as expected). The physical reason for this increase in oscillation frequency is not yet clear.



---

*Source: 197020-2009-10-22.xml*