#### 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 .

#### 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, , or , 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 () 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

When constant, then the -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]:

In the cases above, the lower terminal is set equal to 0, and it is assumed that for . 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 :

Further definitions are obtained by taking first-order derivatives of the integrals defined in (2.2)β(2.4)

where . Caputo-type VO operators are defined by taking the derivative of the function under the integrals

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:

where .

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 :

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, 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

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 such that for any , definitions (2.9) and (2.10) return the same result. However, if is a true constant, such that , 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:

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:

for and [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: , 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 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)**

**(a)**

**(b)**

**(c)**

As expected, the Caputo-based VO operators have values of 0 at , similar to the fractional derivative case. They do not return the zeroth-order derivative at 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 and is only a function of . The Riemann-Liouville and Caputo type operators with the argument , also are shown to be equivalent because . This is analogous to the similarity of the Riemann-Liouville and Caputo fractional derivatives of functions when for [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 , (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 (), the operator in (2.15) returns zero for all values of . However, if is not continuous between and , 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 as a representation of dynamic equilibrium. A nonzero initial condition is treated as a Heaviside function at , 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 where . When , the system is an oscillator with viscous damping, and for , the system is said to be characterized by viscoelastic damping. The equation is recast in dimensionless form using the following scaled parameters:

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

where is the dimensionless amplitude of the forcing, and is the damping ratio and is

Since we are concerned only with stationary behavior, and are not functions of time, so we look for a relationship such as The stationary VO derivative of is

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 yields

We now equate the real and imaginary parts of the above equation to arrive at

The stationary solution for (3.5) can be written as where the amplitude and phase shift are

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

Similarly, the expression for is

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 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 , 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)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

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 , corresponding to the case where the driving frequency is the same as the natural frequency of the system. When , then the systems with higher damping ratios have a higher value for . Once then the behavior is switched with the systems with higher damping ratio having lower values for the order. For all cases, when , 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 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 , 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: () 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, () the VO operator must be effectively evaluated numerically, and () 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 () and the order of position (). The normalization constant is thus , 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.