A computational approach to exponential-type variable-order fractional differential equations

We investigate the properties of some recently developed variable-order differential operators involving order transition functions of exponential type. Since the characterisation of such operators is performed in the Laplace domain it is necessary to resort to accurate numerical methods to derive the corresponding behaviours in the time domain. In this regard, we develop a computational procedure to solve variable-order fractional differential equations of this novel class. Furthermore, we provide some numerical experiments to show the effectiveness of the proposed techniques.


Introduction
Constant-order (CO) fractional operators are powerful tools used to describe various phenomena displaying non-local properties or persistent memory (i.e., non-locality in time). However, in certain scenarios the nature of non-localities can itself vary with respect to time and/or space and one is therefore forced to resort to variable-order (VO) operators (see, e.g., [2,7,8]).
Over the years, several proposals for generalizing CO operators to VO have appeared in the literature. Often these attempts just consist in a simple replacement of the constant fractional order with a function of the independent variable. It is worth noting, however, that even this simple replacement can be performed in more than one way (see, e.g., [42]) and, as pointed out in [24,28], these procedures may come with some shortcomings. For VO derivatives obtained through a naive replacement, as mentioned above, it is indeed often difficult to determine the corresponding VO fractional integrals so that the operators satisfies a generalised fundamental theorem of calculus, thus providing a means to simplify theoretical and numerical analysis.
Recently, a notion of VO fractional derivative, based on the pioneering work of Giambattista Scarpi [29,30] dating back to the early seventies, has been revived. In particular, in [6] this definition of VO derivative has been reviewed and framed in terms of Luchko's theory of generalised fractional derivatives and integrals [17][18][19][20].
Scarpi's formulation of VO integrals and derivatives is obtained as an extension of CO operators in the Laplace domain, as opposed to time domain. Defining and studying operators in the Laplace domain seem to simplify the theoretical treatment of these mathematical objects. Specifically, it turns out that Scarpi's fractional integrals and derivatives are, in the time domain, expressed in terms of Volterra-type integro-differential operators with weakly-singular kernels and that these kernels satisfy the Sonine condition [32], i.e., one of the defining requirements for the validity of a generalised fundamental theorem of calculus [4,12,25,26].
The numerical treatment of integral and differential equations with these VO operators turns out to be more challenging. This is due to the lack of an explicit time-domain representation for Scarpi's operators. In fact, while for VO operators defined in the time domain there exists well-established numerical techniques to solve the corresponding VO fractional differential equations (VO-FDEs) (see, for instance, [5,33,[38][39][40][41]43]) the same does not apply to our case. It is therefore necessary to develop some specific approaches to tackle VO-FDEs in this new framework.
The aim of this work is to provide a further characterisation of VO operators defined in the LT domain. In particular, we will expand on some preliminary results in [6], originally presented only as numerical experiments, providing precise statements and rigorous proofs. Moreover, we will investigate suitable numerical tools to solve VO-FDEs involving these operators. Specifically, we will discuss convolution quadrature rules based on Lubich's fundamental works [14][15][16]. A detailed error analysis will also be presented in order to evaluate weights of the convolution quadrature rule accurately.
Although several VO functions α(t) can be considered, in this work we will focus solely on those of exponential type, thus expanding and building upon previous research [6]. The rationale behind this choice is the fact that such a behaviour provides a smooth transition between an initial order α 1 to a final order α 2 and the dependence of α(t) on the three parameters (initial and final orders, and the transition ratio) is flexible enough to find potential applications in a wide variety of situations. Nonetheless, the general analysis presented here can be applied, albeit with some technical differences, to other VO functions.
This work is organised as follows. In Section 2 we provide a brief review of Scarpi's VO fractional calculus. Section 3 describes an order transition of exponential type and analyses the main properties of the corresponding VO integral and derivative, and their implications for the VO relaxation equation. In Section 4 we describe a numerical approach to solve VO-FDEs and we discuss some technical details concerning the accurate computation of weights in the proposed scheme. Lastly, in Section 5 we present some numerical experiments in order to test and validate the proposed computational approach. In Section 6 we offer some concluding remarks and an outlook on future research.

