Lyapunov Stability of a Fractionally Damped Oscillator with Linear (Anti-)Damping

1 Introduction

This paper is concerned with the development of a Lyapunov stability framework, including Lyapunov’s direct method, for the analysis of stability properties of mechanical systems with fractional damping. The scope of the paper is limited to single degree-of-freedom oscillators with both viscous and fractional damping.

The term fractional refers to fractional calculus, which is a mathematical theory dealing with derivatives and integrals of arbitrary (non-integer) order [1, 2] with a variety of applications in science and engineering [3]. Particularly in mechanics, fractional damping may arise through the modeling of mechanical systems with viscoelastic components. Complex rheological models for viscoelastic materials are often described through an array of classical Kelvin or Maxwell elements, inevitably resulting in a model with a large number of parameters. It has been shown that the viscoelastic behavior of complex materials is in many applications well represented by fractional order elements with only a few parameters [4, 5]. Furthermore, a description using fractional order force-displacements relationships may have much better extrapolation properties on long time-scales. The introduction of fractional calculus in mechanics leads to the concept of a springpot element, being a force law reacting linearly on a fractional derivative of its elongation. In a more general setting, springpot elements may also arise through fractional-order control laws [6].

Many problems in industrial applications originate from (dynamic) instability phenomena, e. g. stick-slip vibrations in oilwell drillstrings, flutter of airfoils, shimmy of vehicles and feedback instabilities in control systems. Methods to rigorously prove stability of linear and nonlinear systems are therefore quintessential. The Lyapunov stability framework, which encompasses the method of Lyapunov functions, forms a central element in the research fields Nonlinear Dynamics and Control Theory [7]. The introduction of springpot elements in (controlled) mechanical systems asks for an extension of the Lyapunov stability framework to non-integer order derivatives. A major complication arises through the non-local character of fractional derivatives, i. e. the force in a springpot element depends on the total history of the elongation. A system with springpot elements has therefore an infinite state which asks for the use of Lyapunov functionals instead of Lyapunov functions in Lyapunov’s direct method, which are introduced in the theory of functional differential equations (FDEs) [89101112]. Special Lyapunov functionals for FDEs with fractional derivatives are introduced in [13, 14], which have been shown to represent the potential energy of an infinite arrangement of springs and dashpots [151617]. The energy expressions for springpots are based on the infinite state or diffusive representation of fractional integrators, which were introduced by Montseny [18], Matignon [19] and have been elaborated by Trigeassou et al. [13, 14, 20, 21]. Beyond that, a lot of work has been done on stability conditions [22, 23] and Lyapunov theory [242526] for fractional differential equations, which cannot directly be used for mechanical systems containing springpots, as the differentiation order is in general irrational for such systems. Furthermore, the Laplace transform method has been used to prove stability of equilibria of mechanical systems containing springpots [27, 28].

The aim of this paper is to give a complete stability analysis of a linear single degree-of-freedom mass-spring-dashpot-springpot system. The analysis encompasses the following results/tasks:

1. The total mechanical energy of the system is derived through the use of the infinite state representation of the springpot element (Section 3).

2. The system is put in the form of an FDE and, based on this, definitions of Lyapunov stability and attractivity of fractionally damped mechanical systems are given. Furthermore, a Lyapunov-Krasovskii theorem is presented for this class of systems (Section 4.2).

3. It is shown that the total mechanical energy is a Lyapunov functional for the system with positive viscous and fractional damping with which stability of the equilibrium can be proven (Section 4.3). Furthermore, a strict Lyapunov functional for positive viscous and fractional damping is derived which rigorously proves asymptotic stability (Section 4.4). However, this Lyapunov functional fails to give a stability result in the case of anti-damping, in which the viscous damping is negative.

4. An extensive eigenvalue analysis is given and, based on that, an expression for the general solution is derived. The eigenvalue analysis reveals that the equilibrium can still be asymptotically stable in the presence of anti-damping (Section 4.5.1).

5. A strict Lyapunov functional for the case of anti-damping is derived in Section 4.5.2. This leads to a Lyapunov stability condition for systems having negative viscous damping but a sufficient amount of positive fractional damping. Moreover, this result opens the way to study global asymptotic stability of nonlinear systems with fractional damping. The theory is demonstrated on a stick-slip oscillator with Stribeck friction law leading to an effective negative viscous damping (Section 4.5.3).

The scope of the present paper is limited to a single degree-of-freedom oscillator, being of course a first step into the direction of multi degree-of-freedom systems. A brief outlook on how these results can be extended to more degrees of freedom will be given in the conclusion section of the paper.

2 Fractional calculus and infinite state representation

We will consider the fractional derivative of Caputo type that is based on the fractional Riemann–Liouville integral, which is defined for an integrable function x = x(t) with t ≥ t₀ and a scalar value α > 0 as

where Γ(α) is the Gamma function. For α = 0 we set It0+0x:=x and it can be seen directly that the choice α = 1 leads to the classical integral. In this paper, we will describe this integral operator by the infinite state representation [21]

