On the construction of m -step methods for FDEs
aa r X i v : . [ m a t h . NA ] M a y On the construction of m -step methods for FDEs Lidia Aceto ∗ , Cecilia Magherini † , Paolo Novati ‡ Abstract
In this paper we consider the numerical solution of Fractional Differ-ential Equations by means of m -step recursions. The construction of suchformulas can be obtained in many ways. Here we study a technique basedon the rational approximation of the generating functions of FractionalBackward Differentiation Formulas (FBDFs). Accurate approximationsallow to define methods which simulate the theoretical properties of theunderlying FBDF with important computational advantages. Numericalexperiments are presented. MSC : 65L06, 65F60, 65D32, 26A33, 34A08.
This paper deals with the solution of Fractional Differential Equations (FDEs)of the type t D αt y ( t ) = g ( t, y ( t )) , t < t ≤ T, < α < , (1)where t D αt denotes the Caputo’s fractional derivative operator (see e.g. [16]for an overview) defined as t D αt y ( t ) = 1Γ(1 − α ) Z tt y ′ ( u )( t − u ) α du, (2)in which Γ denotes the Gamma function. As well known, the use of the Caputo’sdefinition for the fractional derivative allows to treat the initial conditions at t for FDEs in the same manner as for integer order differential equations. Setting y ( t ) = y the solution of (1) exists and is unique under the hypothesis that g iscontinuous and fulfils a Lipschitz condition with respect to the second variable(see e.g. [3] for a proof).As for the integer order case α = 1, a classical approach for solving (1) isbased on the discretization of the fractional derivative (2), which generalizesthe well known Grunwald-Letnikov discretization (see [16, § ∗ Dipartimento di Matematica, Universit`a di Pisa, Italy, [email protected] † Dipartimento di Matematica, Universit`a di Pisa, Italy, [email protected] ‡ Dipartimento di Matematica, Universit`a di Padova, Italy, [email protected] t , t , . . . , t N = T of the time domain with stepsize h = ( T − t ) /N, FBDFs are based on the full-term recursion X nj =0 ω ( p ) n − j y j = h α g ( t n , y n ) , p ≤ n ≤ N, (3)where y j ≈ y ( t j ) and ω ( p ) n − j are the Taylor coefficients of the generating function ω ( α ) p ( ζ ) = ( a + a ζ + ... + a p ζ p ) α (4)= X ∞ i =0 ω ( p ) i ζ i , for 1 ≤ p ≤ , (5)being { a , a , . . . , a p } the coefficients of the underlying BDF. In [12] it is alsoshown that the order p of the BDF is preserved.We remember that for this kind of equations, there is generally an intrinsiclack of regularity of the solution in a neighborhood of the starting point, that is,depending on the function g , one may have y ( t ) ∼ ( t − t ) α as t → t . For thisreason, in order to preserve the theoretical order p of the numerical method,formula (3) is generally corrected as X Mj =0 w n,j y j + X nj =0 ω ( p ) n − j y j = h α g ( t n , y n ) , (6)where the sum P Mj =0 w n,j y j is the so-called starting term, in which M and theweights w n,j depend on α and p (see [2, Chapter 6] for a discussion).Denoting by Π m the set of polynomials of degree not exceeding m , our basicidea is to design methods based on rational approximations of (4), i.e., R m ( ζ ) ≈ ω ( α ) p ( ζ ) , R m ( ζ ) = p m ( ζ ) q m ( ζ ) , p m , q m ∈ Π m . (7)Writing p m ( ζ ) = P mj =0 α j ζ j and q m ( ζ ) = P mj =0 β j ζ j , the above approximationnaturally leads to implicit m -step recursions of the type X nj = n − m α n − j y j = h α X nj = n − m β n − j g ( t j , y j ) , n ≥ m. (8)While the order of the FBDF is lost, we shall see that if the approximation (7)is rather accurate, then (8) is able to produce reliable approximations to thesolution. Starting from the initial data y ( t ) = y , the first m − y , . . . , y m − can be generated by the underlying FBDFs or even consideringlower degree rational approximations.A formula of type (8) generalizes in some sense the methods based on theShort Memory Principle in which the truncated Taylor expansion of (4) is con-sidered (see [16, § α = 1 /
2, the use of a starting formula as in (6), that theoreti-cally should ensure the order of the FBDF, in practice may introduce substantialerrors, causing unreliable numerical solutions. For high-order formulas, this isdue to the severe ill-conditioning of the Vandermonde type systems one has tosolve at each integration step to generate the weights w n,j of the starting term.We also remark that in a typical application α , y and possibly also the function g may be only known up to a certain accuracy (see [3] for a discussion), so thatone may only be interested in having a rather good approximation of the truesolution.For the construction of formulas of type (8), in this paper we present a tech-nique based on the rational approximation of the fractional derivative operator(cf. [14]). After considering a BDF discretization of order p of the first deriva-tive operator, which can be represented by a N × N triangular banded Toeplitzmatrix A p , we approximate Caputo’s fractional differential operator t D αt bycalculating A αp . This computation is performed by means of the contour in-tegral approximation, which leads to a rational approximation of A αp whosecoefficients can be used to define (8). This technique is based on the fact thatthe first column of A αp contains the first N coefficients of the Taylor expansionof ω ( α ) p ( ζ ), so exploiting the equivalence between the approximation of A αp and ω ( α ) p ( ζ ).The outline of the paper is the following. As in [5], the contour integral isevaluated by means of the Gauss-Jacobi rule in Section 2. An error analysisof this approach is outlined in Section 3, together with some numerical experi-ments that confirm its effectiveness. In Section 4 we investigate the reliabilityof this approach for the solution of fractional differential equations. In particu-lar, we present some results concerning the consistency and the linear stability.Finally, in Section 5 we consider the results of the method when applied to thediscretization of two well-known models of fractional diffusion. Denoting by a , a , . . . , a p the p + 1 coefficients of a Backward DifferentiationFormula (BDF) of order p , with 1 ≤ p ≤
6, which discretizes the derivativeoperator (see [8, Chapter III.1] for a background), we consider lower triangularbanded Toeplitz matrices of the type A p = a a a p ... . . . 00 . . . . . . 00 a p · · · a ∈ R ( N +1) × ( N +1) . (9)3n this setting, A αp e , e = (1 , , . . . , T , contains the whole set of coefficients ofthe corresponding FBDF for approximating the solution of (1) in t , t , . . . , t N ,that is e Tj +1 A αp e = ω ( p ) j , ≤ j ≤ N, (10)(cf. (5)). The constraint p ≤ p > . From the theory of matrix functions (see [10] for a survey), we know thatthe fractional power of matrix can be written as a contour integral A α = A πi Z Γ z α − ( zI − A ) − dz, where Γ is a suitable closed contour enclosing the spectrum of A , σ ( A ), in itsinterior. The following known result (see, e.g., [1]) expresses A α in terms of areal integral. Proposition 1
Let A ∈ R N × N be such that σ ( A ) ⊂ C \ ( −∞ , . For < α < the following representation holds A α = A sin( απ ) απ Z ∞ ( ρ /α I + A ) − dρ. (11)Of course the above result holds also in our case since σ ( A p ) = { a } and a > ≤ p ≤ . At this point, for each suitable change of variable for ρ , a k -point quadrature rule for the computation of the integral in (11) yields arational approximation of the type A αp ≈ A p e R k ( A p ) := A p X kj =1 γ j ( η j I + A p ) − , (12)where the coefficients γ j and η j depend on the substitution and the quadratureformula. This technique has been used in [14], where the author applies theGauss-Legendre rule to (11) after the substitution ρ = a α (cos θ ) − α/ (1 − α ) sin θ, which generalizes the one presented in [9] for the case α = 1 / . In order to remove the dependence of α inside the integral we consider thechange of variable ρ /α = τ − t t , τ > , (13)yielding 1 α Z ∞ ( ρ /α I + A p ) − dρ = 2 Z − (cid:18) τ − t t (cid:19) α − (cid:18) τ − t t I + A p (cid:19) − τ (1 + t ) dt = 2 τ α Z − (1 − t ) α − (1 + t ) − α ( τ (1 − t ) I + (1 + t ) A p ) − dt, A αp = 2 sin( απ ) τ α π A p Z − (1 − t ) α − (1 + t ) − α ( τ (1 − t ) I + (1 + t ) A p ) − dt. (14)The above formula naturally leads to the use of a k -point Gauss-Jacobi rule forthe approximation of A αp e and hence to a rational approximation (12).The following result can be proved by direct computation. Proposition 2
Let A p ∈ R N × N be a matrix of type (9), and let A p = a A p .Then the components of ( ξI + A p ) − e , ξ = − , are given by υ ( p )1 ( ξ ) = 1 ξ + 1 ,υ ( p ) j ( ξ ) = c ( p )2 ,j ( ξ + 1) + . . . + c ( p ) j,j ( ξ + 1) j , ≤ j ≤ N, where the coefficients c ( p ) i,j depend on the order p . For p = 1 we simply have { a , a } = { , − } , and hence υ (1) j ( ξ ) = 1( ξ + 1) j , ≤ j ≤ N. The above proposition shows that the components of( τ (1 − t ) I + (1 + t ) A p ) − e are analytic functions in a suitable open set containing [ − ,
1] in its interior,since they are sum of functions of the type(1 + t ) l − ( τ (1 − t ) + a (1 + t )) l , l ≥ , (15)whose pole lies outside [ − ,
1] for τ > a > ≤ p ≤ α, is completely absorbed by the Jacobi weight function so that the Gauss-Jacobirule yields a very efficient tool for the computation of A αp .Increasing k the approximation (12) can be used to approximate the wholeset of coefficients of the FBDFs. We remark that the computation of the vectors( η j I + A p ) − e does not constitute a problem because of the structure of A p (see (9)). We also point out that since our aim is to construct reliable formulasof type (8) we actually do not need to evaluate (12). Indeed we just need toknow the scalars γ j and η j , and then, using an algorithm to transform partialfractions to polynomial quotient, we obtain the approximation z α ≈ z e R k ( z ) = z e p k − ( z ) e q k ( z ) , z = a + a ζ + . . . + a p ζ p , (16)5here e p k − ∈ Π k − and e q k ∈ Π k . This finally leads to the approximation (7)with p m ( ζ ) = ( a + a ζ + . . . + a p ζ p ) e p k − ( a + a ζ + . . . + a p ζ p ) , (17) q m ( ζ ) = e q k ( a + a ζ + . . . + a p ζ p ) , (18)in which m = kp . We remark that whenever the procedure for the definitionof the coefficients γ j and η j has been set for a given α, one can compute thecorresponding coefficients in the m -step formula (8) once and for all. Denoting by J k ( A p ) the result of the Gauss-Jacobi rule for the approximationof J ( A p ) = Z − (1 − t ) α − (1 + t ) − α ( τ (1 − t ) I + (1 + t ) A p ) − dt, by (14) the corresponding approximation to A αp is given by A αp ≈ A p e R k ( A p ) , e R k ( A p ) = 2 sin( απ ) τ α π J k ( A p ) . (19)In this section we analyze the error term componentwise, that is, (see (10), E j := ω ( p ) j − e Tj +1 A p e R k ( A p ) e = e Tj +1 A p (cid:16) A α − p − e R k ( A p ) (cid:17) e (20)= 2 sin( απ ) τ α π e Tj +1 A p ( J ( A p ) − J k ( A p )) e , ≤ j ≤ N, which is the error in the computation of the j -th coefficient of the Taylor expan-sion of ω ( α ) p ( ζ ). Numerically one observes that the quality of the approximationtends to deteriorate when the dimension of the problem N grows. In this sensewe are particularly interested in observing the dependence of the error term on j and k for j ≫ k and to find a strategy to define the parameter τ of the sub-stitution (13) in this situation. As we shall see in the remainder of the paper,this parameter plays a crucial role for the quality of the approximation.We restrict our analysis to the case of p = 1 for which a = − a = 1 . In thissituation, defining the vector r := ( J ( A ) − J k ( A )) e , (21)we have that E j = 2 sin( απ ) τ α π e Tj +1 A r, and therefore | E j | ≤ απ ) τ α π ( | r j | + | r j − | ) . (22)6he analysis thus reduces to the study of the components of the vector (21).By Proposition 2, see also (15), the j -th component of the vector( τ (1 − t ) I + (1 + t ) A ) − e is given by f j ( t ) = (1 + t ) j − ( τ (1 − t ) + 1 + t ) j , (23)so that r j is the error term of the k -point Gauss-Jacobi formula applied to thecomputation of Z − (1 − t ) α − (1 + t ) − α f j ( t ) dt. (24)We start with the following known result, [11]. Theorem 1
The error term of the k -point Gauss-Jacobi formula applied to thecomputation of Z − (1 − t ) α (1 + t ) α g ( t ) dt, α , α > − , g ∈ C k ([ − , , is given by C k,α ,α g (2 k ) ( ξ ) , − < ξ < , where C k,α ,α := Γ( k + α + 1)Γ( k + α + 1)Γ( k + α + α + 1) k !(2 k + α + α + 1) [Γ(2 k + α + α + 1)] (2 k )! 2 k + α + α +1 . (25) Lemma 1
Let < α < and C k,α := C k,α − , − α . Then C k,α ∼ π − k (2 k )! . Proof.
By (25) we easily obtain C k,α = Γ( k + α )Γ( k − α + 1) [Γ( k )] [Γ(2 k )] (2 k )! 2 k − . Using the Legendre formulaΓ (cid:18) k + 12 (cid:19) Γ( k ) = √ π Γ(2 k )2 − k , we have that C k, / = π − k (2 k )! . a, b ∈ (0 , k b − a Γ( k + a )Γ( k + b ) = 1 + O (cid:18) k (cid:19) , we have that Γ( k + α )Γ( k − α +1) = (cid:2) Γ( k + ) (cid:3) (1+ O ( k − )) and, consequently, C k,α → C k, / as k → ∞ . Remark 1
If we set τ = 1 in (13) we obtain r j = 0 for each j = 1 , , . . . , k since f j ∈ Π j − , see (23). From (16), with p = 1 , and (19) one therefore gets (1 − ζ ) α − − e R k (1 − ζ ) = (1 − ζ ) α − − e p k − (1 − ζ ) e q k (1 − ζ ) = O ( ζ k ) , so that e R k (1 − ζ ) is the ( k − , k ) Pad´e approximant of (1 − ζ ) α − with ex-pansion point ζ = 0 . More generally, if τ ∈ (0 , then the resulting rationalapproximation coincides with the same Pad´e approximant with expansion point ζ = 1 − τ (cf. [5, Lemma 4.4]). Numerically, it is quite clear that the best results are obtained for τ strictlyless than 1, so that in what follows we always assume to work with τ ∈ (0 , j. By Theorem 1, weneed to bound (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) in the interval [ − ,
1] in order to bound the error term E j . We start with the following result.
Proposition 3
Let < τ < and a := 1 + τ − τ . (26) For each j and k max [ − , (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! √ a ( √ a − k +2 (cid:18) a + 12 √ a (cid:19) j . Proof.
The function f j can be written as f j ( t ) = (1 + t ) j − ( a + t ) j (cid:18) a + 12 (cid:19) j , a > . Using the Cauchy integral formula we have f (2 k ) j ( t ) = (2 k )!2 πi Z Γ f j ( w )( w − t ) k +1 dw, (27)where Γ is a contour surrounding t but not the pole − a < −
1. We take Γ asthe circle centered at the origin and radius ρ such that 1 < ρ < a , that is, weuse the substitution w = ρe iθ . We obtainmax [ − , (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! ρ ( ρ − k +1 max [0 , π ] (cid:12)(cid:12) f j ( ρe iθ ) (cid:12)(cid:12) . (28)8aking ρ = √ a we have (cid:12)(cid:12) f j ( ρe iθ ) (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)
11 + ρe iθ (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) ρe iθ a + ρe iθ (cid:12)(cid:12)(cid:12)(cid:12) j (cid:18) a + 12 (cid:19) j = (cid:12)(cid:12)(cid:12)(cid:12)
11 + √ ae iθ (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) √ ae iθ a + √ ae iθ (cid:12)(cid:12)(cid:12)(cid:12) j (cid:18) a + 12 (cid:19) j = (cid:12)(cid:12)(cid:12)(cid:12)
11 + √ ae iθ (cid:12)(cid:12)(cid:12)(cid:12) (cid:18) a + 12 √ a (cid:19) j ≤ √ a − (cid:18) a + 12 √ a (cid:19) j . By (28) we immediately achieve the result.The above result is rather accurate only for small values of j . Since f j ( t ) isgrowing in the interval [ − , t , in order to balance this effect. Proposition 4
Let < τ < and a as in (26). For j ≥ k + 2 (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! √ a + 1 + √ a − t ) j − − k ( a + t ) j + k (cid:18) a + 12 (cid:19) j , t ∈ [ − , . (29) Moreover max [ − , (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! (cid:16) √ a + 1 + √ (cid:17) (cid:16) j − k − k +1 (cid:17) j − − k (cid:16) j +2 k k +1 (cid:17) j + k (cid:18) a + 12 (cid:19) j ( a − − k − (30) for j such that a ≤ j + 6 k + 1 j − k − , (31) and max [ − , (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! √ a + 1 + √ a − a + 1) j − k j +12 + k (32) otherwise. Proof.
In (27) we take Γ as the circle centered at t and of radius ρ suchthat 1 + t < ρ < t + a , that is, we use the substitution w = t + ρe iθ . We obtain (cid:12)(cid:12)(cid:12) f (2 k ) j ( t ) (cid:12)(cid:12)(cid:12) ≤ (2 k )! 1 ρ k max [0 , π ] (cid:12)(cid:12) f j ( t + ρe iθ ) (cid:12)(cid:12) . (33)9or t > − ρ = p ( t + a )(1 + t ), so that (cid:12)(cid:12) f j ( t + ρe iθ ) (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)
11 + t + ρe iθ (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) t + ρe iθ a + t + ρe iθ (cid:12)(cid:12)(cid:12)(cid:12) j (cid:18) a + 12 (cid:19) j = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
11 + t + p ( t + a )(1 + t ) e iθ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) × (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + p ( t + a )(1 + t ) e iθ a + t + p ( t + a )(1 + t ) e iθ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) j (cid:18) a + 12 (cid:19) j = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
11 + t + p ( t + a )(1 + t ) e iθ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:18) ta + t (cid:19) j (cid:18) a + 12 (cid:19) j , and hencemax [0 , π ] (cid:12)(cid:12) f j ( t + ρe iθ ) (cid:12)(cid:12) = 1 p ( t + a )(1 + t ) − (1 + t ) (cid:18) ta + t (cid:19) j (cid:18) a + 12 (cid:19) j = p ( t + a ) + p (1 + t ) p (1 + t )( a − (cid:18) ta + t (cid:19) j (cid:18) a + 12 (cid:19) j . By (33) we obtain the bound (29) for each j and k , for t > −
1. By continuity,(29) holds for t ∈ [ − ,
1] if j ≥ k + 2.Now, we observe that the maximum with respect to t of the function(1 + t ) j − − k ( a + t ) j + k is attained at t ∗ = a ( j − k − − ( j + 2 k )4 k + 1 ≥ − . Moreover t ∗ ≤ a verifying (31). Substituting t ∗ in (29) leads to (30). If t ∗ > t = 1 and hence we obtain (32) . In order to derive (30) and (32) we have assumed τ to be a priori fixed.Now, using these bounds, we look for the value of τ which minimize (cid:13)(cid:13)(cid:13) f (2 k ) j (cid:13)(cid:13)(cid:13) fora given j . For j ≥ k + 3, the minimization of the term( a + 1) j ( a − − k − in (30) leads to a (1) = 2 j + 4 k + 32 j − k − , τ (1) = 4 k + 32 j . (34)On the other side, for the same j ≥ k + 3 , the minimization of the term √ a + 1 + √ a − a + 1) j − k in (32) leads to a value a ∗ ≤ j − k + 2 j − k − , which also satisfies (31) and hence does not fulfill the requirement of (32). For a > a ∗ the bound (32) is growing with a and consequently its minimum isattained just for a (2) = j + 6 k + 1 j − k − . Using this value in (32) leads to a bound that is coarser than the one obtained byreplacing a (1) in (30). For this reason, with respect to our estimates, τ given by(34) represents the optimal value for the computation of (24) and consequentlyof ω ( p ) j .The following theorem summarizes the results obtained. The proof followsstraightfully from (22), Theorem 1, Lemma 1, and Propositions 3, 4. Theorem 2
Let < τ < and a := 1 + τ − τ . Then | E j | ≤ − k sin( απ ) τ α Ψ( a, j, k ) . where Ψ( a, j, k ) := √ a ( √ a − ) k +2 (cid:16) a +12 √ a (cid:17) j , j ≤ k + 1 , (cid:0) √ a + 1 + √ (cid:1) ( j − k − k +1 ) j − − k ( j +2 k k +1 ) j k (cid:0) a +12 (cid:1) j ( a − − k − , k + 2 ≤ j ≤ k +1+ a (2 k +1) a − , √ a +1+ √ a − a +1) j − k j +12 + k , j ≥ max (cid:16) k + 2 , k +1+ a (2 k +1) a − (cid:17) . For τ = τ (1) the corresponding expression of Ψ( a, j, k ) is minimized for j ≥ k + 3 . .1 Numerical experiments As already mentioned, the aim of the whole analysis was to have indicationsabout the choice of the parameter τ with respect to the degree k of the formulaand the dimension of the problem N . Unfortunately, the definition of τ as in(34) depends on j , while we need a value which is as good as possible for each1 ≤ j ≤ N . In this sense, the idea, confirmed by the forthcoming experiments,is to use τ (1) with j = N/
2, that is, focusing the attention on the middle of theinterval [0 , N ]. This leads to a choice of τ around the value 4 k/N . We remarkthat the previous analysis was restricted to the case of p = 1 , because of thedifficulties in dealing with the functions f j for p > τ ∗ = arg min τ E ( τ ) , E ( τ ) := (cid:13)(cid:13)(cid:13) A α − p − e R k ( A p ) (cid:13)(cid:13)(cid:13) ∞ , then, in principle, τ ∗ = τ ∗ ( α, k, N, p ) . However, the numerical experiments doneby using the
Matlab optimization routine fminsearch indicate that the depen-dence on α is negligible with respect to the others. In particular, there isnumerical evidence that E ( τ ∗ ) ≈ E (ˆ τ ) whereˆ τ = (7 + p )2 N k. (35)In Figures 1-2 we report the values of E ( τ ∗ ) and E (ˆ τ ) for p = 1 , , respectively.We recall that the corresponding sets of coefficients { a , a , . . . , a p } in (9) aregiven by p = 1 : { , − } ,p = 3 : { / , − , / , − / } . As one can see, all the curves are approximatively overlapped. A “quasi” optimalapproximation of A α − p can be therefore obtained by using the very simpleformula in (35) for choosing τ. Moreover, it is important to remark that suchapproximations are surely satisfactory even with k ≪ N. The previous results are all related to the overall error in the approximationof A αp . Considering that our final goal is the use of such approximation for thesolution of FDEs, it is important to inspect also the componentwise error. Asan example, in Figure 3, we report such errors, i.e., the values of E j defined in(20), in the case of N = 400 , p = 1 , τ = ˆ τ for α = 0 . , . , .
7, and differentvalues of k . We also consider the componentwise errors of the polynomial ap-proximation of the generating function obtained by truncating its Taylor series,with memory length equal to 16. Obviously this is equivalent to approximatewith 0 the coefficients ω ( p ) i of (5), for i >
16, so that the error is just (cid:12)(cid:12)(cid:12) ω ( p ) i (cid:12)(cid:12)(cid:12) . Thecompetitiveness of the rational approach is undeniable.12
10 1510 −10 −5 α = 1/3, N = 300k 5 10 1510 −10 −5 α = 2/3, N = 300k5 10 1510 −10 −5 α = 1/3, N = 600k 5 10 1510 −10 −5 α = 2/3, N = 600k Figure 1: Error behavior of the Gauss-Jacobi rule for the approximation of A α − for τ = τ ∗ (dashed line) and τ = ˆ τ = 4 k/N (solid line). In this section we discuss the use of the described approximation of A α forgetting a k -step method that simulates the FBDF of order 1. The discreteproblem provided by the latter method applied for solving (1) can be writtenin matrix form as follows( A α ⊗ I s ) ( Y − ⊗ y ) = h α G ( Y ) , (36)where s is the dimension of the FDE, I s is the identity matrix of order s, y ∈ R s represents the initial value, h = ( T − t ) /N is the stepsize, = (1 , , . . . , T ∈ R N ,Y = y y ... y N ≈ y ( t ) y ( t )... y ( t N ) , G ( Y ) = g ( t , y ) g ( t , y )... g ( t N , y N ) ≡ g g ... g N . As described in Section 2, the use of a k -point Gauss-Jacobi rule for approxi-mating (14) leads to A α ≈ β β β k ... . . . 00 . . . . . . 00 β k · · · β − α α α k ... . . . 00 . . . . . . 00 α k · · · α ≡ B − C.
10 1510 −6 −4 −2 α = 1/3, N = 300 k 5 10 1510 −6 −4 −2 α = 2/3, N = 300 k5 10 1510 −6 −4 −2 α = 1/3, N = 600 k 5 10 1510 −6 −4 −2 α = 2/3, N = 600 k Figure 2: Error behavior of the Gauss-Jacobi rule for the approximation of A α − for τ = τ ∗ (dashed line) and τ = ˆ τ = 5 k/N (solid line). −15 −10 −5 α =0.3 0 200 40010 −15 −10 −5 α =0.5 0 200 40010 −15 −10 −5 α =0.7 Figure 3: Componentwise error of the Gauss-Jacobi rational approximation withmemory length k = 6 , ,
12, and the polynomial approximation with k = 16.Here, the coefficients { α j } kj =0 and { β j } kj =0 are related to the rational approxi-mation through the formulas (16)–(18), with m = k since p = 1 ,p k ( ζ ) = (1 − ζ ) e p k − (1 − ζ ) = k X j =0 α j ζ j , q k ( ζ ) = e q k (1 − ζ ) = k X j =0 β j ζ j . (37)If we replace A α by B − C in (36) and we multiply both side of the resultingequation from the left by B ⊗ I s , we obtain( C ⊗ I s ) Y − C ⊗ y = h α ( B ⊗ I s ) G ( Y ) , (38)where now Y represents the numerical solution provided by the k -step method.In fact, considering that ( C ) n = 0 for each n = k + 1 , . . . , N, since p k (1) = k X j =0 α j = 0 , (39)14he discrete problem (38) simplifies to n − X j =0 α j ( y n − j − y ) = h α n − X j =0 β j g n − j , n = 1 , . . . , k, (40) k X j =0 α j y n − j = h α k X j =0 β j g n − j , n = k + 1 , . . . , N. (41)Indeed, the equations in (40) allow to get an approximation of the solution overthe first k meshpoints which are then used as starting values for the k -steprecursion in (41). Remark 2
From (39)-(41) follows that the method reproduces exactly constantsolutions, i.e. it is exact if g ( t, y ( t )) ≡ . As it happens in the case of ODEs, a localization of the zeros of the charac-teristic polynomials of the k -step method in (37) is required in order to study itsstability properties. Clearly, such polynomials depend on the parameter τ, i.e. p k ( ζ ) ≡ p k ( ζ ; τ ) and q k ( ζ ) ≡ q k ( ζ ; τ ) since this dependence occurs in e p k − , e q k . The method is therefore based on the following rational approximation(1 − ζ ) α − ≈ e p k − (1 − ζ ; τ ) e q k (1 − ζ ; τ ) ≡ e R k (1 − ζ ; τ ) . (42) Theorem 3
For each τ ∈ (0 , , the adjoint of the characteristic polynomialsof the k -step method, i.e. ζ k p k ( ζ − ; τ ) and ζ k q k ( ζ − ; τ ) , are a Von Neumannand a Schur polynomial, respectively. Proof.
From Remark 1, one obtains e R k (1 − ζ ; τ ) = τ α − e R k (cid:18) − ζτ ; 1 (cid:19) , (43)since d l dζ l (1 − ζ ) α − (cid:12)(cid:12)(cid:12)(cid:12) ζ =1 − τ = τ α − d l dζ l e R k (cid:18) − ζτ ; 1 (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ζ =1 − τ , l = 0 , , . . . , k − . In addition, using the Gauss hypergeometric functions, in [7, Theorem 4.1] ithas been proved that e R k (cid:18) − ζτ ; 1 (cid:19) = F (1 − k, − α − k ; 1 − k ; ( τ − ζ ) /τ ) F ( − k, α − k ; 1 − k ; ( τ − ζ ) /τ ) , or equivalently, by denoting with P ( γ,β ) r the Jacobi polynomial of degree r andby using [15, eq. 142, p. 464] and the symmetry of such polynomials, e R k (cid:18) − ζτ ; 1 (cid:19) = τ ( τ − ζ ) k − P (1 − α,α ) k − (2 τ / ( τ − ζ ) − τ − ζ ) k P ( α − , − α ) k (2 τ / ( τ − ζ ) − . (44)15rom (37) and (42)–(44), one therefore gets p k ( ζ ; τ ) = (1 − ζ ) τ α ( τ − ζ ) k − P (1 − α,α ) k − (2 τ / ( τ − ζ ) − ,q k ( ζ ; τ ) = ( τ − ζ ) k P ( α − , − α ) k (2 τ / ( τ − ζ ) − . It follows that, if we denote with θ i the i th root of P (1 − α,α ) k − then the roots of p k ( ζ ; τ ) are given by ζ i = 1 + τ − θ i θ i > , i = 1 , . . . , k − , ζ k = 1 , (45)where the inequality follows from the fact that the roots of the Jacobi polynomi-als belong to ( − , . Similarly, by denoting with ϑ i the i th root of P ( α − , − α ) k , one deduces that the roots of q k ( ζ ; τ ) read ζ i = 1 + τ − ϑ i ϑ i > , i = 1 , . . . , k. (46)From (45)-(46) the statement follows immediately.An important consequence of the previous result is that the finite recurrencescheme is always 0-stable independently of the stepnumber k and τ ∈ (0 , . More precisely, in the case of g ≡ In this section we examine the consistency of the method. While, theoretically,it is only exact for constant solutions (see Remark 2), numerically one observesthat the consistency is rather well simulated if k is large enough. The analysiswill also provide some hints about the choice of the memory length m . Werestrict our consideration to the case p = 1 ( m = k ) but the generalization isimmediate.For a given y ( t ), the FBDF of order 1 yields the approximation D αt y ( t ) = 1 h α N X j =0 ( − j (cid:18) αj (cid:19) ( y ( t − jh ) − y (0)) + O ( h ) , t = N h.
Let ∆ αh y ( t ) := N X j =0 ( − j (cid:18) αj (cid:19) ( y ( t − jh ) − y (0)) . Writing a rational approximation of degree k to ω ( α )1 ( ζ ) = (1 − ζ ) α as R k ( ζ ) = ∞ X j =0 γ j ζ j , D αt y ( t ) ≈ h α N X j =0 γ j ( y ( t − jh ) − y (0)) . Denoting by R αk,h y ( t ) := N X j =0 γ j ( y ( t − jh ) − y (0)) , we obtain D αt y ( t ) − h α R αk,h y ( t ) == D αt y ( t ) − h α ∆ αh y ( t ) + 1 h α ∆ αh y ( t ) − h α R αk,h y ( t )= O ( h ) + 1 h α X Nj =0 (cid:20) ( − j (cid:18) αj (cid:19) − γ j (cid:21) ( y ( t − jh ) − y (0)) . The consistency of the method is ensured if1 h α X Nj =0 (cid:20) ( − j (cid:18) αj (cid:19) − γ j (cid:21) ( y ( t − jh ) − y (0)) → h → k < ∞ , in what followswe show that numerically, i.e., for h ≥ h >
0, the consistency is well simulatedif k is large enough and if the rational approximation to A αp is reliable.As pointed out in [12], a certain method for FDEs with generating function ω ( α ) ( ζ ) is consistent of order p if h − α ω ( α ) ( e − h ) = 1 + O ( h p ) . In this setting, in order to understand the numerical consistence of our method,we consider the above relation by replacing ω ( α ) ( e − h ) with ω ( α )1 ( e − h ) and R k ( e − h ).In particular, if we set q k ( h ) = log h (cid:16) h − α (cid:12)(cid:12)(cid:12) R k ( e − h ) − ω ( α )1 ( e − h ) (cid:12)(cid:12)(cid:12)(cid:17) , (48)then we obtain h − α (cid:12)(cid:12) R k ( e − h ) (cid:12)(cid:12) ≤ h − α (cid:12)(cid:12)(cid:12) ω ( α )1 ( e − h ) (cid:12)(cid:12)(cid:12) + h − α (cid:12)(cid:12)(cid:12) R k ( e − h ) − ω ( α )1 ( e − h ) (cid:12)(cid:12)(cid:12) . = 1 + O ( h ) + h q k ( h ) This implies that the consistency of the FBDF of the first order is well simulatedas long as q k ( h ) is larger than 1 . In Figure 4, we plot such function for α = 1 / ,τ = 1 /
10, and different values of k.
100 200 300 400 500 6000123 N = h −1 k = 6k = 8k = 10k = 12 Figure 4: Plot of the function q k ( h ) in (48) for α = 1 / , τ = 1 /
10 and differentvalues of k. The previous experiment does not take care of the perturbation introducedin the approximation of the fractional derivative of fractional powers of the in-dependent variable which may be present in the solution of the FDE. In orderto control such perturbations, we therefore consider the following second ex-periment. Going back to formula (47), we let N = 1 /h and y ( t ) = E α ( − t α )where E α ( x ) denotes the one-parameter Mittag-Leffler function (see e.g. [16,Chapter 1]) E α ( x ) = ∞ X k =0 x k Γ( kα + 1) . (49)In Figure 5, we then consider the behavior of the function e q k ( h ) = log h (cid:18) h − α (cid:12)(cid:12)(cid:12)(cid:12)X Nj =0 (cid:18) ( − j (cid:18) αj (cid:19) − γ j,k (cid:19) ( y ( t N − j ) − y (0)) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (50)which, similarly to q k ( h ) , has to be compared with 1 . The values of y ( t ) havebeen computed using the Matlab function mlf from [17] that implements theMittag-Leffler function together with the Schur-Parlett algorithm.We conclude this section by considering what happens with the general as-sumption | y ( t ) | ≤ M . Using this bound, by (47) we consider the function q k ( h ) = log h (cid:18) h − α X Nj =0 (cid:12)(cid:12)(cid:12)(cid:12) ( − j (cid:18) αj (cid:19) − γ j,k (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (51)whose behavior is reported in Figure 6.As already mentioned, the numerical analysis reported in this section canalso be used to select a proper value for k for a fixed time stepping h or viceversa.Figures 4 and 6 are in fact independent of the problem and can be used easilyto this aim. 18
100 200 300 400 500 600012345 N = h −1 k = 6k = 8k = 10k = 12 Figure 5: Plot of the function e q k ( h ) in (50) for α = 1 / , τ = 1 /
10 and differentvalues of k. −1 k = 8k = 10k = 12k = 16 Figure 6: Plot of the function ¯ q k ( h ) in (51), for α = 1 / , τ = 1 /
10 and differentvalues of k. .2 Linear stability For what concerns the linear stability, taking g ( t, y ( t )) = λy ( t ) in (1), we havethat y ( t ) = E α ( λt α ) → | arg( λ ) − π | < (cid:16) − α (cid:17) π, (see (49) and [13]). The absolute stability region of a FBDF is given by thecomplement of n ω ( α ) p ( ζ ) : | ζ | ≤ o so that a good approximation of the gener-ating function should lead to similar stability domains and hence good stabilityproperties. We consider the behavior of methods based on the Gauss-Jacobirule whose corresponding stability regions are given by, see (17)-(18), C \ (cid:26) p m ( ζ ) q m ( ζ ) : | ζ | ≤ (cid:27) . From a theoretical point of view, from Theorem 3 one deduces that for p = 1such regions are always unbounded for each m = k and τ ∈ (0 , . Indeed, asshown in Figure 7, the methods simulate the behavior of the FBDFs rapidly,i.e. already for k and therefore m small. In particular, the stability domain ofthe method of degree k = m = 12 in the left frame of Figure 7 is very close tothe one of the FBDF of order 1 . Figure 7: Left: boundary of the stability domains of the methods based on theGauss-Jacobi rule with k = 6 (dashed line) and k = 12 (solid line) for p = 1 ,τ = 1 / , and α = 1 /
2. Right: boundary of the stability domains of the samemethods of degree k = 8 for p = 1 (inner) to 4 (outer), τ = 1 /
10 and α = 3 / Numerical examples
As first example, we consider the one-dimensional Nigmatullin’s type equation D αt u ( x, t ) = ∂ u ( x, t ) ∂x , t > , x ∈ (0 , π ) ,u (0 , t ) = u ( π, t ) = 0 ,u ( x,
0) = sin x. If we discretize the spatial derivative by applying the classical central differenceson a uniform mesh of meshsize δ = π/ ( s + 1) , we obtain the s -dimensional FDE D αt y ( t ) = Ly ( t ) , y (0) = y , (52)where L = δ − · tridiag(1 , − , y is the sine function evaluated at theinterior grid points. It is known that y is the eigenvector of L corresponding toits largest eigenvalue λ = − ( δ/ /δ . This implies that the exact solutionof (52) is given by, see (49), y ( t ) = E α ( t α λ ) y . In Figure 8 some results are reported. We compare the maximum norm of theerror at each step of the FBDF of order 1 (FBDF1) and the method based onthe Gauss-Jacobi rule for some values of k and α. The initial values for the k -step schemes are defined according to the strategy described in Section 4. Thereference solutions have been computed using the already mentioned Matlab function mlf from [17]. The dimension of the problem is s = 50, and we con-sider a uniform time step h = 1 /N with N = 250 so that h ≈ δ . As one cansee, if we set τ = 1, i.e. if we use the classical Pad`e approximation of (1 − ζ ) α − (see Remark 1), the k -step methods simulate quite well the FBDF1 initially andan improvement of the results can be obtained by slightly increasing (consid-ering to the total number of integration steps) the stepnumber k. A noticeableimprovement can be obtained by choosing a different value of τ. In particular,if we set τ = ˆ τ = 4 k/N (see (35)) then the 6-step method provides a numericalsolution with the same accuracy of the one provided by the FBDF1 over theentire integration interval.As second example we consider the following nonlinear problem D αt u ( x, t ) = ∂ ( p ( x ) u ( x, t )) ∂x + K α ∂ u ( x, t ) ∂x + ru ( x, t ) (cid:18) − u ( x, t ) K (cid:19) ,u (0 , t ) = u (5 , t ) = 0 , t ∈ [0 , ,u ( x,
0) = x (5 − x ) , x ∈ [0 , . This is a particular instance of the time fractional Fokker-Planck equation with anonlinear source term [18]. In population biology, its solution u ( x, t ) representsthe population density at location x and time t and the nonlinear source termin the equation is known as Fisher’s growth term.21 α = 1/2, τ = 1 α = 1/3, τ = 1 α = 2/3, τ = 1 α = 1/2, τ = 1/2 α = 1/3, τ = 1/2 α = 2/3, τ = 1/2 α = 1/2, τ = 4k/N α = 1/3, τ = 4k/N α = 2/3, τ = 4k/N Figure 8: Step by step error (in logarithmic scale) for the solution of (52) forthe FBDF of order 1 (dashed line) and the method based on the Gauss-Jacobirule with k = 6 (solid line) and k = 12 (dash-dotted line).The application of the classical second order semi-discretization in spacewith stepsize δ = 5 / ( s + 1) leads to the following initial value problem D αt y ( t ) = Jy ( t ) + g ( y ( t )) , t ∈ (0 , , y (0) = y , (53)where, for each i = 1 , . . . , s, ( y ( t )) i ≡ y i ( t ) ≈ u ( iδ, t ) , y i (0) = u ( iδ, , ( g ( y )) i = ry i (1 − y i /K ) , and J is a tridiagonal matrix whose significant entries are J ii = p ′ ( x i ) − K α δ , i = 1 , . . . , s,J i,i − = − p ( x i )2 δ + K α δ , J i − ,i = p ( x i − )2 δ + K α δ , i = 1 , . . . , s − . In our experiment, we set α = 0 . , p ( x ) = − , r = 0 . , K α = K = 1 (see [18,Example 5.4]) and s = 90 . We solved (53) over a uniform meshgrid with stepsize h = 1 /
256 by using the FBDF1 and the 6-step method with τ = 24 / . Theso-obtained numerical solutions have the same qualitative behavior as shown inFigure 9 for different times t = 1 / , / , . This is confirmed by the step by stepmaximum norm of the difference between them reported in Figure 10.22 x FBDF1 x k = 6, τ = 24/N t=1/8t=1/2t=1 Figure 9: Numerical solution of (53) with α = 0 . k = 6 at t = 1 / , / , . In this paper we have presented a new approach for the construction of m -stepformulas for the solution of FDEs. The method shows encouraging results inthe discrete approximation of the FDE solution especially if we consider thecomputational saving with respect to the attainable accuracy. Indeed goodresults are attainable with short memory length. Theoretically the method is0-stable and the consistency is well simulated. The linear stability is preserved.We finally remark that even if the paper only deals with the approximationof FBDFs, the ideas can easily be extended to other approaches such as theFractional Adams type methods. It is just necessary to detect the generatingfunction or the corresponding Toeplitz matrix and then apply the techniquepresented in the paper. Acknowledgements
This work was supported by the GNCS-INdAM 2014 project “Metodi numericiper modelli di propagazione di onde elettromagnetiche in tessuti biologici”.
References [1] D.A. Bini, N.J. Higham, B. Meini, Algorithms for the matrix pth root,Numer. Algorithms 39 (2005), 349–378.23 time