On the Selection and Meaning of Variable Order Operators for Dynamic Modeling

Abstract

We review the application of differential operators of noninteger order to the modeling of dynamic systems. We compare all the definitions of Variable Order (VO) operators recently proposed in literature and select the VO operator that has the desirable property of continuous transition between integer and non-integer order derivatives. We use the selected VO operator to connect the meaning of functional order to the dynamic properties of a viscoelastic oscillator. We conclude that the order of differentiation of a single VO operator that represents the dynamics of a viscoelastic oscillator in stationary motion is a normalized phase shift. The normalization constant is found by taking the difference between the order of the inertial term (2) and the order of the spring term (0) and dividing this difference by the angular phase shift between acceleration and position in radians (𝜋), so that the normalization constant is simply 2/𝜋.

1. Introduction

The integer order differential operators of classical calculus (such as the first or second order derivatives) are familiar to anyone who has an active interest in understanding dynamic systems. These differential operators are used to formulate models that accurately describe the majority of physical phenomena and are ubiquitous in the mathematical description of dynamic behavior. However effective these integer order differential operators are in general, there are more complex systems that are better characterized by dynamic behavior that lies in between the normal integer order description. A case in point is the so-called “viscoelastic” behavior, which has characteristics of both elastic (order zero) and viscous (order one) elements. It is thus natural to assume that differential operators of noninteger order, such as a 0.25, 0.50, or 0.75 would provide a convenient mathematical description to analyze these intermediate behaviors. The study of these noninteger differential operators falls under the general subject of what became known as Fractional Calculus, though the orders studied are not strictly limited to rational numbers.

A further generalization of the concept of noninteger order derivatives that is applicable to more complex systems is that of a derivative of varying order. One can find systems where the order of dynamics associated with each element is a function of time or frequency or position or any derivative of the position vector [1–3]. The objective of this work is to identify the most appropriate definition of a variable-order (VO) operator for modeling dynamic systems and to assign to the order of the derivative a physical meaning that will facilitate the understanding of its use in problems of vibration and control. First, we compare all VO operator definitions recently proposed in literature in order to select a definition that is better suited for modeling purposes. We then use the familiar example of viscoelastic harmonic oscillators to connect the order of the derivative to the dynamic properties of the oscillators.

Unlike ordinary derivatives, noninteger order derivatives are integrodifferential operators with either a power-law (in the case of fractional derivatives) or a variable-exponent (for VO derivatives) kernel. Thus, noninteger derivatives are nonlocal by definition and are ideally suited for modeling systems characterized by nonlocal (or memory-laden) behavior. Multiple definitions for a fractional derivative have been proposed, and there is no straightforward geometric or physical interpretation for the meaning of a fractional derivative, although a few interpretations have been proposed [4]. The subject of Fractional Calculus has developed rapidly, especially since in the last four decades, with quite a few recent books dedicated exclusively to the subject (see, e.g., [4–8]). Some applications that involve fractional derivatives include particle motion at small but finite Reynolds numbers [9–12], viscoelastic constitutive equations [13–16], transport dynamics with anomalous diffusion [17], and fractional order control systems [4, 18].

The concept of a VO operator is a much more recent development and is less widely-known. The order of the derivative/integral is allowed to vary over the domain of interest, thus they are suitable for modeling systems with evolving dynamics. Such systems include deformation of viscoelastic materials [19–22] and the mechanics of variable viscoelastic oscillators [1, 3]. Similar to the state of affairs in Fractional Calculus, multiple definitions of a VO derivative have been suggested [1, 2, 19, 23, 24], each preferred for different reasons. Since the kernel of the VO operators have a variable-exponent, analytical solutions to VO differential equations (VODEs) are more difficult to obtain, and have not been the focus of much attention. However, numerical solutions for the VODEs that have been formulated with the various VO operators have been developed [1–3]. In addition, rather than seeking a constant or multiple constants for the order of the derivative in the model that would accurately represent the data, in the VO case a function needs to be determined through mathematical and/or physical arguments. Previous applications of VO derivatives have either explicitly chosen the function [1–3] or found an approximation numerically [19, 20, 22]. This complication can be lessened if the functional form of the order is determined through physical arguments as in [21], where a model based on statistical mechanics was developed to describe the compression of a viscoelastic material.

