Variable-order fractional calculus: a change of perspective
aa r X i v : . [ m a t h . C A ] F e b VARIABLE-ORDER FRACTIONAL CALCULUS: A CHANGE OFPERSPECTIVE
ROBERTO GARRAPPA ∗ , ANDREA GIUSTI † , AND
FRANCESCO MAINARDI ‡ Abstract.
Several approaches to the formulation of a fractional theory of calculus of “variableorder” have appeared in the recent literature to describe physical systems showing special memoryeffects with features changing over time. Unfortunately, most of these proposals lack a rigorous math-ematical framework. Here, we consider an alternative view on the problem, originally proposed byG. Scarpi in the early seventies, based on a naive modification of the Laplace domain representationof standard kernels functions involved in (constant-order) fractional calculus. We discuss how thevariable-order Scarpi derivative can be framed within Kochubei’s General Fractional Calculus, thuspin pointing it is possible to frame variable-order Scarpi’s derivatives and integrals in the contextof a more general and robust mathematical theory. Then, taking advantage of powerful numeri-cal methods for the inversion of Laplace transforms, we discuss some practical applications of thevariable-order Scarpi integral and derivative to describe some kinds of transitions in the features ofthe system memory.
1. Introduction.
Derivatives and integrals of fractional (i.e., non-integer) orderare among the most fashionable tools for modeling phenomena featuring persistentmemory effects (i.e., non-localities in time). Since many physical systems are char-acterized by dynamics involving memory effects whose behaviour changes over time,even transitioning from a fractional order to another, the interest for fractional op-erators soon moved to their variable-order counterparts. Needless to say that thecompelling practical implications of these variable-order objects come at the price ofa more involved mathematical characterization.A naturally looking variable-order generalization of standard fractional derivativesis obtained by replacing the constant order α with a function α : [0 , T ] ⊂ R + → (0 , i.e. ,(1.1) I α ( t )0 f ( t ) = 1Γ (cid:0) α ( t ) (cid:1) Z t ( t − τ ) α ( t ) − f ( τ )d τ, possibly coupled to the RL-like variable-order derivative(1.2) D α ( t )0 f ( t ) = 1Γ (cid:0) − α ( t ) (cid:1) dd t Z t ( t − τ ) − α ( t ) f ( τ )d τ, where we restricted 0 < α ( t ) < D α ( t )0 does not necessarily act as the left-inverse for I α ( t )0 [40, 43]. This propertyis, nonetheless, recovered in the constant-order limit of the theory.Over the years, several proposals for fractional variable-order operators have ap-peared in the literature; for instance, we recall the works by Bohannan [2], Checkinet al. [3], Coimbra [5], Ingman and Suzdalnitsky [21], Kobelev et al. [23, 24], Lorenzo ∗ † Physics and Astronomy Department, Bishop’s University, 2600 College Street, QC J1M 1Z Sher-brooke, Canada ([email protected]), ‡ Physics and Astronomy Department, University of Bologna, Via Irnerio 46, 40126 Bologna, Italy([email protected]). 1
R. GARRAPPA, A. GIUSTI, AND R. MAINARDI and Hartley [28], Pedro et al. [36], Sierociuk et al. [47], Sun et al. [51]. For a com-prehensive review of the literature we refer the interested reader to Samko [39] andSun et al. [50]. However, despite showing seemingly useful for physical applications,these definitions face the conceptual mathematical problems discussed by Samko andRoss [40, 43].It is worth noting that, despite the complications discussed above, variable-ordermethods were employed by Zheng et al. [61, 62] in an attempt of overcoming theinconsistencies of some families of regular-kernel operators (in this regard, see [9] or[20] for a detailed discussion).From a numerical perspective, methods for solving variable-order fractional differ-ential equations (FDEs) have also been analysed to some extent; as a general referencewe mention here the works by Chen et al. [4], Tavares et al. [53], Zhuang et al. [63]and some of the several papers by Karniadakis et al. [58, 59, 60].Viscoelasticity is the perfect playground for variable-order operators, as certainknown scenarios display peculiar transitions from an order to another as a functiontime (see, e.g. , [38] or [6, 15, 17]). Further, in recent years variable-order fractionalcalculus has found some applications also in control theory, see for instance [1, 34].For a review of some of the latest applications of variable-order fractional operatorsin natural sciences we refer the interested reader to [35].To the best of our knowledge, the Italian engineer Giambattista Scarpi was thefirst to propose [44, 45, 46], in the early seventies, the use of time-fractional derivativeswith a time-dependent order. G. Scarpi’s work was inspired by an early model by Smit& de Vries [48] which was aimed at providing a theoretical framework for materialsshowing features intermediate between solids and liquids. Notably, the approachproposed by G. Scarpi was not based on a naive replacement, in the kernel of somefractional derivative, of the constant order α with a variable-order function α ( t ). Theprocedure proposed by Scarpi, instead, acts in a more subtle way at the level of theLaplace transform (LT) domain (on a different basis, however, with respect to theoperator proposed by Coimbra in [5]) and constitutes an interesting novelty withrespect to more traditional approaches.Despite the boom that the active research on fractional calculus has been ex-periencing for the last decade, so far Scarpi’s approach has been mostly overlooked(except for a very recent contribution by Cuesta and Kirane [7] of which we are awarebecause of a private communication).If, on the one hand, Scarpi’s works were the first to introduce this peculiarapproach to variable-order theories, on the other hand, they are solely focused onphysical properties and implications of the proposed methods. In other words, themathematical foundations supporting these object were not analyzed in details. Ad-ditionally, the operators proposed by Scarpi require reliable numerical techniques forhandling the inversion of the LT, which were not available at the time of publicationof Scarpi’s seminal works.Kochubei’s general fractional calculus [25, 27, 26] (see also [32]) offers a self-consistent approach to the generalization of fractional calculus to an enhanced classof kernels. The key aspect of Kochubei’s approach consists in the fact that it relayson a precise characterization of these generalized operators in the Laplace domain.This makes general fractional calculus the perfect tool for providing a precise androbust mathematical framework for Scarpi’s theory. On the numerical side, the severaladvancements in the field of the numerical inversion of the LT, among which we recallthe contribution by Weidemann and Trefethen [57], provide us with the machineryneeded to implement Scarpi’s ideas to their fullest. ARIABLE-ORDER FRACTIONAL CALCULUS α ( t ), must satisfy. In Section 4 we consider some in-structive examples for α ( t ), describing the transition of the fractional order for thecorresponding kernels in the, including transition from elasticity to viscosity. Section5 is devoted the solution of the relaxation equation with the Scarpi derivative and, fi-nally, in Section 6 we provide some concluding remarks. Note that the method used toinvert numerically the LT, allowing the investigation of Scarpi’s fractional operators,is discussed in Appendix A.
2. Scarpi’s variable-order fractional calculus.
In order to introduce, andfurther develop, Scarpi’s ideas on variable-order derivatives we preliminary recall somebackground materials on fractional integrals and derivatives.
In this work we consider functions which are absolutelycontinuous on some interval [0 , T ], i.e. f ∈ AC [0 , T ]. This is a not particularlyrestrictive assumption and, basically, it means that f ( t ) = f (0) + Z t g ( s )d s , where 0 ≤ t ≤ T , for some g ∈ L [0 , T ], with L [0 , T ] the usual space of Lebesgue-integrable functions on [0 , T ].The standard Dzhrbashyan-Caputo notion of fractional derivative of order 0 <α <
1, commonly referred to simply as Caputo derivative, is defined in terms of theweakly-singular Volterra-type integro-differential operator given by(2.1) C D α f ( t ) = Z t φ ( t − τ ) f ′ ( τ )d τ, φ ( t ) = t − α Γ(1 − α ) . The defining property of C D α is that it acts as the left-inverse of the RL integral(2.2) RL I α f ( t ) = Z t ψ ( t − τ ) f ( τ )d τ, ψ ( t ) = 1Γ( α ) t α − , see e.g. [8, 22]. In other words, one has that C D α RL I α f ( t ) = f ( t ) and RL I α C D α f ( t ) = f ( t ) − f (0) , thus implementing a sort of fundamental theorem of fractional calculus [31]. Basically,the Caputo derivative C D α was introduced to provide a regularization of the RL one(2.3) RL D α f ( t ) = dd t Z t φ ( t − τ ) f ( τ )d τ, φ ( t ) = t − α Γ(1 − α ) . Indeed, the Caputo derivative allows to write fractional differential equations (FDEs)of order 0 < α < i.e. , The function ψ ( t ) is known as the Gel’fand-Shilov kernel [14, 19, 33]. R. GARRAPPA, A. GIUSTI, AND R. MAINARDI involving just integer order derivatives. Since these initial value problem have amore straightforward physical interpretation, in this work we focus on regularizedCaputo-like derivatives, though this is done without loss of generality since recastingthe arguments presented here in the RL framework does not involve any particularcomplication.Before moving on with the analysis of the Scarpi derivative it is important torecall some important properties of standard fractional operators. Specifically, werecall that taking advantage of the well-known LTs(2.4) Φ( s ) := L (cid:16) φ ( t ) ; s (cid:17) = s α − , Ψ( s ) := L (cid:16) ψ ( t ) ; s (cid:17) = 1 s α one finds that, assuming that the function f ( t ) admits the LT F ( s ), the LTs of (2.1)and (2.2) yield(2.5) L (cid:16) C D α f ( t ) ; s (cid:17) = s α F ( s ) − s α − f (0) , L (cid:16) RL I α f ( t ) ; s (cid:17) = 1 s α F ( s ) . In order to provide a variable-order generalization of (2.1) we consider a function α ( t ) : [0 , T ] → (0 , , T ]. The restriction on the image of α ( t ) to(0 ,
1) is done to avoid further technical complications.The main idea by Scarpi presented in the pioneering works [44, 45, 46] was todefine a fractional derivative of variable order α ( t ) by generalizing its representationof (2.1) in the Laplace domain rather than in the time domain.If one considers the constant function α ( t ) ≡ α , t >
0, its LT is A ( s ) = α/s andhence one can trivially infer that Φ( s ) and Ψ( s ) in (2.4), can be recast in terms of A ( s ) as Φ( s ) = s sA ( s ) − Ψ( s ) = s − sA ( s ) . Thus, Scarpi’s idea consists in extending this simple argument to any non-constantlocally integrable function α ( t ) with LT A ( s ) := L (cid:16) α ( t ) ; s (cid:17) = Z t e − st α ( t )d t, and define a variable-order derivative by means of a convolution similar to those in(2.1) and (2.3). Definition
Let α : [0 , T ] → (0 , be a locally integrable function, let A ( s ) denote the LT of α ( t ) , and let f ∈ AC [0 , T ] . The Scarpi fractional derivative S D α ( t )0 of variable order α ( t ) is defined as (2.6) S D α ( t )0 f ( t ) := Z t φ α ( t − τ ) f ′ ( τ )d τ, t ∈ [0 , T ] , where the kernel function φ α ( t ) is the inverse LT of Φ α ( s ) = s sA ( s ) − . Clearly, the Scarpi derivative reduces to the standard Caputo one when α ( t )becomes constant. ARIABLE-ORDER FRACTIONAL CALCULUS L (cid:16) S D α ( t )0 f ( t ) ; s (cid:17) = s sA ( s ) F ( s ) − s sA ( s ) − f (0) . However, finding an explicit representation of the kernel φ α ( t ) is not always pos-sible and in the following sections we will explore some computational approaches tothis problem. Remark
It is of interest,especially for solving differential equations, to find an integral operator S I α ( t )0 of convo-lution type, with some kernel ψ α ( t ), such that the fundamental theorem of fractionalcalculus holds also for the Scarpi derivative, namely(2.8) S D α ( t )0 S I α ( t )0 f ( t ) = f ( t ) , S I α ( t )0 S D α ( t )0 f ( t ) = f ( t ) − f (0) . For this to be true the two kernels φ α ( t ) and ψ α ( t ) must satisfy a Sonine equation[25, 41, 42](2.9) Z t φ α ( t − τ ) ψ α ( τ ) = 1 , t > , and, in this case, φ α ( t ) and ψ α ( t ) are said to form a Sonine pair. Sonine pairs havebeen extensively studied in the literature, see e.g. , [42, 41].Given a generic function φ α ( t ), finding the corresponding ψ α ( t ) such that thetwo functions form a Sonine pair is not trivial. However, this problem simplifiessubstantially when working in the Laplace domain. Proposition
Let α : [0 , T ] → (0 , be a locally integrable function, let A ( s ) denote the LT of α ( t ) , and let f ∈ L [0 , T ] . The integral operator (2.10) S I α ( t )0 f ( t ) = Z t ψ α ( t − τ ) f ( τ )d τ, where ψ α ( t ) is such that its LT yields Ψ α ( s ) = s − sA ( s ) , admits S D α ( t )0 as its left-inverse . In other words, the couple { S D α ( t )0 , S I α ( t )0 } satisfy the conditions in (2.8).Proof. If Ψ α ( s ) = s − sA ( s ) , then φ α ( t ) and ψ α ( t ) form a Sonine pair. Indeed, theSonine equation (2.9) in the Laplace domain reads(2.11) Φ α ( s )Ψ α ( s ) = 1 s which is trivially satisfied because of the definition of the Scarpi derivative that re-quires Φ α ( s ) = s sA ( s ) − .
3. Some necessary assumptions.
Not any transition function α ( t ) with val-ues in (0 ,
1) is a suitable candidate to generate a pair of Scarpi-type variable-orderfractional operators. In other words, not all α ( t ) are such that { S D α ( t )0 , S I α ( t )0 } satisfythe fundamental theorem of calculus (2.8). R. GARRAPPA, A. GIUSTI, AND R. MAINARDI
Hanyga in [20] has shown that to ensure that two functions φ ( t ) and ψ ( t ) forma Sonine pair (without moving to the realm of distributions) they must have anintegrable singularity at the origin. This is further supported by the analysis in[9] which showed that operators based on regular kernels satisfy the fundamentaltheorem of calculus (2.8) when they are restricted to spaces of functions with severe(and somewhat artificial) constraints (see also [49]).In the context of Scarpi’s theory it is convenient to follow the approach proposedby Kochubei in [25, 27], and recently reviewed by Luchko and Yamamoto in [32], inorder to embed Scarpi’s ideas in general fractional calculus. This approach indeedrelies on the LT of the kernel of a general fractional derivative, which is actually themain information available for Scarpi’s operators.The theory of general fractional calculus formulated in [27, 32, 25], demands thatthe LT Φ α ( s ) of the derivative kernel φ α ( t ) satisfies the following three conditions C1 : Φ α ( s ) is a Stieltjes function, C2 : Φ α ( s ) → s Φ α ( s ) = s sA ( s ) → ∞ as s → ∞ , C3 : Φ α ( s ) → ∞ and s Φ α ( s ) = s sA ( s ) → s → α ( s ) is a Stieltjes function if it admits the integralrepresentation Φ α ( s ) = Z ∞ K Φ α ( r ) s + r d r , (3.1)with K Φ α ( r ) ≥ C2 and C3 ensure the existence of a function Ψ α ( s ) so that it inverse LT ψ α ( t )forms a Sonine pair with φ α ( t ). In other words, given the Scarpi derivative these con-ditions guarantee the existence of S I α ( t )0 so that the fundamental theorem of calculus(2.8) is satisfied.The condition C1 is instead much stronger. It ensures that φ α ( t ) and ψ α ( t ) arecompletely monotone functions. A consequence of this property is that the solutionof fractional relaxation equations involving these operators are given by completelymonotone functions. This property plays a crucial role in linear viscoelasticity, see e.g. , [33]. However, when considering transitions from an order α to an order α it is generally no longer possible to ensure a completely monotone character for thesolution of the relaxation equation. This means, as discussed in the following Sectionusing some examples, that even very simple transition functions α ( t ) yield a derivativekernel φ α ( t ) violating C1 . Hence, in order to establish the conditions on α ( t ) we willrelax the assumptions of Kochubei’s generalized fractional calculus requiring only thevalidity of C2 and C3 .
4. Physically relevant examples of transition functions.
In this Sectionwe present some examples of variable-order functions α ( t ). We confine to potentiallyphysically interesting scenarios where α ( t ) shows a monotone transition from an initialorder α to a final order α , where the latter is only reached asymptotically as t → ∞ .Various expressions for α ( t ) are presented here and for each of them we show theemerging kernels φ α ( t ) and ψ α ( t ) associated to the corresponding S D α ( t )0 and S I α ( t )0 defined respectively in (2.6) and (2.10).Even when α ( t ) and its LT A ( s ) are given by simple expressions, in general itis not possible to provide an explicit representation of the kernels φ α ( t ) and ψ α ( t ).Therefore, they need to be evaluated numerically by means of the inversion of theirLTs Φ α ( s ) and Ψ α ( s ). On the one hand, this complication constituted the main ARIABLE-ORDER FRACTIONAL CALCULUS
For 0 ≤ α < α ≤ c >
0, we consider the function α ( t ) = α + ( α − α )e − ct describing a variable-order transition from α to α according to an exponential lawwith rate − c . It is simple to evaluate the LT of α ( t ) as A ( s ) = Z ∞ e − st α ( t )d t = α c + α ss ( c + s )and, hence,Φ α ( s ) = s sA ( s ) − = s ( α − c +( α − sc + s , Ψ α ( s ) = s − sA ( s ) = s − α c + α sc + s . It is also easy to verify the conditions C2 and C3 , indeed one has thatlim s →∞ Φ α ( s ) = lim s →∞ s α − = 0 and lim s →∞ s Φ α ( s ) = lim s →∞ s α = + ∞ and lim s → Φ α ( s ) = lim s → s α − = ∞ and lim s → s Φ α ( s ) = lim s → s α = 0 . Furthermore, the spectral distribution K Φ α ( r ) that yields the integral representation(3.1) of Φ α ( s ) can be evaluated by means of the Titchmars [54] inversion formula K Φ α ( r ) = ∓ π Im h Φ α ( s ) (cid:12)(cid:12) s = r e ± i π i that implies K Φ α ( r ) = − π r ( α − c − ( α − rc − r sin h ( α − c − ( α − rc − r π i . Note that K Φ α ( r ) ≥ r ≥ c = 0 or α = α , namely when α ( t ) assumes constant values α or α . This in turn jeopardizes the completely mono-tonic character of φ α ( t ) (and of the solution of the associated relaxation equation) inthe presence of a such a time-dependent order. Transition from α = 0 to α = 1 . By setting for simplicity c = 1, we obtainthe transition function α ( t ) = 1 − e − t , shown in the left panel of Figure 1, and for which one has that A ( s ) = 1 s ( s + 1) , Φ α ( s ) = s − ss +1 , Ψ α ( s ) = s − s +1 . The effects produced by the corresponding kernel φ α ( t ), presented in the rightpanel of Figure 1, in the relaxation equation with the derivative obtained with thiskernel will be better illustrated in the following section. R. GARRAPPA, A. GIUSTI, AND R. MAINARDI
Fig. 1 . Plot of the transition function α ( t ) = 1 − e − t (left plot) and of the corresponding kernel φ α ( t ) (right plot). Generic order transition from α to α . More generally, is more realistic toconsider transitions from two values α and α in (0 , α = 0 . α = 0 .
8. Although it is reasonable to assume, for some physicalmodels, that α and α are close values, we shall consider distant enough values forthese parameters in order to be able to graphically present the asymptotic behaviourof φ α ( t ) and ψ α ( t ) in a nice way.As one can see from Figure 3, the resulting kernels φ α ( t ) and ψ α ( t ) start as thecorresponding kernels of the standard fractional operators of order α and asymptot-ically converge to the kernels of the operators of order α . This behaviour can bebetter appreciated in Figure 4 where φ α ( t ) and ψ α ( t ) are plotted in logarithmic scale. Fig. 2 . Plot of α ( t ) for variable-order transition of exponential type ( c = 2 . ) from α = 0 . to α = 0 . . Fig. 3 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for variable-order transition ofexponential type ( c = 2 . ) from α = 0 . to α = 0 . . ARIABLE-ORDER FRACTIONAL CALCULUS -3 -2 -1 -2 -1 -3 -2 -1 Fig. 4 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for variable-order transition ofexponential type ( c = 2 . ) from α = 0 . to α = 0 . (logarithmic scale). The previousexample can be generalized by replacing the exponential with the Mittag-Leffler (ML)function, i.e. , α ( t ) = α + ( α − α ) E β ( − ct β ) , where E β ( z ) = ∞ X k =0 z k Γ( αk + β )is the one parameter ML function (see, for instance [18]). This procedure gives abetter control on the transition from α to α thanks to the additional parameter β .The representation of this variable-order function α ( t ) is provided in the left panelof Figure 5 for β = 0 . c = 2 .
0. Clearly, the transition presented in Section 4.1 isjust a particular case of the one presented here since e x = E ( x ).It is now fairly easy to compute the LT of α ( t ), that reads [18] A ( s ) = Z ∞ e − st α ( t )d t = α c + α s β s ( c + s β )and, hence,Φ α ( s ) = s sA ( s ) − = s ( α − c +( α − sβc + sβ , Ψ α ( s ) = s − sA ( s ) = s − α c + α sβc + sβ . Conditions C2 and C3 can then be easily checked as discussed for the previouscase. Furthermore, also in this case the corresponding kernels φ α ( t ) and ψ α ( t ) matchthe kernels of the standard fractional operators of order α and α in the limitig casesof the model, as shown in Figures 6 and 7.Note that the parameter β , similarly to the parameter c in the previous case,alters the way in which this transition happens without affecting the initial and finalvalues of the order.In Figure 8 we compare the behaviour of α ( t ) and ψ α ( t ) for the decay of ML-type as we vary the parameter β . Observe that the case β = 1 . erf type. Consider now for 0 < α <α < c > α ( t ) = α + ( α − α ) erf( √ ct )0 R. GARRAPPA, A. GIUSTI, AND R. MAINARDI
Fig. 5 . Plot of α ( t ) for order transition of ML type ( c = 2 . and β = 0 . ) from α = 0 . to α = 0 . . Fig. 6 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for order transition of ML type( c = 2 . and β = 0 . ) from α = 0 . to α = 0 . . -4 -2 -2 -4 -2 Fig. 7 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for order transition of ML type( c = 2 . and β = 0 . ) from α = 0 . to α = 0 . (logarithmic scale). Fig. 8 . Plot of α ( t ) (left plot) and ψ α ( t ) (right plot) for the variable order with decay of MLtype with α = 0 . , α = 0 . , c = 2 . and different values of β . ARIABLE-ORDER FRACTIONAL CALCULUS α to α as shown in Fig-ure 9. Observe that this function can be considered, in some sense, as a furthergeneralization of the variable-order function α ( t ) based on the ML function sinceerf( √ ct ) = √ ct E , ( − tc )with E γα,β ( z ) the three-parameter ML function, also known as Prabhakar function(see, e.g. , [10, 16, 18, 37]).The Laplace transform of α ( t ) is A ( s ) = α s + ( α − α ) √ cs √ s + c = α √ c + α (cid:0) √ s + c − √ c (cid:1) s √ s + c and the corresponding function φ α ( t ) and ψ α ( t ) are depicted in Figures 10 and 11. Fig. 9 . Plot of α ( t ) for order transition of erf type ( c = 2 . ) from α = 0 . to α = 0 . . Fig. 10 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for order transition of erf type( c = 2 . ) from α = 0 . to α = 0 . .
5. Fractional relaxation equation with Scarpi derivative.
The aim of thisSection is to provide a preliminary investigation of the variable-order fractional relax-ation equation(5.1) (cid:26) S D α ( t )0 y ( t ) = − λy ( t ) y (0) = y , where S D α ( t )0 is the Scarpi variable-order fractional derivative introduced in Definition2.1 and λ > R. GARRAPPA, A. GIUSTI, AND R. MAINARDI -4 -2 -2 -4 -2 Fig. 11 . Plot of kernels φ α ( t ) (left plot) and ψ α ( t ) (right plot) for order transition of erf type( c = 2 . ) from α = 0 . to α = 0 . (logarithmic scale). Finding analytical solutions for the initial value problem (5.1) does not seem ingeneral possible since the absence of an explicit representation of the kernel φ α F T t ( t )of S D α ( t )0 . Therefore, tackling this problem from a numerical perspective becomesunavoidable and necessary.Since the linear nature of (5.1), a simple approach consists in exploiting the LTand its numerical inversion. Indeed, by applying the LT to both sides of (5.1), andrecalling Eq. (2.7), one finds s sA ( s ) Y ( s ) − s sA ( s ) − y = − λY ( s ) , with Y ( s ) the LT of the solution y ( t ). Therefore, an algebraic manipulation leads to(5.2) Y ( s ) = y s (cid:0) λ Ψ α ( s ) (cid:1) and hence it is possible to evaluate the solution y ( t ) = L − (cid:0) Y ( s ) ; t (cid:1) in the timedomain by applying again one of the methods for the numerical inversion of the LTas the one described in the Appendix A.To this end we first discuss the solution of the relaxation equation (5.1) for α ( t ) =1 − e − t . Fig. 12 . Plot of the solution y ( t ) of the relaxation equation (5.1), with λ = 1 and y =1 , for variable-order transition α ( t ) = 1 − e − t and comparison with solutions y ( t ) and y ( t ) ofcorresponding standard Hooke’s and Newton’s laws. ARIABLE-ORDER FRACTIONAL CALCULUS y ( t ) obtained by numerical inversion ofthe LT (5.2) instantaneously achieves the value y ( t ) ≡ y / (1 + λ ), and, as expected,asymptotically converges to the solution y ( t ) of the standard relaxation equation.This first result provides an excellent proof of the ability of the Scarpi derivative toproperly describe a transition between scenarios described by operators of differentorder.In the following we present the solutions y ( t ) of the relaxation equation (5.1)with the other transition functions α ( t ) introduced in Section 4. In the various plots,together with the solution y ( t ), we also offer a comparison of y ( t ) with the solutions y ( t ) and y ( t ) of the same relaxation equation with the standard Caputo derivativeof order α and α , respectively.In the first case, see Figure 13, the exponential transition α ( t ) = α +( α − α )e − ct (with α = 0 . α = 0 . c = 2) is considered. Fig. 13 . Plot of the solution y ( t ) of the relaxation equation (5.1), with λ = 1 and y = 1 ,for variable-order transition α ( t ) = α + ( α − α )e − ct , with α = 0 . , α = 0 . and c = 2 , andcomparison with solutions y ( t ) and y ( t ) of the standard fractional relaxation equations of order α and α . The numerical results show how well the solution with the Scarpi derivativematches the solution of the Caputo relaxation equation of order α close to the originand of the Caputo relaxation equation of order α for large t . The box in each figureoffers a closer look of the solutions near to the origin.Similar results are obtained with the transition function α ( t ) = α + ( α − α ) E β ( − ct β ) (with α = 0 . α = 0 . c = 2 . β = 0 .
7) shown in Figure 14, aswell as with the transition function α ( t ) = α + ( α − α ) erf( √ ct ) (with α = 0 . α = 0 . c = 2 .
0) depicted in Figure 15.Alternatively, one can solve the initial value problem in (5.1) by using the integralformulation of the problem(5.3) y ( t ) = y − λ S I α ( t )0 y ( t ) = y − λ Z t ψ α ( t − τ ) y ( τ )d τ, and then apply the convolution quadrature rules devised and studied by Lubich in[29, 30]. These rules have the great advantage of providing accurate approximationsof convolution integrals like the one in (5.3) for which the kernel φ ( t ) is known onlythrough its LT Φ( s ), as it is for the Scarpi integral. Hence, this scheme looks rather4 R. GARRAPPA, A. GIUSTI, AND R. MAINARDI
Fig. 14 . Plot of the solution y ( t ) of the relaxation equation (5.1), with λ = 1 and y = 1 ,for variable-order transition α ( t ) = α + ( α − α ) E β ( − ct β ) , with α = 0 . , α = 0 . , c = 2 . and β = 0 . , and comparison with solutions y ( t ) and y ( t ) of the standard fractional relaxationequations of order α and α . Fig. 15 . Plot of the solution y ( t ) of the relaxation equation (5.1), with λ = 1 and y = 1 , forvariable-order transition α ( t ) = α + ( α − α ) erf( √ ct ) , with α = 0 . , α = 0 . , c = 2 . , andcomparison with solutions y ( t ) and y ( t ) of the standard fractional relaxation equations of order α and α . promising for handling general fractional differential equations, in special way of non-linear type, involving the Scarpi derivative .
6. Concluding remarks.
This paper aims at making the first step toward re-viving Scarpi’s ides on variable-order fractional calculus. We have framed these ideasin terms of Kochubei’s general fractional calculus and we have shown one of the pos-sible numerical approaches needed for handling these derivatives and related initialvalue problems.There are still many open questions which require some care. One of these,for instance, concerns conditions under which a law of exponents S I α ( t )0 S I β ( t )0 f ( t ) = S I α ( t )+ β ( t )0 f ( t ) might hold for these operators. Another example consists in designingnumerical methods to solve more general FDEs involving the Scarpi derivative. Fur-thermore, the potential for applications of this scheme in physics, engineering, andother natural sciences is unbounded, opening a path to a completely uncharted new ARIABLE-ORDER FRACTIONAL CALCULUS
Appendix A. A numerical method for the inversion of the LT.
In this Appendix we provide a detailed description of the method employed inthe previous Sections to numerically invert the LT of the kernels φ α ( t ) and ψ α ( t ) andof the solution of the relaxation equation (5.1).The method is based on the main idea by Talbot [52] consisting in deforming, inthe formula for the inversion of the LT F ( s ) of a function f ( t )(A.1) f ( t ) = 12 π i Z σ +i ∞ σ − i ∞ e st F ( s )d s, the Bromwich line ( σ − i ∞ , σ + i ∞ ) into a different contour C beginning and end-ing in the left complex half-plane. In this way it is possible to obtain an accurateapproximation of the function f ( t ) after applying a suitably chosen quadrature rulealong C , since the strong oscillations of the exponential, and the resulting numericalinstability, are avoided.This approach was successively refined by Weidemann and Trefethen [57] whoprovided a detailed error analysis allowing to properly select the geometry of thecontour C and the quadrature parameters in order to achieve any prescribed accuracy ε > F ( s ) withone or more singularities scattered in the complex plane. However. since in ourexamples we are faced with LTs F ( s ) having just singularities at the origin or on thebranch-cut, the original algorithm introduced in [57] turns out to be good enough.One of the most useful contours used to replace the Bromwich line in (A.1) is aparabolic-shaped contour described by the equation z ( u ) = µ (i u + 1) , −∞ < u < ∞ , where µ > z ( u ), suitably cho-sen to encompass any possible singularity of F ( s ), one obtains the equivalent formu-lation(A.2) f ( t ) = 12 π i Z + ∞−∞ e z ( u ) t F ( z ( u )) z ′ ( u )d u. Hence, the application of a trapezoidal rule with step-size h on a sufficiently largetruncated interval [ − hN, hN ] leads to the approximation(A.3) f h,N ( t ) = h π i N X k = − N e z ( u k ) t F ( z ( u k )) z ′ ( u k ) , u k = hk. The choice of the three parameters µ , h and N is essential to achieve an accurateapproximation of f ( t ) and it is driven by a detailed analysis of the error | f ( t ) − f h,N ( t ) | .6 R. GARRAPPA, A. GIUSTI, AND R. MAINARDI
This in turn consists of two main components: the discretization error (DE) and thetruncation error (TE). By following the analysis in [57], in absence of singularities of F ( s ) (except for the branch-point singularity at the origin and the branch-cut placed,for convenience, on the negative real semi-axis) one can find that | DE | = O (cid:16) e − π/h (cid:17) + O (cid:16) e − π / ( µth )+2 π/h (cid:17) , h → , | T E | = O (cid:16) e µt (cid:0) − ( hN ) (cid:1)(cid:17) , h →
0A more accurate analysis takes into account the round-off error (RE) as well, forwhich (after exploiting | z ′ ( u ) | = 2 √ µ p | z ( u )) | ) the following estimates hold [56] | RE | ≤ ǫhπ e µt N X k =0 | F ( z ( u k )) || z ′ ( u k ) | = 2 ǫ √ µhπ e µt N X k =0 | ˆ F ( z ( u k )) |≈ ǫ e µt √ µπ Z Nh | ˆ F ( s ) | d s, where ǫ is the precision machine and ˆ F ( s ) = F ( s ) s . Obviously, the analysis needsto be customized according to the specific LT F ( s ) which must be inverted. If ˆ F ( s ) isassumed to have a moderate growth and N h is in general not large (in practice veryoften it is h = O (cid:0) N − (cid:1) one can neglect the integral in the estimate of RE (as well asthe 2 √ µ/π term) and just assume | RE | ≈ ǫ e µt .Optimal parameters µ , h and N can be now obtained after balancing the threedifferent errors and imposing that they are proportional to a given prescribed accuracywhich, to simplify the analysis and at the same time ensure accurate results, we selectat the same level of the precision machine ǫ ≈ . × − . Therefore, after imposingthat | DE | ≈ | T E | ≈ | RE | ≈ ε asymptotically as h →
0, and denoting L = − log ǫ ,the balancing of the three errors leads to N = 4 L π , µ = L tπ N , h = 2 πL + L πN . Remark
A.1. In the above analysis we have assumed a moderate growth of ˆ F ( s )as s → ∞ . With respect to the transition α ( t ) considered in our examples thisassumption is truly reasonable in order to compute ψ α ( t ) or the solution y ( t ) of therelaxation equation (5.1) but could be too much optimistic for the evaluation of φ α ( t )which is expected to have a more sustained growth. Although we have obtainedreasonable results the same, we think that a more detailed analysis is necessary if oneaims to compute φ α ( t ) with high accuracy.In the following we report the few lines of a Matlab code for the numerical in-version of the LT F ( s ) on a vector of points t . The code is optimized to evaluatejust functions f ( t ) with real values. The LT F ( s ) is assumed not to have singularitiesexcept a possible one at the origin. L = - log ( eps ) ;N = ceil (4* L /3/ pi ) ;h = 2* pi / L + L /2/ pi / N ^2 ;p = L ^3/4/ pi ^2/ N ^2 ;u = (0: N ) * h ;f = zeros ( size ( t ) ) ;
ARIABLE-ORDER FRACTIONAL CALCULUS for n = 1 : length ( t )mu = p / t ( n ) ;z = mu *( u *1 i +1) .^2 ; z1 = 2* mu *(1i - u ) ;G = exp ( z * t ( n ) ) .* F ( z ) .* z1 ;f ( n ) = ( imag ( G (1) ) /2+ sum ( imag ( G (2: N +1) ) ) ) * h / pi ;end Acknowledgments.
The work of R.Garrappa is supported by INdAM undera GNCS-Project 2020. The work of A.Giusti is supported by the Natural Sciencesand Engineering Research Council of Canada (Grant No. 2016-03803 to V. Faraoni)and by Bishop’s University. The work of A.Giusti and F.Mainardi has been carriedout in the framework of the activities of the Italian National Group for MathematicalPhysics [Gruppo Nazionale per la Fisica Matematica (GNFM), Istituto Nazionale diAlta Matematica (INdAM)].
REFERENCES[1]
G. M. Bahaa , Fractional optimal control problem for variable-order differential systems , Fract.Calc. Appl. Anal., 20 (2017), pp. 1447–1470.[2]
G. W. Bohannan , Comments on time-varying fractional order , Nonlinear Dyn., 90 (2017),p. 2137–2143.[3]
A. Chechkin, R. Gorenflo, and I. Sokolov , Fractional diffusion in inhomogeneous media ,Journal of Physics A: Mathematical and General, 38 (2005), pp. L679–L684.[4]
S. Chen, F. Liu, and K. Burrage , Numerical simulation of a new two-dimensional variable-order fractional percolation equation in non-homogeneous porous media , Comput. Math.Appl., 68 (2014), pp. 2133–2141.[5]
C. Coimbra , Mechanics with variable-order differential operators , Annalen der Physik, 12(2003), pp. 692–703.[6]
I. Colombaro, A. Giusti, and F. Mainardi , A class of linear viscoelastic models based onbessel functions , Meccanica, 52 (2017), pp. 825–832.[7]
E. Cuesta and M. Kirane , On the sub–diffusion fractional initial value problem with timevarying order . submitted, 2020.[8]
K. Diethelm , The Analysis of Fractional Differential Equations , vol. 2004 of Lecture Notes inMathematics, Springer-Verlag, Berlin, 2010.[9]
K. Diethlem, R. Garrappa, A. Giusti, and M. Stynes , Why fractional derivatives withnonsingular kernels should not be used , Fract. Calc. Appl. Anal., 23 (2020), pp. 610–634.[10]
R. Garra and R. Garrappa , The Prabhakar or three parameter Mittag-Leffler function:theory and application , Commun. Nonlinear Sci. Numer. Simul., 56 (2018), pp. 314–329.[11]
R. Garrappa , Numerical evaluation of two and three parameter Mittag-Leffler functions ,SIAM J. Numer. Anal., 53 (2015), pp. 1350–1369.[12]
R. Garrappa and M. Popolizio , Evaluation of generalized Mittag–Leffler functions on thereal line , Adv. Comput. Math., 39 (2013), pp. 205–225.[13]
R. Garrappa and M. Popolizio , Computing the matrix Mittag-Leffler function with applica-tions to fractional calculus , J. Sci. Comput., 77 (2018), pp. 129–153.[14]
I. M. Gel’fand and G. E. Shilov , Generalized Functions. Vol. 1 , AMS Chelsea Publishing,Providence, RI, 2016.[15]
A. Giusti , On infinite order differential operators in fractional viscoelasticity , Frac. Calc. App.Anal., 20 (2017), pp. 854–867.[16]
A. Giusti, I. Colombaro, R. Garra, R. Garrappa, F. Polito, M. Popolizio, andF. Mainardi , A practical guide to Prabhakar fractional calculus , Fract. Calc. Appl. Anal.,23 (2020), pp. 9–54.[17]
A. Giusti and F. Mainardi , A dynamic viscoelastic analogy for fluid-filled elastic tubes , Mec-canica, 51 (2016), pp. 2321–2330.[18]
R. Gorenflo, A. A. Kilbas, F. Mainardi, and S. V. Rogosin , Mittag-Leffler Functions,Related Topics and Applications , Springer Monographs in Mathematics, Springer-VerlagBerlin Heidelberg, 2020.[19]
R. Gorenflo and F. Mainardi , Fractional calculus: integral and differential equations offractional order , in Fractals and fractional calculus in continuum mechanics (Udine, 1996), R. GARRAPPA, A. GIUSTI, AND R. MAINARDIvol. 378 of CISM Courses and Lect., Springer, Vienna, 1997, pp. 223–276. E-printhttp://arxiv.org/abs/0805.3823.[20]
A. Hanyga , A comment on a controversial issue: a generalized fractional derivative cannothave a regular kernel , Fract. Calc. Appl. Anal., 23 (2020), pp. 211–223.[21]
D. Ingman and J. Suzdalnitsky , Control of damping oscillations by fractional differentialoperator with time-dependent order , Comput. Methods Appl. Mech. Engrg., 193 (2004),pp. 5585–5595.[22]
A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo , Theory and Applications of FractionalDifferential Equations , vol. 204 of North-Holland Mathematics Studies, Elsevier ScienceB.V., Amsterdam, 2006.[23]
Y. L. Kobelev, L. Y. Kobelev, and Y. L. Klimontovich , Anomalous diffusion with memorythat depends on time and coordinates , Dokl. Akad. Nauk, 390 (2003), pp. 605–609.[24]
Y. L. Kobelev, L. Y. Kobelev, and Y. L. Klimontovich , Statistical physics of dynamicalsystems with variable memory , Dokl. Akad. Nauk, 390 (2003), pp. 758–762.[25]
A. N. Kochubei , General fractional calculus, evolution equations, and renewal processes , In-tegral Equations Operator Theory, 71 (2011), pp. 583–600.[26]
A. N. Kochubei , Equations with general fractional time derivatives—Cauchy problem , inHandbook of fractional calculus with applications. Vol. 2, De Gruyter, Berlin, 2019,pp. 223–234.[27]
A. N. Kochubei , General fractional calculus , in Handbook of fractional calculus with applica-tions. Vol. 1, De Gruyter, Berlin, 2019, pp. 111–126.[28]
C. F. Lorenzo and T. T. Hartley , Variable order and distributed order fractional operators ,Nonlinear Dynam., 29 (2002), pp. 57–98.[29]
C. Lubich , Convolution quadrature and discretized operational calculus. I , Numer. Math., 52(1988), pp. 129–145.[30]
C. Lubich , Convolution quadrature and discretized operational calculus. II , Numer. Math., 52(1988), pp. 413–425.[31]
Y. Luchko , Fractional derivatives and the fundamental theorem of fractional calculus , Fract.Calc. Appl. Anal., 23 (2020), pp. 939–966.[32]
Y. Luchko and M. Yamamoto , The general fractional derivative and related fractional dif-ferential equations , Mathematics, 8 (2020), p. 2115.[33]
F. Mainardi , Fractional Calculus and Waves in Linear Viscoelasticity , Imperial College Press,London, 2010. An introduction to mathematical models.[34]
P. W. Ostalczyk, P. Duch, D. W. Brzezi´nski, and D. Sankowski , Order functions selectionin the variable-, fractional-order pid controller , in Advances in Modelling and Control ofNon-integer-Order Systems, K. J. Latawiec, M. Lukaniszyn, and R. Stanis lawski, eds.,Cham, 2015, Springer International Publishing, pp. 159–170.[35]
S. Patnaik, J. P. Hollkamp, and F. Semperlotti , Applications of variable-order fractionaloperators: a review , Proc. A., 476 (2020), pp. 20190498, 32.[36]
H. Pedro, M. Kobayashi, J. Pereira, and C. Coimbra , Variable order modeling of diffusive-convective effects on the oscillatory flow past a sphere , Journal of Vibration and Control,14 (2008), pp. 1659–1672.[37]
T. R. Prabhakar , A singular integral equation with a generalized Mittag Leffler function inthe kernel , Yokohama Math. J., 19 (1971), pp. 7–15.[38]
L. Ramirez and C. Coimbra , A variable order constitutive relation for viscoelasticity , Annalender Physik, 16 (2007), pp. 543–552.[39]
S. Samko , Fractional integration and differentiation of variable order: an overview , NonlinearDynam., 71 (2013), pp. 653–662.[40]
S. G. Samko , Fractional integration and differentiation of variable order , Anal. Math., 21(1995), pp. 213–236.[41]
S. G. Samko and R. P. Cardoso , Integral equations of the first kind of Sonine type , Int. J.Math. Math. Sci., (2003), pp. 3609–3632.[42]
S. G. Samko and R. P. Cardoso , Sonine integral equations of the first kind in L p (0 , b ), Fract.Calc. Appl. Anal., 6 (2003), pp. 235–258.[43] S. G. Samko and B. Ross , Integration and differentiation to a variable fractional order , Inte-gral Transform. Spec. Funct., 1 (1993), pp. 277–300.[44]
G. Scarpi , Sopra il moto laminare di liquidi a viscosist`a variabile nel tempo , Atti. Accademiadelle Scienze, Isitituto di Bologna, Rendiconti (Ser. XII), 9 (1972), pp. 54–68.[45]
G. Scarpi , Sulla possibilit`a di un modello reologico intermedio di tipo evolutivo , Atti Accad.Naz. Lincei Rend. Cl. Sci. Fis. Mat. Nat. (8), 52 (1972), pp. 912–917 (1973).[46]
G. Scarpi , Sui modelli reologici intermedi per liquidi viscoelastici , Atti Accad. Naz. LinceiRend. Cl. Sci. Fis. Mat. Nat. (8), 107 (1973), pp. 239–243.ARIABLE-ORDER FRACTIONAL CALCULUS [47] D. Sierociuk, W. Malesza, and M. Macias , Derivation, interpretation, and analog modellingof fractional variable order derivative definition , Appl. Math. Model., 39 (2015), pp. 3876–3888.[48]
W. Smit and H. de Vries , Rheological models containing fractional derivatives , Rheol. Acta,9 (1970), pp. 525–534.[49]
M. Stynes , Fractional-order derivatives defined by continuous kernels are too restrictive , Appl.Math. Lett., 85 (2018), pp. 22–26.[50]
H. Sun, A. Chang, Y. Zhang, and W. Chen , A review on variable-order fractional dif-ferential equations: mathematical foundations, physical models, numerical methods andapplications , Fractional Calculus and Applied Analysis, 22 (2019), pp. 27–59.[51]
H. Sun, W. Chen, and Y. Chen , Variable-order fractional differential operators in anomalousdiffusion modeling , Physica A: Statistical Mechanics and its Applications, 388 (2009),pp. 4586 – 4592.[52]
A. Talbot , The accurate numerical inversion of Laplace transforms , J. Inst. Math. Appl., 23(1979), pp. 97–120.[53]
D. Tavares, R. Almeida, and D. F. M. Torres , Caputo derivatives of fractional variable or-der: numerical approximations , Commun. Nonlinear Sci. Numer. Simul., 35 (2016), pp. 69–87.[54]
E. C. Titchmarsh , Introduction to the Theory of Fourier Integrals , Chelsea Publishing Co.,New York, third ed., 1986.[55]
L. Trefethen, J. Weideman, and T. Schmelzer , Talbot quadratures and rational approxi-mations , BIT, 46 (2006), pp. 653–670.[56]
J. A. C. Weideman , Improved contour integral methods for parabolic PDEs , IMA J. Numer.Anal., 30 (2010), pp. 334–350.[57]
J. A. C. Weideman and L. N. Trefethen , Parabolic and hyperbolic contours for computingthe Bromwich integral , Math. Comp., 76 (2007), pp. 1341–1356.[58]
M. Zayernouri and G. E. Karniadakis , Fractional spectral collocation methods for linearand nonlinear variable order FPDEs , J. Comput. Phys., 293 (2015), pp. 312–338.[59]
F. Zeng, Z. Zhang, and G. E. Karniadakis , A generalized spectral collocation method withtunable accuracy for variable-order fractional differential equations , SIAM J. Sci. Comput.,37 (2015), pp. A2710–A2732.[60]
X. Zhao, Z.-Z. Sun, and G. E. Karniadakis , Second-order approximations for variable orderfractional derivatives: algorithms and applications , J. Comput. Phys., 293 (2015), pp. 184–200.[61]
X. Zheng, H. Wang, and H. Fu , Well-posedness of fractional differential equations withvariable-order Caputo-Fabrizio derivative , Chaos Solitons Fract., 138 (2020), pp. 109966,7.[62]
X. Zheng, H. Wang, and H. Fu , Analysis of a physically-relevant variable-order time-fractional reaction-diffusion model with Mittag-Leffler kernel , Appl. Math. Lett., 112(2021), pp. 106804, 7.[63]