Solving time-fractional differential equation via rational approximation
SSOLVING TIME-FRACTIONAL DIFFERENTIAL EQUATION VIARATIONAL APPROXIMATION ∗ USTIM KHRISTENKO † AND
BARBARA WOHLMUTH † Abstract.
Fractional differential equations (FDEs) describe subdiffusion behavior of dynamicalsystems. Its non-local structure requires to take into account the whole evolution history during thetime integration, which then possibly causes additional memory use to store the history, growingin time. An alternative to a quadrature of the history integral is to approximate the fractionalkernel with the sum of exponentials, which is equivalent to consider the FDE solution as a sumof solutions to a system of ODEs. One possibility to construct this system is to approximate theLaplace spectrum of the fractional kernel with a rational function. In this paper, we use the adaptiveAntoulas–Anderson (AAA) algorithm for the rational approximation of the kernel spectrum whichyields only a small number of real valued poles. We propose a numerical scheme based on thisidea and study its stability and convergence properties. Moreover, we apply the algorithm to atime-fractional Cahn-Hilliard problem.
Key words. time-fractional differential equations, rational approximation, AAA algorithm
AMS subject classifications.
1. Introduction.
Fractional differential equations became a more and morepopular research direction in the last decade, proposing new mathematical tools fordescribing complicated physical phenomena in case where the local differential mod-els face their limitations. In particular, the time-fractional differential equations canbe found in application such as, e.g., tumor growth models and in visco-elasticity.In the general case, the analytical solution of such problems is not available. Thus,the development and study of computational methods and algorithms for FDEs is ofhigh interest and is the object of many research works. The numerical solution ofa FDE is much more complicated and more computationally expensive than for theclassical integer-order problem. In fact, the fractional ordinary differential equations(FODEs) are integral equations involving convolution with singular kernels. Thus,one of the related challenges is the non-local structure of the operators, which takesinto account the whole evolution history during the time integration. This possiblycauses additional memory use to store the history, which grows in time.The most common strategy is directly based on the quadrature schemes for theconvolution integral. The classical one is the so-called L1 scheme based on a fi-nite difference formula [49, 33]. An important class constitutes the fractional lin-ear multi-step methods (FLMM). Pioneering works in this direction has been doneby Lubich [40, 41, 42]. The FLMM include methods of Adams–Moulton/Bashforthtype [20, 21, 60, 62]. Though these methods are conceptually simple and have a highconvergence order, they are shown to experience troubles for certain values of thefractional power [19]. Moreover, the quadrature-based approach is mostly affected bythe curse of non-locality. Precisely, the integration of N time steps has algorithmiccomplexity of order O ( N ) and requires O ( N ) solutions to store. The non-localityproblem can be tackled with various memory-saving techniques including short mem- ∗ Submitted February 11, 2021.
Funding:
This work was funded by the German Research Foundation by grants WO671/11-1and the European Union’s Horizon 2020 research and innovation programme under grant agreementNo 800898. † Lehrstuhl f¨ur Numerische Mathematik, Fakult¨at f¨ur Mathematik (M2), Technische Universit¨atM¨unchen, Garching bei M¨unchen ([email protected], [email protected])1 a r X i v : . [ m a t h . NA ] F e b U. KHRISTENKO, AND B. WOHLMUTH ory principle [50, 14], logarithmic grid [29, 22] and parallel computations [18] for thehistory integral. Consequently, the complexity is possibly reduced to O ( N log N )with memory size O (log N ).Another strategy for the numerical solution of FODEs is based on the approxima-tion of the integral kernel. The crucial role in the approach plays an integral represen-tation of the kernel. The first steps in this direction have been done in [43, 52, 39]. Inthe so-called kernel compression techniques, the kernel is approximated with a sum-of-exponentials [10, 45], leading to a family of ODEs, which can be usually solved inparallel. Therefore, this approach can be seen as a decomposition of the FODE solu-tion into the sum of different modes, each governed by a corresponding local evolutionlaw. This phenomenon clearly illustrates the nature of FODEs, where the computa-tional complexity can be interpreted in terms of a hidden extra-dimension, which isalso reflected in the methods for fractional partial differential equations (compare,e.g., with [12, 8, 11, 56, 31]). The advantage of this methodology is that the historyterm is approximated with a linear combination of a fixed number m of modes, andthus, no special treatment of the memory is needed with the constant additional stor-age of the size O ( m ) and the computational complexity O ( N m ). In many works, thesum-of-exponentials is obtained using a quadrature for the integral representation ofthe kernel, see [37, 46, 32, 61, 2, 7]. An alternative approach for the approximation ofthe kernel consists in polynomial (multi-pole) interpolation of its spectrum [3]. Ourapproach belongs to this class. For the detailed overview of the existing numericalmethods for FODE, we refer in particular to [6, 16, 23].In this work, we propose a new approach based on the approximation of theLaplace spectrum of the convolution kernel with a rational function. We use theadaptive Antoulas–Anderson (AAA) algorithm [48], leading to a multi-pole approxi-mation of the spectrum with real non-negative poles, which transforms to the sum-of-exponentials kernel with an additional singular term. In our numerical settings, wesee that a small number of modes ( m ∼
20) is sufficient to achieve excellent accuracyeven in complex non-linear PDE test cases. Discretizing the obtained system of ODEsfor the modes, we propose new numerical schemes of the Implicit Euler and Crank-Nicolson type. Though we do not discretize explicitly the local integration term, thelocal term coefficients in the numerical schemes, being defined by the rational approx-imation coefficients, naturally reproduce the fractional Adams–Moulton coefficientsarising in linear multi-step methods.The paper is structured as follows. In Section 2, we bring in the necessary defini-tions and formulations and provide some preliminary results. In Section 3, we discussthe rational approximation of the spectrum of the fractional kernel. In Section 4,we analyze the associated approximation error. In Section 5, we suggest numericalschemes based on the discretization of the modal ODE system and discuss their ac-curacy and stability. Finally, in Section 6, we illustrate the performance of the newlyintroduced schemes numerically. In particular, in Subsection 6.1, we consider a lin-ear time-fractional heat equation in 1D with known analytical solution to show theconvergence rate. In the second numerical example in Subsection 6.2, the proposedscheme is applied to a non-linear time-fractional Cahn-Hilliard equation in 2D.
2. Preliminaries.
Let us first introduce some basic definitions of the fractionalderivative and the fractional integral. For a more detailed introduction to the theoryof fractional differential equations and fractional calculus, the reader is referred to [6,
OLVING TIME-FDE VIA RATIONAL APPROXIMATION
34, 17, 36, 44, 50, 47, 51]. The Riemann–Liouville fractional integral is defined as(2.1) I α u ( t ) = K α ∗ u ( t ) = (cid:90) t K α ( t − s ) u ( s ) d s, where the kernel is defined by K α ( t ) = t α − / Γ( α ) , and Γ( α ) is the gamma function. Then, the Riemann–Liouville fractional derivativeis given by ∂ αt u ( t ) = ∂ t (cid:2) I − α u ( t ) (cid:3) . If u ( t ) is sufficiently smooth, then we have ∂ αt [ u ( t ) − u (0)] = I − α [ ∂ t u ( t )] . The right-hand side is the classical Caputo fractional derivative. The formulationon the left hand side, which expresses the Caputo fractional derivative in terms ofRiemann–Liouville fractional derivative, has the advantage that it requires less regu-larity of u ( t ) than the classical definition.Let H be a Hilbert space with inner product (cid:104)· , ·(cid:105) and associated norm (cid:107)·(cid:107) H . Letus consider the associated Bochner space B = L ([0 , T ]; H ) with the norm defined by (cid:107) u (cid:107) B = (cid:82) T (cid:107) u ( t ) (cid:107) H d t . In the current work, we focus on the solution u ∈ B of thefollowing non-linear fractional Cauchy problem: ∂ αt u ( t ) = F [ u ]( t ) , (2.2) u (0) = u , where α ∈ (0; 1], and F : B → ˆ B , with ˆ B = L ([0 , T ]; H (cid:48) ), is a continuous possiblynon-linear operator. The solution of Equation (2.2) can be formally written in termsof (2.1) as(2.3) u ( t ) = u + K α ∗ F [ u ] . Let us denote by L the Laplace transform operator. Then, for any function f ∈ B ,we denote its Laplace transform by ˆ f ( s ) := L { f } ( s ) = (cid:82) T f ( t ) e − st d t , s ∈ C . Inparticular, the Laplace transform of (2.3) writes L { u − u } ( s ) = s − α L { F [ u ] } ( s ) . We are interested in the construction of an approximation of u ( t ) that we intro-duce through(2.4) ˜ u ( t ) = u + ˜ K α ∗ F [˜ u ] , where ˜ K α ( t ) is an approximation of the kernel K α ( t ).
3. Rational approximation of the kernel.
We want to construct a ker-nel ˜ K α ( t ) such that (2.4) yields a good approximation to the solution (2.3) of (2.2),which can be numerically found at a lower computational cost. Let us consider the U. KHRISTENKO, AND B. WOHLMUTH
Laplace transforms of the kernels K α ( t ) and ˜ K α ( t ). We look for ˆ˜ K α ( z ), z ∈ C , inthe form of a rational function, more precisely as a ratio of two polynomials. Thatis, we want to construct ˆ˜ K α ( z ) as a rational approximation of ˆ K α ( z ) = z − α . Since α ∈ (0 , K α ( z ) as polynomials ofthe same degree m . Under the assumption that the polynomials have only simpleroots, the partial fractions decomposition of ˆ˜ K α ( z ) writes:(3.1) ˆ˜ K α ( z ) = m (cid:88) k =1 c k z + d k + c ∞ , with c k ≥ d k ≥
0. Then, the kernel ˜ K α ( t ) is given by a typical sum-of-exponentialsand a singular term: ˜ K α ( t ) = m (cid:88) k =1 c k e − d k t + c ∞ δ ( t ) , where δ ( t ) denotes the Dirac δ -distribution. Hence, the approximation ˜ u ( t ), definedby (2.4), reads as(3.2) ˜ u ( t ) = u + m (cid:88) k =1 u k ( t ) + u ∞ ( t ) , where the modes u k ( t ) are given by u k ( t ) = c k (cid:90) t e − d k ( t − ξ ) F [˜ u ]( ξ ) d ξ, k = 1 , . . . , m,u ∞ ( t ) = c ∞ F [˜ u ]( t ) . That is, the modes u k ( t ), k = 1 , . . . , m , satisfy the following ordinary differentialequation: ∂ t u k ( t ) + d k u k ( t ) − c k F [˜ u ] = 0 , (3.3) u k (0) = 0 . Remark that the uni-modal ( m = 1, c ∞ = 0) setting with c = 1 and d = 0 yieldsthe trivial case α = 1. On the other hand, the situation with only the ”infinity” mode c ∞ = 1 ( m = 0) corresponds to the case α = 0.Note that this expansion (3.2) is similar to the representations in [37, 59, 15],however we obtain the weights and exponents from the rational approximation of thekernel spectrum and not from the integral quadrature. The expansion coefficients c k and d k can be computed using various rational approximation algorithms, e.g., Pad´eapproximant [5], Best Uniform Rational Approximation [53], barycentric rational in-terpolation [9], etc.. For more information on the rational interpolation methods, werefer to [55, 13]. In this work, we employ the adaptive Antoulas–Anderson (AAA)algorithm [48, Fig. 4.1], which in our test cases has demonstrated particular efficiencyand robustness.In order to compute the coefficients c k and the poles − d k in the expansion (3.1),we implement the adaptive Antoulas–Anderson (AAA) algorithm [48, Fig. 4.1]. Letus consider the target function f ( z ) = z α on the interval [ h, T ]. Without loss of OLVING TIME-FDE VIA RATIONAL APPROXIMATION T = 1. Applying the AAA algorithm with the support pointsfrom [ h, T ], using a logarithmic grid, we approximate the target function with a ratioof two polynomials of degree m : P ( z ) = m (cid:88) k =0 p k z k and Q ( z ) = m (cid:88) k =0 q k z k . Then, the rational function r ( z ) = P ( z − ) /Q ( z − ) = (cid:80) mk =0 p k z m − k / (cid:80) mk =0 q k z m − k approximates f ( z − ) = z − α on the interval [ T , h ]. Hence, the partial fractions de-composition of ˆ˜ K α ( z ) = r ( z ) yields the required multi-pole form (3.1).Since the Laplace transform is a compact operator from L ([0 , ∞ ]) to L ([0 , ∞ ]),the problem of its inversion on the real line is, generally speaking, ill-posed [26].Nevertheless, we observe that barycentric rational approximation on the real lineleads to approximation on the vertical line in C . However, the error on the verticalline is amplified with respect to the initial tolerance. In return, approximation onthe real line provides real polynomial coefficients. Moreover, the coefficients c k ≥ d k ≥ E ra := (cid:13)(cid:13)(cid:13) K α − ˜ K Expα (cid:13)(cid:13)(cid:13) L ([ h,T ]) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) h (cid:16) K α ( s ) − ˜ K Expα ( s ) (cid:17) d s − c ∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , where we denote the sum-of-exponentials part of the kernel approximant by˜ K Expα := m (cid:88) k =1 c k e − d k t . In Figure 1 on the left, we show the convergence of this error with respect to thetolerance of the AAA algorithm for T = 1 and h = 10 − . The integral in (3.4)is computed using scipy package [34]. Let us note that the above algorithm in ourimplementation faced its limitations with further decreasing both α and tolerance,when spurious poles appear.In Figure 1 on the right, the number of modes m with the AAA-tolerance 10 − isshown as function of the fractional power α for T = 1 and different h . We can remark asmall number of modes in comparison with other multi-pole approximations, e.g., [3].We can also observe that the number of modes decreases at the limits of the interval(0 , α = 1, we obtain m = 1 with the single mode c = 1 and d = 0. For α = 0, we have m = 0, when only the ”infinity” mode c ∞ = 1 remains.Note that the coefficient c ∞ takes values between 0 and 1, where the extremitiescorrespond respectively to α = 1 and 0.
4. Error analysis.
In this section, we estimate the error E a ( T ) = (cid:107) u − ˜ u (cid:107) B (Theorem 4.3 below). To this end, we will need the following auxiliary lemmas. Lemma
Let g ∈ L ([0 , T ]) and let f ( t ) ∈ L ([0 , T ]) be uniformly Lipschitzcontinuous on [0 , T ] , i.e., there exists C > such that for all t , t ∈ [0 , T ] it holds | f ( t ) − f ( t ) | ≤ C | t − t | . And let c > be fixed. Then, for all t ∈ [ h, T ] , the U. KHRISTENKO, AND B. WOHLMUTH − − − − − − − − − − − − − AAA-tolerance E r a α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . . . . . α N u m b e r o f m o d e s , m h = 10 − h = 10 − h = 10 − h = 10 − h = 10 − Fig. 1:
Left:
Dependence on the AAA algorithm tolerance of the rational approxi-mation error (3.4) for T = 1 and h = 10 − . Right:
Dependence on α of the numberof modes m with the AAA-tolerance 10 − , T = 1 and different values of h . following inequality holds | g ∗ f ( t ) − c f ( t ) | ≤ C h (cid:107) g (cid:107) L ([0 ,h ]) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h (cid:90) g ( s ) d s − c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) · | f ( t ) | + g [ h,T ] ∗ f ( t ) , Proof.
Adding 0 = c f ( t ) (cid:82) h g ( s ) d s − c f ( t ) (cid:82) h g ( s ) d s , we can write t (cid:90) t − h g ( t − s ) f ( s ) d s − cf ( t ) = t (cid:90) t − h g ( t − s ) [ f ( s ) − f ( t )] d s + h (cid:90) g ( s ) d s − c f ( t ) . Note that by our continuity assumption, it holds for the first term that (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) tt − h g ( t − s ) [ f ( s ) − f ( t )] d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ C h (cid:107) g (cid:107) L ([0 ,h ]) . Hence follows the statement.
Lemma
Let v ∈ B , ε ∈ L ([0 , T ]) and a constant C > such that (4.1) (cid:107) v ( t ) (cid:107) H ≤ | ε ( t ) | + C ( K α ∗ (cid:107) v (cid:107) H ) ( t ) . Then, (4.2) (cid:107) v (cid:107) B ≤ (cid:107) ε (cid:107) L ([0 ,T ]) · E α, [ C T α ] , where E α,β [ x ] denotes Mittag-Leffler function [35]: (4.3) E α,β [ x ] = ∞ (cid:88) k =0 x k Γ( αk + β ) . OLVING TIME-FDE VIA RATIONAL APPROXIMATION Proof.
According to [58, Theorem 1], it follows from (4.1) (cid:107) v ( t ) (cid:107) H ≤ | ε ( t ) | + (cid:90) t ∞ (cid:88) k =1 C k ( t − s ) αk − Γ( αk ) | ε ( s ) | d s. Hence, computing the L -norm on [0 , T ] and applying Young’s convolution inequalityto the last term, we obtain (cid:107) v (cid:107) B ≤ (cid:107) ε (cid:107) L ([0 ,T ]) · (cid:32) ∞ (cid:88) k =1 C k T αk Γ( αk + 1) (cid:33) , which corresponds to (4.2). Theorem
Let F : B → ˆ B be a uniformly Lipschitz continuous operator, i.e.,there exist C > and C > such that for any v , v ∈ B , v ∈ H and t ∈ [0 , T ] , (4.4) |(cid:104) F [ v ]( t ) − F [ v ]( t ) , v (cid:105) H (cid:48) , H | ≤ C (cid:107) v ( t ) − v ( t ) (cid:107) H · (cid:107) v (cid:107) H . and for h > , |(cid:104) F [ v ]( t + h ) − F [ v ]( t ) , v (cid:105) H (cid:48) , H | ≤ C h (cid:107) v (cid:107) H . And let u ( t ) = u + K α ∗ F [ u ] , ˜ u ( t ) = u + ˜ K α ∗ F [˜ u ] , such that for h > , there exists C α > that (cid:13)(cid:13)(cid:13) K α − ˜ K Expα (cid:13)(cid:13)(cid:13) L ([0 ,h ]) ≤ C α h α and (4.5) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) h (cid:16) K α ( s ) − ˜ K Expα ( s ) (cid:17) d s − c ∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:13)(cid:13)(cid:13) K α − ˜ K Expα (cid:13)(cid:13)(cid:13) L ([ h,T ]) ≤ C α h α . Then, the following error estimate holds: (4.6) (cid:107) u − ˜ u (cid:107) B ≤ h α · C α E α, [ C T α ] · (cid:0) C + (cid:107) ∂ αt u (cid:107) ˆ B (cid:1) . Proof.
Let us introduce w ( t ) = u + K α ∗ F [˜ u ]( t ). Then, adding 0 = w − w , wefind for the norm(4.7) (cid:107) u − ˜ u (cid:107) H ≤ |(cid:104) u − w, u − ˜ u (cid:105) H | + |(cid:104) w − ˜ u, u − ˜ u (cid:105) H | . We denote by g = K α − ˜ K Expα the difference of the kernels, and by g [ a,b ] a functionwhich coincides with g on [ a, b ] and vanishes elsewhere. Note that we can writeformally ( K α − ˜ K α ) ∗ F [ u ] = (cid:90) tt − h g ( t − s ) F [ u ]( s ) d s − c ∞ F [ u ]( t ) + g [ h,T ] ∗ F [ u ] . Then, by Lemma 4.1, the first term in (4.7) can be bounded as(4.8) |(cid:104) u − w, u − ˜ u (cid:105) H | = (cid:12)(cid:12)(cid:12) (cid:104) ( K α − ˜ K α ) ∗ F [˜ u ] , u − ˜ u (cid:105) H (cid:12)(cid:12)(cid:12) ≤ ε ( t ) · (cid:107) u − ˜ u (cid:107) H , U. KHRISTENKO, AND B. WOHLMUTH where the error ε ( t ) is defined as(4.9) ε ( t ) := C h (cid:107) g (cid:107) L ([0 ,h ]) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) h g ( s ) d s − c ∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) · (cid:107) F [˜ u ] (cid:107) H (cid:48) + (cid:12)(cid:12) g [ h,T ] (cid:12)(cid:12) ∗ (cid:107) F [˜ u ] (cid:107) H (cid:48) . By the continuity assumption (4.4), the following upper bound holds for the secondterm in (4.7): |(cid:104) w − ˜ u, u − ˜ u (cid:105) H | = |(cid:104) K α ∗ ( F [ u ] − F [˜ u ]) , u − ˜ u (cid:105) H | (4.10) ≤ C ( K α ∗ (cid:107) u − ˜ u (cid:107) H ) · (cid:107) u − ˜ u (cid:107) H . Substituting (4.8) and (4.10) to (4.7), we obtain (cid:107) u − ˜ u (cid:107) H ≤ ε ( t ) + C ( K α ∗ (cid:107) u − ˜ u (cid:107) H ) . Hence, by Lemma 4.2, we have(4.11) (cid:107) u − ˜ u (cid:107) B ≤ (cid:107) ε (cid:107) L ([0 ,T ]) · E α, [ C T α ] . Using Young’s convolution inequality and (4.5), we obtain from (4.9): (cid:107) ε (cid:107) L ([0 ,T ]) ≤ C h (cid:107) g (cid:107) L ([0 ,h ]) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h (cid:90) g ( s ) d s − c ∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:107) g (cid:107) L ([ h,T ]) · (cid:107) F [˜ u ] (cid:107) ˆ B (4.12) ≤ C α h α ( C + (cid:107) F [˜ u ] (cid:107) ˆ B ) . Moreover, by the continuity (4.4), it holds that(4.13) (cid:107) F [˜ u ] (cid:107) ˆ B ≤ (cid:107) F [ u ] (cid:107) ˆ B + (cid:107) F [ u ] − F [˜ u ] (cid:107) ˆ B ≤ (cid:107) ∂ αt u (cid:107) ˆ B + C (cid:107) u − ˜ u (cid:107) ˆ B . Thus, under condition that C C α h α ≤ /
2, combining (4.11), (4.12) and (4.13), weeventually obtain (4.6).
5. Numerical schemes.
Let us denote by ˜ u h the numerical approximation of ˜ u ,defined by (3.2), and by u h,k the approximations of the modes u k , k = 1 , . . . , m ,in (3.3). The discretized solution of the system of equations (3.3) can be numericallycomputed with any known numerical scheme. The simplest case of the so-called θ -scheme, including Euler and Crank-Nicolson time-integration, is introduced in thefollowing proposition. Proposition
Applying a standard θ -scheme to the modal system (3.3) , itstraightforwardly yields the following scheme for (3.2) : (5.1) ˜ u n +1 h − β F [˜ u n +1 h ] = β F [˜ u nh ] + u + m (cid:88) k =1 γ k u nh,k , where the discrete modes u nh,k , k = 1 , . . . , m , are updated by (5.2) u n +1 h,k = γ k u nh,k + β k F [˜ u n +1 h ] + β k F [˜ u nh ] , u h,k = 0 , with coefficients (5.3) β k = c k θh θd k h , β k = c k (1 − θ ) h θd k h , γ k = 1 − (1 − θ ) d k h θd k h ,β = m (cid:88) k =1 β k + c ∞ , β = m (cid:88) k =1 β k . OLVING TIME-FDE VIA RATIONAL APPROXIMATION Remark β corresponds to the expansion (3.1) and thus ap-proximates ( θh ) α . Moreover, we have β = ( θh ) α if z = θh is a support point in theAAA algorithm. Remark α = 1, there is only one mode c = 1, d = 0, therefore, theabove scheme reduces to the classical integer-order θ -scheme.Let us note that the proposed scheme does not require the solution of a largecoupled system of equations, but consists in alternating updates of the full solution ˜ u h ,Eq. (5.1), and updates of the modes u h,k , Eq. (5.2). A graphical illustration of thealgorithm is suggested in Figure 2. The modes updates (5.3) are completely decoupledand can be computed in parallel. Besides, they are linear. A non-linear equation ofthe original size has to be solved only once per time-step in (5.1), using any preferrednon-linear solver (e.g., Newton-Raphson or Fixed-point). In particular, in the PDEcase, when F involves a spacial differential operator, the PDE system is solved onlyin (5.1). Moreover, for the updates (5.3), one does not even have to solve a mass matrixsystem. Indeed, instead of computing the modes u nh,k , k = 1 , . . . , m , themselves, onecan proceed with numerical integration computing only M u nh,k , where M stands forthe formal mass matrix. And no explicit computation of the modes is necessary forcomputing the full solution ˜ u nh . ˜ u n (cid:80) Solve non-linear system (5.1):˜ u n , { u nk } mk =1 → ˜ u n +1 ˜ u n +1 (cid:80) { u nk } mk =1 { u n +1 k } mk =1 Update the modes with (5.2)
Fig. 2: Scheme of the n -th time iteration, representing the system (5.1)-(5.2). Boththe integral approximation ˜ u n and the family of modes { u nk } mk =0 have to be updated ineach time step. However, each step requires only one non-linear system of the originalsize to solve.The proposed θ -scheme is simple, however, it does not guarantee unconditionalstability for an arbitrary operator F . In particular, unconditionally stable schemesare usually based on the splitting of the operator [27, 28]. So, let us consider F in theform F [ v ] = F − [ v ] + F + [ v ], where the monotonous operators F − and F + correspondsto the decreasing and strictly increasing parts of F , respectively. That is, for all v , v ∈ B , it holds (cid:104) F − [ v ] − F − [ v ] , v − v (cid:105) H ≤ , (cid:104) F + [ v ] − F + [ v ] , v − v (cid:105) H > . In addition, we rewrite the continuity condition with some C F > (cid:107) F − [ v ] − F − [ v ] (cid:107) H (cid:48) + (cid:107) F + [ v ] − F + [ v ] (cid:107) H (cid:48) ≤ C F (cid:107) v − v (cid:107) H . In the following lemmas, we propose numerical schemes for such F and estimate theassociated discretization error (cid:15) nh := (cid:107) ˜ u ( t n ) − ˜ u nh (cid:107) H . We also introduce the modaldiscretization errors (cid:15) nh,k := (cid:13)(cid:13)(cid:13) u k ( t n ) − u nh,k (cid:13)(cid:13)(cid:13) H .Let us also briefly remark that the solution of a FDE usually exhibits a weaksingularity at the initial time. So, the approximation error near the origin has to be0 U. KHRISTENKO, AND B. WOHLMUTH treated in a special way, see, e.g., [41, 62]. For the following error estimates, we focuson the evolution far enough from the singularity and thus assume the initial state tobe smooth enough to ignore this effect.
Lemma
Let ˜ u nh be defined by the following time-steppingscheme: (5.4) ˜ u n +1 h − β F − [˜ u n +1 h ] = β F + [˜ u nh ] + u + m (cid:88) k =1 γ k u nh,k , where the discrete modes u nh,k , k = 1 , . . . , m , are updated by (5.5) u n +1 h,k = γ k u nh,k + β k F [˜ u n +1 h ] , u h,k = 0 , with coefficients given by γ k = 11 + d k h , β k = c k h d k h , β = m (cid:88) k =1 β k + c ∞ . Then, ˜ u nh approximates ˜ u ( t ) , provided the discretization error of order h : (cid:15) nh = (cid:107) ˜ u ( t n ) − ˜ u nh (cid:107) H ≤ O ( mh ) · e C F [( t n ) α + ε ra ] , where ε ra := max s ∈ [ T , h ] (cid:12)(cid:12)(cid:12) ˆ K ( s ) − ˆ˜ K ( s ) (cid:12)(cid:12)(cid:12) stands for the rational approximation error.Proof. Taylor expansion of u k ( t n ) at the point t n +1 , with the first derivative givenby (3.3), yields u k ( t n +1 ) = u k ( t n ) − h (cid:0) d k u k ( t n +1 ) − c k F [˜ u ]( t n +1 ) (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) ∂ t u k ( t n +1 ) + O ( h ) . Hence, we express u k ( t n +1 ):(5.6) u k ( t n +1 ) = γ k u k ( t n ) + β k F [˜ u ]( t n +1 ) + γ k O ( h ) . Summing up the modes, u and u ∞ , we obtain˜ u ( t n +1 ) − β F [˜ u ]( t n +1 ) = u + m (cid:88) k =1 γ k u k ( t n ) + O ( mh ) . Note that β = O ( h α ). Then, using Taylor expansion of F + [˜ u ]( t n +1 ) at the point t n ,we can write(5.7) ˜ u ( t n +1 ) − β F − [˜ u ]( t n +1 ) = β F + [˜ u ]( t n ) + u + m (cid:88) k =1 γ k u k ( t n ) + O ( mh + h α ) . Recall that F is Lipschitz continuous and F − monotonously decreases, which implies(5.8) |(cid:104) F [˜ u ]( t n ) − F [˜ u h ] , ˜ u ( t n ) − ˜ u nh (cid:105) H | ≤ C F | (cid:15) nh | , (cid:104) F − [˜ u ]( t n ) − F − [˜ u h ] , ˜ u ( t n ) − ˜ u nh (cid:105) H ≤ . OLVING TIME-FDE VIA RATIONAL APPROXIMATION (cid:15) n +1 h ≤ β C F (cid:15) nh + m (cid:88) k =1 γ k (cid:15) nh,k + O ( mh ( h + h α /m )) , and(5.10) (cid:15) n +1 h,k ≤ γ k (cid:15) nh,k + β k C F (cid:15) n +1 h + O ( h ) , (cid:15) h,k = 0 . Recursive substitution in (5.10) yields to (cid:15) nh,k ≤ β k C F n (cid:88) j =1 γ n − jk (cid:15) jh + O ( h ) . And substituting this to (5.9), we end up with (cid:15) n +1 h ≤ β C F (cid:15) nh + C F m (cid:88) k =1 n (cid:88) j =1 β k γ n +1 − jk (cid:15) jh + O ( mh ) , (cid:15) h = 0 . Hence, the discrete Gr¨onwall inequality completes the proof: (cid:15) nh ≤ O ( mh ) · e C F [ (cid:80) mk =1 (cid:80) n − j =0 β k γ jk + c ∞ ] ≤ O ( mh ) · e C F [( t n ) α + ε ra ] , where we used the following bound for the exponent: n − (cid:88) j =0 β k γ jk = β k − γ nk − γ k = c k − γ nk d k ≤ c k d k (cid:18) −
11 + nd k h (cid:19) ≤ c k nh d k nh . Remark θ = 1 / − ≤ γ k ≤ θ = 1 /
2, moreover, γ k approaches − d k , giving rise to a stiffproblem. Thus, the higher modes produce undesired oscillation of the solution. Sincemax d k grows when α decays, the oscillations become more dominant the smaller α is.An example can be found in the next section (Figure 4). This observation motivates usto introduce in the following lemma a modified scheme which expresses more stabilitybut preserves the order. Lemma
Let ˜ u nh be defined by the following time-stepping scheme: (5.11) ˜ u n +1 h − β F − [˜ u n +1 h ] = β F + [˜ u nh ] + u + m (cid:88) k =1 γ k u nh,k , where the discrete modes u nk , k = 1 , . . . , m , are updated by (5.12) u n +1 h,k = γ k u nh,k + β k F [˜ u n +1 h ] + β k F [˜ u nh ] , u h,k = 0 , with coefficients (5.13) γ k = e − d k h , β k = c k γ k − (1 − d k h ) d k h , β k = c k − (1 + d k h ) γ k d k h ,β = m (cid:88) k =1 ( β k + β k ) + c ∞ = m (cid:88) k =1 c k d k (1 − γ k ) + c ∞ . U. KHRISTENKO, AND B. WOHLMUTH
Then, ˜ u nh approximates ˜ u ( t ) , provided the discretization error of order h α : (cid:15) nh = (cid:107) ˜ u ( t n ) − ˜ u nh (cid:107) H ≤ O ( h α ) · e C F [( t n ) α + ε ra ] . where ε ra := max s ∈ [ T , h ] (cid:12)(cid:12)(cid:12) ˆ K ( s ) − ˆ˜ K ( s ) (cid:12)(cid:12)(cid:12) stands for the rational approximation error.Proof. From (3.3), the modes satisfy the recurrence relation(5.14) u k ( t n +1 ) = u k ( t n ) e − d k h + c k t n +1 (cid:90) t n e − d k ( t − s ) F ( s, ˜ u ( s )) d s. Let us consider the following quadrature rule for an arbitrary function f ( t ) andscalar λ ≥ t n +1 (cid:90) t n e − λ ( t − s ) f ( s ) d s = a f ( t n +1 ) + f ( t n )2 + a f ( t n +1 ) − f ( t n ) h + O ( a h )= (cid:16) a a h (cid:17) f ( t n +1 ) + (cid:16) a − a h (cid:17) f ( t n ) + O ( a h )with coefficients given as a = a ( λ ) = h (cid:90) e − λs d s ≤ h λh , a = a ( λ ) = h (cid:90) e − λs (cid:18) h − s (cid:19) d s Applying the above quadrature rule to (5.14), given β k := c k (cid:16) a ( d k )2 + a ( d k ) h (cid:17) and β k = c k (cid:16) a ( d k )2 − a ( d k ) h (cid:17) , we obtain(5.15) u k ( t n +1 ) = γ k u k ( t n ) + β k F [˜ u ]( t n +1 ) + β k F [˜ u ]( t n ) + O ( c k a ( d k ) h ) . Note that 0 ≤ a ( d k ) ≤ h d k h and thus 0 ≤ (cid:80) mk =1 c k a ( d k ) ≤ h α − c ∞ ≤ h α .Summing up the modes, u and u ∞ , we obtain˜ u ( t n +1 ) − β F [˜ u ]( t n +1 ) = β F [˜ u ]( t n ) + u + m (cid:88) k =1 γ k u k ( t n ) + O ( h α ) , where(5.16) β := m (cid:88) k =0 β k + c ∞ and β := m (cid:88) k =0 β k . Note that β k + β k ≤ c k h d k h and thus β + β = O ( h α ). Then, using Taylor expansionof F + [˜ u ]( t n +1 ) at the point t n and respectively F − [˜ u ]( t n ) at the point t n +1 , we write(5.17) ˜ u ( t n +1 ) − β F − [˜ u ]( t n +1 ) = β F + [˜ u ]( t n ) + u + m (cid:88) k =1 γ k u k ( t n ) + O ( h α ) . OLVING TIME-FDE VIA RATIONAL APPROXIMATION (cid:15) n +1 h ≤ β C F (cid:15) nh + m (cid:88) k =1 γ k (cid:15) nh,k + O ( h α ) , and(5.19) (cid:15) n +1 h,k ≤ γ k (cid:15) nh,k + C F (cid:0) β k (cid:15) n +1 h + β k (cid:15) nh (cid:1) + O (cid:18) c k h d k h (cid:19) , (cid:15) h,k = 0 . where we used continuity of F and monotonicity of F − as in (5.8). Recursive substi-tution in (5.19), yields to (cid:15) nh,k ≤ C F n (cid:88) j =1 ( β k + β k ) γ n − jk (cid:15) jh + O (cid:18) c k h d k h (cid:19) . Substituting this to (5.18), we thus write the estimation (cid:15) n +1 h ≤ C F β (cid:15) nh + m (cid:88) k =1 n (cid:88) j =1 ( β k + β k ) γ n +1 − jk (cid:15) jh + O ( h α ) , (cid:15) h = 0 . Hence, by the discrete Gr¨onwall inequality, we finally obtain (cid:15) nh ≤ O ( h α ) · e C F [ (cid:80) mk =1 (cid:80) n − j =0 ( β k + β k ) γ jk + c ∞ ] ≤ O ( h α ) · e C F [( t n ) α + ε ra ] , using the following bound for the exponent: n − (cid:88) j =0 ( β k + β k ) γ jk = c k d k (1 − γ k ) 1 − γ nk − γ k = c k t n (cid:90) e − d k t d t ≤ c k t n d k t n . Remark
Remark u ∞ which does not however enter to the history part). Moreover, let us remarkthat the coefficients β and β in (5.16) approximate the 2nd order fractional Adams–Moulton coefficients [60], β AM = h α / Γ( α + 2) and β AM = αh α / Γ( α + 2), respectively.The values of the coefficients for α = 0 . α = 0 . U. KHRISTENKO, AND B. WOHLMUTH − − − − − − h r e l. e rr o r α = 0 . β β β AM β AM α = 0 . β β β AM β AM Fig. 3: Comparison of the coefficients β and β in (5.13) for the modified Crank-Nicolson scheme (Lemma 5.6) with the 2nd order fractional Adams–Moulton coeffi-cients β AM = h α / Γ( α + 2) and β AM = αh α / Γ( α + 2), respectively.
6. Numerical examples.
In this section, we illustrate the proposed scheme inapplication to the two following examples. We first consider a simple linear case, moreprecisely, the one-dimensional fractional heat equation, where the analytical solutionis known and given by a Mittag-Leffler function, so that we can study the accuracy andconvergence rate. Then, the scheme is applied to the more complex non-linear Cahn-Hilliard equation and compared to a classical fractional time-stepping scheme. Bothproblems are discretized in space with Finite Elements using the FEniCS package [1].
Let us consider one-dimensional fractionalheat equation with homogeneous Dirichlet boundary conditions: ∂ αt u ( t, x ) − ∂ x u ( t, x ) = 0 , t > , x ∈ (0 , ,u (0 , x ) = sin( π x ) ,u ( t,
0) = u ( t,
1) = 0 . Its analytical solution is given by u ( t, x ) = E α, (cid:2) − π t α (cid:3) sin( π x ), see, e.g., [35],where E α,β [ x ] is the Mittag-Leffler function (4.3). Note that in this example, we have F [ v ] = F − [ v ] = ∂ x v , i.e. F + [ v ] ≡
0. And thus, the scheme in Lemma 5.4 coincideswith the implicit Euler scheme in Proposition 5.1 ( θ = 1).The solution of a FDE usually has a weak singularity at the initial time andthus pollutes the error at the beginning of the time interval. Thus, the optimalconvergence rate can then not be observed globally. Therefore, the error in the originhas to be treated in a special way, see, e.g., [41, 62]. Here we just ignore this effectand consider the error on the shortened interval [ T / , T ], where the effect of theinitial singularity has been faded out. So, we are going to study the relative error E r := (cid:107) u − ˜ u h (cid:107) / (cid:107) u − u (cid:107) , where the norm is defined via (cid:96) -norm in time on theshortened interval [ T / , T ] and L -norm in space (using FE). We use 1000 P -elementsfor the space discretization and the AAA-tolerance 10 − with 100 support points forthe rational approximation.We start with a comparison of the Crank-Nicolson scheme with its stable mod-ification (Lemma 5.6). The corresponding solutions (norms) with the time step size h = 10 − for the cases α = 0 .
1, 0 .
3, 0 . OLVING TIME-FDE VIA RATIONAL APPROXIMATION α decreases. .
02 0 .
04 0 .
06 0 .
08 0 . t (cid:107) u (cid:107) α = 0 . α = 0 . α = 0 . .
02 0 .
04 0 .
06 0 .
08 0 . t (cid:107) u (cid:107) α = 0 . α = 0 . α = 0 . Fig. 4: Illustration of instability effects of the Crank-Nicolson scheme (left) for α = 0 . .
3, 0 . h = 10 − . The effect is stronger for smaller α . However, there is nooscillatory effect in case of the modified scheme (5.11)-(5.12) (right).Then, we study the convergence rate of the relative error E r with respect to thetime step size h for the Implicit Euler scheme (5.4)-(5.5) and the modified Crank-Nicolson scheme (5.11)-(5.12). The results for α = 1, 0 .
8, 0 .
5, 0 .
3, 0 .
1, 0 .
03, 0 .
01 areplotted in Figure 5 and confirm the theoretical error bounds suggested in Lemmas 5.4and 5.6, respectively. In particular, for the first scheme, we clearly observe linearconvergence rate for all α , while the second scheme presents convergence of order h α . − − − − − − − − h E r IE α = 1 α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . E r :Slope h (left) or h α (right): − − − − − − − − h E r CN Fig. 5: Convergence rate of the Implicit Euler scheme (5.4)-(5.5) (left) and the mod-ified Crank-Nicolson scheme (5.11)-(5.12) (right).6
U. KHRISTENKO, AND B. WOHLMUTH
Let Ω = (0 , be a unit squaredomain and B = L ([0 , T ]; H ), where H = H (Ω) × H (Ω). Moreover, let (cid:104)· , ·(cid:105) and (cid:107)·(cid:107) denote the scalar product and the norm in L (Ω), respectively. Then, weformulate the following non-linear Cahn-Hilliard problem: find ( u, µ ) ∈ B satisfyingfor all ( v, w ) ∈ H and all t ∈ (0 , T ],(6.1) (cid:104) ∂ αt u, v (cid:105) = −(cid:104) M ∇ µ, ∇ v (cid:105) , (cid:104) µ, w (cid:105) = (cid:104) ψ ( u ) , w (cid:105) + ε (cid:104)∇ u, ∇ w (cid:105) , provided homogeneous Neumann boundary conditions, the initial state u (0) = u , theconstant mobility M and the surface parameter ε . The non-linear function ψ ( u ) isdefined as derivative of the potential Φ( u ) = ( u − : ψ ( u ) = Φ (cid:48) ( u ) = u − u. The Ginzburg–Landau free energy of the system is defined as E ( u ) = (cid:90) Ω (cid:20) Φ( u ) + ε |∇ u | (cid:21) dΩ . Note that the function ψ ( u ) is Lipschitz continuous on the interval [ − , u ). Therefore, if the initial conditions are contained in theinterval, the problem (6.1) satisfies the conditions of Theorem 4.3.Due to the ”double-well” structure of the potential, presenting both convex andconcave parts, stability of time-schemes for Cahn-Hilliard equation is a sophisticatedquestion and is a subject of numerous works. The fully implicit time-schemes are onlyconditionally stable [25]. Unconditionally stable schemes include so-called gradientstability, providing monotone decay of the discretized Ginzburg–Landau energy, e.g.,splitting to implicit convex and explicit concave parts [27, 57] or others [24, 30]. Thesituation becomes more complicated in the case of fractional derivative, since evenfor the analytical solution, the associated Ginzburg–Landau energy is not proved tobe monotone [54]. Note that the schemes presented in Lemmas 5.4 and 5.6 naturallyallow splitting techniques, ensuring stability. Splitting the potential into convex andconcave parts, we consider ψ ( u ) = ψ + ( u ) + ψ − ( u ), where the functions ψ + ( u ) and ψ − ( u ) are monotonously increasing and decreasing in [ − , ψ + ( u ) = 2 u, ψ − ( u ) = u − u. (6.2)Let H be an appropriate finite elements space. Implementing the modified Crank-Nicolson scheme (5.11)-(5.13) for discretization of the problem (6.1), we compute ateach time step n the discrete solution pair ( u n +1 h , µ n +1 h ) ∈ H satisfying (cid:104) u n +1 h − H n , v (cid:105) = − β (cid:104) M ∇ ˜ µ n +1 h , ∇ v (cid:105) , (6.3) (cid:104) ˜ µ n +1 h , w (cid:105) = (cid:104) ψ + ( u n +1 h ) + ψ − ( u nh ) , w (cid:105) + ε (cid:104)∇ u n +1 h , ∇ w (cid:105) , for all ( v, w ) ∈ H , with the history term defined as H n = u + (cid:80) mk =1 γ k u nh,k , wherethe modes u nk , k = 1 , . . . , m , are updated as follows:(6.4) (cid:104) u n +1 h,k , v (cid:105) = γ k (cid:104) u nh,k , v (cid:105) − (cid:104) M ∇ (cid:2) β k µ n +1 + β k µ n (cid:3) , ∇ v (cid:105) , u h,k = 0 , (cid:104) µ n +1 h , w (cid:105) = (cid:104) ψ ( u n +1 h ) , w (cid:105) + ε (cid:104)∇ u n +1 h , ∇ w (cid:105) , OLVING TIME-FDE VIA RATIONAL APPROXIMATION β , β k , β k and γ k given in (5.13). Remark that a linear choice ofthe increasing part in (6.2) leads to a linear implicit part in the time-scheme. That is,though the problem (6.1) is non-linear, the numerical solution of (6.3) requires onlya linear solver at each time step.For our simulation, we fix the constant mobility M = 0 .
05, the surface parame-ter ε = 0 .
03 and the final time T = 4. For discretization in space, we use ( Q , Q )elements on a 64 ×
64 quadrilateral mesh. For the rational approximation, we use theAAA-tolerance 10 − with 100 support points. The initial state is u ( xxx ) = (cid:88) i =1 tanh r − | xxx − xxx i |√ ε + 3 , xxx ∈ Ω , with xxx = (0 . , . xxx = (0 . , . xxx = (0 . , . xxx = (0 . , .
3) and r = 0 . r centered at xxx i . Due to the surfacetension, the bubbles tend to coalesce in time [38]. However, the process proceeds withdifferent speed for different values of the fractional order α . In particular, for a small α , the coalescence accelerates in the beginning but then slows down with respect tolarger values of α . This effect can be observed in Figure 6, where different states(computed with h = 2 − ) are shown for α = 0 .
1, 0 .
3, 0 .
5, 0 . t = 0, 0 .
4, 3 . E r := (cid:107) u − u h (cid:107) / (cid:107) u − u (cid:107) ,where the norm is again defined via (cid:96) -norm in time on the shortened interval [ T / , T ]and L -norm in H . In Figure 7 on the right, there are plotted the convergence ratesof the error with respect to the time step size h for the same values of α . We canobserve that the error convergence respects the theoretical bounds.
7. Conclusion.
In this work, we proposed a new numerical method for solvingfractional in time differential equation. The method is based on the approximationof the Laplace spectrum of the fractional convolution kernel with a rational function,more precisely, a multi-pole series with an additional constant term. To this end,we used the barycentric rational interpolation with the adaptive Antoulas–Anderson(AAA) algorithm. This leads to the approximation of the kernel itself with a sum-of-exponentials with an additional singular term. Thus, the solution of the FODEis represented as a sum of a small number of modes which solve a system of ODEsand can be updated in parallel. Since the number of modes is fixed and definedonly by the fractional kernel, the computational complexity and the memory load isconstant for each time step. We proposed the associated variants of the implicit Eulerand Crank-Nicolson numerical schemes with convergence orders O ( h ) and O ( h α ),respectively. The accuracy of the schemes is illustrated using a linear problem withknown analytical solution. The method has been also applied to a non-linear fractionalCahn-Hilliard problem in 2D. REFERENCES[1] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring,M. E. Rognes, and G. N. Wells. The fenics project version 1.5.
Archive of NumericalSoftware , 3(100), 2015.[2] D. Baffet. A gauss–jacobi kernel compression scheme for fractional differential equations.
Jour-nal of Scientific Computing , 79(1):227–248, 2019.[3] D. Baffet and J. S. Hesthaven. A kernel compression scheme for fractional differential equations.
SIAM Journal on Numerical Analysis , 55(2):496–520, 2017. U. KHRISTENKO, AND B. WOHLMUTH α = 0 . α = 0 . α = 0 . α = 0 . t = t = . t = . t = Fig. 6: Fractional non-linear Cahn-Hilliard equation with α = 0 .
1, 0 .
3, 0 .
5, 0 .
9. Fourbubbles coalesce with different speed: for smaller α , the coalescence accelerates in thebeginning but then slows down with respect to larger values of α . . . . . . t E α = 0 . α = 0 . α = 0 . α = 0 . − . − − . − − − − − h E r E r h α α = 0 . α = 0 . α = 0 . α = 0 . Fig. 7: Fractional non-linear Cahn-Hilliard equation with α = 0 .
1, 0 .
3, 0 .
5, 0 . Left:
Evolution in time of the Ginzburg–Landau energy.
Right:
Convergence rate for therational approximation based scheme (6.3)-(6.4).
OLVING TIME-FDE VIA RATIONAL APPROXIMATION [4] E. G. Bajlekova et al. Fractional evolution equations in Banach spaces . Citeseer, 2001.[5] G. A. Baker, G. A. Baker Jr, G. Baker, P. Graves-Morris, and S. S. Baker.
Pade Approximants:Encyclopedia of Mathematics and Its Applications, Vol. 59 George A. Baker, Jr., PeterGraves-Morris , volume 59. Cambridge University Press, 1996.[6] D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo.
Fractional calculus: models and numer-ical methods , volume 3. World Scientific, 2012.[7] L. Banjai and M. L´opez-Fern´andez. Efficient high order algorithms for fractional integrals andfractional differential equations.
Numerische Mathematik , 141(2):289–317, 2019.[8] L. Banjai, J. M. Melenk, R. H. Nochetto, E. Otarola, A. J. Salgado, and C. Schwab. Tensor femfor spectral fractional diffusion.
Foundations of Computational Mathematics , 19(4):901–962, 2019.[9] J.-P. Berrut, R. Baltensperger, and H. D. Mittelmann. Recent developments in barycentricrational interpolation. In
Trends and applications in constructive approximation , pages27–51. Springer, 2005.[10] G. Beylkin and L. Monz´on. On approximation of functions by exponential sums.
Applied andComputational Harmonic Analysis , 19(1):17–48, 2005.[11] A. Bonito and J. Pasciak. Numerical approximation of fractional powers of elliptic operators.
Mathematics of Computation , 84(295):2083–2110, 2015.[12] L. Caffarelli and L. Silvestre. An extension problem related to the fractional laplacian.
Com-munications in partial differential equations , 32(8):1245–1260, 2007.[13] O. S. Celis and A. Cuyt.
Practical rational interpolation of exact and inexact data: theory andalgorithms . Universiteit Antwerpen, Faculteit Wetenschappen, Departement Wiskunde . . . ,2008.[14] W. Deng. Short memory principle and a predictor–corrector approach for fractional differentialequations.
Journal of Computational and Applied Mathematics , 206(1):174–188, 2007.[15] K. Diethelm. An investigation of some nonclassical methods for the numerical approximationof Caputo-type fractional derivatives.
Numerical Algorithms , 47(4):361–390, 2008.[16] K. Diethelm. An investigation of some nonclassical methods for the numerical approximationof caputo-type fractional derivatives.
Numerical Algorithms , 47(4):361–390, 2008.[17] K. Diethelm.
The analysis of fractional differential equations: An application-oriented expo-sition using differential operators of Caputo type . Springer Science & Business Media,2010.[18] K. Diethelm. An efficient parallel algorithm for the numerical solution of fractional differentialequations.
Fractional Calculus and Applied Analysis , 14(3):475–490, 2011.[19] K. Diethelm, J. M. Ford, N. J. Ford, and M. Weilbeer. Pitfalls in fast numerical solversfor fractional differential equations.
Journal of computational and applied mathematics ,186(2):482–503, 2006.[20] K. Diethelm, N. J. Ford, and A. D. Freed. A predictor-corrector approach for the numericalsolution of fractional differential equations.
Nonlinear Dynamics , 29(1-4):3–22, 2002.[21] K. Diethelm, N. J. Ford, and A. D. Freed. Detailed error analysis for a fractional adams method.
Numerical algorithms , 36(1):31–52, 2004.[22] K. Diethelm and A. D. Freed. An efficient algorithm for the evaluation of convolution integrals.
Computers & Mathematics with Applications , 51(1):51–72, 2006.[23] K. Diethelm, R. Garrappa, and M. Stynes. Good (and not so good) practices in computationalmethods for fractional calculus.
Mathematics , 8(3):324, 2020.[24] Q. Du and R. A. Nicolaides. Numerical analysis of a continuum model of phase transition.
SIAM Journal on Numerical Analysis , 28(5):1310–1322, 1991.[25] C. M. Elliott. The cahn-hilliard model for the kinetics of phase separation. In
Mathematicalmodels for phase change problems , pages 35–73. Springer, 1989.[26] C. L. Epstein and J. Schotland. The bad truth about laplace’s transform.
SIAM review ,50(3):504–520, 2008.[27] D. J. Eyre. Unconditionally gradient stable time marching the cahn-hilliard equation. In
Materials Research Society Symposium Proceedings , volume 529, pages 39–46. MaterialsResearch Society, 1998.[28] D. J. Eyre. An unconditionally stable one-step scheme for gradient systems.
Unpublishedarticle , pages 1–15, 1998.[29] N. J. Ford and A. C. Simpson. The numerical solution of fractional differential equations: speedversus accuracy.
Numerical Algorithms , 26(4):333–346, 2001.[30] H. Gomez and T. J. Hughes. Provably unconditionally stable, second-order time-accurate,mixed variational methods for phase-field models.
Journal of Computational Physics ,230(13):5310–5327, 2011.[31] S. Harizanov and S. Margenov.
Positive approximations of the inverse of fractional powers of U. KHRISTENKO, AND B. WOHLMUTH
SPD M-matrices , pages 147–163. Springer, 2018.[32] S. Jiang, J. Zhang, Q. Zhang, and Z. Zhang. Fast evaluation of the caputo fractional derivativeand its applications to fractional diffusion equations.
Communications in ComputationalPhysics , 21(3):650–678, 2017.[33] B. Jin, R. Lazarov, and Z. Zhou. An analysis of the l1 scheme for the subdiffusion equationwith nonsmooth data.
IMA Journal of Numerical Analysis , 36(1):197–221, 2016.[34] E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001–.[35] L. Kexue and P. Jigen. Laplace transform and fractional differential equations.
Applied Math-ematics Letters , 24(12):2019–2023, 2011.[36] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo.
Theory and applications of fractionaldifferential equations , volume 204. elsevier, 2006.[37] J.-R. Li. A Fast Time Stepping Method for Evaluating Fractional Integrals.
SIAM Journal onScientific Computing , 31(6):4696–4714, 2010.[38] H. Liu, A. Cheng, H. Wang, and J. Zhao. Time-fractional allen–cahn and cahn–hilliard phase-field models and their numerical investigation.
Computers & Mathematics with Applica-tions , 76(8):1876–1892, 2018.[39] M. L´opez-Fern´andez, C. Lubich, and A. Sch¨adle. Adaptive, fast, and oblivious convolution inevolution equations with memory.
SIAM Journal on Scientific Computing , 30(2):1015–1037, 2008.[40] C. Lubich. On the stability of linear multistep methods for volterra convolution equations.
IMA Journal of Numerical Analysis , 3(4):439–465, 1983.[41] C. Lubich. Discretized fractional calculus.
SIAM Journal on Mathematical Analysis , 17(3):704–719, 1986.[42] C. Lubich. Convolution quadrature and discretized operational calculus. i.
Numerische Math-ematik , 52(2):129–145, 1988.[43] C. Lubich and A. Sch¨adle. Fast convolution for nonreflecting boundary conditions.
SIAMJournal on Scientific Computing , 24(1):161–182, 2002.[44] F. Mainardi.
Fractional calculus and waves in linear viscoelasticity: an introduction to math-ematical models . World Scientific, 2010.[45] W. McLean. Exponential sum approximations for t- β . In Contemporary ComputationalMathematics-A Celebration of the 80th Birthday of Ian Sloan , pages 911–930. Springer,2018.[46] W. McLean, I. H. Sloan, and V. Thom´ee. Time discretization via laplace transformation of anintegro-differential equation of parabolic type.
Numerische Mathematik , 102(3):497–522,2006.[47] K. S. Miller and B. Ross.
An introduction to the fractional calculus and fractional differentialequations . Wiley, 1993.[48] Y. Nakatsukasa, O. S`ete, and L. N. Trefethen. The aaa algorithm for rational approximation.
SIAM Journal on Scientific Computing , 40(3):A1494–A1522, 2018.[49] K. Oldham and J. Spanier.
The fractional calculus theory and applications of differentiationand integration to arbitrary order . Elsevier, 1974.[50] I. Podlubny.
Fractional differential equations: an introduction to fractional derivatives, frac-tional differential equations, to methods of their solution and some of their applications .Elsevier, 1998.[51] S. G. Samko, A. A. Kilbas, O. I. Marichev, et al.
Fractional integrals and derivatives , volume 1.Gordon and Breach Science Publishers, Yverdon Yverdon-les-Bains, Switzerland, 1993.[52] A. Sch¨adle, M. L´opez-Fern´andez, and C. Lubich. Fast and oblivious convolution quadrature.
SIAM Journal on Scientific Computing , 28(2):421–438, 2006.[53] H. R. Stahl. Best uniform rational approximation ofx α on [0, 1]. Acta mathematica , 190(2):241–306, 2003.[54] T. Tang, H. Yu, and T. Zhou. On energy dissipation theory and numerical stability for time-fractional phase-field equations.
SIAM Journal on Scientific Computing , 41(6):A3757–A3778, 2019.[55] L. N. Trefethen.
Approximation theory and approximation practice , volume 164. Siam, 2019.[56] P. N. Vabishchevich. Numerically solving an equation for fractional powers of elliptic operators.
Journal of Computational Physics , 282:289–302, 2015.[57] X. Wu, G. Van Zwieten, and K. Van der Zee. Stabilized second-order convex splitting schemesfor cahn–hilliard models with application to diffuse-interface tumor-growth models.
Inter-national journal for numerical methods in biomedical engineering , 30(2):180–203, 2014.[58] H. Ye, J. Gao, and Y. Ding. A generalized gronwall inequality and its application to a fractionaldifferential equation.
Journal of Mathematical Analysis and Applications , 328(2):1075–1081, 2007.OLVING TIME-FDE VIA RATIONAL APPROXIMATION [59] L. Yuan and O. P. Agrawal. A numerical scheme for dynamic systems containing fractionalderivatives. Journal of Vibration and Acoustics, Transactions of the ASME , 124(2):321–324, 2002.[60] M. Zayernouri and A. Matzavinos. Fractional adams–bashforth/moulton methods: an applica-tion to the fractional keller–segel chemotaxis system.
Journal of Computational Physics ,317:1–14, 2016.[61] F. Zeng, I. Turner, and K. Burrage. A stable fast time-stepping method for fractional integraland derivative operators.
Journal of Scientific Computing , 77(1):283–307, 2018.[62] Y. Zhou, J. L. Suzuki, C. Zhang, and M. Zayernouri. Fast imex time integration of nonlinearstiff fractional differential equations. arXiv preprint arXiv:1909.04132arXiv preprint arXiv:1909.04132