The infinite state z(ω,t) fulfills the above differential equation ∀ω≥0 and the fractional integral is obtained by integrating all contributions z(ω,t) weighted by the function

The equivalence of eqs. (1) and (2) is derived in [15]. The representation (2) is very useful to give a mechanical interpretation of a fractional derivative (Section 3) and to formulate Lyapunov functionals for fractionally damped mechanical systems (Section 4). Using the variation of constants formula, we can formulate eq. (2) as

Finally we introduce the fractional Caputo derivative for an absolutely continuous function x = x(t) and 0 < α < 1 as

3 Springpot: Mechanical representation and potential energy

In this section, we will briefly introduce the Caputo springpot (Figure 1) as an abstract mechanical element and discuss its mechanical representation to gain a potential energy expression which may be used for Lyapunov stability considerations of mechanical systems containing springpots. Again, details may be found in [15]. A springpot is defined by its constitutive equation

where f is the force acting on the springpot which results in an elongation q depending on the coefficient c > 0, initialization time t₀ < 0 and differentiation order α∈(0,1). The time interval [t₀,0] represents the entire significant history of the springpot, i. e. for earlier time-instants t ≤ t₀ we assume q(t) = 0 and f(t) = 0.

Figure 1:

Force acting on a springpot.

Together with eqs. (2) and (5) we derive the infinite state representation of a Caputo springpot as

The infinite state y in eq. (7) may, similar as in eq. (4), be expressed by the variation-of-constants formula

which will be useful for stability considerations later on. The above representation leads to a mechanical analogue model of a springpot, which is a parallel arrangement of an infinite number of Maxwell elements (Figure 2) as in [16], where the forces g(ω,t)dω of the Maxwell elements are integrated to the resulting force

Figure 2:

Schematic mechanical representation of a Caputo springpot.

on the system. The springs of the Maxwell elements are characterized by their elongation q_s(ω,t) and spring constant k(ω)dω and the dashpots by the elongation q_d(ω,t) and constant d(ω)dω such that the elongation of the system q(t) appears as

with the incremental internal force

Differentiation of eq. (10) and substitution of eq. (11) leads to

Comparison of eqs. (9) and (12) to (7) results in the identification

Furthermore, we obtain an interpretation of the infinite state y of the Caputo springpot as

Finally, we consider the energy of the mechanical equivalent system which is the potential energy stored in the springs of the Maxwell elements, i. e.

which can be reformulated with eqs. (13) and (14) as

4 Stability

4.1 Introduction

The mechanical representation and the potential energy expressed in terms of the infinite state y were derived more detailed in [15]. The energy expression in eq. (16) was used to prove stability of the equilibrium of a mass-spring-springpot system. In the following we want to extend this approach to consider stability of the same system when linear (anti-)damping is introduced, i. e. we regard the system (Figure 3)

(17) mq¨(t)=−dq˙(t)−cCDt0+αq(t)−kq(t),t≥0

with mass m, elongation q(t), spring coefficient k, springpot coefficient c, damping coefficient d and differentiation order α∈(0,1) and given initial functions φ1,φ2∈CB((−∞,0];R) such that

(18) q(t)=φ1(t),t≤0,q˙(t)=φ2(t),t≤0,

Figure 3:

Mass-spring-dashpot-springpot system.

where φ_i(t) = 0 for t ≤ t₀, i = 1,2. Again, with the help of eqs. (7) and (8) we reformulate eq. (17) as an FDE