The aim of this work is to compare the Variable Order operator definitions that have been proposed and to select the operator with the characteristics that are critical for the success of a dynamic model. We compare the various VO operators based on a very simple criteria: the VO operator must return the correct fractional derivative that corresponds to the argument of the functional order. In other words, if the argument 𝑞(𝑡) is equal to, say, 0, 0.50 or 1, then the function itself (zero-order derivative), the half-derivative, or the first derivative must be returned as the output of the operator. We will see that only two definitions previously proposed satisfy this elementary requirement. Of the two operator definitions that satisfy this property, one is more efficient from the numerical standpoint, and is therefore adopted in the remainder of this work. The appropriate operator then is used to study the somewhat familiar problem of a harmonically forced oscillator with viscoelastic damping. The goal of this second part of this work is to illustrate how a familiar problem in dynamics can be used to understand the meaning of a VO operator, and to understand how the dynamics in this familiar problem is affected by the physical parameters of the system using a VO analysis.

The next section presents an overview of the various VO operator definitions and a brief comparison of the VO operators applied to a harmonic and other bounded function. Subsequent to selecting the operator, we propose a VO model for the harmonically forced oscillator with viscoelastic damping of order 𝑝 (0<𝑝<1) and conduct a stationary analysis that yields a very concrete meaning to the order of the operator.

2. Variable Order Operators

The VO operator definitions that have been proposed are either direct extensions of the fractional calculus definitions or generalizations that arise from Laplace or Fourier transformations. In the direct extension approach, the constant exponent in the fractional operator is replaced with a function. For example, a VO integral is defined in [23] as

𝒟^{−𝑞(𝑡)𝑐+}𝑓(𝑡)=1Γ[𝑞(𝑡)]^𝑡𝑐(𝑡−𝜎)^{𝑞(𝑡)−1}𝑓(𝜎)d𝜎.(2.1) When 𝑞(𝑡)=𝛼= constant, then the 𝛼th-order fractional integral is recovered. Other definitions can be formulated by changing the form of the argument of the exponent to be 𝑞=𝑞(𝜎) or 𝑞=𝑞(𝑡−𝜎)and considering the Gamma function under the integral sign [24, 25]:

¹D^{0−𝑞(𝑡)}𝑓(𝑡)=^𝑡0(𝑡−𝜎)^{𝑞(𝑡)−1}Γ[𝑞(𝑡)]𝑓(𝜎)d𝜎,(2.2)₂D^{0−𝑞(𝑡)}𝑓(𝑡)=^𝑡0(𝑡−𝜎)^{𝑞(𝜎)−1}Γ[𝑞(𝜎)]𝑓(𝜎)d𝜎,(2.3)₃D^{0−𝑞(𝑡)}𝑓(𝑡)=^𝑡0(𝑡−𝜎)^{𝑞(𝑡−𝜎)−1}Γ[𝑞(𝑡−𝜎)]𝑓(𝜎)d𝜎.(2.4) In the cases above, the lower terminal is set equal to 0, and it is assumed that 𝑓(0)=0 for 𝑡<0. Since (2.3) and (2.4) involve the variable of integration within the exponent, then this implies memory in the order, with the past states having a stronger effect on the order for definition (2.4) [24]. Also, the full convolution form of (2.4) enables use of the convolution properties to study the operator.

Similarly, a VO derivative definition can be obtained by directly substituting 𝑞=𝑞(𝑡) in the Riemann-Liouville fractional derivative definition [23] valid for 0<𝑞<1:

D^{𝑞(𝑡)𝑐+}𝑓(𝑡)=1Γ[1−𝑞(𝑡)]dd𝑡^𝑡𝑐𝑓(𝜎)(𝑡−𝜎)^𝑞(𝑡)d𝜎.(2.5) Further definitions are obtained by taking 𝑚 first-order derivatives of the integrals defined in (2.2)–(2.4)

