Exponentially Convergent Numerical Method for Abstract Cauchy Problem with Fractional Derivative of Caputo Type

We present an exponentially convergent numerical method to approximate the solution of the Cauchy problem for the inhomogeneous fractional differential equation with an unbounded operator coefficient and Caputo fractional derivative in time. The numerical method is based on the newly obtained solution formula that consolidates the mild solution representations of sub-parabolic, parabolic and sub-hyperbolic equations with sectorial operator coefficient $A$ and non-zero initial data. The involved integral operators are approximated using the sinc-quadrature formulas that are tailored to the spectral parameters of $A$, fractional order $\alpha$ and the smoothness of the first initial condition, as well as to the properties of the equation's right-hand side $f(t)$. The resulting method possesses exponential convergence for positive sectorial $A$, any finite $t$, including $t = 0$ and the whole range $\alpha \in (0,2)$. It is suitable for a practically important case, when no knowledge of $f(t)$ is available outside the considered interval $t \in [0, T]$. The algorithm of the method is capable of multi-level parallelism. We provide numerical examples that confirm the theoretical error estimates.

1. Problem Formulation and Introduction.In this paper, we consider a Cauchy problem for the following fractional order differential equation: (1) Here, ∂ α t denotes the Caputo fractional derivative of order α with respect to t where u (n) (s) is the usual integer order derivative, n = ⌈α⌉ is the smallest integer greater or equal to α and Γ(•) is Euler's Gamma function.The operator ∂ α t provides a generalization of the classical differential operator ∂ ∂t = ∂ 1 t .For non-integer α, the action of Caputo fractional derivative is essentially nonlocal in time.In addition to that, the memory kernel from ∂ α t , α < 1 has a mild singularity at 0. These two facts have a profound impact on the analytical and numerical properties of solutions to fractional differential equation (1).If α < 1, this equation is called sub-parabolic.Similarly, when α > 1, the equation is called sub-hyperbolic.We direct the reader to [33] for a more concise introduction into the subject of fractional derivatives and the theory of associated ordinary differential equations.
The numbers ρ s > 0 and ϕ s < π/2 are called spectral parameters (characteristics) of A. In addition to the assumptions on the location of spectrum, we suppose that the resolvent of A: R (z, A) ≡ (zI − A) −1 satisfies the bound (3) (zI − A) −1 ≤ M 1 + |z| outside the sector Σ and on its boundary Γ Σ .Following the established convention [23], we will call such sectorial operators strongly positive.We accompany equation (1) with the usual initial condition (4a) u(0) = u 0 , for the solution and additional condition for its derivative, when 1 < α < 2: The theory of fractional Cauchy problems for differential operators was developed in the works [45,49,12].The abstract setting, considered here, has been theoretically studied in [34,3,31] for α ∈ (0, 1) then in [4] for α ∈ [1, 2) and, most recently, in [54].In the current work, we focus on the numerical evaluation of the mild solution to problem (1), ( 4) that is given by the following result.
The bulk of the existing research is devoted to the particular cases of (1) when A is specified as a strongly elliptic linear partial differential or, more generally, pseudodifferential operator with the domain D(A) that is dense in X [49,13,56].These cases also include the fractional powers of elliptic operators are encompassed by the class of strongly positive operators [16] considered in Theorem 1 and below.In this regard, the shape of Sp(A) justifies the choice of the range (0, 2) for α, as a maximally possible for the considered class of A (see [54] for a more detailed discussion).
There exists a considerable body of work devoted to numerical methods for evolution fractional differential equations (see [20,10,9,11] and the references therein).Philosophically, it can be subdivided into methods that directly approximate the components of (1), or its integral analogue, and those that make use of more elaborate solution approximations.The methods from the first class are sequential in nature and have algebraic convergence order that typically does not exceed 2, even for the multi-step methods [19], due to the intrinsic fractional-kernel singularity [52].In addition, at each time-step, these methods need to query the entire solution history in order to evaluate ∂ α t or J α , numerically.In the consequence of that, they are computationally costly and memory constrained.Nonetheless, the methods from this class are popular due to their simplicity [53], numerical stability [19] and the ability to handle non-smooth initial data [28].The second class of numerical methods is represented by the works [7,2,14,26,32], to name a few.These methods are based on the clever solution approximations that result in a time-stepping scheme requiring only a small number of previous solution states for the next state evaluation.With some exceptions (e.g., [26]), these methods are also O(h p ).
Spectral methods from [43,44,6,57] deserve a separate mention.Although formally belonging to the second class, they make use of the exponentially convergent contour-based propagator approximation, which permits to evaluate the transient component of the solution to the linear problem without time-stepping.The authors of these works, however, do not apply it to (1), (4) directly.Instead, they consider a special proxy problem ∂ t u + I 1−α Au = g where I α is a nonlocal operator equal to ∂ α t , if α < 1, or to J α , otherwise.It was shown in [44], that the existing methodology for parabolic problems [21,60,59] can be transferred to the mild solution of such proxy problem with all important numerical features of the solution algorithms preserved, including uniform exponential convergence for t ∈ [0, T ] and the capacity for multi-level parallelism.Despite being simple and efficient, the proxy-problem idea has certain ramifications when applied to (1), (4).Firstly, there is no easy way to incorporate the initial condition from (4b) into the proxy problem formulation, so all existing works consider u ′ (0) = 0. Secondly, the methods from [43,44,46,6] operating on the Laplace transform image of the right-hand side g are prone to errors when the original f from (1) is not given in the closed form.Hence, they are unsuitable for many applications.Meanwhile, formula (5), which serves as a base for our numerical method, does not require any extra knowledge about the right-hand side f ∈ W 1,1 ([0, T ], X) beside the values f (0) and f ′ (t), t ∈ [0, T ].In addition to that, the rigorous analysis from [44,57] addresses a version of the proxy problem where ∂ α is a Riemann-Liouville (RL) fractional derivative.Cauchy problems with RL derivative are simpler in the sense of propagator representation [54], but they are compatible with (1), (4) only under some additional assumptions.
It is fair to point out that the majority of the mentioned methods are designed to handle the nonlinear fractional differential equation, more general than (1).With the view of similar nonlinear extensions in mind, in this work we would like to prioritize those properties of the solution method for (1), (4), which will make such extensions possible.Let us for the moment assume that f = f (t, u).Then, representation (5) can be used as a base for the sequential time-stepping scheme [27,36] or as the fine propagator in a more parallelization-friendly ParaExp-type scheme [17].In both cases, the method will be free of the issues with approximating ∂ α t u in the vicinity of t = 0, provided that the proposed approximation of (5) converges uniformly.Such application scheme also justifies the use of a moderate in size final time T .If, more generally, we assume that Au = A(t, u), then the problem in question can be reduced to (1), (4) using collocation [5,23] or a similar in nature time-stepping scheme, inspired by [25].In such scenario, A(t, u) is approximated by A(t k , u k ) having spectral characteristics that may vary drastically with k (see Cahn-Hilliard equation from [15], for instance), and the right-hand side in the form A(t k , u k )−A(t, u), which makes sense only locally.Thus, the solution method should be able to reliably handle operators with arbitrary spectral parameters and right-hand sides that are unknown a priori.
Taking the aforementioned properties into account, below we devise an exponentially convergent approximation for (5) by building upon a well-established technique [42,55,61,21] that involves the application of a trapezoidal quadrature rule to the parametrized contour integral from (7).
In Section 2 we study a question regarding the choice of the suitable integration contour for such parametrization.The proposed time-independent hyperbolic contour Γ = Γ I is valid for the wide class of sectorial operators A with fixed ϕ s < π min 1  2 , 1 − α 2 and arbitrary ρ s > 0. The parameters of Γ I are derived using the set of constraints that utilizes all available analyticity of the propagator, and therefore, maximize theoretical convergence speed of the sinc-quadrature applied to S α (t).Section 3 is devoted to the development and justification of the numerical method.Using the moderate smoothness assumption u 0 ∈ D(A γ ), γ > 0, in Section 3.2 we propose an exponentially convergent approximation of S α (t)u 0 , that does not degrade for small t like the similar methods from [14,47].Additionally, the approximation is numerically stable for sectorial operators with the spectrum arbitrary close to the origin.This new result is made possible by extending the idea of resolvent correction, originally introduced for S 1 (t) in [21], to the class of abstract integrands with a scalar-part singularity; see Lemma 3, below.In Section 3.3 and 3.4, we apply the developed approximation of S α (t) to turn solution representation (5) into the exponentially convergent numerical method.A priori error estimates given by Theorems 7 and 12 characterize the method's convergence in terms of the smoothness of u 0 , f ′ (t), values α, ϕ s and the size of T .
The implementation details are provided by Algorithms 1 and 2 which are capable of multilevel parallelism: at the level of solution evaluation for each of the desired t's; at the level of evaluating resolvents for the set of different quadrature points z m and at the level of solving stationary problem that pertains to the resolvent evaluation for the fixed z m .
The mentioned numerical properties are experimentally verified in Example 1 and Example 2, for the homogeneous and inhomogeneous part of the solution, respectively.Both examples consider the negative Laplacian with tunable spectral characteristics in place of A and a conventional eigenfunction-based initial data.Such restriction on the form of initial data permits us to evaluate the space component of solution explicitly, thus removing its contribution to the overall error.The restriction is relaxed in Example 3, which is devoted to the experimental analysis of a fully discretized numerical scheme based on the combination of the developed method with a finitedifference stationary solver.In all three examples, a stable numerical behavior of the approximated solution is observed for α ∈ [0.1, 1.9] and T ≤ 5.

