EFFECT OF PERTURBATION IN THE NUMERICAL SOLUTION OF FRACTIONAL DIFFERENTIAL EQUATIONS

. The equations describing engineering and real-life models are usually derived in an approximated way. Thus, in most cases it is necessary to deal with equations containing some kind of perturbation. In this paper we consider fractional diﬀerential equations and study the eﬀects on the continuous and numerical solution, of perturbations on the given function, over long-time intervals. Some bounds on the global error are also determined.


1.
Introduction. Although the history of fractional calculus dates back to the 17th century, the extensive analysis of differential equations of non integer order is quite recent since the interest for their practical application has emerged only during the last few decades.
The increasing interest towards the applications of fractional operators is motivated by several reasons, in particular by the possibility of taking into account the hereditary effects commonly observed in real-life systems when the action of some external stimulation does not have just instantaneous effects but depends on the whole past dynamic of the system.
From a more practical point of view, it is the presence of a further parameter (namely the fractional order) which has raised most of the interest in contexts such as control theory, electro-magnetism and other industrial and engineering applications. Most of the models are indeed derived by matching experimental data into a theoretical model, but the fitting of the observed data can be usually made just in an approximated way; the possibility of suitably tuning the fractional order clearly helps to device more accurate models.
Nevertheless, although fractional calculus offers the possibility of improving the reliability of several models, it cannot be disregarded that the main parameters and functions involved in any model can be only estimated, often in a quite rough way. It is therefore necessary to take into account that the derivation and construction of the models lead to equations affected by some kind of perturbation; an important task is therefore to study the effects of these perturbations on the behavior of the solution and their consequent impact in numerical simulations.
The analysis of perturbations for systems of integer order has been extensively investigated in the past and several works are available today in the literature (see e.g. [4,Sec. 7.1], [22] and the bibliography therein). This is not the case however of FDEs for which, as far as we know, this topic has not been yet adequately studied. The aim of this work is therefore to fill this gap and perform an analysis of perturbations for the solution of FDEs.
Some results on data perturbations can be found in the Diethlem's book [9] where, however, the analysis is confined to bounded intervals. Since our main interest concerns what happens over long time intervals, our starting point is a study carried out by de Hoog and Anderssen in [8] on non singular Volterra integral equations.
Thus, inspired by the analysis presented in [8], we will study the effects of perturbations on the true solution and we will investigate how they affect the solution obtained numerically by the product-integration trapezoidal rule, a method which is widely used for the numerical solution of FDEs.
The paper is therefore organized as follows. In Section 2 we recall some introductory material on fractional differential equations and we introduce the class of equations which will be the subject of our investigation. Thus we first analyze, in Section 3, the effects of perturbations on the exact solution and hence we move to study, in Section 4, their influence on the numerical solution and we provide some bounds for the global error. Section 5 is devoted to present some numerical experiments with the aim of validating the theoretical findings and some concluding remarks are discussed in Section 6 together with some possible extensions of the results presented in this paper. 2. Fractional differential equations. In this paper we will consider the initial value problem for the linear FDE where 0 < α < 1 is the fractional (i.e. non integer) order and 0 D α t denotes the fractional derivative according to the Caputo's definition

EFFECT OF PERTURBATION IN THE NUMERICAL SOLUTION OF FDES 2681
The equation (1) can be rewritten as the corresponding Volterra integral equation (VIE) where F (t) denotes the contribution of the sourcing term which can be explicitly evaluated as From several results available in the literature (e.g., see [2,9]) we know that the continuity on [0, +∞) of k(t) and f (t) is sufficient to assure the existence of a unique solution to (2).
In this paper we want to investigate about the sensitiveness of the exact solution of (1) and of its numerical approximation obtained by trapezoidal product integration method, with respect to certain classes of perturbations on the function k. Therefore, we will consider the problem (1) in which the kernel k(t) is replaced bỹ k(t), and δ(k(t)) = k(t) −k(t) is some perturbation that we will assume sufficiently small. The corresponding perturbed VIE (2) becomes The first task in our analysis will be devoted to derive some bounds for the error evaluated as the difference between the exact solution of the original problem and the exact solutions of the perturbed problem, namely and only after we will consider numerical solutions for the exact and the perturbed problems and we will study the corresponding errors.
3. Influence of perturbations on the exact solution. We start the analysis by studying the boundedness of the solution of the problem under investigation. To this purpose, after denoting and assuming the main hypothesis we are able to show that the exact solution y(t) of equations (1) and (2) turns out to be bounded (this result has already been proved in [2]; we report here a simple alternative proof for the sake of completeness).
and η must be replaced byη in (6). In order to analyze the effect of the perturbation on the analytical solution, under hypotheses (4) and (7), we consider the perturbation error equation. From (2) and (3) the error δy(t) = y(t) −ỹ(t) satisfies Theorem 3.2. Assume that (4) and (7) hold for equations (2) and (3) respectively, then the perturbation error δy in the solution satisfies Proof. Since for Theorem 3.1 y(t) is bounded, the term 1 t 0 (t−s) α−1 δ(k(s))y(s)ds can be seen as the bounded forcing function of a weakly singular VIE. Thus, The proof comes straightforward by applying the result in (6) to (10).