D^0𝑞(𝑡)𝑓(𝑡)=d^𝑚d𝑡^𝑚^𝑡0(𝑡−𝜎)^{𝑚−𝑞(𝑡,𝜎)−1}Γ[𝑞(𝑡,𝜎)]𝑓(𝜎)d𝜎,(2.6) where 𝑚−1<𝑞(𝑡)<𝑚. Caputo-type VO operators are defined by taking the derivative of the function under the integrals

D^0𝑞(𝑡)𝑓(𝑡)=^𝑡0(𝑡−𝜎)^{𝑚−𝑞(𝑡,𝜎)−1}Γ[𝑞(𝑡,𝜎)]𝑓^(𝑚)(𝜎)d𝜎,(2.7) where 𝑓^(𝑚) denotes the 𝑚th integer order derivative of 𝑓(𝑡). VO operators based on other fractional derivative definition forms have also been proposed. Samko and Ross [23] introduce a VO operator based on the Marchaud fractional derivative:

D^{𝑐𝑞(𝑡)}𝑓(𝑡)=𝑓(𝑡)Γ[1−𝑞(𝑡)](𝑡−𝑐)^𝑞(𝑡)+𝑞(𝑡)Γ[1−𝑞(𝑡)]^𝑡𝑐𝑓(𝑡)−𝑓(𝜎)(𝑡−𝜎)^1+𝑞(𝑡)d𝜎,(2.8) where 0<Re𝑞(𝑡)<1.

Using a different approach, Coimbra [1] begins with the Laplace transform of the Caputo operator to obtain a VO operator definition. Treating 𝑞(𝑡) as a running parameter and inverting back into the time domain yields the following definition valid for 𝑞(𝑡)<1:

𝐃^𝑞(𝑡)𝑓(𝑡)=1Γ[1−𝑞(𝑡)]^𝑡0+(𝑡−𝜎)^{−𝑞(𝑡)}𝑓⁽¹⁾(𝜎)d𝜎+(𝑓(0+)−𝑓(0−))𝑡^{−𝑞(𝑡)}Γ[1−𝑞(𝑡)].(2.9) Extension of definition (2.9) to values of 𝑞(𝑡) larger than unity is possible as long as the higher-order derivatives are defined and integrable. Although this definition began with the Laplace transform of the Caputo operator, it is not strictly a Caputo generalization because of the addition of the initial condition term. As a result, D^𝑞(0)𝑓(0)≠0 for any function 𝑓(𝑡) as is the the case for a Caputo-based operator. Soon et al. [3] show that a properly weighted sum of fractional order derivatives terms approximates this single term VO operator when a large number of terms are used, which implies convergence of both the operator and the numerical method.

Although definitions (2.8) and (2.9) were defined independently through different methods, they are similar. After integrating (2.8) by parts and simplifying, we arrive at

1Γ[1−𝑞(𝑡)]^𝑡𝑐𝑓⁽¹⁾(𝜎)(𝑡−𝜎)^𝑞(𝑡)d𝜎+𝑓(𝑐)Γ[1−𝑞(𝑡)](𝑡−𝑐)^𝑞(𝑡).(2.10) The difference between (2.9) and (2.10) lies in the term that is evaluated at the lower terminal. For situations in which the system is assumed to be in dynamical equilibrium for 𝑡<0 such that 𝐷²𝑓(𝑡<0)=0 for any 𝑓(𝑡), definitions (2.9) and (2.10) return the same result. However, if 𝑓(0) is a true constant, such that 𝑓(0⁻)=𝑓(0⁺)=𝑎, then (2.9) would return 0 for the derivative, whereas (2.10) will not.

In summary, there are a total of nine VO operator definitions to be compared:

¹D^0𝑞(𝑡)𝑓(𝑡)=1Γ[1−𝑞(𝑡)]dd𝑡^𝑡0𝑓(𝜎)(𝑡−𝜎)^{−𝑞(𝑡)}d𝜎,(2.11)₂D^0𝑞(𝑡)𝑓(𝑡)=𝑓(𝑡)Γ[1−𝑞(𝑡)](𝑡)^𝑞(𝑡)+𝑞(𝑡)Γ[1−𝑞(𝑡)]^𝑡0𝑓(𝑡)−𝑓(𝜎)(𝑡−𝜎)^1+𝑞(𝑡)d𝜎,(2.12)₃₋₅D^0𝑞(𝑡)𝑓(𝑡)=dd𝑡^𝑡0(𝑡−𝜎)^{𝑞(𝑡,𝜎)−1}Γ[𝑞(𝑡,𝜎)]𝑓(𝜎)d𝜎,(2.13)₆₋₈D^0𝑞(𝑡)𝑓(𝑡)=^𝑡0(𝑡−𝜎)^{𝑞(𝑡,𝜎)−1}Γ[𝑞(𝑡,𝜎)]𝑓⁽¹⁾(𝜎)d𝜎,(2.14)₉D^{0−𝑞(𝑡)}𝑓(𝑡)=1Γ[1−𝑞(𝑡)]^𝑡0+(𝑡−𝜎)^{−𝑞(𝑡)}𝑓⁽¹⁾(𝜎)d𝜎+(𝑓(0+)−𝑓(0−))𝑡^{−𝑞(𝑡)}Γ[1−𝑞(𝑡)],(2.15) where 𝑞(𝑡,𝜎) in definitions (2.13) and (2.14) signify the three arguments: 𝑞(𝑡,𝜎)=𝑞(𝑡), 𝑞(𝑡,𝜎)=𝑞(𝜎), and 𝑞(𝑡,𝜎)=𝑞(𝑡−𝜎). Each of the above definitions is defined for real derivative orders between 0 and 1. The value at all the lower terminals is set to 0, since we are interested in applying the operators to physical processes not necessarily at steady state (we examine a stationary problem in the next section). As is the case with fractional derivatives, there is no single VO derivative (or integral) definition that is widely considered to be the “correct" definition. Samko and Ross prefer definition (2.12) because the operator retains the symmetry on power functions that is found in the case of constant orders, that is:

D^{𝑞(𝑡)𝑐+}+(𝑡−𝑐)^𝛼=Γ[𝛼+1]Γ[𝛼+1−𝑞(𝑡)](𝑡−𝑐)^{𝛼−𝑞(𝑡)}(2.16) for 0<Re𝑞(𝑡)<1 and 𝛼>−1 [6]. Lorenzo and Hartley prefer the full convolution VO integral definition (2.4) because it satisfies the index rule for certain functions [25] and is time-invariant [24].

The approach chosen here is to determine which operator when acting upon a function returns the fractional derivative of the function at the corresponding time. This is an important characteristic from the aspect of physical modeling because it signifies that the operator yields a continuous transition of all orders of differentiation between integer orders. Thus, a smooth transition from zero-order dynamics to first-order dynamics is possible. For a qualitative comparison of the the VO operators, we look at the 𝑞(𝑡)=𝑡 derivative of two bounded functions: sin(2𝜋𝑡), and erfc(𝑡). The 𝑞(𝑡)=𝑡 derivative of both functions is computed numerically using a product trapezoidal rule to evaluate the convolution integrals [3, 26]. Plots of the 𝑡 derivative of sin(2𝜋𝑡) and erfc(𝑡) are shown in Figures 1-2. The 0, 0.25, 0.50, 0.75, and 1st-order derivatives of the sinusoidal function are shown for comparison since both the Riemann-Liouville and Caputo fractional derivatives return the same result. We show only the 0- and 1st-order derivatives of the erfc function since the different fractional derivative definitions yield different results. Also, note that Samko's Marchaud-based definition and Coimbra's definition coincide for these conditions, and therefore only definition (2.15) is plotted.

(a)

(b)

(c)

Figure 1 