Variable-order fractional integrals and derivatives
Let us begin by reviewing the key ideas of Scarpi's approach to VO calculus. To this end, we will follow [6], to which we refer the interested reader for a more detailed discussion on the subject. 1 The (CO) Riemann-Liouville integral of order α > 0 is defined in terms of a convolution integral as Analogously, the (CO) Caputo fractional derivative of order 0 < α < 1 reads .
These two operators satisfy the so-called fundamental theorem of fractional calculus, i.e., and a generalisation of the fundamental theorem of calculus, i.e., The validity of these properties can be directly linked to the fact that the kernels ψ(t) and ϕ(t) satisfy the Sonine condition [32], i.e., Furthermore, it is worth pointing out that the Laplace transforms (LT) of ψ(t) and ϕ(t), appearing in Eq. (2.1) and Eq. (2.2), read Ψ(s) := L ψ(t) ; s = s −α , Φ(s) := L ϕ(t) ; s = s α−1 . (2.5) Now, in order to generalise the definitions in Eq. (2.1) and Eq. (2.2) let be a function such that: A2. the analytical form of the LT of α, i.e., A(s) = L α(t) ; s , is known; Remark. While A1 guarantees the existence of the LT of α, A2 is not strictly necessary but significantly simplifies both the theoretical and numerical analyses. A3 is instead essential to define integral kernels satisfying the Sonine condition (see [6]).
Scarpi's proposal, later reformulated and extended in [6], sparks from the simple observation that in the constant order case α(t) = α ∈ (0, 1) for all t ∈ [0, T ], then A(s) = α/s and, therefore, the LTs of the kernels ψ and ϕ in Eqs. (2.5) can be rewritten as Ψ(s) = s −sA(s) and Φ(s) = s sA(s)−1 . Following this line of thought, considering now a function α(t) that satisfies the conditions A1-A3, one can define the operators with kernels such that In other words, Scarpi's operators are defined in terms of kernel functions specified in the Laplace (complex) domain rather than on the time domain.
We observe that analytical expressions for ψ α (t) and ϕ α (t) are, in general, not known in this approach. Hence, the treatment of (2.6) or (2.7) turns out to be numerical for the most part. Thus, developing numerical techniques capable of addressing these sort of problems is, in fact, among the main objectives of this work.
It is now easy to see that ψ α (t) and ϕ α (t) naturally satisfy the Sonine equation since, in fact, in the Laplace domain it holds that as one can easily check from the expressions in Eq. (2.8). This means that ψ α (t) and ϕ α (t) form a Sonine pair and therefore one finds that and Scarpi's operators satisfy the fundamental theorem of fractional calculus.
3 Variable-order transitions of exponential type The general framework described in Section 2 applies to any function α(t) that satisfies the assumptions A1, A2 and A3. In this section we take a closer look at VO integrals and derivatives by selecting a specific class of functions for α(t). The reason for this choice consists in the fact that the investigation of qualitative properties of both operator kernels and solutions of FDEs with a general VO function α(t) would represent a rather unfeasible task unless one fixes a specific function. Moreover, selecting a specific VO function α(t) significantly simplifies the development of computational procedures, as we shall discuss in a forthcoming section. We focus therefore on the case of transition functions α(t), defied on R + , describing a smooth transition from the order 0 < α 1 < 1, at t = 0, to 0 < α 2 < 1 as t → +∞. More precisely, we consider a family of exponential transitions parametrised by a transition rate c > 0 and given by which satisfies assumption A1, A2, and A3 and whose LT reads The objective of this work is to investigate VO operators with α(t) as in (3.10) from a purely mathematical perspective. We observe, however, that α(t) can be viewed as a solution of a relaxation equation, which is ubiquitous in physics.
For the sake of convenience, and to emphasise the dependence on the asymptotic orders α 1 and α 2 in the exponential transitions, we adopt a slightly different notation compared to the more general one presented in (2.8). Specifically, we denote the LTs of the generalised VO kernels by Ψ α 1 ,α 2 (s) and Φ α 1 ,α 2 (s), i.e., Similarly, we denote the corresponding kernels (of S I For the sake of a lighter notation we omit to mention explicitly the dependence of the kernels on the transition rate c.
Furthermore, we will refer to the usual kernels of the fractional integral and derivative of constant order α 1 by denoting the corresponding LTs. Similar expressions hold for the constant order α 2 . The convergence region of the LT A(s), as well as that of Ψ α 1 (s), Ψ α 2 (s), Φ α 1 (s), and Φ α 2 (s), includes the whole positive complex half-plane. Moreover, A(s) has two poles at s = 0 and s = −c. Additionally, in order to make Ψ α 1 ,α 2 (s) and Φ α 1 ,α 2 (s) single-valued one has to choose a branch-cut, that we set on the negative real axis.
Analytical expressions for ψ α 1 ,α 2 (t) and ϕ α 1 ,α 2 (t) are not known. Nonetheless, a lot can be learned about these functions from their LTs. For this reason we shall first briefly review some basic properties of the LT.

Basic properties of the Laplace transform
Let f (t) be a piece-wise continuous function of exponential order a, i.e., there exists an M > 0 such that for some t 0 ≥ 0 it is |f (t)| ≤ M e at for all t ≥ t 0 . Then, for any s ∈ C, with ℜ(s) > a the LT exists (i.e., the above integral converges). Moreover, one also finds that The following standard theorems will play an important role in the derivation of the main results presented in the following.
Theorem 1 (Initial value (IV) theorem for the LT (see [13])). Let f be continuous for t > 0, of exponential order and with f ′ (t) possessing at worst an integrable singularity at t = 0. Then if the two limits exist.
Theorem 2 (Final value (FV) theorem for the LT (see [3,11])). Let f (t) : [0, T ] → R and assume that its LT F (s) does not have any singularities in the right half-plane H c := s ∈ C | ℜ(s) ≥ 0 , except for possibly a simple pole at the origin. Then, if the two limits exist.
For a more detailed discussion of these results and of the Laplace transform method, in general, we refer the interested reader to [23].

Connecting VO and constant-order kernels
In [6], based on some numerical experiments, it was observed that ψ α 1 ,α 2 (t) and ϕ α 1 ,α 2 (t) tend to approach, asymptotically in t, the behaviour of constant-order kernels. This intuition is further supported by the following analytical results.
Proposition 3. Let 0 < α 1 , α 2 < 1 and c > 0. Then Proof. We begin by observing that which implies, in light of Theorem 1, that Similarly, observing that and therefore, in light of Theorem 2, one has which concludes the proof.
A similar result can also be proven for the kernel ϕ α 1 ,α 2 (t) of the VO fractional derivative.
Proof. Analogous to the proof of Proposition 3.

Relaxation equation
Let us now consider the fractional relaxation equation the VO fractional derivative realizing the exponential order transition (3.10), λ > 0 a real parameter and y 0 any positive initial value.
In the Laplace domain Eq. (3.12) can be formulated as The exact solutions of the corresponding relaxation equations with standard Caputo fractional derivative (2.2) of constant order 0 < α 1 < 1 and 0 < α 2 < 1 are known to be respectively .
Finding a general result concerning exact solution of the VO-FDE (3.12) appears to be a rather hard task. Nonetheless, we can still point out a relationship between a solution of (3.12) and the one of the corresponding CO counterparts at early and late times.
There exist a function f (t) such that the solution y(t) of the relaxation equation (3.12) satisfies Moreover, if H(s) does not have singularities on H c , there exists a function g(t) such that y(t) = E α 2 (−t α 2 λ)y 0 + g(t)y 0 , lim t→∞ g(t) = 0.

Proof. Recall that
As in the proof of Proposition 3 we write Then, one finds that .
This result is particularly important since it shows that solutions of the relaxation VO-FDE (3.12) converge to solutions of the CO relaxation FDEs of order α 1 and α 2 respectively for t → 0 + and t → ∞. This result had been already observed in [6], albiet only at the numerical level. In Figure 2 we show the difference between solutions of VO and CO relaxation equations where, for simplicity we considered y 0 = 1.

A numerical approach to exponential-type variable-order differential equations
Let us now introduce, for the exponential VO transition α(t) described by (3.10), the VO-FDE which, after applying S I α(t) 0 to both sides, and in view of (2.9) (and with the notation introduced in Section 3), can be equivalently reformulated as the integral equation y(t) = y 0 + S I α(t) 0 f (t, y(t)), i.e.
The main difficulty to numerically approximate solutions of (4.14) lies in the absence of an analytical formulation of ψ α 1 ,α 2 (t). Since just its LT Ψ α 1 ,α 2 (s) is known, it is natural to exploit convolution quadrature rules (CQRs) introduced by Lubich in his pioneering works [14,15], and later discussed in [16]. The main feature of these rules is that they allow to approximate the convolution integral in (4.15) without requiring the explicit knowledge of the convolution kernel, but rather just of its LT. The application to (4.15) is however not straightforward. Not only because Lubich's theory requires the fulfillment of some assumptions to be valid, but also in view of some difficulties for the computation of convolution weights.
In this section, after briefly review CQRs in Lubich's framework, we will study their application in the context of VO operators and, specifically, we will devise a strategy for the accurate computation of convolution weights supported by a detailed error analysis.
A discrete CQR extending the (ρ, σ) LMM to solve (4.15), and preserving the same convergence order, say p, of the LMM, reads as where ω n are the coefficients in the expansion of The second sum in (4.17) is introduced to deal with a possible lack of smoothness of the solution at the origin (which is typical of solutions of FDEs and it is reasonably expected in this VO setting since the result in Proposition 3).
Since the method is expected to inherit the convergence order p = 1 of the backward Euler method [14,Theorem 3.1], there is no need of introducing starting weights in (4.17) and therefore its extension to the VO integral equation (4.15) simply reads as where the convolution weights ω n are the coefficients in the expansion of Remark. It can be of interest to observe that in the CO case the CQR obtained by exploiting the generating function of the backward Euler method has the same weights of the Grünwald-Letnikov (GL) discretization scheme (we refer, for instance, to [3,9,27,31] for more details about this scheme).
The method proposed here can be in some sense considered as a sort of generalization of the GL scheme to VO operators (2.6) and (2.7).

Computation of weights
Computing convolution weights ω n in (4.19) is a challenging task which must be accomplished with high precision to avoid loss of accuracy and, possibly, in a fast way. From (4.20) we observe that ω n are coefficients in the Taylor expansion of Ψ [h] α 1 ,α 1 (ξ) around the origin. Therefore, they can be expressed in terms of Cauchy integrals over contours C encompassing ξ = 0 in the region of analiticity of Ψ [h] α 1 ,α 2 (ξ) is analytic the whole complex plane C except the real line [1, +∞). Hence, we select C as a circle with center at the origin and radius 0 < ρ < 1.
a formula which lends itself to being calculated by means of an efficient fast Fourier transform algorithm. The quadrature rule (4.22) depends on two parameters: the contour radius ρ and the number L of quadrature nodes. They are strictly related since their ratio determines the step-size of the quadrature rule and, hence, the discretization error. However, in the finite precision arithmetic of computers, the parameter ρ plays a further important role: it indeed determines the distance of the integration contour from singularities of Ψ [h] α 1 ,α 2 (z) and choosing a contour too close to the singularities may leads to strong round-off errors. It is therefore necessary to separately study discretization and round-off errors and take both of them into account for selecting parameter ρ and L. The aim is minimizing the error and taking it below a given tolerance τ > 0 (possibly close to the precision machine), namely to obtain |ω n − ω (ρ,L) n | < τ, n = 0, 1, . . . N.

Discretization error
The main error source for computation of weights ω n is the discretization error in the quadrature rule (4.22). In this respect we are able to provide the following result. Proposition 6. Let 0 < ρ < r < 1 and a = log(r/ρ). Then for any n ∈ N and L > 0 it is Proof. Consider the 2π-periodic function v n (θ) = Ψ [h] α 1 ,α 2 ρe iθ e −inθ , such that ω n = ρ −n 2π 2π 0 v n (θ)dθ.

Round-off errors
The discretization error (4.23) is not the only error affecting the accuracy of computed weights ω n . When high accuracy and/or a large number of weights is requested, round-off errors due the finite-precision arithmetic of computers must be taken into account.
To introduce round-off errors in our analysis, we represent weights which are actually computed as with |ϵ k | ≈ ϵ and ϵ the precision machine. Thus, if we consider since 0 < ρ < 1, and in view of Lemma 8 below, we obtain the following rough estimation holding for any n = 0, 1, . . . , N and for a reasonably small step-size h < 1 − r.
The overall error is hence , and the goal is to keep the round-off error smaller than the discretization error τ , thus to possibly neglect it. We can therefore fix a safety factory 0 < F s < 1, for instance F s ≈ 0.01 ÷ 0.1, and impose the round-off error to be proportional to F s τ , thus to select (4.24)

A strategy for parameters selection
The proposed strategy to select parameters ρ, r and L is to first fix a value for r < 1, which for convenience must be chosen close to 1, and determine by (4.24) a suitable contour radius ρ < r on the basis of the number N of nodes which must be computed. When (4.24) leads to a value of ρ ≥ r, the value of r must be suitably increased. The tolerance τ should be chosen as small as possible to consider weights as computed in a virtually exact way, but not too strict to impose a huge number of quadrature nodes to achieve it. Therefore, τ is selected in the range 10 −12 ÷ 10 −13 .
Once ρ is selected, it is necessary to provide an estimation for the value of M r,ρ in (4.23). It seems not possible to provide a theoretical estimation for M r,ρ and therefore it will be necessary to sample Ψ [h] α 1 ,α 2 z on a sufficiently large number of points z ∈ C, with Im z ≥ − log(r/ρ) and determine its maximum.
Finally, it is simple to use (4.23) to obtain the number L of nodes in the trapezoidal quadrature rule (4.22) to to achieve the accuracy |ω n − ω (ρ,L) n | < τ by selecting

Two technical results
In this subsection we collect two technical results which has been used for the estimation of round-off error in the previous subsection. We confine their presentation here with the only aim of lightening the reading of the paper.
where θ ξ = arg(1 − ξ) and A(x, y) = Proof. First rewrite Ψ [h] and observe that A(x, y) and B(x, y) are, respectively, the real and the imaginary part of the first term of the argument in the above exponential, i.e.
Hence, after evaluating the proof follows by means of elementary manipulations.

Numerical experiments
To validate the numerical scheme discussed in the previous section, we perform here some tests. Although the scheme can be applied, in principle, to any linear or nonlinear problem, we restrict experiments to just some linear problems since the nonlinear case could require a more detailed analysis.
It is quite natural to start by testing the method with the relaxation equation (3.12) for which a reference solution can be approximated, with high accuracy, by numerical inversion of the LT as explained, and already exploited, in Subsection 3.3 (it is worth mentioning that no exact solutions are in general available with VO-FDEs).
For different sets of values α 1 , α 2 , c and λ, in Figure 4 we show solutions y ⋆ (t) obtained by inversion of the LT and used as reference solutions.
In Table 1 we collect, for a decreasing sequence of the step-size h, errors E(h) at T = 4 between numerical solutions y h (T ) and references y ⋆ (T ), together with the Estimated Order of Convergence (EOC) evaluated as As we can clearly observe, the test confirms that, as predicted from theory, the scheme preserves the first order convergence in the VO case as well.   We introduce now in (3.12) a forcing term to obtain the VO-FDE The reference solutions, presented in Figure 5, are now evaluated by the same method proposed in this work but with the smaller step-size h = 2 −10 . Also with the test problem (5.26) the comparison between reference and numerical solutions allows to confirm the first-order convergence as one can clearly observe from Table 2.   We now build the 2-dimensional linear system of VO-FDEs by selecting p ± = cos θ ± √ cos 2 θ + 1 and choosing the parameter θ so that when the same system has CO derivatives of order α 1 and α 2 , it shows a stable behavior. In particular, θ is selected in order to eigenvalues of the matrix system of (5.27) lie in the sector | arg(z)| > σ, σ = max{α 1 , α 2 }π/2, and we used the value θ = 1.1 · σ in all experiments.
The reference solutions are presented in Figure 6. We have presented here a comparison with the solutions of the same system with CO derivatives of order α 1 and α 2 to show how the behaviour for small t (the box in each plot) and for large t follows the behaviour of solutions of the CO system with order α 1 and α 2 respectively.
Also with the linear system (5.27) the comparison between reference and numerical solutions allows to confirm the first-order convergence of the numerical method under investigation, as shown in Table 3.

Concluding remarks
In this work we have discussed in detail both theory and computation of some VO-FDEs with order transition functions of exponential type. Specifically, first we have investigated the relationships between VO operators and CO one at early and late times. We have presented a numerical method, based on Lubich's CQRs, that extends the backward Euler scheme to the VO case. In particular, we have developed an algorithm for the computation of weights based on a detailed error analysis. Numerical experiments have, further, shown the robustness of the proposed approach.
In order to obtain more accurate numerical solutions, the natural continuation of the proposed research will entail the development of higher-order methods, for which a preliminary analysis of solution regularity at the origin is necessary to properly select regularization weights [14][15][16]. Some seeds for this analysis are in Section 3 at least for what concerns relaxation equations. However this matter deserves a deeper discussion which is left to a future study.
It will be, moreover, of particular interest to extend the methodology presented in this study to other classes of transition functions α(t). In particular, to address the problem of VO operators with transition functions α(t) for which the analytic expression of their LT is not known (thus to avoid condition A2 in Section 2). Clearly, such a problem can only be tackled at the numerical. It is worth mentioning that this research line is mostly uncharted territory to date [36], although it has the potential of becoming a stimulating and prolific research topic in the near future.