Contour of Integration.
It is well known that the choice of integration contour Γ in ( 7) is critical to the performance of the numerical evaluation of opera-tor function based on the contour integral representation [59,61,21,37].Judicious contour selection involves the analysis of the interplay between the shape of the integration contour, analytical properties of the parametrized integrand and their impact on the performance of a quadrature rule, that is used to evaluate the resulting integral numerically.The authors of [23] showed that the hyperbolic contour is the most convenient choice for the quadrature-based numerical evaluation of abstract functions with sectorial operator argument.Below, we extend their analysis to the case of fractional propagator S α (t).
Let us consider the following hyperbolic contour: (8) with the parameters a 0 , a I , b I that are called shift, first and second semi-axes, respectively.Admissible range of values for these parameters is determined from Theorem 1 that enforces the integration contour Γ = Γ I to encircle the singularities of the integrand in (7) for , because in such case the norm of integrand on Γ will decay faster than the exponential.This observation transforms into the condition a I > 0 for the first semi-axis of hyperbola from (8).The condition b I > 0 for the second semi-axis is induced by the orientation of z α (ξ).We also have to make sure that this curve does not intersect the spectrum of −A.It is worth noting that, for any ϕ ∈ [0, π], the function z α maps the sector Σ(0, ϕ/α) into the sector Σ(0, ϕ).Such mappings can be associated with the Dunford-Cauchy representation of the fractional powers of A [1,23].They are often studied in the theory of fractional resolvent families [35] and associated Cauchy problems [40].
For non-negative a I , b I , the hyperbolic contour Γ I is contained within the region Σ(a 0 − a I , ϕ I ) \ Σ(a 0 , ϕ I ).Here, ϕ I is the angle between positive real semi-axes and asymptotes of Γ I : a 0 + ρe ±iϕI depicted in Figure 1 (b), i.e., tan ϕ I = − bI aI .Consequently, the pair of positive contour parameters a I , b I is admissible if z(ξ) ∈ Σ(0, π−ϕs α ) \ Σ a, π 2 , for some a > a 0 , i.e., Next, we move on to derive exact formulas for a 0 , a I , b I .Let us assume that the chosen set of parameters satisfies (9).The substitution of z(ξ) from ( 8) into (7) yields (10) where z ′ (ξ) = −a I sinh(ξ) + ib I cosh(ξ).The illustration provided by Figure 1 shows that both scalar and operator parts of the parametrized integrand F α (t, ξ), t ∈ [0, T ] remain analytic when ξ is extended into a certain complex neighborhood D of R.
According to the general theory of numerical integration [8], an accuracy of quadrature formula is characterized by a norm of the error-term in the Hardy space H p (D) of functions, defined on the domain D ⊂ C. The shape of D depends on the chosen type of quadrature.For the reasons that are soon to be understood, we approximate integral (10) by the sinc-quadrature formula [51,21]: with the discretization parameter N ∈ N and the step-size h = h(N, F α ).Then, D is formed by an infinite horizontal strip D d of the half-height d: The detailed error analysis of (11) will be presented in Section 3.2.For now, it is sufficient to say that the error of sinc-quadrature decays as O(e −πd/h ) if the integrand is exponentially decaying and belongs to H p (D d ) [51].Thus, in order to achieve a faster convergence rate of quadrature (11), we need to maximize the height of the strip D d , where F α remains analytic, by tuning the parameters of Γ I .
Let us consider the family of curves Γ(ν) = {a 0 −a I cosh (ξ + iv)+ib I sinh (ξ + iv) : ξ ∈ (−∞, ∞)}, which extends the definition of Γ I = Γ(0) to the arguments with nonzero imaginary part ν.Observe, that for a fixed ν > 0, the curve Γ(ν) is also a hyperbola, albeit with different semi-axes a(ν), b(ν): Hence, the mapping w → z(w) transforms D d into the region of complex plane bounded by two hyperbolas z(ξ + id), z(ξ − id), which will be denoted as Γ s and Γ c , correspondingly.We choose parameters a 0 , a I , b I , so that Γ s has the vertex at zero and its asymptotes form the angle φ s ≡ min π, π−ϕs α with R + , as shown in Figure 1 (c).In addition, we require that the asymptotes of Γ c form the angle φ c ∈ π 2 , φ s with R + , which will be called a critical angle; see Figure 1 (b).The above requirements for Γ I , Γ s , Γ c are codified in the following system of equations: which is sufficient to ensure (9) and will lead to the maximal possible d, when φ c = π/2.The system composed from the first two equations is linear with respect to a I , b I ; thus, By that means, the left-hand side of the third equation is transformed as which, after back-substitution, implies tan (φ s − 2d) = tan φ c .Due to the constraints on d, φ c , φ s we are interested only in the following solution of the last equation: For φ c = π/2, and an arbitrary fixed a 0 > 0 we obtain (15) Here, α ∈ (0, 2) is the order of fractional derivative from (1), ϕ s is the spectral angle parameter defined in (2) and a 0 ∈ R + is given.
3. Numerical Method.To begin with the description of the numerical scheme, let us introduce some notation.We rewrite formula (5) in the form Here, u h (t) denotes the solution to the homogeneous part (f (t) ≡ 0) of the given problem (1), (4) and u ih (t) the solution to the inhomogeneous part (u 0 = u 1 ≡ 0): 3.1.Alternative Propagator Representation.We consider the representation of the solution to the homogeneous part u h (t) first.In the seminal paper [21], Gavrylyuk and Makarov showed that the numerical method for S 1 (t) = e −At naively obtained from representation ( 7) is unsuitable for small values of t because its accuracy degrades when t approaches 0. They traced back the root cause of this behavior to the fact that the considered representation of e −At is, formally speaking, divergent at t = 0, which result in the unremovable error of the quadrature-based numerical method for such t.It turns out that propagator representation (7) poses the same adverse feature for any fractional α.One could learn more about its impact on a numerical solution of (1) by analyzing the results of works [14,47].
In order to get around the divergence issue, we propose an alternative formula for S α,1 (t), constructed in the vein of [22,21].It is based on the following proposition, which can be regarded as a generalization of Lemma 3.3 from [44].
Proposition 2. Let A be the sectorial operator satisfying the conditions of Theorem with some constant K > 0 and M defined by (3).
−1 remains analytic and bounded for any z α / ∈ Sp(A) ∪ {0}, so its Neumann series converges unconditionally.Therefore, the last transformation is justified by the fact that R(z, A) and A 1−γ commutes.Target estimate (17) follows directly from the above formula, after we apply inequality (2.30) from [23] with z = z α .
It is worth noting that, if the argument x posses certain spatial regularity x ∈ D(A γ ), γ > 0, estimate (17) guarantees a faster decay of the corrected term's norm The next result defines an improper integral representation for the components of u h (t) from ( 16) and shows the way how the aforementioned correction is incorporated into the formula for S α (t).
Lemma 3. Assume that the given A and α satisfy the conditions of Theorem 1.For any u 0 ∈ D(A γ ), γ > 0 and u 1 ∈ X the operator functions S α (t)u 0 , S α,2 (t)u 1 admit the following representation: where ) and a 0 , a I , b I are specified by (15).Moreover, for arbitrary finite t ≥ 0 integrals in (18) and ( 19) are uniformly convergent.
Proof.Assume Γ is a contour fulfilling the conditions of Theorem 1. Due to the estimate 1+|z| α dz from [54], the integral representation of operator function S α,2 (t) is uniformly convergent for any bounded non-negative t.Formula ( 19) is obtained as a result of the parametrization of ( 7) on the contour Γ I defined by (8).In order to prove (18), we apply the identity e zt /z = 1 to rewrite (7) with β = 1 in the following manner: (20) Proposition 2 and the inequality |e zt | < max {e ℜ(z)t , 1}, t ∈ [0, T ], T > 0 guarantee that the last integral converges uniformly, so we are permitted to parameterize it on the contour Γ = Γ I defined by (8).This yields representation (18).
Lemma 3 is essential for all remaining analysis.Unlike (7) or (10), the new representation of S α,1 (t) by formula ( 18) remains convergent at t = 0.For that matter, it can be used as vehicle for the uniformly convergent numerical method.We shall use the term "corrected propagator representation" as a reference to (18).