(19) { q˙(t)=v(t),v˙(t)=−kmq(t)−dmq˙(t)           −cm∫t0−t0∫0∞μ1−α(ω)eωsdω vt(s)ds,

where

vt(s)=v(t+s),s∈(−∞,0]

and use the associated stability theory. For the cases d > 0 (damping) and d < 0 (anti-damping) we use different methods to prove Lyapunov stability of the equilibrium of eq. (17).

4.5 Linear anti-damping

For the case d < 0, whose physical interpretation is explained and motivated in Section 4.5.3, we expect the equilibrium of eq. (17) to remain stable only for certain values of d and it appears to be much more difficult to find Lyapunov functionals that yield stability criteria. Therefore, we start with the Laplace transform method to find sufficient conditions for stability and try to find similar conditions with the help of a Lyapunov functional. Finally, we compare the results.

4.5.1 Laplace transform method

Before we introduce the Laplace transform of eq. (17), we derive the Laplace transform of the fractional derivative. Therefore, consider the Laplace transform of (7)

(39) sLy(ω,t)(s)−y(ω,0)=−ωLy(ω,t)(s)+sQ(s)−q0,L CDt0+αq(t)(s)=∫0∞μ1−α(ω)Ly(ω,t)(s)dω.

with Laplace transform Q(s) of q(t) and initial value q(0) = q₀. Substitution of the first equation of (39) in the second results in

(40) L CDt0+αq(t)(s)=∫0∞μ1−α(ω)sQ(s)−q0+y(ω,0)ω+sdω.

We reformulate (40) with the help of the following relation.

Proposition 4.6

(41) ∫0∞μα(ω)ω+sdω=s−α,s∈C∖R−, α∈(0,1)

Proof.

Due to the relation for the Laplace transform of e−ωt

L{e−ωt}(s)=∫0∞e−ωte−stdt=∫0∞e−(ω+s)tdt=−1ω+se−(ω+s)t0∞=1ω+s

we obtain eq. (41) using the formula

Γ(α)Γ(1−α)=πsin(απ)

and Fubini’s Theorem as

∫0∞μα(ω)ω+sdω=sin(απ)π∫0∞ω−α∫0∞e−(ω+s)tdtdω=sin(απ)π∫0∞e−st∫0∞ω−αe−ωtdωdt=sin(απ)π∫0∞e−stΓ(1−α)tα−1dt=sin(απ)πΓ(1−α)Γ(α)s−α=s−α.

   □

This leads to the Laplace transform of the fractional derivative

(42) L CDt0+αq(t)(s)=sαQ(s)−sα−1q0+∫0∞μ1−α(ω)y(ω,0)ω+sdω.

Hence, we obtain the Laplace transform of eq. (17) as

(43) ms2Q(s)−sq0−v0=−kQ(s)−c(sαQ(s)−sα−1q0+∫0∞μ1−α(ω)y(ω,0)ω+sdω)−d(sQ(s)−q0)

with Laplace transform Q(s) of q(t) and initial values q(0) = q₀, q˙(0)=v0. Solving equation (43) for Q(s) leads to

(44) Q(s)=m(sq0+v0)ms2+ds+csα+k+csα−1q0−∫0∞μ1−α(ω)y(ω,0)ω+sdωms2+ds+csα+k+dq0ms2+ds+csα+k.

The inverse Laplace transform may be obtained integrating along a Hankel contour and using the residue theorem, similar as described in [27, 28]. We will accomplish the entire derivation later and see, that, similar to the classical case, stability of the equilibrium depends on the real part of the poles of the right-hand side of eq. (44), i. e. we consider the equation

(45) ms2+ds+csα+k=0.

Let s=reiθ, we obtain real and imaginary part of eq. (45) as

(46) mr2cos(2θ)+drcos(θ)+crαcos(αθ)+k=0,mr2sin(2θ)+drsin(θ)+crαsin(αθ)=0.

From (46) we want to derive conditions, such that the roots of eq. (45) are located in the left-half complex plane. Therefore, we first consider the critical case for stability θ=±π2, which turns (46) into

(47) −mr2+crαcosαπ2+k=0,dr+crαsinαπ2=0.

For fixed m, c, k and α the first equation of (47) has a unique solution r = r_* > 0, which may be inserted in the second equation to compute a critical d < 0 for stability

(48) dcrit=−csinαπ2r∗α−1.

By numerical solution of eq. (47), we obtain the critical negative damping parameter dcrit depending on the value of α∈(0,1) and the parameters m, c, k, see Figure 4. The value |dcrit| is a measure for the damping capability of the springpot. As expected, it holds that dcrit→0 for α→0, as in this case the springpot degenerates to a spring, which stores energy and dcrit→−c for α→1, as the springpot becomes a dashpot. The dependency of dcrit on α for α∈(0,1) may change drastically for different parameters and it is quite interesting that dcrit<−c can be achieved for certain values of α, i. e. a springpot can induce higher damping than a dashpot with the same coefficient. An example for this phenomenon is given in Section 4.5.3, where the knowledge (and tuning) of the system parameters m, k, c and α lead to the characterization of a Stribeck friction law. From the critical case for stability we now derive the following inequality conditions on r, such that a solution s=reiθ of (45) is located in the left-half complex plane. More specifically, we obtain the following proposition.

Figure 4:

Critical negative damping parameter depending on α∈(0,1) for different parameters.

Proposition 4.7.

Let the inequalities

(49) −mr2+crαcosαπ2+k≤0,

(50) dr+crαsinαπ2>0

have a non-empty solution set for r > 0. Then there exists a pair of complex conjugate roots s=reiθ, sˉ=re−iθ of (45) such that π2<θ<π2−α. Furthermore, there exists no solution outside the sectors θ∈π2,π and θ∈−π,−π2.

Remark 4.8.

To depict the solution set of inequalities (49) and (50), consider Figure 5 below, where r∗<R∗ has to hold.

Figure 5:

Representation of the solution set of inequalities (49) and (50).

Proof of Proposition 4.7.

From (46) we see that for each root s=reiθ of equation (45), its complex conjugate sˉ = re−iθ is another root. Therefore, we only consider 0 ≤ θ ≤ π. We examine the following cases.

1. Case 1: θ = 0In this case, eq. (45) degenerates to one equation

   (51) mr2+dr+crα+k=0.

   Using eqs. (49) and (50), we can estimate

   mr2+dr+crα+k>crα1+cosαπ2−sinαπ2+2k.

   As the function

   h(α):=1+cosαπ2−sinαπ2

   fulfills

   h(0)=2,h(1)=0,h′(α)=−π2sinαπ2+cosαπ2<0∀α∈(0,1),

   we obtain

   mr2+dr+crα+k>0

   and there exists no solution of eq. (51).

2. Case 2: θ = πThe second equation of (46) in this case reads as

   crαsin(απ)=0,

   which has no solution for α∈(0,1) except r = 0, which does not solve the first equation of (46)

   mr2−dr+crαcosαπ+k=0.

3. Case 3: 0<θ<π2Multiplying the first equation of (46) by cos(θ) and the second by sin(θ) sums up as

   (52) mr2cos(θ)+dr+crαcos((1−α)θ)+kcos(θ)=0.

   The left-hand side of eq. (52) may be estimated with (49) and (50) as

   mr2cos(θ)+dr+crαcos((1−α)θ)+kcos(θ)>crα(cos(θ)cosαπ2+cos((1−α)θ)−sinαπ2)+2kcos(θ)>0,

   because

   cos(θ)cosαπ2>0

   and

   cos((1−α)θ)−sinαπ2=cos((1−α)θ)−cos(1−α)π2>0.

   Hence, there is no solution of eq. (45) for θ∈0,π2 and α∈(0,1).

4. Case 4: π2<θ<πMultiplying the second equation of (46) by cos(αθ) and subtracting the first equation multiplied by sin(αθ) leads to

   mr2sin((2−α)θ)+drsin((1−α)θ)−ksin(αθ)=0,

   which may be solved for r > 0 as

   (53) r(θ)=−d2msin((1−α)θ)sin((2−α)θ)+d2msin((1−α)θ)sin((2−α)θ)2+kmsin(αθ)sin((2−α)θ).

   Furthermore, multiplying the first equation of (46) by sin(2θ) and subtracting the second equation multiplied by cos(2θ) leads to

   (54) drsin(θ)+crαsin((2−α)θ)+ksin(2θ)=0.

   As the first and the last term on the left-hand side of eq. (54) are negative for θ∈π2,π, the second term has to be positive to solve the equation, i. e.

   sin((2−α)θ)>0⇒π2<θ<π2−α.

   Now, we consider the left-hand side of eq. (54) as a function of θ

   g(θ):=dr(θ)sin(θ)+crα(θ)sin((2−α)θ)+ksin(2θ),θ∈π2,π2−α

   with r(θ) given by eq. (53). The function g is continuous for θ∈π2,π2−α and it holds that

   gπ2=dr+crαsin(2−α)π2=dr+crαsinαπ2>0

   as follows from (50). Furthermore, it can be seen from (53), that there exists a constant C > 0, such that

   limθ→π2−αr(θ)=limθ→π2−αCsin((2−α)θ)=∞,

   so that

   limθ→π2−αg(θ)=limθ→π2−α[ dCsin((2−α)θ)sin(π2−α)                    + c(Csin((2−α)θ))αsin((2−α)θ)                    + ksin(2π2−α) ]=−∞.

   Therefore, there is at least one root of g, i. e. one pair of complex conjugate solutions of eq. (45) such that θ∈π2,π2−α.

   □

As we have found conditions for solutions of the characteristic eq. (45) to be in the left-half complex plane, we want to prove asymptotic stability of the trivial solution by inverse Laplace transform using fundamental ideas of complex analysis. Therefore, we reformulate eq. (44) as

(55) Q(s)=ms+d+csα−1ms2+ds+csα+kq0+mms2+ds+csα+kv0−cms2+ds+csα+k∫0∞μ1−α(ω)y0(ω)ω+sdω.

Similar as in [27] we consider the function

(56) Ξ(s):=ms+d+csα−1ms2+ds+csα+k

and we compute the inverse Laplace transform ξ(t) of Ξ(s)=L{ξ(t)}(s). As

ξ(0)=lims→∞sΞ(s)=1,

we obtain

L{ξ˙(t)}(s)=sΞ(s)−ξ(0)=ms2+ds+csαms2+ds+csα+k−1=−kms2+ds+csα+k,

which, together with eq. (55) leads to the solution

(57) q(t)=q0ξ(t)−mkv0ξ˙(t)+ck∫0∞μ1−α(ω)y(ω,0)∫0te−ω(t−τ)ξ˙(τ)dτdω

of eq. (17) and we examine the asymptotic behavior of q from ξ and ξ˙. Therefore, we determine the inverse Laplace transform

ξ(t)=12πi∫σ−i∞σ+i∞Ξ(s)estds,Re(σ)>0

with the help of the residue theorem

(58) 12πi∫ΣΞ(s)estds=∑jResΞ(s)est,sj,

where s_j are the roots of eq. (45) and the closed curve Σ (Figure 6) is split up in six parts, such that

ξ(t)=∑jResΞ(s)est,sj−12πilimR→∞ε→0∫II−VIΞ(s)estds.

First, we compute the residues for a pair of complex conjugate roots s1, s2=s1ˉ of eq. (45). As s_1/2 are simple poles of Ξ, we obtain the residue by derivation of the denominator as

ResΞ(s)est,s1+ResΞ(s)est,s2=ms1+d+cs1α−12ms1+d+cαs1α−1es1t+ms2+d+cs2α−12ms2+d+cαs2α−1es2t.

As the two addends are conjugate, we obtain with s1=a+ib=reiθ

ResΞ(s)est,s1+ResΞ(s)est,s2=2Rems1+d+cs1α−12ms1+d+cαs1α−1es1t=2eatcos(bt)f1(r,θ)f3(r,θ)+2eatsin(bt)f2(r,θ)f3(r,θ).

with

Figure 6:

Curve Σ used for integration to apply the residue theorem.

(59) f1(r,θ)=2m2r2+d2+3mrdcos(θ)+(1+α)cdrα−1cos((1−α)θ)+(2+α)mcrαcos((2−α)θ)+c2αr2(α−1),f2(r,θ)=(2−α)mcrαsin((2−α)θ)+mdrsin(θ)+(1−α)cdrα−1sin((1−α)θ),f3(r,θ)=4m2r2+d2+4mdrcos(θ)+2cdαrα−1cos((1−α)θ)+4mcαrαcos((2−α)θ)+c2α2r2(α−1).

We continue considering the contribution of the integral along the paths II−VI to the value of ξ. There is no contribution of II, because

12πi∫IIΞ(s)estds=12πi∫π2πΞσ+ReiϕeσteRcos(ϕ)teiRsin(ϕ)tiReiϕdϕ≤12πC1(t)e−C2RtR⋅π2R⟶R→∞0,

with C1,C2>0 as Ξ(s)→0 for s→∞ and cos(ϕ) < 0 for ϕ∈(π2,π). The same argumentation holds for path VI. For path IV, we obtain

12πi∫IVΞ(s)estds=−12πi∫−ππΞεeiϕeεt(cos(ϕ)+isin(ϕ))iεeiϕdϕ⟶ε→00,

as sΞ(s)→0 for s→0. Finally, for III and V we obtain a contribution

12πi∫III,VΞ(s)estds=12πi∫0∞Ξωeiπ−Ξωe−iπe−ωtdω=1π∫0∞ImΞωeiπe−ωtdω

with

1πImΞωeiπ=−1πkcωα−1sin(απ)(mω2−dω+k)2+2cωαcos(απ)(mω2−dω+k)+c2ω2α=−μ1−α(ω)kc(mω2−dω+k)2+2cωαcos(απ)(mω2−dω+k)+c2ω2α.

This leads us to the inverse Laplace transform of Ξ

(60) ξ(t)=∑jodd(2eajtcos(bjt)f1(rj,θj)f3(rj,θj)+2eajtsin(bjt)f2(rj,θj)f3(rj,θj))+∫0∞μ1−α(ω)Z(ω)e−ωtdω

for roots sj/j+1=aj±ibj=rje±iθj of eq. (45) where

Z(ω)=kc(mω2−dω+k)2+2cωαcos(απ)(mω2−dω+k)+c2ω2α.

The asymptotic behavior of ξ is determined by the exponential functions in the first addends of eq. (60), which decay, as aj<0∀j, if the inequalities (49) and (50) have a non-empty solution set. Furthermore, the asymptotic behavior of the last term in eq. (60) may be estimated as follows. It holds that

(mω2−dω+k)2+2cωαcos(απ)(mω2−dω+k)+c2ω2α>(mω2−dω+k−cωα)2≥0.

Hence, Z is continuous and bounded in [0,∞) and by the mean value theorem, there exists C₃ > 0 such that

∫0∞μ1−α(ω)Z(ω)e−ωtdω=C3∫0∞μ1−α(ω)e−ωtdω=C3t−αΓ(1−α),

which leads to algebraic decay of order α for the last term in eq. (60) for t→∞. For ξ˙, we obtain the expression

(62) ξ˙(t)=∑jodd2eajtf3(rj,θj)((ajf1(rj,θj)+bjf2(rj,θj))cos(bjt)+(ajf2(rj,θj)−bjf1(rj,θj))sin(bjt))−∫0∞μ1−α(ω)Z(ω)ωe−ωtdω,

where the first terms again describe an exponentially decaying oscillation and the last term fulfills

−∫0∞μ1−α(ω)Z(ω)ωe−ωtdω=−αC3t−α−1Γ(1−α),

which again implies algebraic decay, this time of order 1 + α for t→∞. To conclude asymptotic stability of the trivial solution of eq. (17) from the asymptotic behavior of ξ, we still have to consider the last term in eq. (57). Therefore, we recall from eq. (8), that the initial infinite state y(ω,0) has the form

(63) y(ω,0)=∫t00eωτq˙(τ)dτ,

which may be estimated as

(64) |y(ω,0)|=∫t00eωτq˙(τ)dτ=eωτq(τ)t00−ω∫t00eωτq(τ)dτ≤2∥q∥∞+ω∫−∞0eωτdτ∥q∥∞=3∥q∥∞,

which is bounded if q∈CB((−∞,0];R). Furthermore, consider the reformulation

(65) ∫0te−ω(t−τ)ξ˙(τ)dτ=ddt∫0te−ω(t−τ)ξ(τ)dτ−ξ(0)e−ωt

of the inner integral in the last term of eq. (57). The last term in eq. (65) results in a term

∫0∞μ1−α(ω)y(ω,0)e−ωtdω≤3∥q∥∞t−αΓ(1−α)

in eq. (57). Substitution of the exponential terms of ξ in the last term of eq. (57) using eq. (65) leads to the estimation

∫0∞μ1−α(ω)y(ω,0)ddt∫0te−ω(t−τ)esjτdτdω≤3∥q∥∞∫0∞μ1−α(ω)1ω+sjsjesjt+ωe−ωtdω≤3∥q∥∞(sj∫0∞μ1−α(ω)ω+sjdωeRe(sj)t+∫0∞μ1−α(ω)e−ωtdω)=3∥q∥∞sjαeRe(sj)t+t−αΓ(1−α).

with roots s_j of eq. (45). For the algebraic decay part in ξ we obtain, again using the mean value theorem a constant C₄ > 0 and the term

∫0∞μ1−α(ω)y(ω,0)×ddt∫0te−ω(t−τ)∫0∞μ1−α(η)Z(η)e−ητdηdτdω≤C4ddt∫0∞μ1−α(ω)∫0te−ω(t−τ)τ−αΓ(1−α)dτdω=C4ddtCD0+αt1−αΓ(2−α)=C4|1−2α|t−2αΓ(2−2α).

in the last term of eq. (57). In summary, we obtain sufficient conditions (49) and (50) for global asymptotic stability of the equilibrium of eq. (17), from which we retrieve a Lyapunov functional in Section 4.5.2.

Remark 4.9.

It is even possible to obtain a purely exponential solution of eq. (17) without algebraic decay. Choose the initial function q(τ)=esjτ for τ∈(−∞,0] (which is not in CB((−∞,0];R)) for a root s_j of eq. (45). This leads to initial conditions

q˙(τ)=sjesjt,y(ω,0)=∫−∞0eωτsjesjτdτ=sjω+sj.

Using this function in the Laplace transform (55) leads to

Q(s)=ms+d+csα−1ms2+ds+csα+k+msjms2+ds+csα+k−csjms2+ds+csα+k∫0∞μ1−α(ω)(ω+s)(ω+sj)dω=ms+d+csα−1+msjms2+ds+csα+k−csj∫0∞μ1−α(ω)ω+sjdω−∫0∞μ1−α(ω)ω+sdω(s−sj)(ms2+ds+csα+k)=ms+d+csα−1+msj(s−sj)(s−sj)(ms2+ds+csα+k)−csjsjα−1−sα−1(s−sj)(ms2+ds+csα+k)=1s−sj.

Hence, the solution is

q(t)=esjt,∀t,

which shows that the integral term in eq. (42) should in general not be omitted.

4.5.2 Lyapunov functional

To formulate a Lyapunov functional, we will use the identities of the next proposition.

Proposition 4.10.

For α∈(0,1) and r > 0, the identities

(66) ∫0∞μ1−α(ω)ω2+r2dω=cos(απ2)rα−2

(67) ∫0∞μ1−α(ω)ωω2+r2dω =sin(απ2)rα−1

hold.

Proof.

Substitute η = ω² and dη=2ωdω in the integral and obtain

∫0∞μ1−α(ω)ω2+r2dω=sin(απ)π∫0∞ωα−1ω2+r2dω=sin(απ)2π∫0∞ηα2−1η+r2dη=sin(απ)2sinαπ2∫0∞μ1−α2(η)η+r2dη.

Using sine-double-angle formula and eq. (41), we directly obtain eq. (66). The proof of eq. (67) is analogous.   □

In the following, let r_* be the solution of the first equation of (47). We reformulate eq. (17) using eqs. (7), (36) and Proposition 4.10 as

(68) mq¨(t)=−kq(t)−dq˙(t)−cr∗2∫0∞μ1−α(ω)ω2+r∗2y(ω,t)dω−c∫0∞μ1−α(ω)ωω2+r∗2ωy(ω,t)dω=−k+ccosαπ2r∗αq(t)−d+csinαπ2r∗α−1q˙(t)+cr∗2∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω+c∫0∞μ1−α(ω)ω2+r∗2ωy˙(ω,t)dω.

Introducing the quantities

(69) k˜=k+ccosαπ2r∗α,d˜=d+csinαπ2r∗α−1,

and new coordinates

(70) q˜(t)=q(t)−ck˜r∗2∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω=kk˜q(t)+ck˜r∗2∫0∞μ1−α(ω)ω2+r∗2y(ω,t)dω

and

(71) v˜(t)=q˙(t)−cm∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)dω