Plots of the 𝑡 derivatives of sin(2𝜋𝑡). The points are the 𝛼=0,0.25,0.50,0.75, and 1st-order fractional derivatives at 𝑡=𝛼. (a) Definition (2.13) with 𝑞(𝑡,𝜎)=𝑞(𝑡) (thin line), 𝑞(𝑡,𝜎)=𝑞(𝜎) (medium line), and 𝑞(𝑡,𝜎)=𝑞(𝑡−𝜎) (thick line). Note that none of the definitions match the 1st-order derivative at 𝑡=1. (b) Caputo-type VO operators (2.14) with 𝑞(𝑡,𝜎)=𝑞(𝑡) (thin line), 𝑞(𝑡,𝜎)=𝑞(𝜎) (medium line), and 𝑞(𝑡,𝜎)=𝑞(𝑡−𝜎) (thick line). Note that the variant of (2.14) with argument 𝑡 matches the corresponding 𝛼 fractional derivatives at 𝑡=𝛼. (c) Definition (2.11) (thick line) and Coimbra's operator (2.15) (thin line). The 𝑡-derivative defined by (2.15) is equivalent to the corresponding fractional derivatives at all the points.

(a)

(b)

(c)

Figure 2 

Plots of the 𝑡 derivatives of erfc(𝑡). The points are the 0- and 1st-order derivatives at 𝑡=0 and 𝑡=1. Also shown for comparison is the 𝑡 derivative from Coimbra's operator (2.15) (dashed line). (a) Riemann-Liouville type operator (2.13) with 𝑞(𝑡,𝜎)=𝑞(𝑡) (thin line), 𝑞(𝜎) (medium line), 𝑞(𝑡−𝜎) (thick line). (b) Caputo-type definition (2.14) with 𝑞(𝑡,𝜎)=𝑞(𝑡) (thin line), 𝑞(𝜎) (medium line), 𝑞(𝑡−𝜎) (thick line). In this case the operator (2.14) with exponent 𝑞(𝑡) does not match Coimbra's operator, but is equivalent to the first derivative at 𝑡=1. (c) Definition (2.11) (thin line).

As expected, the Caputo-based VO operators have values of 0 at 𝑡=0, similar to the fractional derivative case. They do not return the zeroth-order derivative at 𝑡=0 for any functions that deviate from this (see Figure 2). Definition (2.14) is equivalent to (2.15) for the sine function, since 𝑓(𝑡) is continuous at 𝑡=0 and 𝑞(𝑡) is only a function of 𝑡. The Riemann-Liouville and Caputo type operators with the argument (𝑡−𝜎), also are shown to be equivalent because 𝑓(0)=0. This is analogous to the similarity of the Riemann-Liouville and Caputo fractional derivatives of functions when 𝑓⁽⁰⁾(0)=0 for 0≤𝑞<1 [4]. Definitions (2.12) and (2.15) are the only operators that have the desirable property of returning the corresponding 𝑞𝑡ℎ order fractional derivative of 𝑥(𝑡) when 𝑞(𝑡)=𝑝 for both the sine and erfc funcitons. However, the convergence of (2.12) to that of (2.15) is slower due to the stronger singularity that must be evaluated in the convolution integral. Also, in the case of a true constant function 𝑓(𝑡) where 𝑓(0⁻)=𝑓(0⁺), (2.12) would not return 0 as the derivative, so we conclude that definition (2.15) is preferable for modeling dynamic systems.

Note that the operator defined in (2.15) is dynamically consistent with the causal behavior of the initial conditions. In other words, when 𝑥(𝑡) is a true constant from −∞ to the initial time (𝑡=0⁺), the operator in (2.15) returns zero for all values of 𝑞(𝑡). However, if 𝑓(𝑡) is not continuous between 𝑡=0⁻ and 𝑡=0⁺, the operator returns the appropriate Heaviside contribution to the integral value of 𝒟^𝑞(𝑡)𝑓(𝑡). In accordance with this causal definition, we take the value of the physical variable 𝑓(𝑡) to be identically null from −∞ to 0⁻ as a representation of dynamic equilibrium. A nonzero initial condition is treated as a Heaviside function at 𝑡=0, and therefore included in the second term of the definition of the operator (2.15).