Propagator Approximation.
As we can see from Lemma 3, the task of approximating the homogeneous part u h (t) of the mild solution to (1), (4), defined by (16), is reduced to the task of numerically evaluating improper integrals (18) and (19).In this part, we describe how this is achieved using the trapezoidal quadrature rule.Then, we proceed to study the accuracy of the obtained approximation using the theory of sinc-quadrature [51] along with its generalizations to propagator approximations [21].In what follows, symbols S N α,β (t) are used to denote the operators that approximate S α,β (t).
For some h > 0 and N ∈ N, let The functions F α,2 (ξ), F α,2 (ξ), z(ξ) and the parameter γ in the above formulas for propagator approximations have the meaning prescribed by Lemma 3. Similarly to S α (t), we use S N α (t) to denote S N α,1 (t), where appropriate in the sequel.Recall that Γ I is symmetric with respect to the real axis; hence, one can further reduce the number of summands in (21) using the following argument [23].
Remark 4. Let z denote the complex conjugate of z.If the operator A is defined in such a way that R(z, A) + R(z, A) = 2R(ℜz, A), for any z ∈ C \ Sp(A) and x β is defined over the field of real numbers, then and the number of resolvent evaluations for S N α,β (t) in formula ( 21) can be reduced from 2N + 1 to N + 1.
The error of ( 21) admits the following decomposition: where • is the norm of X, as before.This two-term representation of the error is common in the analysis of the accuracy of sinc-quadrature (see Section 3.2 in [51]).
The contribution from the first term is responsible for the replacement of integrals of F α,β (t, ξ) from ( 18) and ( 19) by the infinite series S ∞ α,β (t) of discrete function values F α,β (t, kh).As such, it is commonly called the discretization error.To determine the value of h, one needs to balance it with the contribution from a so-called truncationerror term, that comes second in the formula above.
Let H 1 (D d ) be a family of all functions F : C → X, which are analytic in the strip D d , equipped with the norm The discretization errors of (21) satisfy the estimate [51,21]: Thus, in order to bound this term, one needs to obtain estimates for the H 1 (D d ) norms of the functions F α,β (t, z), β = 1, 2. These are provided by the next lemma.
Lemma 5. Let A be a sectorial operator satisfying the conditions of Theorem 1.For any t ≥ 0, α ∈ (0, 2), x 1 ∈ D(A γ ), x 2 ∈ X, γ > 0 and arbitrary small δ > 0 (21), split out the scalar part in each norm and use bounds ( 17), (3) for the operatordependent parts, correspondingly.As a result, we obtain When the variable ξ is extended from the real line into the strip, the integration hyperbola z(ξ), adopted here from Lemma 3, transforms into the parametric family of hyperbolas with a(ν), b(ν) are defined by (12).Consider the function determines the behavior of η 1 (s, b 0 ) for the values of b 0 that belong to the interval (0, a 0 ), induced by the identity η 1 tanh ξ, a0 . Therefore, for any b 0 , the maximum of η 1 (s, b 0 ) is attained at s = 0, whence The norm F α,1 (t, w)x can be further estimated as ( 26) Here, r(ξ, ν) is a strictly positive bounded function that is defined by the equality 1 + |z(w)| α = r(ξ, ν) cosh α ξ.Solving it for r(ξ, ν) gives us (25).Now, we turn our attention to . By inspecting the derivative we learn that its sign is also determined by a sign of the numerator with only one real root s = 0. Two other roots of η ′ 2 (s, b 0 ) are non-real because the quadratic function The obtained estimates for F α,β (t, w) , β = 1, 2 demonstrate that these norms are exponentially decaying for any t ≥ 0 as ξ → ∞.Consequently, the integral terms from F α,β (t, w) H 1 (D d ) over the vertical parts of ∂D d−δ (ǫ) vanish in the limit ǫ → 0 and we end up with the following expression: After estimating the last integral using the bounds obtained above, we remove the dependence of the integrands on r(ξ, ν) by bounding its value from below with a positive function r 0 (ν) < r(ξ, ν) and, subsequently, evaluate the obtained integrals explicitly.This yields the pair of objective estimates from (23), with constants K 1 = (1 + M )K2 αγ+1 and K 2 = 4M .The lemma is proved.
The proof of this Lemma relies on the established estimates for F α,β (t, w) and is analogous to the proof of the respective part in Theorem 3.1.7from [51].For brevity, we omit it here.

