An L1 approximation for a fractional reaction-diffusion equation, a second-order error analysis over time-graded meshes
aa r X i v : . [ m a t h . NA ] S e p AN L APPROXIMATION FOR A FRACTIONALREACTION-DIFFUSION EQUATION, A SECOND-ORDER ERRORANALYSIS OVER TIME-GRADED MESHES ∗ KASSEM MUSTAPHA † Abstract.
A time-stepping L L L Key words.
Fractional diffusion, L
1. Introduction.
Consider the following time-fractional diffusion equation,(1) ∂ t u ( x, t ) + ∂ − αt A u ( x, t ) = f ( x, t ) , for x ∈ Ω and 0 < t < T ,with initial condition u ( x,
0) = u ( x ), where ∂ t = ∂/∂t , Ω is a convex polyhedraldomain in R d ( d ≥ A u ( x, t ) = −∇ · ( κ α ( x ) ∇ u ( x, t )) + d ( x ) u ( x, t ) . The diffusivity coefficient c ≤ κ α ≤ c on Ω for some positive constants c and c , and the reaction coefficient d is such that the bilinear form associated with the ellipticoperator A (see (5)) is positive definite on the Sobolev space H (Ω). That is, it issufficient (but not necessary) to impose that d ≥ κ and d are assumedto be sufficiently regular functions.The fractional exponent is restricted to the range 0 < α < ∂ − αt u = ∂ t I α u , wherethe fractional integration operator I α is defined by I α v ( t ) = Z t ω α ( t − s ) v ( s ) ds, ω α ( t ) = t α − Γ( α ) . We impose a homogeneous Dirichlet boundary condition,(2) u ( x, t ) = 0 for x ∈ ∂ Ω and 0 < t < T .Over the last decade, various time-stepping numerical methods were invesitigatedfor solving the fractional diffusion equation (1), see for example [24, 26] and relatedrefreences therein. The motivation of this paper is propose and analyze a second-orderaccurate time-stepping L t = 0 [18, 22]. Such graded meshes (6) were originally used in the ∗ This work was supported by the KFUPM. † Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals,Dhahran 31261, Saudi Arabia ([email protected]).1
KASSEM MUSTAPHA context of Volterra integral equations with weakly singular kernels, see for example[2, 3, 4], and see also [25] for a recent concrete superconvergence error analysis. Lateron, time-graded meshes were successfully used to improve the performance of differentnumerical methods applied to fractional diffusion and fractional wave equations, seefor instance [19, 23, 24, 27]. Nonuniform meshes are flexible and reasonably convenientfor practical implementation, however they can significantly complicate the numericalerror analysis of schemes. The time-graded mesh properties are carefully used in ourerror analysis to achieve optimal-order convergence rates. The designed approach isnovel and concise, some innovative ideas are employed to estimate efficiently certaincandidates. For completeness, we discretize in space using the standard Galerkinfinite elements, where the error analysis is also examined.To the best of our knowledge, we are not aware of any work that showed asecond-order error bounds of the popular time-stepping L κ to be constant and the reaction coefficient d to be zero):(3) I − α ∂ t u ( x, t ) + A u ( x, t ) = f ( x, t ) , for x ∈ Ω and 0 < t < T ,various types of L − α was proved. Furthermore, the singularity of thecontinuous solution u near t = 0 was taken into account in a few papers only, how-ever the rest frequently ignored this fact. In contrast, a time-stepping discontinuousPetrov-Galerkin method using piecewise polynomials of degree m was introduced andanalyzed in [27] for solving problem (3). For the special case m = 1, this methodreduces to a second-order accurate time-stepping L L u of problem (1) satisfies the following regularity properties:(4) k u ( t ) k ≤ M and k u ′ ( t ) k + t − α/ k u ′′ ( t ) k + t − α/ k u ′′′ ( t ) k ≤ M t σ − , for some positive constants M and σ . In (4), ′ denotes the time partial derivative and k · k ℓ is the norm on the usual Sobolev space H ℓ (Ω) which reduces to the L (Ω)-normwhen ℓ = 0 denoted by k·k . As an example, when f ( t ) ≡ u ∈ H (Ω) ∩ H . − (Ω),these assumptions hold true for σ = α − , see [18, 22] for more details.At each time level t n , an optimal O ( τ t σ + α − /γn )-rate of convergence is proved inTheorem 3.5, assuming that the time mesh exponent γ > max { / ( σ + α/ , / ( σ +3 α/ − / } (see (6) for the definition of the time-graded mesh). Noting that, for1 / ≤ α < σ + α/ ≤ σ + 3 α/ − /
2, and so it is sufficient to assume γ > / ( σ + α/ t n . At the preliminarystage, our error analysis makes use of the inequality in the next lemma [21, Lemma2.3] which will eventually enable us to establish pointwise estimates for certain terms.
N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL Lemma
Let < α ≤ . If the function φ is in the space W (cid:0) (0 , t ); L (Ω) (cid:1) satisfies φ (0) = I α φ ′ (0) = 0 , then k φ ( t ) k ≤ ω − α ( t ) Z t hI α φ ′ ( t ) , φ ′ ( t ) i dt, where h u, v i is the L -inner product on the spatial domain Ω.Although the main scope of this paper is on the optimal error analysis of thetime-stepping L O ( h ) ( h is themaximum spatial mesh element size) is derived assuming that σ > (1 − α ) /
2, seeTheorem 4.1. Numerically, it is observed that this condition is not necessary. In thispart of our error analysis, the next lemma [28, Lemma 3.1] is used.
Lemma
If the functions φ and ψ are in the space L (cid:0) (0 , t ); L (Ω) (cid:1) , then for < α < and for ǫ > , (cid:12)(cid:12)(cid:12)(cid:12)Z t h φ, I α ψ i ds (cid:12)(cid:12)(cid:12)(cid:12) ≤ ǫ (1 − α ) Z t h φ, I α φ i ds + ǫ Z t h ψ, I α ψ i ds . Unfortunately, the coefficient − α ) blows up as α approaches 1 − , and consequently,the error bounds blowup. Such a blowup phenomenon, which was highlighted andinvestigated recently in [6], occurs in the error analysis (but not in numerical experi-ments [20]) of various numerical methods applied to different time-fractional diffusionmodels, see for example [11, 12, 13, 14, 15, 20, 26, 27, 30]. This blowup behavior ap-pears to be an artifact of the method of proof, see Remark 4.2 where the blowup coef-ficient is controlled assuming that σ > (1 − α ) /
2. To validate this, some compatibilityconditions on the initial data u (for example, u ∈ H /α (Ω) with u , A u ∈ H (Ω))and also on the source term f are needed. Noting that, in the limiting case, α → − ,problem (1) reduces to the classical equation (10) and our fully-discrete scheme in(28) amounts to the sandard time-stepping Crank-Nicolson (see (11)) combined withthe (linear) spatial standard continuous Galerkin method. A straightforwrd analysisleads to an optimal time-space second-order convergence rate [31].Finally, in section 5, a second-order convergence of the L O ( τ min { γ ( σ + α ) , } )-rates for different choicesof the time-graded mesh exponent γ and the fractional exponent α . These resultsindicate that the condition γ > max { / ( σ + α/ , / ( σ + 3 α/ − / } in Theorem3.5 is pessimistic. Practically, it is enough to choose γ = 2 / ( σ + α ) to guaranteean O ( τ ) accuracy. Furthermore, the numerical results in Table 4 showed O ( h )-rates of convergence in space for different values of α even though the assumption σ > (1 − α ) / A ( · , · ) : H (Ω) × H (Ω) → R denotes the bilinear form associatedwith the elliptic operator A , which is symmetric and positive definite, defined by(5) A ( v, w ) = h κ α ∇ v, ∇ w i + h d v, w i . Throughout the paper, C is a generic constant which may depend on the parameters M , σ , T , Ω, and γ , but is independent of τ and h . KASSEM MUSTAPHA
2. Numerical method.
This section is devoted to introduce our semi-discretetime-stepping L t i = ( i τ ) γ , for 0 ≤ i ≤ N, for γ ≥ , with τ = T /γ /N, where N is the number of subintervals. Denote by τ n = t n − t n − the length of the n th subinterval I n = ( t n − , t n ), for 1 ≤ n ≤ N . It is not hard to show that such atime-graded mesh has the following properties [19]: for n ≥ , (7) t n ≤ γ t n − , γτ t − /γn − ≤ τ n ≤ γτ t − /γn , τ n − τ n − ≤ C γ τ min(1 , t − /γn ) . For a given function v defined on the time interval [0 , T ] , let v n = v ( t n ) for0 ≤ n ≤ N. With this grid function, we associate the backward difference, ∂v n = v n − v n − τ n . To define our time-stepping numerical scheme, integrating problem (1) over the timeinterval I n ,(8) Z t n t n − u ′ ( t ) dt + Z t n t n − ∂ − αt A u ( t ) dt = Z t n t n − f ( t ) dt . Our L U , which is a continuous linear polynomial in the timevariable on each closed subinterval [ t n − , t n ] , is defined by replacing u with U in (8),(9) U n − U n − + Z t n t n − ∂ − αt A U ( t ) dt = Z t n t n − f ( t ) dt, for 1 ≤ n ≤ N, with U = u . As α → − , the fractional model problem (1) amounts to the classicalreaction-diffusion equation:(10) u ′ ( x, t ) + A u ( x, t ) = f ( x, t ) , for x ∈ Ω and 0 < t < T ,and the time-stepping L U n − U n − + τ n A ( U n + U n − ) / Z t n t n − f ( t ) dt, which is the time-stepping Crank-Nicolson method for problem (10). Motivated bythis, a generalized Crank-Nicolson scheme for the fractional reaction-diffusion equa-tion (1), defined by(12) U n − U n − + Z t n t n − ∂ − αt A U ( t ) dt = Z t n t n − f ( t ) dt, with U = u , was developed in [24], where U ( t ) = ( U j + U j − ) / t ∈ I j . Therein, the theoreticaland numerical convergence results confirmed O ( τ α )-rates in time over sufficientlytime-graded meshes. Both schemes (9) and (12) are computationally similar, howeverthe theoretical and numerical results show better convergence rates of the L N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL ω nj = Z t j t j − ω α ( t n − s ) ds and b ω nj = Z t j t j − Z t j s ω α ( t n − q ) dq ds, for j ≤ n. Hence Z t n t n − ∂ − αt U ( t ) dt = ( I α U )( t n ) − ( I α U )( t n − ) , with( I α U )( t n ) = n X j =1 Z t j t j − ω α ( t n − s ) (cid:16) U j − +( s − t j − ) ∂U j (cid:17) ds = n X j =1 ( ω nj U j − + b ω nj ∂U j ) . Then, the numerical scheme in (9) is equivalent to(13) U n + τ αn Γ( α + 2) A U n = U n − − ατ αn Γ( α + 2) A U n − − n − X j =1 (cid:16) ( ω nj − ω n − ,j ) A U j − + ( b ω nj − b ω n − ,j ) A ∂U j (cid:17) + Z t n t n − f ( t ) dt.
3. Error analysis.
In this section, we study the error bounds from the time-stepping scheme (9). A preliminary estimate will be derived in the next lemma. Forconvenience, putting(14) η ( t ) = η n = 1 τ n Z t n t n − ∂ − αt ( u − ˇ u )( s ) ds, for t ∈ I n , where the piecewise linear polynomial ˇ u interpolates u at the time nodes, that is,ˇ u ( t ) = u j − + ( t − t j − ) ∂u j for t j − ≤ t ≤ t j with 1 ≤ j ≤ N . Lemma
For ≤ n ≤ N, we have k U n − u ( t n ) k ≤ Ct − αn n X j =1 τ j k η j k . Proof.
From equations (8) and (9),( U j − u ( t j )) − ( U j − − u ( t j − )) + Z t j t j − A ∂ − αt ( U − ˇ u )( t ) dt = τ j A η j . Taking the inner product with v := R t j t j − ∂ − αt ( U − ˇ u ) dt = R t j t j − I α ( U − ˇ u ) ′ dt (because U − ˇ u (0) = U − u = 0), and using A (cid:16) Z t j t j − I α ( U − ˇ u ) ′ dt, Z t j t j − I α ( U − ˇ u ) ′ dt (cid:17) ≥ β (cid:13)(cid:13)(cid:13) Z t j t j − I α ( U − ˇ u ) ′ dt (cid:13)(cid:13)(cid:13) , for some positive constant β depends on Ω (due to the Poincar´e inequality), τ j Z t j t j − h ( U − ˇ u ) ′ , I α ( U − ˇ u ) ′ i dt + β (cid:13)(cid:13)(cid:13) Z t j t j − I α ( U − ˇ u ) ′ dt (cid:13)(cid:13)(cid:13) ≤ τ j (cid:28) A η j , Z t j t j − I α ( U − ˇ u ) ′ dt (cid:29) . KASSEM MUSTAPHA
An application of the Cauchy-Schwarz inequality leads to τ j (cid:28) A η j , Z t j t j − I α ( U − ˇ u ) ′ dt (cid:29) ≤ β τ j k η j k + β (cid:13)(cid:13)(cid:13) Z t j t j − I α ( U − ˇ u ) ′ dt (cid:13)(cid:13)(cid:13) , and consequently,(15) Z t j t j − h ( U − ˇ u ) ′ , I α ( U − ˇ u ) ′ i dt ≤ Cτ j k η j k . Summing over the variable j , Z t n h ( U − ˇ u ) ′ , I α ( U − ˇ u ) ′ i dt ≤ C n X j =1 τ j k η j k . Hence, using Lemma 1.1, the desired bound is obtained.The current task is to estimate the candidate k η j k which is very delicate. Theapproach is novel and some new ideas are used. Starting from the fact that ( u − ˇ u )(0) =0 (since ˇ u interpolates u at the time nodes t n for 0 ≤ n ≤ N ), we observe ∂ − αt ( u − ˇ u )( t ) = I α ( u − ˇ u ) ′ ( t ) = j X i =1 Z min { t i ,t } t i − ω α ( t − s ) (cid:16) u ′ ( s ) − ∂u i (cid:17) ds, for t ∈ I j . However, u ( t i ) − u ( t i − ) = τ i u ′ ( t i ) − τ i u ′′ ( t i ) + 12 Z t i t i − ( q − t i − ) u ′′′ ( q ) dq, and hence, after some manipulations, one can show that u ′ ( s ) − ∂u i = e ( s ) + e ( s ) for s ∈ I i , where e ( s ) = Z t i s ( q − s ) u ′′′ ( q ) dq − τ i Z t i t i − ( q − t i − ) u ′′′ ( q ) dq,e ( s ) = (cid:0) s − t i − / (cid:1) u ′′ ( t i ) . Thus, splitting η j as: η j = τ − j (cid:16) η j + η j (cid:17) , where(16) η j = Z t j t j − I α e ( t ) dt and η j = Z t j t j − I α e ( t ) dt. Estimating k η j k will be the topic of the next lemma. Lemma
For ≤ j ≤ N, we have k η j k ≤ Cτ j τ t α/ σ − − /γj , with γ > / ( σ + α/ . N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL Proof.
Expanding η j as η j = j X i =1 Z t j t j − Z min { t i ,t } t i − ω α ( t − s ) e ( s ) ds dt. From the definition of e and the regularity assumption (4), for s ∈ I , (17) k e ( s ) k ≤ Z t s ( q − s ) k u ′′′ ( q ) k dq + 12 t Z t q k u ′′′ ( q ) k dq ≤ M Z t s q q σ + α/ − dq + M t Z t q q σ + α/ − dq ≤ C max { s σ + α/ − , t σ + α/ − } . Hence, for j = 1 , (18) k η k ≤ C Z t Z t ω α ( t − s ) max { s σ + α/ − , t σ + α/ − } ds dt = C Z t max n Γ( σ + α/ σ + 3 α/ t σ +3 α/ − , ω α +1 ( t ) t σ + α/ − o dt = C max n Γ( σ + α/ σ + 3 α/ t σ +3 α/ , ω α +2 ( t ) t σ + α/ − o ≤ Cτ α/ σ . For the case j ≥ , noting first that k η j k ≤ Z t j t j − (cid:16) Z t ω α ( t − s ) k e ( s ) k ds + j X i =2 Z min { t i ,t } t i − ω α ( t − s ) k e ( s ) k ds (cid:17) dt. For estimating the first term on the right-hand side in the above inequality, using t − s ≥ t j − − s = ( j − γ (cid:16) t − s ( j − γ (cid:17) ≥ ( j − γ (cid:16) t − s (cid:17) ≥ ( j/ γ (cid:16) t − s (cid:17) , and the achieved bound in (17), Z t ω α ( t − s ) k e ( s ) k ds ≤ Cj γ ( α − Z t ( t − s ) α − max { s σ + α/ − , t σ + α/ − } ds = Ct α − j τ γ (1 − α ) Z t ( t − s ) α − max { s σ + α/ − , t σ + α/ − } ds ≤ C t α − j τ γ (1 − α ) t α/ σ − = C τ t α − j t α/ σ − /γ ≤ C τ t α/ σ − /γ − j , for γ ≥ / ( σ + α/ . On the other hand, for s ∈ I i with i ≥ , from the definition of e , the regularityassumption (4), and the time mesh property (7), we have k e ( s ) k ≤ Cτ i Z t i t i − k u ′′′ ( q ) k dq ≤ Cτ i Z t i t i − q σ + α/ − dq ≤ Cτ i t σ + α/ − i ≤ Cτ t σ + α/ − − /γi ≤ Cτ s σ + α/ − − /γ . KASSEM MUSTAPHA
Thus, j X i =2 Z min { t i ,t } t i − ω α ( t − s ) k e ( s ) k ds ≤ Cτ Z tt ( t − s ) α − s σ + α/ − − /γ ds ≤ Cτ Z t ( t − s ) α − s σ + α/ − − /γ ds ≤ Cτ t α/ σ − − /γ , for γ > / ( σ + α/ . Gathering the above contribution and using again the time mesh property (7), k η j k ≤ Cτ α/ σ + Cτ j τ t α/ σ − − /γj , for j ≥ γ > / ( σ + α/ . From this bound, and the achieved estimate in (18), the desired result is obtained.It remains to estimate k η j k . The technical result in the next lemma is needed. Lemma
For γ ≥ , and for a given positive sequence { a i } , we have j X i =2 a i (cid:12)(cid:12)(cid:12) L i − ,j − − τ i − τ i L i,j (cid:12)(cid:12)(cid:12) ≤ Cτ t α − /γj − j max i =2 (cid:16) a i τ i − (cid:17) , for 2 ≤ i ≤ j, with L i,j := 12 Z t i t i − ( s − t i − )( t i − s ) ω α ( t j − s ) ds, for 1 ≤ i ≤ j ≤ N. Proof.
For γ = 1 , L i − ,j − − τ i − τ i L i,j = 0, and so, we have nothing to show. For γ > , from the substitution s = τ − i (cid:16) ( t i − q ) t i − + ( q − t i − ) t i − (cid:17) , we observe(19) L i − ,j − = 12 τ i − τ i Z t i t i − ( q − t i − )( t i − q ) ω α ( t j − − s ) dq. Since t j − − s ≤ t j − q , L i − ,j − ≥ τ i − τ i L i,j . This leads to j X i =2 a i (cid:12)(cid:12)(cid:12) L i − ,j − − τ i − τ i L i,j (cid:12)(cid:12)(cid:12) = 12 j X i =2 a i τ i − τ i Z t i t i − ( q − t i − )( t i − q )[ ω α ( t j − − s ) − ω α ( t j − q )] dq ≤ j max i =2 (cid:16) a i τ i − τ i (cid:17) j X i =2 Z t i t i − [ ω α ( t j − − s ) − ω α ( t j − q )] dq ≤ j max i =2 (cid:16) a i τ i − (cid:17) j X i =2 (cid:18) τ i τ i − Z t i − t i − ω α ( t j − − v ) dv − Z t i t i − ω α ( t j − q ) dq (cid:19) . A simple manipulation shows that j X i =2 Z t i t i − ω α ( t j − q ) dq = Z t j t ω α ( t j − q ) dq = ω α +1 ( t j − t ) ≥ ω α +1 ( t j − ) = j X i =2 Z t i − t i − ω α ( t j − − v ) dv, N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL j X i =2 a i (cid:12)(cid:12)(cid:12) L i − ,j − − τ i − τ i L i,j (cid:12)(cid:12)(cid:12) ≤ j max i =2 (cid:16) a i τ i − (cid:17) j X i =2 (cid:16) τ i τ i − − (cid:17) Z t i − t i − ω α ( t j − − q ) dq. By the mesh properties in (7), τ i τ i − − τ i − τ i − ) τ − i − ≤ Cτ t − /γi τ − t /γ − i − ≤ Cτ t − /γi − , for i ≥ , and therefore, using γ > , j X i =2 (cid:16) τ i τ i − − (cid:17) Z t i − t i − ω α ( t j − − q ) dq ≤ Cτ j X i =2 t − /γi − Z t i − t i − ω α ( t j − − q ) dq ≤ Cτ Z t j − q − /γ ω α ( t j − − q ) dq ≤ Cτ t α − /γj − . This completes the proof.Now, we are ready to bound k η j k . Lemma
For η j defined as in (16) with j ≥ , we have k η j k ≤ Cτ τ j t α/ σ − − /γj , for γ ≥ / ( σ + α/ . Proof.
Splitting η j follows by reversing the order of integration then integratingby parts, η j = j X i =1 u ′′ ( t i ) Z t j t j − Z min { t i ,t } t i − ( s − t i − / ) ω α ( t − s ) ds dt = j − X i =1 u ′′ ( t i ) Z t i t i − ( s − t i − / )[ ω α +1 ( t j − s ) − ω α +1 ( t j − − s )] ds + u ′′ ( t j ) Z t j t j − ( s − t j − / ) ω α +1 ( t j − s ) ds = j − X i =1 u ′′ ( t i )[ L i,j − − L i,j ] − u ′′ ( t j ) L j,j . Thus, η j can be decomposed as: η j = − η j , − η j , − η j , , with η j , = j X i =2 [ u ′′ ( t i ) − u ′′ ( t i − )] L i,j , (20) η j , = j X i =2 u ′′ ( t i − ) (cid:18) − τ i − τ i (cid:19) L i,j , (21) η j , = u ′′ ( t ) L ,j − j X i =2 u ′′ ( t i − ) (cid:18) L i − ,j − − τ i − τ i L i,j (cid:19) . (22)0 KASSEM MUSTAPHA
By the regularity assumption in (4) and the time mesh properties in (7),(23) k η j , k ≤ C j X i =2 τ i Z t i t i − k u ′′′ ( q ) k dq Z t i t i − ( t j − s ) α − ds ≤ C j X i =2 τ i t σ + α/ − i Z t i t i − ( t j − s ) α − ds ≤ Cτ j X i =2 t σ + α/ − /γi Z t i t i − ( t j − s ) α − ds ≤ Cτ Z t j t s σ + α/ − /γ ( t j − s ) α − ds ≤ Cτ t σ +3 α/ − /γj , for γ ≥ / ( σ + α/ . The next task is to estimate η j , . Seeing that1 − τ i − /τ i ≤ τ − i ( τ i − τ i − ) ≤ Cτ τ − i t − /γi (the third time mesh property in (7) is used here), yields(24) k η j , k ≤ Cτ j X i =2 τ − i t − /γi k u ′′ ( t i − ) k τ i Z t i t i − ( t j − s ) α − ds ≤ Cτ j X i =2 Z t i t i − ( t j − s ) α − s σ + α/ − /γ ds ≤ Cτ Z t j t ( t j − s ) α − s σ + α/ − /γ ds ≤ Cτ t σ +3 α/ − /γj , for γ ≥ / ( σ + α/ . To estimate η j , , we use again the regularity assumption in (4) and then apply Lemma3.3 with t σ − i − in place of a i , and get j X i =2 k u ′′ ( t i − ) k (cid:12)(cid:12)(cid:12) L i − ,j − − τ i − τ i L i,j (cid:12)(cid:12)(cid:12) ≤ C j X i =2 t σ − i − (cid:12)(cid:12)(cid:12) L i − ,j − − τ i − τ i L i,j (cid:12)(cid:12)(cid:12) ≤ Cτ t α − /γj − j max i =2 (cid:16) t σ + α/ − i − τ i − (cid:17) ≤ Cτ t α − /γj − j max i =2 (cid:16) t σ + α/ − /γi − (cid:17) ≤ Cτ t α/ σ − /γj − , for γ ≥ / ( σ + α/ . By the definition of η j , and the above contribution, we obtain(25) k η j , k ≤ k u ′′ ( t ) k L ,j + Cτ t α/ σ − /γj − ≤ Cτ σ + α/ Z t ω α ( t j − s ) ds + Cτ t α/ σ − /γj − ≤ Cτ α/ σ + Cτ t α/ σ − /γj − , for γ ≥ / ( σ + α/ . Therefore, to complete the proof, combining the estimates from (23), (24) and (25),and use the inequality τ ≤ Cτ j t /γ − j (follows from the second mesh property in (7)). N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL γ is not sharp. Theorem
Let U be the time-stepping solution defined by (9) and let u bethe solution of the fractional reaction-diffusion problem (1) . Assume that u satisfiesthe regularity assumptions in (4) . If the time mesh exponent γ is greater than themaximum of { / ( σ + α/ , / ( σ + 3 α/ − / } with σ + 3 α/ − / > , then k U n − u ( t n ) k ≤ Cτ t σ + α − /γn ≤ Cτ , for 1 ≤ n ≤ N .
Proof.
From the decomposition η j = τ − j (cid:0) η j + η j (cid:1) , and the established boundsof η j and η j in Lemma 3.2 and Lemma 3.4, respectively,(26) n X j =1 τ j k η j k ≤ Cτ n X j =1 τ j t α/ σ − − /γ ) j ≤ Cτ Z t n t t α/ σ − − /γ ) dt ≤ Cτ max { t σ +3 α − /γ − , t σ +3 α − /γ − n } , for γ > / ( σ + α/ . Inserting this in the achieved bound in Lemma 3.1 and using γ ≥ / ( σ + 3 α/ − /
4. Fully-discrete solution.
To compute our numerical solution, we thereforeseek a fully-discrete solution U h by discretizing (9) in space via the standard Galerkinfinite element method. To this end, let T h be a family of regular (conforming) tri-angulation of the domain Ω and let h = max K ∈T h (diam K ) , where h K denotes thediameter of the element K. Let V h ⊂ H (Ω) denote the usual space of continuous,piecewise-linear functions on T h that vanish on ∂ Ω. Let W ( V h ) ⊂ C ([0 , T ]; V h ) denotethe space of linear polynomials on [ t n − , t n ] for 1 ≤ n ≤ N , with coefficients in V h . Taking the inner of (9) with a test function χ ∈ H (Ω), and apply the first Greenidentity. Then, the semi-discrete L U satisfies(27) h U n − U n − , χ i + Z t n t n − A ( ∂ − αt U ( t ) , χ ) dt = Z t n t n − h f ( t ) , χ i dt, with U = u , Motivated by this, our fully-discrete computational solution U h ∈ W ( V h ) is definedas: for 1 ≤ n ≤ N, (28) h U nh − U n − h , v h i + Z t n t n − A ( ∂ − αt U h ( t ) , v h ) dt = Z t n t n − h f ( t ) , v h i dt ∀ v h ∈ V h , with U h = R h u , where R h : H (Ω) → V h is the Ritz projection defined by A ( R h w, v h ) = A ( w, v h ) , ∀ v h ∈ V h . Following the derivation used to obtain (13), the scheme in (28) is equivalent to h U nh , v h i + τ αn Γ( α + 2) A ( U nh , v h ) = h U n − h , v h i − ατ αn Γ( α + 2) A ( U n − h , v h ) − A (cid:18) n − X j =1 (cid:16) ( ω nj − ω n − ,j ) U j − h + ( b ω nj − b ω n − ,j ) ∂U jh (cid:17) , v h (cid:19) + Z t n t n − h f ( t ) , v h i dt . KASSEM MUSTAPHA
For 1 ≤ p ≤ d h := dim V h , let φ p ∈ V h denote the p th basis function associated withthe p th interior node ~x p , so that φ p ( ~x q ) = δ pq and U nh ( ~x ) = d h X p =1 u nh ( ~x p ) φ p ( ~x ) . We define d h × d h matrices: M = [ h φ q , φ p i ] , G = [ A ( φ q , φ p )] , and the d h -dimensionalcolumn vectors U nh and F n with components U nh ( ~x p ) and R t n t n − h f ( t ) , φ p i dt, respec-tively. Therefore, the fully-discrete scheme (28) has the following matrix representa-tions: (cid:16) M + τ αn Γ( α + 2) G (cid:17) U nh = (cid:16) M − ατ αn Γ( α + 2) G (cid:17) U n − h − n − X j =1 (cid:16)(cid:0) ω nj − ω n − ,j (cid:1) G U j − h + ( b ω nj − b ω n − ,j ) G ∂ U jh (cid:17) + F n . Therefore, at each time level t n , the numerical scheme (28) reduces to a finitesquare linear system, and so the existences of U nh follows from its uniqueness. Thelatter follows from the fact that both matrices M and G are positive definite.Turning now into the error analysis, we introduce the following notations: θ ( t ) = U h ( t ) − R h ˇ u ( t ) and ρ ( t ) = u ( t ) − R h u ( t ) , and . Since ˇ u interpolates u at the time nodes, θ n = U nh − R h ˇ u ( t n ) = U nh − R h u ( t n ) . Thence,the pointwise time error U nh − u ( t n ) can be decomposed as(29) U nh − u ( t n ) = [ U nh − R h u ( t n )] − [ u ( t n ) − R h u ( t n )] = θ n − ρ n . The estimate of the second term follows easily from the Ritz projector approxi-mation property and the first regularity assumption in (4),(30) k ρ ( t n ) k ≤ Ch k u ( t n ) k ≤ Ch , for 0 ≤ n ≤ N .The next duty is to estimate θ n . From the weak formulation of problem (1); h u ( t j ) − u ( t j − ) , χ i + Z t j t j − A ( ∂ − αt u ( t ) , χ ) dt = Z t j t j − h f ( t ) , χ i dt ∀ χ ∈ H (Ω) , the numerical scheme (28), and the decomposition in (29), we have τ j h ∂θ j , v h i + Z t j t j − A (cid:0) ∂ − αt ( U h − ˇ u )( t ) , v h (cid:1) dt = τ j h ∂ρ n , v h i + Z t j t j − A (cid:0) ∂ − αt ( u − ˇ u )( t ) , v h (cid:1) dt, ∀ v h ∈ V h . From the orthogonality property of the Ritz projection, and the definition of η in (14),(31) τ j h ∂θ j , v h i + Z t j t j − A (cid:0) ∂ − αt θ ( t ) , v h (cid:1) dt = τ j [ h ∂ρ j , v h i + A ( η j , v h )] , ∀ v h ∈ V h . N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL θ = U h − R h ˇ u (0) = R h u − R h u = 0, R t j t j − ∂ − αt θ ( t ) dt = R t j t j − I α θ ′ ( t ) dt. Now,setting v h = R t j t j − I α θ ′ ( t ) dt and applying the Poincar´e inequality, then the secondterm in (31) is ≥ β k R t j t j − I α θ ′ ( t ) dt k , for some positive constant β depends on Ω . This and the fact that ∂θ j = θ ′ ( t ) (constant in time) for t ∈ I j , lead to τ j Z t j t j − h θ ′ , I α θ ′ i dt + β (cid:13)(cid:13)(cid:13) Z t j t j − I α θ ′ dt (cid:13)(cid:13)(cid:13) ≤ τ j (cid:28) ∂ρ j , Z t j t j − I α θ ′ dt (cid:29) + τ j A (cid:16) η j , Z t j t j − I α θ ′ dt (cid:17) . By the Cauchy-Schwarz inequality, the last term is ≤ β τ j k η j k + β (cid:13)(cid:13)(cid:13) Z t j t j − I α θ ′ dt (cid:13)(cid:13)(cid:13) , and consequently,(32) τ j Z t j t j − h θ ′ , I α θ ′ i dt + β (cid:13)(cid:13)(cid:13) Z t j t j − I α θ ′ dt (cid:13)(cid:13)(cid:13) ≤ τ j Z t j t j − h ˇ ρ ′ , I α θ ′ i dt + Cτ j k η j k , where ˇ ρ ( t ) = ρ j − + ( t − t j − ) ∂ρ j for t ∈ I j . Dividing both sides by τ j , and then,summing over the variable j and using the inequality(33) Z t n h ˇ ρ ′ , I α θ ′ i dt ≤ − α ) Z t n h ˇ ρ ′ , I α ˇ ρ ′ i dt + 12 Z t n h θ ′ , I α θ ′ i dt, (Lemma 1.2 is used here), we reach Z t n h θ ′ , I α θ ′ i dt ≤ − α ) Z t n h ˇ ρ ′ , I α ˇ ρ ′ i dt + Cτ j k η j k . Thanks to Lemma 1.1,(34) k θ n k ≤ Ct − αn (cid:16) Z t n |h ˇ ρ ′ , I α ˇ ρ ′ i| dt + n X j =1 τ j k η j k (cid:17) . To estimate R t n |h ˇ ρ ′ , I α ˇ ρ ′ i| dt , split it as (recall that ˇ ρ ′ ( t ) = ∂ρ j on I j )(35) Z t n |h ˇ ρ ′ ( t ) , I α ˇ ρ ′ ( t ) i| dt ≤ k ∂ρ k Z t Z t ω α ( t − s ) dsdt + n X j =2 k ∂ρ j k k ∂ρ k Z t j t j − Z t ω α ( t − s ) dsdt + n X j =2 k ∂ρ j k j X i =2 k ∂ρ i k Z t j t j − Z min { t i ,t } t i − ω α ( t − s ) dsdt . By the definition of the function ρ , the Ritz projection error bound in (30) with u ′ inplace of u , and the regularity assumption (4), we obtain(36) k ∂ρ j k = τ − j (cid:13)(cid:13)(cid:13) Z t j t j − ( R h u ′ − u ′ )( s ) ds (cid:13)(cid:13)(cid:13) ≤ Ch τ − j Z t j t j − k u ′ ( s ) k ds ≤ Ch τ − j Z t j t j − s σ − ds, for j ≥ . KASSEM MUSTAPHA If σ − ≥ k ∂ρ j k ≤ Ch t σ − j . Thus,(37) Z t n |h ˇ ρ ′ , I α ˇ ρ ′ i| dt ≤ Ch t σ − n Z t n Z t ω α ( t − s ) dsdt ≤ Ch t σ + α − n , for σ ≥ . Now, turning into the case σ − <
0, which is probably more interesting. Assumingthat σ > (1 − α ) /
2, a similar bound will be achieved next, see (38). Using (36), thefirst term on the right-hand side of (35) is bounded by Ch τ − (cid:16) Z t s σ − ds (cid:17) ω α +2 ( t ) ≤ Ch t σ + α − . Using (36) and the first time mesh property in (7), the second candidate on theright-hand side of (35) is ≤ Ch n X j =2 t σ − j τ Z t s σ − ds Z t j t j − Z t ω α ( t − s ) ds dt ≤ Ch n X j =2 ( t t j ) σ − Z t j t j − Z t ω α ( t − s ) ds dt ≤ Ch n X j =2 Z t j t j − t σ − Z t s σ − ω α ( t − s ) ds dt ≤ Ch Z t n t t σ − Z t s σ − ω α ( t − s ) ds dt ≤ Ch Z t n t t σ + α − dt, while, the last term in (35) is ≤ Ch n X j =2 t σ − j j X i =2 τ i Z I i s σ − ds Z t j t j − Z min { t i ,t } t i − ω α ( t − s ) ds dt ≤ Ch n X j =2 t σ − j j X i =2 t σ − i Z t j t j − Z min { t i ,t } t i − ω α ( t − s ) ds dt ≤ Ch n X j =2 j X i =1 Z t j t j − t σ − Z min { t i ,t } t i − s σ − ω α ( t − s ) ds dt ≤ Ch Z t n t t σ − Z t s σ − ω α ( t − s ) ds dt ≤ Ch Z t n t t σ + α − dt . Therefore, gathering the above estimates, we conclude that(38) Z t n |h ˇ ρ ′ , I α ˇ ρ ′ i| dt ≤ Ch t σ + α − n , for (1 − α ) / < σ < . From the decomposition (29), the Ritz projection in (30), the inequality in (34),the achieved bounds in (26) and (37), and the above estimate, the error result in thenext convergence theorem holds true. It is claimed that for a sufficiently time-gradedmesh, the proposed fully-discrete scheme is second-order accurate in both time andspace. The numerical results in the forthcoming section confirm that the imposedassumption on the time mesh exponent γ is pessimistic. Furthermore, these results N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL O ( h )-rates of convergence in space, although the imposed condition σ > (1 − α ) / O ( h )-rate of convergence was carriedout without this assumption [13]. Theorem
Let U h be the numerical solution defined by (28) and let u bethe solution of the fractional reaction-diffusion problem (1) . Assume that u satisfiesthe regularity assumptions in (4) with σ > (1 − α ) / . If the time mesh exponent γ > max { / ( σ + α/ , / ( σ + 3 α/ − / } , then k U nh − u ( t n ) k ≤ C ( τ + h ) , for 1 ≤ n ≤ N. We end this section with the following remark.
Remark C in (34)blows up as α → − . To control this phenomena, Lemma 1.2 (and consequently, theinequality in (33)) should be avoided. Since ˇ ρ ′ ( t ) = ∂ρ j for t ∈ I j , an application ofthe Cauchy-Schwarz inequality yields τ j Z t j t j − h ˇ ρ ′ , I α θ ′ i dt = (cid:28) τ j ∂ρ j , Z t j t j − I α θ ′ dt (cid:29) ≤ Cτ j k ∂ρ j k + β (cid:13)(cid:13)(cid:13)(cid:13) Z t j t j − I α θ ′ ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) . Substitute this in (32) gives Z t j t j − h θ ′ , I α θ ′ i dt ≤ Cτ j k ∂ρ j k + Cτ j k η j k . Summing over j , follows by using Lemma 1.1 with θ in place of φ , we notice that k θ n k ≤ Ct − αn (cid:16) n X j =1 τ j k ∂ρ j k + n X j =1 τ j k η j k (cid:17) , where the constant C in the above bound does not blowup as α → − . The remaining exercise is to estimate P nj =1 τ − j k ∂ρ j k . From (36), n X j =1 τ j k ∂ρ j k ≤ Ch n X j =1 τ − j (cid:16) Z t j t j − s σ − ds (cid:17) ≤ Ch n X j =1 Z t j t j − s σ − ds = Ch Z t n s σ − ds ≤ Ch t σ − n , for σ > / . For σ = 1 / , these steps can be slightly adjusted to show an O ( h ) bound of theabove term, but with a logarithmic coefficient log( t n /t ) for n ≥ .
5. Numerical results.
To support the achieved theoretical convergence resultsin Theorems 3.5 and 4.1, this section is devoted to perform some numerical experi-ments (on a typical test problem). In the fractional model problem (1), we choose u ( x, y ) = x (1 − x ), κ α = d = 1, and f = 0 . The time and space domains are chosento be the intervals [0 ,
1] and (0 , , respectively. Separation of variables yields theseries representation of the solution:(39) u ( x, t ) = 8 ∞ X m =0 λ − m sin( λ m x ) E α ( − λ m t α ) , λ m := (2 m + 1) π − , KASSEM MUSTAPHA
M γ = 1 γ = 2 γ = 3 γ = 420 3.40e-02 1.09e-02 2.15e-03 8.02e-0440 2.78e-02 0.292 5.16e-03 1.083 7.05e-04 1.607 2.13e-04 1.91180 2.19e-02 0.346 2.34e-03 1.142 2.46e-04 1.518 5.57e-05 1.935160 1.65e-02 0.402 1.10e-03 1.085 8.66e-05 1.509 1.44e-05 1.951320 1.20e-02 0.458 5.44e-04 1.019 3.05e-05 1.505 3.70e-06 1.962640 8.48e-03 0.506 2.71e-04 1.001 1.08e-05 1.503 9.44e-07 1.972 Table 1
Errors and convergence rates ( r t ) for α = 0 . and for different choices of γ. M γ = 1 γ = 2 γ = 2 . γ = 320 2.51e-02 2.21e-03 8.26e-04 6.50e-0440 1.50e-02 0.748 7.60e-04 1.540 2.12e-04 1.964 1.65e-04 1.97980 8.32e-03 0.846 2.66e-04 1.513 5.44e-05 1.960 4.17e-05 1.982160 4.53e-03 0.879 9.36e-05 1.509 1.46e-05 1.897 1.05e-05 1.987320 2.53e-03 0.840 3.30e-05 1.505 3.93e-06 1.893 2.65e-06 1.991640 1.47e-03 0.779 1.16e-05 1.503 1.06e-06 1.892 6.64e-07 1.995 Table 2
Errors and convergence rates ( r t ) for α = 0 . and for different choices of γ. where E α ( t ) := P ∞ p =0 t p Γ( αp +1) is the Mittag-Leffler function.The initial data u ∈ ˙ H . − (Ω) ∩ H (Ω). Thus, as expected from the regularityanalysis in [18, 22], the regularity properties in (4) hold true for σ = α − / . For the numerical illustration of the convergence rates from the time-stepping L h so that the time errors aredominant. Therefore, by Theorems 3.5 and 4.1, we expect to observe O ( τ )-rates ofconvergence for γ > max { / ( σ + α/ , / ( σ + 3 α/ − / } = max { / (3 α − ) , / (7 α − − } , with σ + 3 α/ − / >
0. However, the results in Tables 1–3 are more optimistic, O ( τ )-rates were observed for γ ≥ / ( σ + α ) = 8 / (5 α − ), for different values of α. Moreover, these results confirm that the condition σ +3 α/ − / > u in (39) of problem (1)by truncating the Fourier series in (39) after 60 terms. To measure the error in thenumerical solution, we computed E N,M := max ≤ n ≤ N k U nh − u ( t n ) k , where N is the number of time subintervals, while M is the number of uniform spacemesh elements. Noting that, the spatial L -norm was evaluated using the two-pointGauss quadrature rule on the finest spatial mesh. The convergence rates r t (in time)and r x (in space) were calculated from the relations r t ≈ log (cid:16) E N,M /E N,M (cid:17) , when h r x ≪ τ r t ,r x ≈ log (cid:16) E N,M /E N, M (cid:17) , when τ r t ≪ h r x . For the graphical interpretation, we fixed N = 160 and M = 1200, so the timeerror is dominant. Figure 1 shows how the error on uniform and nonuniform timemeshes varies with t for various choices of α , using a log scale. N L1 METHOD FOR A FRACTIONAL DIFFUSION MODEL M γ = 1 γ = 1 . γ = 220 9.7410e-03 1.8813e-03 5.4226e-0440 4.2516e-03 1.1961 6.7825e-04 1.4719 1.3563e-04 1.999380 2.0954e-03 1.0207 2.3983e-04 1.4998 3.4030e-05 1.9948160 1.0687e-03 9.7134 8.4803e-05 1.4998 8.5585e-06 1.9914320 5.3632e-04 9.9475 2.9991e-05 1.4996 2.1776e-06 1.9746 Table 3
Errors and convergence rates ( r t ) for α = 0 . and for different choices of γ. t -7 -6 -5 -4 -3 -2 E rr o r s = 1 = 4 t -7 -6 -5 -4 -3 -2 E rr o r s = 1 = 2.7 t -7 -6 -5 -4 -3 E rr o r s = 1 = 2 Fig. 1 . The error k U nh − u ( t n ) k as a function of t n . The fractional exponent α = 0 . in theleft figure, while α = 0 . in the middle one, and α = 0 . in the right figure. To demonstrate the O ( h )-rates from the spatial discretization by Galerkin finiteelements, the time mesh size is refined so that the errors in space are dominant. Theexpected convergence orders are displayed in Table 4 for α = 0 . , . , and 0 . . Theseresults also illustrate that the condition σ > (1 − α ) / α > / σ = α − / , however an O ( h )-rate wasobserved despite α not being greater than 2 / .
6. Concluding remarks. An L L γ can be further relaxed. Due to several difficulties, improving thechoice of γ is beyond the scope of this work, it will be a subject of future research. REFERENCES[1] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput.Phys., 280 (2015), 424–438.[2] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations Meth-ods. Cambridge University Press, Cambridge (2004).[3] H. Brunner, A. Pedas, and Vainikko, The piecewise polynomial collocation method for weaklysingular Volterra integral equations. Math. Comp., 68 (1999), 1079–1095.[4] G. A. Chandler and I. G. Graham, Product integration-collocation methods for noncompactintegral operator equations, Math. Comp., 50 (1988), 125–138. KASSEM MUSTAPHA
M α = 0 . α = 0 . α = 0 .
810 8.4612e-04 8.4612e-04 8.4612e-0420 1.9932e-04 2.0858 1.9932e-04 2.0858 1.9932e-04 2.085840 4.8140e-05 2.0498 4.8140e-05 2.0498 4.8140e-05 2.049880 1.1795e-05 2.0291 1.1449e-05 2.0720 1.1449e-05 2.0720160 3.0940e-06 1.9306 2.9063e-06 1.9780 2.7481e-06 2.0587
Table 4