Linear Fractionally Damped Oscillator

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 is 0≤𝑣≤1. At the lower end (𝑣=0) the equation represents an undamped oscillator and at the upper end (𝑣=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).

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,

𝐷²_𝑡𝑥+𝜆₀𝐷^𝜈_𝑡𝑥+𝜔²𝑥=0. (1.1) 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 for 0<𝜈<1

ℒ(₀^𝑅𝐿𝐷^𝜈_𝑡𝑓(𝑡))=𝑠^𝜈𝐹(𝑠)−₀^𝑅𝐿𝐷_𝑡^𝜈−1𝑓(𝑡)|₀,ℒ(^𝐶₀𝐷^𝜈_𝑡𝑓(𝑡))=𝑠^𝜈𝐹(𝑠)−𝑓(𝑡)|₀. (1.2) 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 become 3/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

𝐷²_𝑡𝑥+𝜆𝐷_𝑡𝑥+𝜔²𝑥=0. (2.1) The constants 𝜆 and 𝜔 are taken to be real and positive. 𝜆 is the damping force per unit mass. 𝜔² is the restoring force per unit mass. In both cases (fractional and Nonfractional) the following initial conditions will be used:

𝑥(0)=𝑥₀,𝐷_𝑡𝑥(0)=𝑥₁. (2.2) Transforming (2.1) together with the initial conditions (2.2) gives the following:

𝑠²𝑋(𝑠)−𝑠𝑥₀−𝑥₁+𝜆(𝑠𝑋(𝑠)−𝑥₀)+𝜔²𝑋(𝑠)=0, (2.3) or

𝑋(𝑠)=𝑠𝑥₀𝑠²+𝜆𝑠+𝜔²+𝑥₁+𝜆𝑥₀𝑠²+𝜆𝑠+𝜔². (2.4) 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 𝑠 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