2
, φ α = min{π, π−ϕs α } and φ c ∈ π 2 , φ α , a 0 > 0 are given.The constants C β from (27) and (28) are independent of t, N .Proof.To obtain error bounds for the approximants S N α,β , we depart from the previously established decomposition and then use the results of Lemmas 5 and 6 to estimate the right-hand sides.This transcribes into Having the aim of balancing the order of error contributions from each term inside the brackets, we make two involving exponential functions asymptotically equal as N → ∞.This yields two independent equations with the solutions described by (29).After that, we substitute these expressions into the previously established error estimates to get the following bounds: which are reduced to ( 27) and (28) after denoting γ, 0).It follows from ( 27) and (28), that the value of the contour parameter a 0 can be used to control the error contribution of the factor e a0t .Throughout the rest of this work, we set a 0 = π/6 to make this factor reasonably bounded: e a0t ≤ e 5π/6 ≤ 14.

Numerical Scheme for Homogeneous Part of Solution.
We approximate the homogeneous part u h (t) of the solution to (1), (4) defined by ( 16) using the numerical methods for propagators approximation constructed in Section 3.2.Then, for every fixed N > 0, the approximation u N h (t) to u h (t) is defined as The error of u N h (t) is characterized by the following corollary, which is an immediate consequence of Theorem 7.
Corollary 8. Assume that the operator A, initial values u 0 , u 1 and the fractional order α satisfy the conditions of Theorem 7 with x 1 = u 0 , x 2 = u 1 .For any given N ∈ N, the approximate solution u N h (t), defined by (30) with N 1 = N , N 2 = αγN , converges to the homogeneous solution u h (t) of (1), (4) and the following error bound is valid: The constant C γ is dependent on A, u 0 , and independent of t, N .
It is important to note that the smoothness assumptions for u 0 , enforced by Corollary 8 and Theorem 7, are compatible with the similar assumptions made in [22,23] for the Cauchy problem with the integer order derivative.For a more concise discussion on the impact of the initial data smoothness on the properties of solution to problem (1), ( 4) we direct the reader to [28].
To compute the approximation u N h (t), we suggest to use Algorithm 1 provided below.In this algorithm, the evaluation of each propagator S N β α,β , β = 1, 2 is decoupled into two cycles.The first cycle is responsible for the evaluation of resolvents (z(mh 1 ) α I + A) −1 at the quadrature points of Γ I .This amounts to the solution of 2N β + 1 linear equations that are all mutually independent and hence can be solved in parallel.
Algorithm 1 Algorithm for computing the homogeneous part approximation u N h (t).
INPUT: α, u0, u1, t k , ϕs, N, γ OUTPUT: u N h (t k ) 1: N1 := N ; N2 := αγN1 2: Calculate aI , bI and h1, h2 by ( 15) and ( 29) for each t k do 16: end for 18: end if If A is the discretization of a certain partial differential operator, every resolvent equation from line 4 of the algorithm is actually a system of linear equations.When this is the case, one can leverage additional level of parallelism here, as long as the size and the solution method of the resolvent equation warrant that and the computing environment permits for such possibility.Furthermore, the total number of resolvent evaluations in Algorithm 1 can be reduced all the way down to N 1 + N 2 + 2 if the initial data and A satisfy the conditions from Remark 4.
Given the solution of resolvent equations obtained in the first cycle of Algorithm 1, the second cycle computes the resulting propagator approximation.As we can see from line 8 of Algorithm 1, for every fixed t = t k this step amounts to calculating the weighted sum of resolvents.Hence, its computation is apparently independent of the computed values of solution at different times and can be performed simultaneously.Such feature of the method alone results in the substantial computational advantage over existing sequential time-discretization methods [38,2,10,32], because the average computation cost per u N h (t k ) for any k ∈ {1, . . ., K} is independent on the value of t k ∈ [0, T ] and, unlike in the case of a sequential method, this cost goes down when K grows.Even in the worst-case scenario of K = 1, t k << 1, our method should still remain competitive with the mentioned sequential methods due to its parallelization capability and the uniform exponential convergence.We postpone a more detailed comparison with existing methods until Example 3, where a fully discretized problem is considered.It is important to point out that the described multi-level parallel evaluation strategy is well suited for the multi-node computing architectures, with each node containing the combination of a central processing unit and multiple hardware accelerators, that are ubiquitous nowadays.
For certain realizations of A (c.f.[43]) and large values of N , the resolvent evaluation steps of Algorithm 1 might lead to the numerical instability when |z| is large.This problem can be alleviated by modifying lines 4-5 and 12-13 as described in [44,Eq. (2.18)].Another noticeable feature of the above algorithm is its use of the resolvent evaluations with the complex arguments.This may require additional attention from the implementation point of view if the resolvent is evaluated numerically, for instance using the finite element method software that does not support complex arithmetic.Alternatively, one could deal with the complex resolvent arguments by redefining A via the embedding of its domain into the real space of higher dimensionality.This is always possible, since the resolvent equations from lines 4 and 12 of the algorithm are linear in z.Such modifications are unnecessary for the numerical experiments conducted below; hence, we do not incorporate them into the algorithms for simplicity.
Example 1.Let us consider the standard example problem in which A is a onedimensional Laplacian accompanied by the Dirichlet boundary conditions on [0, L]: where a > 0 is some predefined constant.The initial values u 0 , u 1 are chosen to be the eigenfunctions of the operator A with indices k 0 , k 1 , correspondingly: The exact solution of fractional Cauchy problem (1), (4) with such A, u 0 , u 1 and f (t) = 0 can be represented as follows (see Section 1.3 in [4]): Here, are the eigenvalues of A and H(•) is the Heaviside function, which is added to make the above solution formula valid for all α ∈ (0, 2).
It is easy to verify that for any z ∈ C \ Sp(−A) and k ∈ N, the resolvent R(z α , −A) sin πkx L admits the following representation Hence, all the resolvent evaluations in Algorithm 1 for such u 0 , u 1 can be conducted explicitly.This allows us to focus on analyzing the error contribution from the numerical method for u h (t), given by (30), in the absence of the error associated with the discretization of spatial operator A. The results, presented below, were obtained using the implementation1 of Algorithm 1 developed in Matlab.The standard double precision IEEE 754 arithmetic (and its extension to complex numbers) is used for computations everywhere in this and other examples.The evaluation of E α,β (z) was performed via the contour method from [18], using the accompanied Matlab implementation.The interested reader may also consider alternative methods from [50,41].
The behavior of the exact solution u(t, x) for the simplest case c = 1, L = 1 is shown in Figure 2, where it is plotted as a function of time for different values of α at x = 0.5.In the sub-parabolic case α ≤ 1, the solution remains positive for positive u 0 and it is monotonously decaying toward zero as t → T .More specifically, for small α (see graphs for α = 0.1, 0.3 in the left plot of Figure 2), |u(t)| has a fast initial decay which tends to be getting slower as t progresses.This effect becomes less noticeable as α goes toward 1, at which point u(t, x) = E 1,1 (−π 2 t) sin πx = e −π 2 t sin πx.In the sub-hyperbolic case α > 1 (see the right plot of Figure 2), the solution exhibits more complex behavior.It is akin to the damped oscillations with the initial amplitude equal to u 0 (E α,2 (0) = 0 by definition) and the amount of damping that decays as α approaches 2. To quantify the error of the numerical solution to Cauchy problem (1), ( 4), ( 32), (33), calculated using Algorithm 1, we define The behavior of E h (t, x) as a function of t for fixed x = 0.5, L = 1, k 0 = 1, k 1 = 4, a = 1, ϕ s = π/60, γ = 1 and different values of α, N is illustrated in Figure 3, using a semi-logarithmic scale for the plots.All graphs clearly illustrate the decay of E h (t, x) on the whole time interval as N increases.In the consequence of Theorem 7, the error E h (t, x) also depends on how much fractional order α deviates from 1.For α ≤ 1, this happens due to the direct presence of α in error bound (27) as a factor.For α > 1, this is explained by the contribution of α to the factor c = √ 2πd, from the same error bound, via (14).
Aside of that, for α ≤ 1, we witness a sharp drop of E h (t, x) in the vicinity t = 0 (see Figure 3 (a)-(c)), which does not seem to be predicted by the error bound.This behavior is attributed to the rather pessimistic estimate |e z(ξ)t | ≤ e (a0−a(ν) cosh ξ)t ≤ e a0t , which was used to account for the contribution of the t-dependent term into both truncation and discretization errors of S N α,1 (t)u 0 (see the proof of Lemmas 5 and 6 above).Similar phenomenon was observed in [44], where a related fractional problem was considered.The influence of factor e a0t became more evident for larger t, as seen from the graphs of Figure 3 (d)-(f ).For fixed N , the amplitude of error oscillations increases when t approaches 5 but remains approximately equal α-wise (visually larger amplitude oscillations for smaller α in Figure 3 are caused by the semi-log nature of the plots).This observation supports the theoretical claim from Theorem 7 that the growth of E h (t, x) in time is not influenced by α or d.In order to analyze the error dependency on the position of Sp(A), we evaluate the sup-norm error err h (N ) for several values of diffusivity constant a = 10 −5 , 0.1, 1, 10 from (32), and a range of α values (see Figure 4).The magnitude of the quantity ρ s = inf z∈Sp(A) ℜz = aπ 2 corresponding to a = 10 −5 in Figure 4 (a) is characteristic for problems with a singularly perturbed A [30] and, in particular, advection (convection)dominated flows [48].Our prior experiments suggest that existing numerical methods [21,42,44,55], with the integration contour which lies entirely in the same halfplane as Sp(A), face certain difficulties in handling problems with such small ρ s .Those are caused by the implicit rescaling of z(ξ) needed to fit z(D d ) between Sp(A) and the origin.In contrast, the current method does not experience any accuracy degradation related to ρ s → 0, because the integration contour Γ I encircles Sp(A) ∪ {0}.In fact, Figure 4 shows that the sup-norm error decays exponentially with the order proportional to √ αN as prescribed by (31), for all analyzed values of a.The convergence results of our method for α ≤ 1 are similar to those obtained in [6] for the specific case of (1), when A comes from the viscoelastic beam model.3.4.Numerical Scheme for the Inhomogeneous Part.In this part, we apply the propagator approximation method from Section 3.2 to obtain an efficient numerical algorithm for the inhomogeneous part u ih (t) of the mild solution to (1), (4), defined by (16).This formula combines the action of S α (t) on a certain vector from X with the subsequent action of the integral operator.The numerical evaluation of such composition amounts to the reevaluation of S α (s) at each quadrature point {s k } K k=1 , needed to approximate the outer integral.As we have learned from the properties of the numerical method developed in Section 3.2, this is not a problem for the first term of u ih (t), where the argument f (0) is fixed, because only 2N + 1 parallel resolvent evaluations are needed.For the second term, however, the numerical evaluation of S α (t − s)J α f ′ (s) for every new value of t requires the reevaluation of resolvents for the entire set of new quadrature points on Γ I .This leads to the solution of up to (2N +1)K additional stationary problems and may require additional storage and inter-process communication, when the parallel computing model is used for evaluation.To reduce the number of required resolvent evaluations, we take advantage of the fact that the operator-dependent part of S α (t) in representation ( 18) is itself a linear operator on X; hence, it can be interchanged with the integral operator acting in t only Here, we used formula (20) under the assumption that J α f ′ (s) ∈ D(A γ ), γ > 0 and, then, relied upon the uniform convergence of the corrected representation of S α (t − s) with respect to s ∈ [0, t], that had been established earlier.As we can see from the newly obtained representation, now the evaluation of the time-dependent part is performed in the resolvent's argument.This reduces the number of parallel resolvent evaluations per t to 2N +1 for this term.The last representation permits us to rewrite the inhomogeneous part of solution u ih (t) in the form (34) t 0 e z(ξ)(t−s) J α f ′ (s) ds dξ.
Next, we address another ingredient essential to the numerical evaluation of (34), which is an efficient quadrature method for the Riemann-Liouville integral J α v(t) defined by (6).While evaluating this integral numerically, it is important to select the quadrature rule that, on the one hand, can handle the endpoint singularity appearing in the integrand when α < 1 and, on the other hand, is able to provide exponentially convergent approximation.Among existing quadrature rules, only sincquadrature on a finite interval satisfies two mentioned properties simultaneously (see [51]).We construct a version of such quadrature rule by transforming J α into the integral over (−∞, ∞) and then applying the chosen sinc-quadrature formula.Let s = te p / (1 + e p ), then The reader may note that the singularity (t − s) α−1 from original definition (6) of J α v(t) is no longer present in the last integral and the new kernel of J α v(t) decays exponentially as p → ∞.More precisely, there exist a constant c > 0, such that (35) e p (1 + e p ) α+1 = e − p α+1 + e αp α+1
Our intent here is to approximate J α v(t) by the time-dependent operator J N α v(t), that takes into account the difference in a speed of kernel's decay as p → ±∞, illustrated by the above bound.To introduce the approximation The constants a, b > 0 will be referred to as the decay orders (or the decay order if a = b).
d , which is a proper subset of D 2 d as long as t < T .Consequently, if the assumptions regarding v(z) are fulfilled, the integrand from (38) belongs to the class of functions L 1,α (D d ) for any t ∈ (0, T ].Then, the results regarding the convergence of (37) to (38), as well as the form of (37) itself, and the error estimate stated in (39) follow from [51,Thm. 4.2.6].
Remark 10.The results of Proposition 9 remain valid if, instead of the boundedness of v(z), we assume that the integrand from (38) belongs to the class L a,b (D d ).
The presence of factor t α in error estimate (39) makes it possible to use fewer terms in (37) as t decreases, if the end goal is to reach the prescribed accuracy uniformly in t.Let us assume that the desired accuracy is achieved for some t 0 by setting N = N 0 , then for t ∈ (0, t 0 ): .
After solving this equation for N , we obtain ( 40) Formula (40) becomes instrumental in the situations where one needs to numerically evaluate J N α v(t) for a range of t-values.This is the case of the inhomogeneous part of solution representation given by (34), whose terms contain the integrals of J α f ′ (s).
Before addressing the question on how to numerically evaluate (34), we would like to consider a prerequisite problem on how to quantify the contribution of the error in the argument v(t) of J N α v(t) to the overall error of approximation to J α v(t).Corollary 11.Assume that functions v, v satisfy the assumptions of Proposition 9.If v(t) − v(t) ≤ κ, for all t ∈ [0, T ], then the error of approximation J N α v(t), satisfies the bound Proof.We rewrite (41) as The first term of this error decomposition is estimated by (39), so we focus on the second term (42) ≥ h was used to cancel out h in the last estimation step.The combination of ( 39), ( 42) and ( 43) completes the proof.
With all the necessary results in place, now we move on to construct the approximation to u ih (t).To achieve that, we apply approximations ( 21), (37) and discretize the remaining time-dependent integrals in a similar fashion as the Riemann-Liouville integral J α (see Proposition 9).The rationale for such integral discretizations will become apparent when we analyze the error below.Meanwhile, let us introduce the proposed approximation u N ih (t) of the inhomogeneous solution u ih (t) from ( 34): (44) where F α,1 (ξ), z(ξ) are defined in Lemma 3 and Theorem 12. Let A be a sectorial operator satisfying the assumptions of Theorem 7. If the function f (t) from (1) admits the analytic extension to the "eye-shaped" domain D 2 d , d ∈ (0, π/2) and ( 45) with some χ > 0, then for any α ∈ (0, 2) and t ∈ [0, T ], the approximation u N ih (t) from (44) converges to the inhomogeneous part u ih (t) of the mild solution to (1), ( 4), defined by (16).Moreover, for any fixed N ∈ N, the following error bound is valid: with c = √ 2πd, provided that the values of N i and h, h i in (44) are determined by the following formulas Here, d and independent of t, N .Proof.We analyze error u ih (t) − u N ih (t) of ( 44) in a term-by-term manner.The first error term is estimated via Corollary 11 and Theorem 7, applied in succession: The error bound for the second term can be decomposed as with η 1 being the quadrature error of the outer integral: stated here after the substitution s = ψ(p) is performed therein, whereas η 2 is the compound error of the discretized Riemann-Liouville operators: It is worth noting that the last series is a specific version of the one from (42), with α = 1.Thus, formula (43) along with the bound from Proposition 9, warranted by the analyticity assumptions on f ′ (z), yield Let us return to η 1 .The aforementioned analyticity of f ′ (z) induces the uniform convergence of the integral for J α f ′ (z) in formula (38)  d permits us to conclude that ψ ′ (p)J α f ′ (tψ(p)) is analytic for p ∈ D d .Due to the form of ψ ′ (p), it is also exponentially decaying as |p| → ∞ in D d ⊆ C, with the decay order 1.Hence, the error of sinc-quadrature η 1 admits the bound established in Theorem 4.2.6 from [51]: We treat the third term from (44) in a similar way as the second term, albeit this time the error decomposition is conducted after the application of (20) in reverse: Here, the quantity η 4 is used to denote the quadrature error for the outer integral: The upper bound for the last integrand is determined by the properties of the norm dξ , which can be estimated using inequality (26).Indeed, setting x 1 = J α f ′ (s) in ( 21) reveals that the above integrand equals to the expression for F α,1 (t − s, ξ).As such, it admits the estimate Here, we used the relation stemming from the equivalence of definitions (38) and (6).The previous chain of estimates leads us to the following bound: (51) Therefore, the integrand from ( 49) is analytic in D d ′ , with the value of d ′ = φs 2 − π 4 , which is equal to d from (15).Furthermore, the integrand's norm is exponentially decaying as a function of p ∈ D d , as a consequence of (51), the boundedness of sup Next, we use (27) in conjunction with (43) and (50) to estimate the propagator approximation error η 3 The remaining summand η 5 represent the effect of discretized Riemann-Liouville operators on the error of the third approximation term in (44).We estimate it as Assumption (45) enables us to estimate the last norm via Proposition 9: This decouples the inner and outer series in the above estimate for η 5 .Thus, with r m = inf p∈R r 0 (p) being strictly greater than zero, due to (25).By combining two previously obtained bounds with (43), we arrive at , here, is independent of N 5 .The bounds derived for the quantities η i , i = 1, . . ., 5 show that they all are exponentially decaying as N i → ∞.We make these error bounds asymptotically equal to the error of the first term from (44), that is decaying on the order of e −c √ αχN , provided that εN 0 = αχN in the error estimate from the beginning of the proof.The resulting equations for N i are as follows: where ε = min {1, α}, as per Proposition 9.The solution of these equations gives us (47).By collecting the derived error bounds for the terms of u N ih (t), we end up with This bound is reduced to (46) by absorbing the individual constants into C χ,f , while retaining the asymptotic behavior with respect to α, χ and t.The derived bound also proves the convergence of approximation (44) to (34) and, therefore, to the original definition for the inhomogeneous part of the solution given by ( 16).
Theorem 12 demonstrates that the proposed numerical method to approximate u ih (t) inherits essential properties of the numerical method for propagator approximation, it is based upon.Firstly, the constructed approximation u N ih (t) is exponentially convergent on the whole interval t ∈ [0, T ].Secondly, bound (46) exhibits, similar to, (27) dependence on the fractional order α and the argument smoothness parameter χ.Thirdly, just like (34), formula (44) permits for an independent evaluation of resolvents R(z α , −A) for the different values z ∈ Γ I .Moreover, the presence of factor t α in (46) guaranties that the approximation u N ih (t) matches the asymptotic behavior of the inhomogeneous part u ih (t), when t → 0+ [54].
Algorithm 2 Algorithm for computing the inhomogeneous part approximation u N ih (t).
The following example is aimed to numerically verify the quality of approximation (44) to the inhomogenous part of solution given by (34) or (16).
Example 2. Let A be defined as in Example 1. Furthermore, let f from (1) be a product of the eigenfunction of A and the polynomial: where m and c i , i = 0, . . ., m are given.For such f , we have The inhomogeneous part of the solution to (1), ( 4) takes the form u(t, x) =c ′ 0 sin which is derived using the fractional propagator representation from Example 1.The integrals from the above formula for u(t) cannot be evaluated explicitly for arbitrary α.Thus, we rely upon the numerical evaluation of u(t) via exponentially convergent quadrature formulas (38) and [51, Theorem 4.2.6] with discretization parameters N J and N I , correspondingly.The analysis conducted in the Proof of Proposition 9 suggests to set N J = N I / min{1, α}.This leaves us with only one discretization parameter N I , which has to be chosen large enough for the error of the approximated u(t) to be negligible with respect to the error of the numerical solution u N ih (t).The latter one is obtained by Algorithm 2 for the data specified in (32), ( 33) and (52), using the explicit resolvent evaluation formula from Example 1 and the software implementation mentioned there.We fix m = 1, c 0 = 1, c 1 = 1, k 0 = 1, k 1 = 4, L = 1 in (52) and, after conducting several numerical experiments, settle with N I = 256.The resulting behavior of u(t) is visualized in Figure 5. the exact solution to problem (1), ( 4) is defined as  x 2 (x − 1) , so the right-hand side of (1) takes the form ( 54) x 2 (x − 1) − 12x 2 + 6x(b + 1) − 2b.
Such f (t) permits us to study one important practical aspect of the developed solution method.Namely, what happens to the accuracy of a fully-discretized solution when f ′ (t) does not formally belong to the domain of A, but the discretization A satisfies A f ′ (t) < ∞?
Let A be m × m matrix obtained by a second-order finite-difference discretization of operator (32) on a grid ∆ d = {(i − 1)/(m − 1)} m i=1 .Then, the discretized righthand side f ′ (t) ∈ (R m , ∞ ) is defined by the projection of f ′ (t) onto ∆ d : f ′ (t) = (f ′ (t, 0), f ′ (t, x 2 ), . . ., f ′ (t, L)) T .We  As we can see, the function f ′ (t), t > 0 does not satisfy the boundary conditions from (32); hence, f ′ (t) / ∈ D(A).Furthermore, when α > 1, this function posses an (5) using the combination of efficient methods for the contour evaluation of the propagators S α,β (t), β = 1, 2 and tailored quadrature rules for the discretization of the Riemann-Liouiville and convolution integral operators from (5).As a result, the numerical evaluation of ( 5) is reduced to the solution of a sequence of independent linear stationary problems.The accuracy estimates established by Theorems 7 and 12 remain valid uniformly in time for the entire range α ∈ (0, 2), under the moderate smoothness assumptions u 0 ∈ D(A γ ), f ′ (z) ∈ D(A χ ), with some γ, χ > 0 and all z ∈ D 2 d , defined by (36).These results recover the previously existing error estimates for parabolic problems [22,21], when α is set to 1 and T < ∞.All the theoretical results are verified experimentally.This includes the results from Corollary 8 and Theorem 12 regarding the approximation of homogeneous and inhomogeneous parts of the solution, which are experimentally considered in Example 1 and 2. Here, we put extra effort to demonstrate that the constructed solution approximation is numerically stable for α ∈ [0.1, 1.9] and practically capable of handling operators A with a broad range of spectral characteristics.It encompasses the class of so-called singularly perturbed operators, that are modeled in Example 1 by the Laplacian with a very small distance between the Sp(A) and the origin (see Figure 4).Additionally, in Example 3, we considered a fully discretized solution scheme for (1), (4) to practically verify the robustness of the constructed approximation toward errors caused by the discretization in space.Naturally, the mentioned benefits of the developed method come with some limitations.Among such, we mention the required analyticity of f (z) in the complex neighborhood D 2 d of time interval (0, T ).On the one hand, this is a considerably stronger assumption on the problem's right-hand side than the assumption f ∈ W 1,1 ([0, T ], X), imposed by the solution existence result from Theorem 1.On the other hand, such analyticity assumptions are typical for the theory of exponentially convergent quadrature [8,51].Moreover, the quadrature rule chosen in Theorem 12 permits for a practically realistic situation when f ′ (t) has an integrable singularity at t = 0.The ability of the method to handle such class of f was experimentally demonstrated in Example 3. Another possible limitation of the current method is its practical viability only for moderate T ≤ 20.Nonetheless, existing numerical evidence suggests that the long-term stability of the method could be improved by some nonessential modifications.Larger values of T ≈ 200 are necessary for certain parameter identification problems [62], which, along with the mentioned nonlinear and nonlocal extensions of the given problem, are going to be considered in the future works.