we obtain the system

(72) { q ˜ . (t)=   v ˜ . (t), v ˜ . (t),=− k ˜ m q ˜ (t)− d ~ m v ˜ (t)− d ~ c m 2 ∫ 0 ∞ μ 1−α (ω) ω 2 + r * 2 ωy(ω,t)dω.

Note that the first equation in (72) holds, as r_* is a solution of the first equation of (47). Moreover, we consider the candidate Lyapunov functional

(73) V3(t,qt,vt)=m2v˜t2(0)+k˜2q˜t2(0)+d˜c2m∫0∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω

and prove inequality (26) for V₃ w.r.t. the functions q_t and v_t. Therefore, consider the split of the integral term in eq. (73)

∫0∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω=∫01μ1−α(ω)ω2+r∗2ωy2(ω,t)dω+∫1∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω

and use the mean value theorem for the first term and the inequality ω ≥ 1 in the second term to find a constant C˜>0, such that

(74) ∫0∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω≥C˜∫0∞μ1−α(ω)ω2+r∗2y2(ω,t)dω.

Moreover, we use Hölder’s inequality to obtain

(75) ∫0∞μ1−α(ω)ω2+r∗2y(ω,t)dω2≤∫0∞μ1−α(ω)ω2+r∗2dω⋅∫0∞μ1−α(ω)ω2+r∗2y2(ω,t)dω

and

(76) ∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)dω2≤∫0∞μ1−α(ω)ωω2+r∗2dω⋅∫0∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω.