Influence of perturbed data on the numerical solution.
To study the effects of perturbations on the numerical solution, we consider a uniform mesh I N = {t n = nh} n=0,...,N , with h a positive step-size and we rewrite equation (2) in a piecewise way on the mesh I N as for n = 0, . . . N . The trapezoidal product integration method (see [12,14,24]) consists in replacing k(s)y(s) in each subinterval by the first-degree interpolant polynomial where y j is the numerical approximation to the solution y(t j ). The discretization scheme therefore reads as for n = 0, . . . , N . This method is commonly (and successfully) applied also to FDEs, mainly in the framework of predictor-corrector methods [10]. The convergence analysis of the discretization method (12) over a finite time interval, for h → 0 and constant N h, can be found in [11,12].
This section is devoted to study the behavior of the numerical method (12), as N → +∞ and h > 0, when applied to problem (2) and to its perturbed version (3). In practice, we are looking for bounds on the global error y(t n ) − y n and on the total error y(t n ) −ỹ n , (withỹ n the numerical solution of the perturbed problem) which are constant with respect to N ∈ N and that hold for all h ∈ (0, h * ] for some h * > 0. From now on we will refer to the following set of assumptions holding for the function k(t) in equation (2): actually, the assumption b) is the same introduced in (4) but it is reported here for clearness.
The following lemma represents the link between the continuous and discrete classes of problems considered here.
Lemma 4.1. Assume that, for the function k(t) in equation (2), the assumptions a) and c) hold and that sup t≥0 K α (t) < +∞, furthermore, let Then, we have Proof.
The last inequality holds since (t n − s) α−1 ≤ h α−1 for any s ∈ [t j , t j+1 ], for 0 ≤ j ≤ n − 2. Thus, for any N ∈ N, Here we have used the fact that s − t n−1 ≤ h and we have bounded the term tn tn−1 (tn−s) α h |k(t n−1 )|ds, by using the same technique as above. The result comes straightforward by taking into account the boundedness of K α (t) and the form of A in (13).
Observe that, in the previous result, if hypothesis b) holds there surely exists h * > 0 such that Let us now consider the local error of the discretization (12) and we provide a bound for the global error Theorem 4.2. Assume that a), b) and c) hold. Then, for any h ∈ (0, h * ], with h * defined in (15), it is Proof. By subtracting (12) from (11), we get, for n = 0, 1, . . . , N, N ∈ N, Therefore we obtain Thanks to the hypothesis b), we choose h * > 0 according to (15), from which the proof follows.
The boundedness of the global error comes directly from the following theorems.
Theorem 4.3. Assume that, for equation (2), it holds Then for any h ∈ (0, h * ] the local error T n (h) satisfies Proof. For 0 ≤ j ≤ n − 1, denote with T n,j (h) the restriction of the local error (16) to the interval [t j , t j+1 ], i.e. Thus, Hence, tn 0 |(k(x)y(x)) |dx, and the proof follows since the well-known identity Γ(α + 1) = αΓ(α). Theorems 4.2 and 4.3 allow us to finally state the following result, which implies a stable behavior of the numerical solution for any fixed h > 0 satisfying (15). For h → 0, a convergence result to the solution of (1) can be easily derived from here.
where the constant C is given in (19).
Proof. According to Theorem 3.1, assumptions b) and (5) imply that the analytical solution to equation (2) is bounded. Therefore, if c) is true, the condition sup t≥0 t 0 |k(x)y (x)|dx < +∞ is sufficient to guarantee that (19) is satisfied. In order to complete the proof, we consider (18) and use the result in Theorem 4.2 on the boundedness of the local error.
When the kernel is affected by perturbations the error analysis of the numerical methods has to take into account an additional term depending on the magnitude of the perturbations.
Let us consider, as in Section 2, the FDE (1) or, equivalently, the corresponding VIE (2) with the perturbed kernelk(t). The numerical approximation obtained by means of the trapezoidal product integration rule is clearlỹ for any n ∈ N.
By assuming that alsok(t) satisfies the same assumptions a), b) and c) made for k(t), namely it is clear that there exists h * * such to guaranteẽ whereÃ has the same form of A in (13) with k and k replaced byk andk respectively.
We are therefore able to present the main result allowing to establish a bound for the error between the exact solution of the original problem and the numerical solution which one actually computes by working on the perturbed problem.
Theorem 4.5. Assume that a), b), and c) hold for k in equation (2) and that the perturbed kernelk satisfies a * ), b * ), and c * ); furthermore assume that (5) holds. Then for any h ∈ (0, h * * ] it is whereT n (h) is the local error of (21).
Proof. The result in Lemma 4.1 applied to equation (3) with a * ), b * ), and c * ) gives (24) In order to obtain a bound for the total error in trapezoidal product integration approximation plus perturbation, we rewrite it as the sum of the following two terms y(t n ) −ỹ n = y(t n ) −ỹ(t n ) + (ỹ(t n ) −ỹ n ). Thus, from (10) and (18)  The sequence of the stepsizes {h n = t n+1 − t n , n = 0, 1, ..} is, in this case, strictly increasing and bounded, with Graded meshes are of great importance for accuracy reasons since they allow to avoid the drop in the order of convergence which usually happens when product integration methods are applied to problems, like (1), whose exact solution has a non smooth behavior at the origin. It is straightforward to see that the results presented in this paper can be proved in this more general context.