Through the direct comparison of the various proposed VO operator definitions, we selected the operator that has fundamental characteristics that are desirable for physical modeling. Definition (2.15) proposed in [1] represents a continuous transition between the integer order derivatives, and returns a zero value for the derivative of a function that is constant from −∞<𝑡<∞; so we select it as the most appropriate definition. Now with the chosen operator, we proceed to connect the behavior of a VO operator with a physical quantity that is characteristic of a selected memory-laden system.

3. Stationary Analysis for Viscoelastic Oscillators

One of the drawbacks of VO modeling is that to date there is no clear physical understanding of what a VO derivative represents. The objective of this section is to illustrate how a variable order differential operator may be used to understand a familiar problem: the stationary analysis of a constant order viscoelastic oscillator. We will use a single VO operator of order 𝑞 to replace multiple terms of constant order differential operators, including the viscoelastic term of constant order 𝑝. We seek an analytical expression for 𝑞 to examine the effects of the parameters of the system on the dynamics of the oscillator. The exact expression for 𝑞(𝜔) is obtained from the stationary analysis of the problem. In order words, we look for two functions 𝑞 that would allows us to replace the multiterm differential equation describing the steady motion of the oscillator with a single-term variable order equation. Note that the objective of this section is not to rehash the analysis of the constant order viscoelastic oscillator, rather we seek to find meaning for a VO operator by comparing its order with the physical parameters in a well-known problem. The reader who is interested in the details of the dynamics of constant order viscoelastic oscillators should consult [4, 27, 28].

The equation of motion for the constant order viscoelastic oscillator is 𝑚𝐷²𝑥(𝑡)+𝑐𝐷^𝑝𝑥(𝑡)+𝑘𝐷⁰𝑥(𝑡)=𝐹(𝑡),(3.1) where 𝐹(𝑡)=𝐹₀cos(Ω𝑡). When 𝑝=1, the system is an oscillator with viscous damping, and for 0<𝑝<1, the system is said to be characterized by viscoelastic damping. The equation is recast in dimensionless form using the following scaled parameters:

𝑥̃=𝑥𝐿_𝑐,𝑡̃=𝑡𝜔_𝑛,Ω=Ω𝜔_𝑛,(3.2) where 𝐿_𝑐 is a characteristic length (or amplitude of the motion) and 𝜔_𝑛=√𝑘/𝑚 is the undamped natural frequency of the system. The dimensionless equation of motion is

𝐷²𝑥̃+𝜉𝐷^𝑝𝑥̃+𝑥̃=𝐹_𝑜cosΩ𝑡̃,(3.3) where 𝐹₀=𝐹₀/𝑚𝜔^2𝑛𝐿_𝑐 is the dimensionless amplitude of the forcing, and 𝜉 is the damping ratio and is

𝜉=𝑐𝑘^𝑝/2−1𝑚^𝑝/2.(3.4)

Since we are concerned only with stationary behavior, 𝜂 and 𝑞 are not functions of time, so we look for a relationship such as 𝜂𝑝,Ω,𝜉𝐷^𝑞𝑝,Ω_,𝜉𝑥̃=𝐷²+𝜉𝐷^𝑝+𝐷⁰𝑥̃=Re𝐹₀expiΩ𝑡̃.(3.5) The stationary VO derivative of 𝑒ⁱΩ𝑡̃ is

𝐷^𝑞(𝑝,Ω_,𝜉)𝐴𝑒ⁱΩ𝑡̃=iΩ^𝑞(𝑝,Ω_,𝜉)𝐴𝑒ⁱΩ𝑡,̃(3.6) where definition (2.15) with a lower terminal of −∞ is used since we are dealing with a stationary problem where the initial conditions are irrelevant. Rewriting (3.6) with (3.3) as the stationary solution exp(iΩ𝑡)̃ yields