Using the three inequalities above and Proposition 4.10, we can estimate (73) as

(77) V3(t,qt,vt)≥m2v˜t2(0)+d˜c4msinαπ2r∗α−1×∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)dω2+k˜2q˜t2(0)+d˜cC˜4mcosαπ2r∗α−2×∫0∞μ1−α(ω)ω2+r∗2y(ω,t)dω2

Finally, applying the general relation

(78) (a+b)2+γb2=γ1+γa2+a1+γ+1+γb2

for a,b,γ∈R, γ>0 on the first two and the last two terms of (77) using eqs. (70) and (71), we obtain inequality (26) for V₃. Furthermore, we compute the rate of V₃ as

V˙3=mv˜˙t(0)v˜t(0)+k˜q˜t(0)q˜˙t(0)+d˜cm∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)y˙(ω,t)dω.

Inserting the dynamics from (72), we obtain

(79) V˙3=−d˜v˜2(t)+d˜c2m2(∫0∞μ1−α(ω)ω2+r*2ωy(ω,t)dω)2 

(80) −d˜cm∫0∞μ1−α(ω)ω2+r*2ω2y2(ω,t)dω.

Again, using Hölder’s inequality leads to

(81) ∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)dω2≤∫0∞μ1−α(ω)ω2+r∗2dω⋅∫0∞μ1−α(ω)ω2+r∗2ω2y2(ω,t)dω