Numerical experiments.
In this section we present some numerical experiments with the aim of validating the theoretical results derived in the previous sections. To this purpose, we consider, for a real parameter λ > 0, the test problem with the kernel function k(t) = exp(−λt) together with an oscillating source term where ω > 0 is the oscillatory frequency.
In order to provide an estimation of the constant η in the main hypothesis (4) we preliminarily observe that, after replacing in K α (t) the exponential with its series expansion and integrating term by term, it it is possible to evaluate where E α,β (z) is the Mittag-Leffler function [16], a function playing a fundamental role in fractional calculus. As we can see from Figure 1 this function has just one point of maximum and converges towards 0 when t → ∞.  The unique extrema of K α (t) can be hence located by finding the root of its derivative and a straightforward derivation enable us to write thus allowing to evaluate the value of η by means of standard numerical techniques.
To this purpose we observe that η appears as a monotonically decreasing function with respect to λ (see Figure 2). Therefore, condition (4) is surely satisfied for sufficiently large values of λ, which obviously depend on α. The solution obtained by means of the trapezoidal product integration method, for α = 0.8, λ = 5 and ω = 1 is presented in Figure 3, where the initial condition y(0) = 1 has been considered. This solution is obtained with a reasonably small step-size h = 2 −12 ≈ 2.4 × 10 −4 in order to obtain a quite accurate numerical approximation that, due to the convergence properties of the method, can be used as reference solution. In the same Figure 3 the bound (6) derived in Theorem 3.1 is also reported for comparison (we have numerically evaluated η ≈ 0.202948). and we consider only step-sizes h satisfying (15) (the corresponding value η + h α A, which must be less than 1, is also reported in the first column of each table).
As we can observe from the above results, the theoretical predictions of Theorem 4.2 are verified in a quite clear way.
To introduce a perturbed problem we consider now, for a small δλ, a different valueλ = λ + δλ for the perturbed kernelk(t) = exp(−λt). Obviously alsoλ is chosen such that (7) holds.
In Figure 4 we compare the difference δy(t) between the exact and perturbed solutions whenλ = 6 together with the bound (9) derived in Theorem 3.2.    For the same valueλ = 6 we present now in Tables 4, 5 and 6 the comparison between the global errors for the perturbed problem with kernelk(t) = e −λt and the bounds derived thanks to Theorem 4.5.     tables indeed show the difference between a reference solution of the original problem and the numerical solutions of the perturbed problem; no kind of convergence can be expected, as h → 0, between the solutions of two different problems and these numerical experiments just aim to show how this difference can be predicted, and in some sense taken under control, by means of Theorem 4.5. As expected, however, the sequence of the discretized differences sup n=0,N |ẽ n | converges, as h → 0, to the continuous differences sup 0 |δy(t)| which for the problems presented in Tables 4, 5 and 6 has been numerically estimated, respectively, as 9.213083 × 10 −2 , 6.318025 × 10 −2 and 4.962170 × 10 −2 .
Remark 2. In our tests we observe that the bounds (18) on the global errors are quite sharp. However, for some problems where, for example, 0 < 1 − η 1, the upper bound in (18) may be considerably large, unless considering unreasonably small stepsizes. These cases may be only qualitatively treated by the theory developed here and no information about the magnitude of the global errors can be deduced. 6. Concluding remarks and extension of the results. In this paper we have studied the effects of perturbations to a class of non-autonomous linear FDEs and some bounds for the errors have been provided both for the exact and the numerical solution obtained by means of the trapezoidal product integration method.
For ease of presentation in this paper we have considered just some simplified cases; however most of the results can be extended to more general situations. One example is, as observed in Remark 1, the case of graded meshes which are often more appropriate than uniform meshes when using product integration methods to approximate (1).
In a more realistic situation, some perturbations occur also on the initial datum y 0 and on the forcing term f (t); the essence of the results in Theorems 3.2 and 4.5 does however not change once we require that sup t≥0 δ(F (t)) = y 0 −ỹ 0 + 1 Γ(α) sup t≥0 t 0 (t − s) α−1 (f (s) −f (s))ds < +∞.