1𝑥(𝑡)=Residue−𝑒2𝜋𝑖^𝑠𝑡𝑋(𝑠)𝑑𝑠. (2.5) Recall that the Bromwich contour begins at 𝛾−𝑖∞ goes vertically up to 𝛾+𝑖∞ (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 𝛾−𝑖∞. 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

𝑠²+𝜆𝑠+𝜔²=0. (2.6) There are three different cases.

(1)𝜆>2𝜔; 2 unequal real roots that are negative

𝑠_1,2=−𝜆±𝜆²−4𝜔²2. (2.7)

(2)𝜆=2𝜔; 2 repeated real roots that are negative

𝑠₃𝜆=−2=−𝜔. (2.8)

(3)𝜆<2𝜔; 2 complex roots whose real parts are negative

𝑠_4,5=−𝜆±𝑖4𝜔²−𝜆²2. (2.9) See Figure 1 for a graphical representation of the location of the roots.

Figure 1 

Note 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 ±𝑖𝜔.

Computing the residue for case one gives

Residue=lim_{𝑠→𝑠1,2}(𝑠−𝑠_1,2)𝑒^𝑠𝑡(𝑠𝑥₀𝑠²+𝜆𝑠+𝜔²+𝑥₁+𝜆𝑥₀𝑠²+𝜆𝑠+𝜔²), (2.10) or

𝑒𝑥(𝑡)=^𝑠1𝑡2𝑠₁+𝜆(𝑠₁𝑥₀+𝑥₁+𝜆𝑥₀𝑒)+^𝑠2𝑡2𝑠₂+𝜆(𝑠₂𝑥₀+𝑥₁+𝜆𝑥₀). (2.11) As 𝑠₁ and 𝑠₂ 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

𝑥(𝑡)=𝑒^−𝛼𝑡(𝑥₀cos(𝜌𝑡)+(2𝑥₁+𝜆𝑥₀)2𝜌sin(𝜌𝑡)), (2.12) where 𝜌=𝜔²−𝜆²/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

lim_{𝑠→−𝜔}𝑑𝑑𝑠((𝑠+𝜔)²𝑒^𝑠𝑡(𝑠𝑥₀+𝑥₁+𝜆𝑥₀𝑠²+𝜆𝑠+𝜔²)). (2.13) Recall that 𝜆=2𝜔, this allows the denominator to be factored. The limit then becomes

lim_{𝑠→−𝜔}𝑑𝑑𝑠(𝑒^𝑠𝑡(𝑠𝑥₀+𝑥₁+2𝜔𝑥₀)), (2.14) or

𝑥(𝑡)=𝑒^−𝜔𝑡(𝑡(𝜔𝑥₀+𝑥₁)+𝑥₀). (2.15) 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

𝑠²𝑋−𝑠𝑥₀−𝑥₁+𝜆(𝑠^𝜈𝑋−𝑠^𝜈−1𝑥₀)+𝜔²𝑋=0, (3.1) or

𝑋=𝑠𝑥₀+𝑥₁+𝜆𝑠^𝜈−1𝑥₀𝑠²+𝜆𝑠^𝜈+𝜔². (3.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 𝑠. Hence, a Hankel contour will be used. This contour starts at 𝛾−𝑖∞ and goes vertically up to 𝛾+𝑖∞ (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 𝛾−𝑖∞.

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:

𝑠²+𝜆𝑠^𝜈+𝜔²=0, (3.3) Which, for an arbitrary 𝜈, is not a trivial problem. To determine if there are solutions, and if so how many, let 𝑠=𝑟𝑒^𝑖𝜃 then (3.3) breaks into 2 equations, a real and an imaginary part

𝑟²cos(2𝜃)+𝜆𝑟^𝜈cos(𝜈𝜃)+𝜔²𝑟=0,₂sin(2𝜃)+𝜆𝑟^𝜈sin(𝜈𝜃)=0. (3.4) 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)

(−𝜆sin(𝜈𝜃))sin(2𝜃)^2/(2−𝜈)cos(2𝜃)+𝜆(−𝜆sin(𝜈𝜃))sin(2𝜃)^{𝜈/(2−𝜈)}cos(𝜈𝜃)+𝜔²=0. (3.5) 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

((sin(𝜈𝜃))^𝜈(sin(2𝜃))²)^1/(2−𝜈)𝜔sin((2−𝜈)𝜃)=(𝜆^1/(2−𝜈))². (3.6)

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

lim_{𝜃→𝜋/2+}((sin(𝜈𝜃))^𝜈(sin(2𝜃))²)^1/(2−𝜈)sin((2−𝜈)𝜃)=∞,lim_{𝜃→𝜋/(2−𝜈)−}((sin(𝜈𝜃))^𝜈(sin(2𝜃))²)^1/(2−𝜈)sin((2−𝜈)𝜃)=0. (3.7) 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,

𝜕𝜕𝜃{((sin(𝜈𝜃))^𝜈(sin(2𝜃))²)^1/(2−𝜈)sin((2−𝜈)𝜃)}<0. (3.8) Computing the derivative, doing some algebra, and throwing away over all factors that are always positive we have

𝜈²sin²(2𝜃)−4𝜈sin(2𝜃)sin(𝜈𝜃)cos((2−𝜈)𝜃)+4sin²(𝜈𝜃)>0. (3.9) On the restricted domain sin(𝜈𝜃)>0, sin(2𝜃)<0, and cos((2−𝜈)𝜃)≤1. This reduces (3.9) to

𝜈²sin²(2𝜃)−4𝜈sin(2𝜃)sin(𝜈𝜃)+4sin²(𝜈𝜃)>0. (3.10) Equation (3.10) can now be factored into a perfect square and prove the assertion made in(3.8)

(𝜈sin(2𝜃)−2sin(𝜈𝜃))²>0. (3.11) 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 2 

Now the question of the poles that has been settled the solution to (1.1) can be generated. Denote the two poles as

𝑠_6,7=𝛽±𝑖𝜎=𝑟𝑒^±𝑖𝜃, (3.12) where 𝛽 and 𝜎 are determined from the r and 𝜃 values that satisfy (3.6) in the usual way, 𝑟=𝛽²+𝜎² and tan(𝜃)=𝜎/𝛽. Note that 𝛽 is negative, the two solutions are in the second and third quadrants and 𝑠₇ is the complex conjugate of 𝑠₆. Note also that 𝛽 plays the role of −𝜆/2from the Nonfractional case. The poles are of order one and the residue is given by

Residue=lim_𝑠→𝑠6(𝑠−𝑠₆)𝑒^𝑠𝑡(𝑠𝑥₀+𝑥₁+𝜆𝑠^𝜈−1𝑥₀𝑠²+𝜆𝑠^𝜈+𝜔²)+lim_𝑠→𝑠7(𝑠−𝑠₇)𝑒^𝑠𝑡(𝑠𝑥₀+𝑥₁+𝜆𝑠^𝜈−1𝑥₀𝑠²+𝜆𝑠^𝜈+𝜔²)=𝑒^𝑠6𝑡(𝑠₆𝑥₀+𝑥₁+𝑥₀𝜆𝑠₆^𝜈−12𝑠₆+𝜈𝜆𝑠₆^𝜈−1)+𝑒^𝑠6𝑡(𝑠₆𝑥₀+𝑥₁+𝑥₀𝜆𝑠₆^𝜈−12𝑠₆+𝜈𝜆𝑠₆^𝜈−1), (3.13) where 𝑠₇ has been replaced by 𝑠₆. After some algebra this can be reduced to

2𝑒^𝛽𝑡𝑥cos(𝜎𝑡)⁰(2𝑟²+𝜈𝜆²𝑟^2𝜈−2+𝜆𝑟^𝜈(𝜈+2)cos(𝜃(𝜈−2)))+𝔄4𝑟²+4𝜈𝜆^𝜈𝑟^𝜈cos((2−𝜈)𝜃)+𝜈²𝜆^2𝜈𝑟^2𝜈−2+2𝑒^𝛽𝑡𝑥sin(𝜎𝑡)⁰(𝜆𝑟^𝜈(𝜈−2)sin((𝜈−2)𝜃))+𝑥₁(2𝑟sin(𝜃)+𝜈𝜆𝑟^𝜈−1sin(𝜃(𝜈−1)))4𝑟²+4𝜈𝜆𝑟^𝜈cos((2−𝜈)𝜃)+𝜈²𝜆²𝑟^2𝜈−2, (3.14) where 𝔄 denote to 𝑥₁(2𝑟cos(𝜃)+𝜈𝜆𝑟^𝜈−1cos(𝜃(𝜈−1))).

For the contour integral the only contributions come from the paths along the negative real axis

𝜆𝜋^∞₀(𝑅𝑥₀−𝑥₁)sin(𝜈𝜋)+(𝑥₀/𝑅)(𝑅²+𝜔²)sin(𝜋(𝜈−1))(𝑅²+𝜔²)²+2𝜆𝑅^𝜈(𝑅²+𝜔²)cos(𝜈𝜋)+(𝜆𝑅^𝜈)²𝑒^−𝑅𝑡𝑅^𝜈𝑑𝑅. (3.15) 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

𝑥(𝑡)=𝐴𝑒^𝛽𝑡cos(𝜎𝑡)+𝐵𝑒^𝛽𝑡sin(𝜎𝑡)−Decayfunction. (3.16) 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 𝑒^𝛽𝑡 is similar to the damping factor for the Nonfractional case, 𝑒^{−𝜆𝑡/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(𝑠₆). 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) 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,

𝜎=0𝜆≥2𝜔,𝜎=4𝜔²−𝜆²2𝜆<2𝜔. (4.2) For 0<𝜈<1 there will always be a nonzero frequency. Note that 0≤4𝜔²−𝜆²/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(𝑠₆), 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

𝑠²+𝜆𝑠^𝜈+𝜔²=0, (4.3) 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 𝑑𝑠/𝑑𝜈 (remember, 𝜆 and 𝜔 are being held fixed)

𝑑𝑠𝑑𝜈=−𝜆𝑠^𝜈ln(𝑠)𝑠2𝑠²+𝜆𝑠^𝜈𝜈. (4.4) The imaginary part of this equation is

𝑑𝜎𝑑𝜈=Im(𝑑𝑠𝑑𝜈). (4.5) Specifically, consider this equation at 𝜈=0

𝑑𝜎|𝑑𝜈_𝜈=0=Im((𝑠²+𝜔²)ln(𝑠)𝑠2𝑠²−𝜈(𝑠²+𝜔²))|_𝜈=0=𝜆(ln(𝜆+𝜔²))4√𝜆+𝜔². (4.6) This gives three initial slopes for the rate of change of 𝜎 with respect to 𝜈.

𝜆+𝜔²>1⇒ The frequency initially increases with increasing damping order.𝜆+𝜔²=1⇒The frequency initially is not changing with increasing damping order.𝜆+𝜔²<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.

Figure 3 is a representative graph of case one, 𝜆+𝜔²>1. The green graph is for 𝜆=𝜔=1, the black graph is for 𝜆=2and𝜔=1, and the red graph is for 𝜆=3and𝜔=1.

Figure 3 

Figure 4 is a representative graph from case two, 𝜆+𝜔²=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 4 

Figure 5 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.