and together with eq. (66), we finally obtain

(82) V˙3≤−d˜v˜2(t)−d˜cm1−cmcosαπ2r∗α−2×∫0∞μ1−α(ω)ω2+r∗2ω2y2(ω,t)dω

where, due to eq. (47)

1−cmcosαπ2r∗α−2>0⇔mr∗2−ccosαπ2r∗α=k>0.

In summary, we found a Lyapunov functional V₃, such that V˙3≤0, which has the form of an energy functional w.r.t. the new coordinates q˜t and v˜t but it has the disadvantage that we cannot prove asymptotic stability of the trivial solution of eq. (17) with the help of V₃ as in the classical case of the damped linear oscillator using the total mechanical energy as a Lyapunov function. Hence, we introduce another Lyapunov functional V₄, which is related to the functional V₂ in Section 4.4, to prove the following proposition.

Proposition 4.11.

Let m,k,c > 0, α∈(0,1) and let r = r_* > 0 be the solution of

−mr2+crαcosαπ2+k=0.

Let d∈R be such that the inequality

dr∗+cr∗αsinαπ2>0

holds. Then the trivial solution of eq. (17) is asymptotically stable.

Proof.

We introduce the Lyapunov functional

(83) V4(t,qt,vt)=12mv˜t2(0)+12k˜q˜t2(0)+d˜24mq˜t2(0)+12d˜q˜t(0)v˜t(0)+d˜c2m∫0∞μ1−α(ω)ω2+r∗2ωy2(ω,t)dω+d˜2c4m2∫0∞μ1−α(ω)ω2+r∗2ω2Y2(ω,t)dω−d˜2c4m2ck˜r∗2∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω2

