A high-order L2 type difference scheme for the time-fractional diffusion equation
aa r X i v : . [ m a t h . NA ] F e b A high-order L2 type difference scheme for the time-fractional diffusionequation ⋆ Anatoly A. Alikhanov a, ∗ , Chengming Huang b a North-Caucasus Federal University, Pushkin str. 1, Stavropol, 355017, Russia b School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China
Abstract
The present paper is devoted to constructing L2 type difference analog of the Caputo fractional derivative.The fundamental features of this difference operator are studied and it is used to construct difference schemesgenerating approximations of the second and fourth order in space and the (3 − α ) th-order in time for thetime fractional diffusion equation with variable coefficients. Stability of the schemes under considerationas well as their convergence with the rate equal to the order of the approximation error are proven. Thereceived results are supported by the numerical computations performed for some test problems. Keywords: fractional diffusion equation, finite difference method, stability, convergence
1. Introduction
A significant growth of the researches’ attention to the fractional differential equations has been noticedlately . It is brought about by many effective applications of fractional calculation to various branches ofscience and engineering [1, 2, 3, 4, 5, 6]. For instance, we cannot dispense with mathematical language offractional derivatives when it comes to the description of the physical process of statistical transfer which,as it is well known, brings us to diffusion equations of fractional orders [7, 8].Let us consider the time fractional diffusion equation with variable coefficients ∂ α t u ( x, t ) = L u ( x, t ) + f ( x, t ) , < x < l, < t ≤ T, (1) u (0 , t ) = 0 , u ( l, t ) = 0 , ≤ t ≤ T, u ( x,
0) = u ( x ) , ≤ x ≤ l, (2)where ∂ α t u ( x, t ) = 1Γ(1 − α ) t Z ∂u ( x, η ) ∂η ( t − η ) − α dη, < α < α , L u ( x, t ) = ∂∂x (cid:18) k ( x, t ) ∂u∂x (cid:19) − q ( x, t ) u,k ( x, t ) ≥ c > q ( x, t ) ≥ f ( x, t ) are given functions.The time fractional diffusion equation constitutes a linear integro - differential equation. Its solutionin many cases cannot be found in an analytical form; as a consequence it is required to apply numerical ⋆ The reported study was jointly funded by RFBR (No. 20-51-53007) and NSFC (No. 12011530058) ∗ Corresponding author
Email addresses: [email protected] (Anatoly A. Alikhanov), [email protected] (Chengming Huang)
Preprint submitted to Journal February 18, 2021 ethods. Nevertheless, in contrast to the classical case, when we numerically approximate a time fractionaldiffusion equation on a certain time layer, we need information about all the previous time layers. That iswhy algorithms for solving the time fractional diffusion equations are rather labour-consuming even in one- dimensional case. When we pass to two - dimensional and three - dimensional problems, their complexitygrows significantly. In this respect constructing stable differential schemes of higher order approximation isa major task.A common difference approximation of fractional derivative (3) is the so-called L ∂ α t j +1 u ( x, t ) = 1Γ(1 − α ) j X s =0 u ( x, t s +1 ) − u ( x, t s ) t s +1 − t s t s +1 Z t s dη ( t j +1 − η ) α + r j +1 , (4)where 0 = t < t < . . . < t j +1 , and r j +1 is the local truncation error. In the case of the uniformgrid, τ = t s +1 − t s , for all s = 0 , , . . . , j + 1, it was proved that r j +1 = O ( τ − α ) [9, 10, 11]. The L L ∂ α t f ( t )of the function f ( t ) is to replace the integrand f ( t ) inside the integral by its piecewise linear interpolatingpolynomial (see [2, 9] ). A simple technique for improving the accuracy of L L − α to the order r + 1 − α , where r ≥ α ∈ (0 , O ( τ − α ), called L − − α ) th-order numericalformula (called the L α ∈ (0 , L − σ formula)to approximate the Caputo fractional derivative ∂ α t f ( t ) at a special points with the numerical accuracy oforder 3 − α was derived. Then some finite difference methods based on the L − σ formula were proposedfor solving the time-fractional diffusion equation. In [23, 24] L − σ formula was generalized and appliedfor solving the multi-term, distributed and variable order time-fractional diffusion equations.Difference schemes of the heightened order of approximation such as the compact difference scheme[13, 16, 17, 18, 15] and spectral method [10, 19, 20] were used to enhance the spatial accuracy of fractionaldiffusion equations.By means of the energy inequality method, a priori estimates for the solution of the Dirichlet, Robin andnon-local boundary value problems for the diffusion-wave equation with the Caputo fractional derivativehave been found in [14, 25, 26].In the present paper we construct L O ( τ − α ) for each α ∈ (0 , − α )th-order in time for the time fractional diffusion equation with variable coefficients are built. By means ofthe method of energy inequalities, the stability and convergence of these schemes are proven. Numericalcomputations of some test problems confirming reliability of the obtained results are implemented. Themethod can be without difficulty expanded to other time fractional partial differential equations with otherboundary conditions. 2 . The L2 type fractional numerical differentiation formula In this section we study a difference analog of the Caputo fractional derivative with the approximationorder O ( τ − α ) and explore its fundamental features.We consider the uniform grid ¯ ω τ = { t j = jτ, j = 0 , , . . . , M ; T = τ M } . For the Caputo fractionalderivative of the order α , 0 < α <
1, of the function u ( t ) ∈ C [0 , T ] at the fixed point t j +1 , j ∈ { , , . . . , M − } the following equalities are valid ∂ α t j +1 u ( t ) = 1Γ(1 − α ) t j +1 Z u ′ ( η ) dη ( t j +1 − η ) α = 1Γ(1 − α ) t Z u ′ ( η ) dη ( t j +1 − η ) α + 1Γ(1 − α ) j X s =2 t s +1 Z t s u ′ ( η ) dη ( t j +1 − η ) α . (5)On each interval [ t s − , t s ] (1 ≤ s ≤ j ), applying the quadratic interpolation Π ,s u ( t ) of u ( t ) that usesthree points ( t s − , u ( t s − )), ( t s , u ( t s )) and ( t s +1 , u ( t s +1 )), we arrive atΠ ,s u ( t ) = u ( t s − ) ( t − t s )( t − t s +1 )2 τ − u ( t s ) ( t − t s − )( t − t s +1 ) τ + u ( t s +1 ) ( t − t s − )( t − t s )2 τ , (Π ,s u ( t )) ′ = u t,s + u ¯ tt,s ( t − t s +1 / ) , (6)and u ( t ) − Π ,s u ( t ) = u ′′′ ( ¯ ξ s )6 ( t − t s − )( t − t s )( t − t s +1 ) , (7)where t ∈ [ t s − , t s +1 ], ¯ ξ s ∈ ( t s − , t s +1 ), t s − / = t s − . τ , u t,s = ( u ( t s +1 ) − u ( t s )) /τ , u ¯ t,s = ( u ( t s ) − u ( t s − )) /τ .In (5), we make use of Π ,s u ( t ) in order to approximate u ( t ) on the interval [ t s − , t s ] (1 ≤ s ≤ j ). Inview of the equality t s Z t s − ( η − t s − / )( t j +1 − η ) − α dη = τ − α − α b ( α ) j − s +1 , ≤ s ≤ j (8)with b ( α ) l = 12 − α (cid:2) ( l + 1) − α − l − α (cid:3) − (cid:2) ( l + 1) − α + l − α (cid:3) , l ≥ , from (5) and (6) we get the difference analog of the Caputo fractional derivative of order α (0 < α <
1) forthe function u ( t ), at the points t j +1 ( j = 1 , , . . . ), in this form: ∂ α t j +1 u ( η ) = 1Γ(1 − α ) t Z u ′ ( η ) dη ( t j +1 − η ) α + 1Γ(1 − α ) j X s =2 t s +1 Z t s u ′ ( η ) dη ( t j +1 − η ) α ≈ − α ) t Z (Π , u ( η )) ′ dη ( t j +1 − η ) α + 1Γ(1 − α ) j X s =2 t s +1 Z t s (Π ,s u ( η )) ′ dη ( t j +1 − η ) α − α ) t Z u t, + u ¯ tt, ( η − t / )( t j +1 − η ) α dη + 1Γ(1 − α ) j X s =2 t s +1 Z t s u t,s + u ¯ tt,s ( η − t s +1 / )( t j +1 − η ) α dη = τ − α Γ(2 − α ) (cid:16) ( a ( α ) j − b ( α ) j − b ( α ) j − ) u t, + ( a ( α ) j − + b ( α ) j − + b ( α ) j ) u t, (cid:17) + τ − α Γ(2 − α ) j X s =2 ( − b ( α ) j − s u t,s − + ( a ( α ) j − s + b ( α ) j − s ) u t,s )= τ − α Γ(2 − α ) (cid:16) ( a ( α ) j − b ( α ) j − b ( α ) j − ) u t, + ( a ( α ) j − + b ( α ) j − + b ( α ) j − b ( α ) j − ) u t, (cid:17) + τ − α Γ(2 − α ) j − X s =2 ( − b ( α ) j − s − + a ( α ) j − s + b ( α ) j − s ) u t,s + ( a ( α )0 + b ( α )0 ) u t,j ! == τ − α Γ(2 − α ) j X s =0 c ( α ) j − s u t,s = ∆ α t j +1 u, (9)where a ( α ) l = ( l + 1) − α − l − α , l ≥ j = 1 c ( α ) s = ( a ( α )0 + b ( α )0 + b ( α )1 , s = 0 ,a ( α )1 − b ( α )1 − b ( α )0 , s = 1 , (10)for j = 2 c ( α ) s = a ( α )0 + b ( α )0 , s = 0 ,a ( α )1 + b ( α )1 + b ( α )2 − b ( α )0 , s = 1 ,a ( α )2 − b ( α )2 − b ( α )1 , s = 2 , (11)and for j ≥ c ( α ) s = a ( α )0 + b ( α )0 , s = 0 ,a ( α ) s + b ( α ) s − b ( α ) s − , ≤ s ≤ j − ,a ( α ) j − + b ( α ) j − + b ( α ) j − b ( α ) j − , s = j − ,a ( α ) j − b ( α ) j − b ( α ) j − , s = j. (12)We name the fractional numerical differentiation formula (9) for the Caputo fractional derivative of order α (0 < α <
1) the L2 formula.
Lemma 2.1.
For any α ∈ (0 , , j = 1 , , . . . , M − and u ( t ) ∈ C [0 , t j +1 ] | ∂ α t j +1 u − ∆ α t j +1 u | = O ( τ − α ) . (13) Proof.
Let ∂ α t j +1 u − ∆ α t j +1 u = R + R j +12 , where R = 1Γ(1 − α ) t Z u ′ ( η ) dη ( t j +1 − η ) α − − α ) t Z (Π , u ( η )) ′ dη ( t j +1 − η ) α − α ) t Z ( u ( η ) − Π , u ( η )) ′ dη ( t j +1 − η ) α = − α Γ(1 − α ) t Z ( u ( η ) − Π , u ( η )) dη ( t j +1 − η ) α +1 = − α − α ) t Z u ′′′ ( ¯ ξ ) η ( η − t )( η − t )( t j +1 − η ) − α − dη,R j +12 = 1Γ(1 − α ) j X s =2 t s +1 Z t s u ′ ( η ) dη ( t j +1 − η ) α − − α ) j X s =2 t s +1 Z t s (Π ,s u ( η )) ′ dη ( t j +1 − η ) α = 1Γ(1 − α ) j X s =2 t s +1 Z t s ( u ( η ) − Π ,s u ( η )) ′ ( t j +1 − η ) − α dη = − α Γ(1 − α ) j X s =2 t s +1 Z t s ( u ( η ) − Π ,s u ( η )) ( t j +1 − η ) − α − dη = − α − α ) j X s =2 t s +1 Z t s u ′′′ ( ¯ ξ s )( η − t s − )( η − t s )( η − t s +1 )( t j +1 − η ) − α − dη, Next we estimate the errors R and R j +12 :For j = 1 we have (cid:12)(cid:12) R (cid:12)(cid:12) = α − α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t Z u ′′′ ( ¯ ξ ) η ( η − t )( η − t )( t − η ) − α − dη (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ αM τ − α ) t Z ( t − η ) − α dη = 2 − α αM − α ) τ − α , For j ≥ (cid:12)(cid:12) R (cid:12)(cid:12) = α − α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t Z u ′′′ ( ¯ ξ ) η ( η − t )( η − t )( t j +1 − η ) − α − dη (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ √ αM τ − α ) t Z ( t j +1 − η ) − α − dη = √ M τ − α − α ) (cid:0) ( j − − α − ( j + 1) − α (cid:1) ≤ √ − − α ) M − α ) τ − α , (cid:12)(cid:12)(cid:12) R j +12 (cid:12)(cid:12)(cid:12) = α − α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) j X s =2 t s +1 Z t s u ′′′ ( ¯ ξ s )( η − t s − )( η − t s )( η − t s +1 )( t j +1 − η ) − α − dη (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ √ αM τ − α ) j − X s =2 t s +1 Z t s ( t j +1 − η ) − α − dη + αM τ − α ) t j +1 Z t j ( t j +1 − η ) − α dη = √ αM τ − α ) t j Z t ( t j +1 − η ) − α − dη + αM τ − α ) t j +1 Z t j ( t j +1 − η ) − α dη √ M τ − α ) (cid:0) τ − α − t − αj − (cid:1) + αM τ − α ) τ − α − α ≤ √
39 + α (1 − α ) ! M − α ) τ − α . Lemma 2.2.
For all α ∈ (0 , and s = 1 , , , . . . − α ( s + 1) α < a s < − αs α , (14) α (1 − α )( s + 2) α +1 < a s − a s +1 < α (1 − α ) s α +1 , (15) α (1 − α )12( s + 1) α +1 < b s < α (1 − α )12 s α +1 , (16) Proof.
The validity of Lemma 2.2 follows from the following equalities: a s = (1 − α ) Z dξ ( s + ξ ) α ,a s − a s +1 = α (1 − α ) Z dη Z dξ ( s + ξ + η ) α +1 ,b s = α (1 − α )2 − α Z ηdη s +1+ η Z s +1 − η dξξ α +1 . For j = 1 we have c ( α )0 = 2 + α α (2 − α ) , c ( α )1 = 2 − α α (2 − α ) , c ( α )0 + 3 c ( α )1 = 2 − α (1 − α )2 − α > . For j ≥
2, the next lemma shows properties of the coefficient c ( α ) s defined in (11) and (12) Lemma 2.3.
For any α ∈ (0 , and c ( α ) s ( ≤ s ≤ j , j ≥ ) the following inequalities are valid · − α ( j + 1) α < c ( α ) j < − αj α , (17) c ( α )0 > c ( α )2 > c ( α )3 > . . . > c ( α ) j − > c ( α ) j − > c ( α ) j , (18) c ( α )0 + 3 c ( α )1 − c ( α )2 > . (19) Proof.
For j ≥ c ( α ) j = a ( α ) j − b ( α ) j − b ( α ) j − < a ( α ) j < − αj α .c ( α ) j = a ( α ) j − b ( α ) j − b ( α ) j − > − α ( j + 1) α − α (1 − α )12 j α +1 − α (1 − α )12( j − α +1 . − α ( j + 1) α (cid:18) − α · (cid:18) j + 1 j (cid:19) α · j − α · (cid:18) j + 1 j − (cid:19) α · j − (cid:19) > − α ( j + 1) α (cid:18) − · · − · (cid:19) = 1116 · − α ( j + 1) α Inequality (17) is proved. Let us prove inequality (18). c ( α )0 − c ( α )2 ≥ a ( α )0 − a ( α )2 + b ( α )0 − b ( α )2 + b ( α )1 − b ( α )3 > , For j ≥
5, 2 ≤ s ≤ j − c ( α ) s − c ( α ) s +1 = a ( α ) s − a ( α ) s +1 − b ( α ) s − + 2 b ( α ) s − b ( α ) s +1 = 12 − α (cid:0) − ( s + 2) − α + 3( s + 1) − α − s − α + ( s − − α (cid:1) − (cid:0) − ( s + 2) − α + 3( s + 1) − α − s − α + ( s − − α (cid:1) = α (1 − α ) Z dz Z dz Z dz ( s − z + z + z ) α +1 − α (1 − α )(1 + α )2 Z dz Z dz Z dz ( s − z + z + z ) α +2 = α (1 − α ) Z dz Z dz Z (cid:16) − α · s − z + z + z (cid:17) ( s − z + z + z ) α +1 dz > α (1 − α )( s + 2) α +1 − α Z dz Z dz Z dz z + z + z . Since Z dz Z dz Z dz z + z + z = 12 (44 ln 2 −
27 ln 3) < ,c ( α ) s − c ( α ) s +1 > α (1 − α )( s + 2) α +1 (cid:18) − α (cid:19) > α (1 − α )2( s + 2) α +1 > . For j ≥ c ( α ) j − − c ( α ) j − = a ( α ) j − − a ( α ) j − − b ( α ) j − + 2 b ( α ) j − − b ( α ) j − − b ( α ) j > α (1 − α )2 j α +1 − b ( α ) j > α (1 − α )2 j α +1 − α (1 − α )12 j α +1 = 5 α (1 − α )12 j α +1 > . For j ≥ c ( α ) j − − c ( α ) j = a ( α ) j − − a ( α ) j − b ( α ) j − + 2 b ( α ) j − + 2 b ( α ) j > a ( α ) j − − a ( α ) j − b ( α ) j − + 2 b ( α ) j − − b ( α ) j > α (1 − α )2( j + 1) α +1 > . Inequality (18) is proved.For j = 2 we get c ( α )0 + 3 c ( α )1 − c ( α )2 = a ( α )0 + 3 a ( α )1 − a ( α )2 − b ( α )0 + 7 b ( α )1 + 7 b ( α )2 a ( α )0 − a ( α )1 − b ( α )0 = 3 − − α − − α . Since, for any function f ( x ) ∈ C [0 , f (0) = 0, f (1) = 0 and f ′′ ( x ) < x ∈ (0 ,
1) then f ( x ) > x ∈ (0 , c ( α )0 + 3 c ( α )1 − c ( α )2 > f ( α ) = 3 − − α − − α > α ∈ (0 , . For j = 3 we get c ( α )0 + 3 c ( α )1 − c ( α )2 = a ( α )0 + 3 a ( α )1 − a ( α )2 − b ( α )0 + 7 b ( α )1 − b ( α )2 − b ( α )3 > a ( α )0 − a ( α )1 − b ( α )0 + 4( a ( α )1 − a ( α )2 − b ( α )3 ) > − − α − − α > . Lemma 2.4.
For any real constants c , c such that c ≥ max { c , − c } , and { v j } j = Mj =0 the following in-equality holds v j +1 ( c v j +1 − ( c − c ) v j − c v j − ) ≥ E j +1 − E j , j = 1 , . . . , M − , (20) where E j = r c − c r c + 3 c ! v j + r c − c v j − r c − c r c + 3 c ! v j − ! , j = 1 , , . . . , M. Proof.
The proof of Lemma 2.4 immediately follows from the next equality v j +1 ( c v j +1 − ( c − c ) v j − c v j − ) − E j +1 + E j = r c − c − r c + 3 c ! v j +1 − r c − c v j + r c − c r c + 3 c ! v j − ! ≥ . Lemma 2.5. [22] If g j +1 j ≥ g j +1 j − ≥ . . . ≥ g j +10 > , j = 0 , , . . . , M − then for any function v ( t ) definedon the grid ω τ one has the inequalities v j +1 g ∆ α t j +1 v ≥ g ∆ α t j +1 ( v ) + 12 g j +1 j (cid:16) g ∆ α t j +1 v (cid:17) , (21) where g ∆ α t j +1 y i = j X s =0 (cid:0) y s +1 i − y si (cid:1) g j +1 s , is a difference analog of the Caputo fractional derivative of the order α ( < α < ). Lemma 2.6.
For any function v ( t ) defined on the grid ω τ one has the inequality v j +1 ∆ α t j +1 v ≥ τ − α Γ(2 − α ) ( E j +1 − E j ) + 12 ¯∆ α t j +1 v = τ − α Γ(2 − α ) ( E j +1 − E j ) − τ − α − α ) ¯ c ( α ) j v , (22) where ¯∆ α t j +1 v = τ − α Γ(2 − α ) j X s =0 ¯ c ( α ) j − s ( v s +1 − v s ) , j = 1 , , . . . , M, ¯ c ( α )0 = c ( α )2 , ¯ c ( α )1 = c ( α )2 , ¯ c ( α ) s = c ( α ) s , s = 2 , , . . . , j, for j = 1 , , , . . . , M, E j = E j ( c ( α )0 − c ( α )2 , c ( α )1 − c ( α )2 ) , E j = E j + 12 j − X s =0 ¯ c ( α ) j − − s v s +1 . roof. For j = 1 we have v ∆ α t v = τ − α Γ(2 − α ) v (cid:16) c ( α )0 ( v − v ) + c ( α )1 ( v − v ) (cid:17) = τ − α Γ(2 − α ) v (cid:16) ( c ( α )0 − c ( α )2 ) v − ( c ( α )0 − c ( α )1 ) v − ( c ( α )1 − c ( α )2 ) v (cid:17) + τ − α Γ(2 − α ) c ( α )2 (cid:0) v − v v (cid:1) ≥ τ − α Γ(2 − α ) ( E − E ) + τ − α − α ) c ( α )2 (cid:0) v − v (cid:1) = τ − α Γ(2 − α ) ( E − E ) + 12 ¯∆ α t v . For j = 2 , , . . . , M − v j +1 ∆ α t j +1 v = τ − α Γ(2 − α ) v j +1 j X s =0 c ( α ) j − s ( v s +1 − v s ) τ − α Γ(2 − α ) v j +1 (cid:16) ( c ( α )0 − c ( α )2 )( v j +1 − v j ) + ( c ( α )1 − c ( α )2 )( v j − v j − ) (cid:17) + v j +1 ¯∆ α t j +1 v = τ − α Γ(2 − α ) v j +1 (cid:16) ( c ( α )0 − c ( α )2 ) v j +1 − ( c ( α )0 − c ( α )1 ) v j − ( c ( α )1 − c ( α )2 ) v j − (cid:17) + v j +1 ¯∆ α t j +1 v ≥ τ − α Γ(2 − α ) ( E j +1 − E j ) + 12 ¯∆ α t j +1 v . In addition, the following equality holds¯∆ α t j +1 v = τ − α Γ(2 − α ) j X s =0 ¯ c ( α ) j − s ( v s +1 − v s ) = τ − α Γ(2 − α ) j X s =0 ¯ c ( α ) j − s v s +1 − j − X s =0 ¯ c ( α ) j − − s v s +1 − ¯ c ( α ) j v ! .
3. A difference scheme for the time fractional diffusion equation
In this section for problem (1)–(2) a difference scheme with the approximation order O ( h + τ − α ) isconstructed. The stability of the constructed difference scheme as well as its convergence in the grid L - norm with the rate equal to the order of the approximation error is proved. The obtained results aresupported with numerical calculations carried out for a test example. Lemma 3.1. [22] For any functions k ( x ) ∈ C x and v ( x ) ∈ C x the following equality holds true: ddx (cid:18) k ( x ) ddx v ( x ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) x = x i = k ( x i +1 / ) v ( x i +1 ) − ( k ( x i +1 / ) + k ( x i − / )) v ( x i ) + k ( x i − / ) v ( x i − ) h + O ( h ) . (23)9et u ( x, t ) ∈ C , x,t be a solution of problem (1)–(2). Then we consider equation (1) for ( x, t ) = ( x i , t j +1 ) ∈ Q T , i = 1 , , . . . , N − j = 1 , , . . . , M − ∂ α t j +1 u = L u ( x, t ) | ( x i ,t j +1 ) + f ( x i , t j +1 ) . (24)On the basis of Lemmas 2.1 and 3.1 we have ∂ α t j +1 u = ∆ α t j + σ u + O ( τ − α ) L u ( x, t ) | ( x i ,t j +1 ) = Λ u ( x i , t j +1 ) + O ( h ) , where the difference operator Λ is defined as follows(Λ y ) i = (( ay ¯ x ) x − dy ) i = a i +1 y i +1 − ( a i +1 + a i ) y i + a i y i − h − d i y i ,y ¯ x,i = y i − y i − h , y x,i = y i +1 − y i h , with the coefficients a j +1 i = k ( x i − / , t j +1 ), d j +1 i = q ( x i , t j +1 ). Let ϕ j +1 i = f ( x i , t j +1 ), then we get thedifference scheme with the approximation order O ( h + τ − α ):∆ α t j +1 y i = Λ y j +1 i + ϕ j +1 i , i = 1 , , . . . , N − , j = 1 , , . . . , M − , (25) y (0 , t ) = 0 , y ( l, t ) = 0 , t ∈ ω τ , y ( x,
0) = u ( x ) , x ∈ ω h , (26) Remark.
We assume that the solution y i is found with the order of accuracy O ( h + τ − α ). For example,we can use L , τ ] with step τ = O ( τ − α − α ). Theorem 3.1.
The difference scheme (25)–(26) is unconditionally stable and its solution meets the followinga priori estimates: M − X j =1 (cid:16) k y j +1 k + k y j +1¯ x ] | (cid:17) τ ≤ M k y k + k y k + M − X j =1 k ϕ j +1 k τ , (27) where k y ] | = N P i =1 y i h , M > is a known number independent of h and τ .Proof. Taking the inner product of the equation (25) with y j +1 , we have (cid:16) y j +1 , ∆ α t j +1 y (cid:17) − (cid:0) y j +1 , Λ y j +1 (cid:1) = (cid:0) y j +1 , ϕ j +1 (cid:1) . (28)Using Lemma 2.6, we obtain (cid:16) y j +1 , ∆ α t j +1 y (cid:17) ≥ τ − α Γ(2 − α ) ( E j +1 − E j ) + 12 ¯∆ α t j +1 k y k = τ − α Γ(2 − α ) ( E j +1 − E j ) − τ − α − α ) ¯ c ( α ) j k y k , , j = 1 , , . . . , M − , where E j = s c ( α )0 − c ( α )1 s c ( α )0 + 3 c ( α )1 − c ( α )2 k y j k (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)s c ( α )0 − c ( α )1 y j − s c ( α )0 − c ( α )1 s c ( α )0 + 3 c ( α )1 − c ( α )2 y j − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) . E j = E j + 12 j − X s =0 ¯ c ( α ) j − − s k y s +1 k . For the difference operator Λ using Green’s first difference formula for the functions vanishing at x = 0and x = l , we get ( − Λ y, y ) ≥ c k y ¯ x ] | .From (28), using that (cid:0) y j +1 , ϕ j +1 (cid:1) ≤ c l k y j +1 k + l c k ϕ j +1 k ≤ c k y j +1¯ x ] | + l c k ϕ j +1 k , one obtains the inequality τ − α Γ(2 − α ) ( E j +1 − E j ) + c k y j +1¯ x ] | ≤ l c k ϕ j +1 k + τ − α − α ) ¯ c ( α ) j k y k . (29)Multiplying inequality (29) by τ and summing the resulting relation over j from 1 to M − Numerical computations are executed for a test problem on the assumption that the function u ( x, t ) = sin( πx ) (cid:0) t α + t + 1 (cid:1) is the exact solution of problem (1)–(2) with the coefficients k ( x, t ) = 2 − cos( xt ), q ( x, t ) = 1 − sin( xt ) and l = 1, T = 1.The errors ( z = y − u ) and convergence order (CO) in the norms k · k and k · k C (¯ ω hτ ) , where k y k C (¯ ω hτ ) =max ( x i ,t j ) ∈ ¯ ω hτ | y | , are given in Table 1. Table 1 demonstrates that as the number of the spatial subintervals and time steps increases, while h = τ − α , then the maximum error decreases, as it is expected and the convergence order of the approximatescheme is O ( h ) = O ( τ − α ), where the convergence order is given by the formula: CO= log h h k z kk z k ( z i isthe error corresponding to h i ). Table 2 shows that if h = 1 / O ( τ − α ), wherethe convergence order is given by the following formula: CO= log τ τ k z kk z k .
4. A compact difference scheme for the time fractional diffusion equation
In this section for problem (1)–(2), we create a compact difference scheme with the approximation order O ( h + τ − α ) in the case when k = k ( t ) and q = q ( t ). The stability and convergence of the constructeddifference scheme in the grid L - norm with the rate equal to the order of the approximation error areproved. The found results are supported by the numerical calculations implemented for a test example.11 able 1: The error and the convergence order in the norms k · k and k · k C (¯ ω hτ ) when decreasing time-grid size for differentvalues of α = 0 .
1; 0 .
5; 0 . τ − α = h . α τ h max ≤ j ≤ M k z j k CO max ≤ j ≤ M k z j k C (¯ ω hτ ) CO max ≤ j ≤ M k z j ¯ x ] | CO Table 2: The error and the convergence order in the norms k · k and k · k C (¯ ω hτ ) when decreasing time-grid size for differentvalues of α = 0 .
3; 0 .
5; 0 . h = 1 / . α τ max ≤ j ≤ M k z j k CO max ≤ j ≤ M k z j k C (¯ ω hτ ) CO max ≤ j ≤ M k z j ¯ x ] | CO Let a difference scheme be put into a correspondence with differential problem (1)–(2) in the case when k = k ( t ) and q = q ( t ):∆ α t j +1 H h y i = a j +1 y j +1¯ xx,i − d j +1 H h y j +1 i + H h ϕ j +1 i , i = 1 , . . . , N − , j = 0 , , . . . , M − , (30) y (0 , t ) = 0 , y ( l, t ) = 0 , t ∈ ω τ , y ( x,
0) = u ( x ) , x ∈ ω h , (31)where H h v i = v i + h v ¯ xx,i / i = 1 , . . . , N − a j +1 = k ( t j +1 ), d j +1 = q ( t j +1 ), ϕ j +1 i = f ( x i , t j +1 ).From Lemma 2.1 it follows that if u ∈ C , x,t , then the difference scheme has the approximation order O ( τ − α + h ). 12 .2. Stability and convergence Theorem 4.1.
The difference scheme (30)–(31) is unconditionally stable and its solution meets the followinga priori estimate: M − X j =1 (cid:16) kH h y j +1 k + k y j +1¯ x ] | (cid:17) τ ≤ M kH h y k + kH h y k + M − X j =1 kH h ϕ j +1 k τ , (32) where M > is a known number independent of h and τ .Proof. Taking the inner product of the equation (30) with H h y j +1 = ( H h y ) j +1 , we have( H h y j +1 , ∆ α t j +1 H h y ) − a j +1 ( H h y j +1 , y j +1¯ xx )+ d j +1 ( H h y j +1 , H h y j +1 ) = ( H h y j +1 , H h ϕ j +1 ) . (33)We transform the terms in identity (33) as( H h y j +1 , ∆ α t j +1 H h y ) ≥ τ − α Γ(2 − α ) ( E j +1 − E j ) + 12 ¯∆ α t j + σ kH h y k == τ − α Γ(2 − α ) ( E j +1 − E j ) − τ − α − α ) ¯ c ( α ) j kH h y k , , j = 1 , , . . . , M − , where E j = s c ( α )0 − c ( α )1 s c ( α )0 + 3 c ( α )1 − c ( α )2 kH h y j k + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)s c ( α )0 − c ( α )1 H h y j − s c ( α )0 − c ( α )1 s c ( α )0 + 3 c ( α )1 − c ( α )2 H h y j − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , E j = E j + 12 j − X s =0 ¯ c ( α ) j − − s kH h y s +1 k . − ( H h y j +1 , y j +1¯ xx ) = − ( y j +1 , y j +1¯ xx ) − h k y j +1¯ xx k = k y j +1¯ x ] | − N − X i =1 ( y j +1¯ x,i +1 − y j +1¯ x,i ) h ≥ k y j +1¯ x ] | − k y j +1¯ x ] | = 23 k y j +1¯ x ] | , ( H h y j +1 , H h ϕ j +1 ) ≤ ε kH h y j +1 k + 14 ε kH h ϕ j +1 k = ε N − X i =1 y j +1 i − + 10 y j +1 i + y j +1 i +1 ! h + 14 ε kH h ϕ j +1 k ≤ ε k y j +1 k + 14 ε kH h ϕ j +1 k ≤ εl k y j +1¯ x ] | + 14 ε kH h ϕ j +1 k . In view of the above-performed transformations, from identity (33) at ε = c l we get the inequality τ − α Γ(2 − α ) ( E j +1 − E j ) + c k y j +1¯ x ] | ≤ l c kH h ϕ j +1 k + τ − α − α ) ¯ c ( α ) j kH h y k . The following process is similar to the proof of Theorem 3.1, and we leave it out.13he norm kH h y k is equivalent to the norm k y k , which follows from the inequalities512 k y k ≤ kH h y k ≤ k y k . Using a priori estimate (32), we obtain the convergence result.
Theorem 4.2.
Let u ( x, t ) ∈ C , x,t be the solution of problem (1)–(2) in the case k = k ( t ) , q = q ( t ) , and let { y ji | ≤ i ≤ N, ≤ j ≤ M } be the solution of difference scheme (30)–(31). Then it holds true that vuut M − X j =1 (cid:16) k z j +1 k + k z j +1¯ x ] | (cid:17) τ ≤ C R (cid:0) τ − α + h (cid:1) , ≤ j ≤ M, where z ji = u ( x i , t j ) − y ji and C R is a positive constant independent of τ and h .4.3. Numerical results Numerical calculations are performed for a test problem when the function u ( x, t ) = sin( πx ) (cid:0) t α + t + 1 (cid:1) is the exact solution of the problem (1)–(2) with the coefficients k ( x, t ) = 2 − cos( t ), q ( x, t ) = 1 − sin( t ) and l = 1, T = 1.The errors ( z = y − u ) and convergence order (CO) in the norms k · k and k · k C (¯ ω hτ ) , where k y k C (¯ ω hτ ) =max ( x i ,t j ) ∈ ¯ ω hτ | y | , are given in Table 1. Table 3 shows that as the number of the spatial subintervals and time steps increases keeping h = τ − α ,the maximum error decreases, as it is expected and the convergence order of the compact difference schemeis O ( h ) = O ( τ − α ), where the convergence order is given by the formula: CO= log h h k z kk z k ( z i is the errorcorresponding to h i ). Table 4 demonstrates that if h = 1 / O ( τ − α ),where the convergence order is given by the following formula: CO= log τ τ k z kk z k .
5. Conclusion
In the current paper we construct a L O ( τ − α ). The fundamental features of this difference operator are studied.New difference schemes of the second and fourth approximation order in space and the 3 − α approximationorder in time for the time fractional diffusion equation with variable coefficients are also constructed. Thestability and convergence of these schemes with the rate equal to the order of the approximation error areproved. The method can be without difficulty expanded to include other time fractional partial differentialequations with other boundary conditions.Numerical tests entirely corroborating the found theoretical results are implemented. In all the calcula-tions Julia v1.5.1 is used. References [1] A. M. Nakhushev, Fractional Calculus and its Application, FIZMATLIT, Moscow, 2003 (in Russian).[2] K. B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.[3] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999.[4] R. Hilfer (Ed.), Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000.[5] A. A. Kilbas, H. M. Srivastava, J. J. Trujillo, Theory and Applications of Fractional Differential Equation, Elsevier,Amsterdam, 2006. able 3: The error and the convergence order in the norms k · k and k · k C (¯ ω hτ ) when decreasing time-grid size for differentvalues of α = 0 .
1; 0 .
5; 0 . τ − α = ( h/ . α τ h max ≤ j ≤ M k z j k CO max ≤ j ≤ M k z j k C (¯ ω hτ ) CO max ≤ j ≤ M k z j ¯ x ] | CO Table 4: The error and the convergence order in the norms k · k and k · k C (¯ ω hτ ) when decreasing time-grid size for differentvalues of α = 0 .
3; 0 .
5; 0 . h = 1 / . α τ max ≤ j ≤ M k z j k CO max ≤ j ≤ M k z j k C (¯ ω hτ ) CO max ≤ j ≤ M k z j ¯ x ] | CO [6] V. V. Uchaikin, Method of Fractional Derivatives, Artishok, Ul’janovsk, 2008 (in Russian).[7] R. R. Nigmatullin, Realization of the generalized transfer equation in a medium with fractal geometry, Physica Status (B):Basic Res.[8] K. V. Chukbar, Stochastic transport and fractional derivatives, Zh. Eksp. Teor. Fiz. 108 (1995), 1875- 1884[9] Z. Z. Sun, X. N. Wu, A fuly discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209.[10] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225(2007) 1553–1552.[11] A.A. Alikhanov, Numerical methods of solutions of boundary value problems for the multi-term variable-distributed orderdiffusion equation, Appl. Math. Comput. 268 (2015) 12–22.[12] M. Kh. Shkhanukov-Lafishev, F.I. Taukenova, Difference methods for solving boundary value problems for fractional ifferential equations, Comput. Math. Math. Phys. 46(10) (2006) 1785–1795.[13] C. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdif-fusion equations, SIAM J. Scien. Comput. 32(4) (2010) 1740–1760.[14] A.A. Alikhanov, Boundary value problems for the diffusion equation of the variable order in differential and differencesettings, Appl. Math. Comput. 219 (2012) 3938–3946.[15] Yuan-Ming Wang, Lei Ren, A high-order L A. A. Samarskii, V. B. Andreev,