Figure 1 .
Figure 1.Schematic plot of the complex neighborhood D ≡ D d of R where the parametrized integrand Fα(t, ξ) remains analytic and exponentially decaying for any t ∈ [0, T ] (a) along with the image of D d under the mapping v → z(v) defined by Γ I (b) and the region z α (v), v ∈ D d (c).The "forbidden" regions of complex plane are indicated by "beige" color).(α = 1.3, ρs = π, ϕs = π/6).

Figure 4 .
Figure 4. Sup-norm error err h of the approximate solution u N h to problem (1), (4) with f (t) = 0, A, u 0 , u 1 being defined by (32), (33) and L = 1, k 0 = 1, k 1 = 4. Graphs from sublots correspond to the different values of diffusivity constant: (a) a = 1 × 10 −5 ; (b) a = 0.1; (c) a = 1; (d) a = 10; J N α v(t) properly, let us to recall the following definition [51, Definition 3.1.5].The function f is said to belong to the class L a,b (D d ) if it is analytic in D d and there exist a constant c > 0 such that for all z ∈ D d :

Proof.
When t = T , the function te z /(1 + e z ) maps the infinite horizontal strip D d of half-height d into the "eye-shaped" region D 2 d (see [51, Example 1.7.5])around the interval [0, T ].For smaller values of t, it maps D d into the region tD 2 with respect to z ∈ D 2 d .Furthermore, for an arbitrary value of p ∈ (−∞, ∞), the function zψ(p) from (38) maps the convex region D 2 d defined by (36), onto itself.By repeating the argument from the Proof of Proposition 9, these two facts and the relation D d ψ − → D 2

z∈D 2 d
A χ f ′ (z) and the convergence of the integral, established before.Hence, similarly to η 1 , the bound from Theorem 4.2.6 of[51] yields
: end for 10: for each t k do