and check the conditions in Theorem 4.3. Once again, we use Hölder’s inequality to obtain

∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω2≤∫0∞μ1−α(ω)ω2+r∗2dω⋅∫0∞μ1−α(ω)ω2+r∗2ω2Y2(ω,t)dω.

Together with (76), estimation (38) and Proposition 4.10 we obtain

V4(t,qt,vt)≥14mv˜t2(0)+12k˜q˜t2(0)+d˜c2msin(απ2)r∗α−1∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)dω2+d˜2c4m2ck˜r∗2k˜ccos(απ 2)r∗α−1×∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω2.

All coefficients in this estimation are positive and, again using eqs. (70), (71) and (78), we obtain relation (26) for V₄. Furthermore, we compute the rate of V₄ as

V˙4=mv˜t(0)v˜˙t(0)+k˜q˜t(0)q˜˙t(0)+d˜2v˜t(0)q˜˙t(0)+d˜2q˜t(0)v˜˙t(0)+d˜22mq˜t(0)q˜˙t(0)+d˜cm∫0∞μ1−α(ω)ω2+r∗2ωy(ω,t)y˙(ω,t)dω+d˜2c2m2∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)ωY˙(ω,t)dω−d˜2c2m2ck˜r∗2∫0∞μ1−α(ω)ω2+r∗2ωY(ω,t)dω×∫0∞μ1−α(ω)ω2+r∗2ωY˙(ω,t)dω