𝜂iΩ^𝑞(𝑝,Ω_,𝜉)=1−Ω²+iΩ^𝑝𝜉.(3.7) We now equate the real and imaginary parts of the above equation to arrive at

𝑞𝑝,Ω,𝜉=tan⁻¹⎛⎜⎜⎝𝜉Ω^𝑝sin((𝜋/2)𝑝)1−Ω²+Ω^𝑝𝜉cos((𝜋/2)𝑝)⎞⎟⎟⎠2𝜂𝑝,Ω,𝜉=1𝜋,Ω^𝑞1−Ω²²+Ω^𝑝𝜉²+21−Ω²Ω^𝑝𝜉cos(𝜋/2)𝑝.(3.8) The stationary solution for (3.5) can be written as 𝑥̃(𝑡)̃=𝐴cos(Ω𝑡̃−𝜙) where the amplitude 𝐴 and phase shift 𝜙 are

𝐴=𝐹⁰𝜙1−=Ωta²n²₋₁+Ω^𝑝𝜉²+21−Ω²Ω^𝑝𝜉cos(𝜋/2)𝑝,⎛⎜⎜⎝𝜉Ω^𝑝sin((𝜋/2)𝑝)1−Ω²+Ω^𝑝𝜉cos((𝜋/2)𝑝)⎞⎟⎟⎠.(3.9) The procedure used to obtain the functional form for 𝑞 and 𝜂 is easily extendable to systems that consist of multiple viscoelastic terms. For 𝑛 terms, the expression for 𝑞 is

𝑞=tan⁻¹⎛⎜⎜⎝∑^𝑛𝑘=1𝜉_𝑘Ω^𝑝𝑘sin(𝜋/2)𝑝_𝑘1−Ω²+∑^𝑛𝑘=1𝜉_𝑘Ω^𝑝𝑘cos(𝜋/2)𝑝_𝑘⎞⎟⎟⎠2𝜋.(3.10) Similarly, the expression for 𝜂 is

𝜂=1Ω^𝑞⎷1−Ω²+^𝑛_𝑘=1𝜉_𝑘Ω^𝑝𝑘cos𝜋2𝑝^𝑘²+^𝑛_𝑘=1𝜉_𝑘Ω^𝑝𝑘sin𝜋2𝑝^𝑘².(3.11)

Equations (3.8)-(3.9) clearly show that 𝑞 represents a scaled phase shift and 𝜂 is the scaled ratio of the amplitude of the forcing to the response. Thus, for the stationary solution of the oscillator the order of the derivative is connected to a physical quantity that is characteristic of the system. The multiterm equation in frequency that represents the stationary motion of a viscoelastic operator can be replaced by a single term parametric operator in frequency, where both the order of the derivative and the scaling function 𝜂 have physical meaning. A phase shift of 𝜋/2 implies that the response is proportional to the velocity and hence to the 1st-order derivative, while a phase shift of 𝜋 implies that it is proportional to the acceleration or to the 2nd-order derivative. The order of the VO operator that captures the whole dynamics of the systems is thus naturally connected with the phase shift between the response and the forcing.

Plots of 𝑞 and 𝜂 versus Ω for damping orders of 𝑝=0.25,0.50,0.75, and 1 and various values of the damping ratio are shown in Figures 3, 4, 5, and 6. Also shown are the maps of 𝑞 and 𝜂.

(a)

(b)

(c)

(d)

Figure 3 