Inserting the dynamics and using eqs. (81) and (74), we obtain

V˙4≤−d˜2v˜2(t)−d˜cm1−cmcosαπ2r∗α−2×∫0∞μ1−α(ω)ω2+r∗2ω2y2(ω,t)dω−k˜d˜2mq˜2(t)−d˜2cC˜2m2∫0∞μ1−α(ω)ω2+r∗2y2(ω,t)dω.

Again, with the help of the estimations (75), (81) and eqs. (70), (71) and (78), one can prove that V˙4 fulfills inequality (29) w.r.t. q and v. This completes the proof.   □

Remark 4.12.

1. As a special case, Proposition 4.11 proves asymptotic stability of the trivial solution of eq. (17) for the case d = 0. Previously in Section 4.3, with the help of the energy functional we could only prove stability but not attractivity of the trivial solution.

2. The conditions for asymptotic stability in Proposition 4.11 are equivalent to the necessary and sufficient conditions obtained by the eigenvalue analysis in Proposition 4.7. In that sense, the choice of the functionals V₃ and V₄ is optimal.

4.5.3 Example: Stick-slip oscillator [293031]

In this section, we describe a model of a mechanical system, where effective negative linear damping occurs in the linearization of the equation of motion around an equilibrium and we give sufficient conditions for local asymptotic stability of the equilibrium using Proposition 4.11. Consider a mass m suspended by a spring with spring coefficient k, and a springpot with coefficient c and differentiation order α∈(0,1), which is sliding on a conveyor belt as in Figure 7. By q we denote the displacement of the mass. The belt moves with a constant velocity vdr>0 in the direction of q and we assume friction between the mass and the belt, which leads to a friction force F_T, such that the equation of motion reads as

(84) mq¨(t)=FT−cCDt0+αq(t)−kq(t).

Figure 7:

Fractionally damped oscillator with dry friction and graph of a set-valued force law describing the Stribeck effect.

For the friction force F_T in the slip phase we consider the force law

(85) FT=−μ(vrel)FNsign(vrel),vrel≠0

where

(86) vrel(t)=q˙(t)−vdr

is the relative velocity between mass and belt,

(87) FN=mg

is the normal force acting on the mass and µ = µ(v) is the friction coefficient depending on the relative velocity, where the function µ increases (at least) for small negative values of vrel, i. e.

(88) μ′(−|v|)>0,|v|≪1,

which is known as the Stribeck effect. We consider the slip equilibrium q^* of eq. (84)

(89) 0=FT−kq*,

(90) −FT=μ(−vdr)FNsign(−vdr)=−μ(−vdr)mg,

which implies

(91) q∗=μ(−vdr)mgk

We introduce a new coordinate

(92) qˉ=q−q∗,

such that we reformulate eq. (84) in terms of qˉ as

(93) mqˉ¨(t)+cCDt0+αqˉ(t)+k(qˉ(t)+q∗)=−μ(qˉ˙(t)−vdr)mgsign(qˉ˙(t)−vdr).

Linearizing the right-hand side of eq. (93) near the equilibrium leads to

(94) mqˉ¨(t)+cCDt0+αqˉ(t)+k(qˉ(t)+q∗)=μ(−vdr)mg+μ′(−vdr)mgqˉ˙(t)+Oqˉ˙2,

which together with eq. (91) leads to the linearized equation

(95) mqˉ¨(t)−μ′(−vdr)mgqˉ˙(t)+cCDt0+αqˉ(t)+kqˉ(t)=Oqˉ˙2

Hence, using Proposition 4.11, we obtain the condition

(96) cr∗α−1sinαπ2>μ′(−vdr)mg

for local asymptotic stability of the slip equilibrium q^*, where again r_* is the solution of

−mr2+crαcosαπ2+k=0.