Plots of 𝑞 and 𝜂 for 𝑝=0.25. (a) 𝑞(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). Note that 𝑞 aysmptotes to 2 as Ω increases. (b) Map of 𝑞(Ω,𝜉) where lower values are represented by darker colors. When Ω=1, 𝑞=0.25 for any value of 𝜉. (c) 𝜂(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (d) Map of 𝜂(Ω,𝜉). The lowest values of 𝜂 occur near Ω=1 and for small damping ratio.

(a)

(b)

(c)

(d)

Figure 4 

Plots of 𝑞 and 𝜂 for 𝑝=0.50. (a) 𝑞(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). The asymptotic behavior is similar to the case when 𝑝=0.25. (b) Map of 𝑞(Ω,𝜉) where lower values are represented by darker colors. When Ω=1, 𝑞=0.50 for any value of 𝜉. (c) 𝜂(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (d) Map of 𝜂(Ω,𝜉).

(a)

(b)

(c)

(d)

Figure 5 

Plots of 𝑞 and 𝜂 for 𝑝=0.75. (a) 𝑞(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (b) Map of 𝑞(Ω,𝜉) where lower values are represented by darker colors. When Ω=1, 𝑞=0.75 for any value of 𝜉. (c) 𝜂(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (d) Map of 𝜂(Ω,𝜉).

(a)

(b)

(c)

(d)

Figure 6 

Plots of 𝑞 and 𝜂 for 𝑝=1, signifying the oscillator with viscous damping. (a) 𝑞(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (b) Map of 𝑞(Ω,𝜉) where lower values are represented by darker colors. When Ω=1, 𝑞=0.75 for any value of 𝜉. (c) 𝜂(Ω,𝜉) versus Ω for 𝜉=0.1,0.5,1,1.5, and 2 (going from light to dark). (d) Map of 𝜂(Ω,𝜉).

The expression for 𝑞 reveals regions in which the three terms of the original equation of motion are dominant. For example, regions in which the order of the derivative is near 2 (such as systems with lower damping ratio and higher forcing frequency) are primarily dominated by the inertial term. Lower values of damping ratio and orders of viscoelastic damping reach the asymptotic value more quickly, suggesting that the change in the dynamics (and also the phase shift) is more sensitive to changes in frequency in those cases. The dependence of the order of the derivative on the damping ratio also changes when Ω=1, corresponding to the case where the driving frequency is the same as the natural frequency of the system. When Ω<1, then the systems with higher damping ratios have a higher value for 𝑞. Once Ω>1 then the behavior is switched with the systems with higher damping ratio having lower values for the order. For all cases, when Ω=1, then 𝑞=𝑝 for any damping ratio.

A scaled behavior identical to 𝑞 is shown in the plots of 𝜂. The normalized amplitude ratio reaches an asymptotic value of 𝜂=1 as Ω increases. Thus, at higher frequencies the amplitude ratio is proportional to Ω^−𝑝, and the order of the damping and the damping ratio do not have any effect on the amplitude of the stationary motion. Similar to the case for 𝑞, 𝜂 approaches the asymptotic value more quickly for small viscoelastic damping orders 𝑝 and lower damping ratios. Also for the cases with damping order 𝑝<1, the peak amplitude response shifts to higher frequencies for increasing damping ratio. For increasing 𝑝 and damping ratio, the peak response begins to flatten out. The maps of 𝜂 show that as the order of the viscoelastic damping increases, the regions where the amplitude of the response is damped increase.

4. Conclusions

This work advances our understanding of the use of variable order (VO) differential operators in dynamics in two substantial ways. First, we compare several definitions proposed in literature, and select the most suitable definition based on a few criteria: (1) the VO operator must be able to return all intermediate values between 0 and 1 that correspond to the argument of the order of differentiation, (2) the VO operator must be effectively evaluated numerically, and (3) all derivatives of a true constant (a function that is constant from −∞ to +∞) must be zero. The operator defined in (2.15) satisfies these criteria for modeling dynamic systems. We then proceed to illustrate the meaning of a variable order of differentiation by analyzing a familiar dynamical problem (the stationary analysis of a viscoelastic oscillator). We determine that the order of differentiation for a single operator describing all dynamic elements in the stationary equation of motion (mass, damping and spring) is equal to the normalized phase shift. The normalization is easily understood as the quantity that transforms the maximum phase shift between acceleration and position (𝜋) into the maximum difference between the order of acceleration (2) and the order of position (0). The normalization constant is thus 2/𝜋, and the variable order differentiation is seen as being just the normalized phase shift in this problem, which gives us a straightforward interpretation of the meaning of variable orders of differentiation in dynamic systems.