Fast second-order evaluation for variable-order Caputo fractional derivative with applications to fractional sub-diffusion equations
aa r X i v : . [ m a t h . NA ] F e b Noname manuscript No. (will be inserted by the editor)
Fast second-order evaluation for variable-order Caputo fractionalderivative with applications to fractional sub-diffusion equations
Jia-li Zhang · Zhi-wei Fang · Hai-wei Sun Received: date / Accepted: date
Abstract
In this paper, we propose a fast second-order approximation to the variable-order (VO)Caputo fractional derivative, which is developed based on L σ formula and the exponential-sum-approximation technique. The fast evaluation method can achieve the second-order accuracyand further reduce the computational cost and the acting memory for the VO Caputo fractionalderivative. This fast algorithm is applied to construct a relevant fast temporal second-order andspatial fourth-order scheme ( F L σ scheme) for the multi-dimensional VO time-fractional sub-diffusion equations. Theoretically, F L σ scheme is proved to fulfill the similar properties of thecoefficients as those of the well-studied L σ scheme. Therefore, F L σ scheme is strictly provedto be unconditionally stable and convergent. A sharp decrease in the computational cost and theacting memory is shown in the numerical examples to demonstrate the efficiency of the proposedmethod. Keywords variable-order Caputo fractional derivative · exponential-sum-approximation method · fast algorithm · time-fractional sub-diffusion equation · stability and convergence Mathematics Subject Classification (2010)
Over the last few decades, the fractional calculus has gained much attention of both mathematicsand physical science due to the non-locality of fractional operators. In fact the fractional calculushas been widely applied in various fields including the biology, the ecology, the diffusion, and the
This work is supported in part by research grants of the Science and Technology Development Fund, Macau SAR (fileno. 0118/2018/A3), and University of Macau (file no. MYRG2018-00015-FST).1Department of Mathematics, University of Macau, MacaoE-mail: [email protected] of Mathematics, University of Macau, MacaoE-mail: [email protected] Author. Department of Mathematics, University of Macau, MacaoE-mail: [email protected] Jia-li Zhang et al. control system [3, 14, 18, 20, 25, 27–29]. Recently, more and more researchers revealed that manyimportant dynamical problems exhibit the fractional order behavior that may vary with time, space,or some other factors, which leads to the concept of the variable-order (VO) fractional operator;see [21, 36, 38]. This fact indicates that the VO fractional calculus is an expected tool to providean effective mathematical framework to characterize the complex problems in fields of science andengineering; for instances, anomalous diffusion [7, 8, 37, 43], viscoelastic mechanics [9, 10, 15, 33], andpetroleum engineering [26].The VO fractional derivative can be regarded as an extension of the constant-order (CO)fractional derivative [39]. An interesting extension of the classical fractional calculus was proposedby Samko and Ross [31] in which they generalized the Riemann-Liouville and Marchaud fractionaloperators in the VO case. Later, Lorenzo and Hartley [21] first introduced the concept of VOoperators. According to the concept, the order of the operator is allowed to vary as a functionof independent variables such as time and space. Afterwards, various VO differential operators withspecific meanings were defined. Coimbra [9] gave a novel definition for the VO differential operator bytaking the Laplace-transform of the Caputo’s definition of the fractional derivative. Soon, Coimbra,and Kobayashi [34] showed that Coimbra’s definition was better suited for modeling physicalproblems due to its meaningful physical interpretations; see also [30]. Moreover, the Coimbra’sVO fractional derivative could be considered as the Caputo-type definition, which is defined asfollows [9, 36] C D α ( t ) t u ( t ) ≡ Γ (1 − α ( t )) Z t u ′ ( τ )( t − τ ) α ( t ) d τ, (1)where α ( t ) ∈ [ α, α ] ⊂ (0 ,
1) is the VO function and Γ ( · ) denotes the Gamma function.Since the problems described by the VO fractional operator are difficult to handle analytically,possible numerical implementations of the VO fractional problems are given. In [42], Zhao, Sun, andKarniadakis derived two second-order approximations for the VO Caputo fractional derivative andprovided the error analysis. For the VO time-fractional sub-diffusion equations, Du, Alikhanov,and Sun [11] proposed L σ scheme, which makes use of the piecewise high-order polynomialinterpolation of the solution. The resulting method was investigated to be unconditionally stableand second-order convergent via the energy method.In consequence of the nonlocal property and historical dependence of the fractional operators,the aforementioned numerical methods always require all previous function values, which leadsto an average O ( n ) storage and computational cost O ( n ), where n is the total number of thetime levels. To overcome this difficulty, many efforts have been made to speed up the evaluationof the CO fractional derivative [2, 4, 13, 16, 17, 22–24, 40]. Nevertheless, the coefficient matrices ofthe numerical schemes for the VO fractional problems lose the Toeplitz-like structure and the VOfractional derivative is no longer a convolution operator. Therefore, those fast methods for the COfractional derivative cannot be directly applied to VO cases. Recently, Fang, Sun, and Wang [12]proposed a fast algorithm for the VO Caputo fractional derivative based on a shifted binary blockpartition and uniform polynomial approximations of degree r . Compared with the general directmethod, the proposed algorithm reduces the memory requirement from O ( n ) to O ( r log n ) and thecomplexity from O ( n ) to O ( rn log n ). Zhang, Fang, and Sun [41] developed an exponential-sum-approximation (ESA) technique to speed up L O (log n ) and the computational cost to O ( n log n ), respectively. Both of those fast methods were proved to achieve the convergence orderof 2 − α with expected parameters.Motivated by the need of fast high-order approaches for the VO Caputo fractional derivative (1),we combine L σ formula [11] with the ESA technique [41] to produce a fast evaluation formula, itle Suppressed Due to Excessive Length 3 which is called F L σ formula later. Compared with L σ formula, the fast algorithm reduces theacting memory from O ( n ) to O (log n ) and flops from O ( n ) to O ( n log n ). Then F L σ formula isapplied to construct a fast temporal second-order and spatial fourth-order difference scheme ( F L σ scheme) for the multi-dimensional VO time-fractional sub-diffusion equations, which significantlyreduces the memory requirement and computational complexity. We present the properties of thecoefficients of F L σ formula, which ensure the stability and the convergence of the proposedscheme. Numerical experiments are provided to show the sharp decrease in the CPU time andstorage of the fast algorithm with the same accuracy as the direct method.The paper is organized as follows. In Section 2, we present F L σ formula approaching theVO Caputo fractional derivative utilizing the ESA technique. Then we construct F L σ schemefor multi-dimensional VO fractional sub-diffusion problems in Section 3 and investigate the stabilityand convergence. In Section 4, reliability and efficiency are confirmed by some numerical examples.Concluding remarks are given in Section 5. We firstly introduce L σ formula [11] for the VO Caputo fractional derivative. Consider theVO Caputo fractional derivative defined by (1) with t ∈ (0 , T ] ( T ≥ n , let ∆t = T /n and t k = k∆t for k = 0 , , . . . , n . We denote t k + σ k = t k + σ ( t k ) ∆t for k = 0 , , . . . , n − σ ( t k ) is determined by the equation [11] F ( σ ) = σ − (cid:16) − α ( t k + σ∆t ) (cid:17) = 0 , which has a unique root σ ( t k ) ∈ ( ,
1) being conveniently calculated by Newton’s method. In thefollowing, we denote σ k = σ ( t k ), α k + σ k = α ( t k + σ k ).From the definition (1), the VO Caputo derivative at the point t k + σ k is presented as C D α ( t ) t u ( t ) | t = t k + σk = 1 Γ (cid:0) − α k + σ k (cid:1) Z t k + σk u ′ ( τ )( t k + σ k − τ ) α k + σk d τ. (2)To discretize the VO Caputo fractional derivative, we denote the quadric interpolation on the timeinterval [ t k − , t k ] with 1 ≤ k ≤ n − t k , t k +1 ]with 0 ≤ k ≤ n − L u ( τ ) = X p = − u ( t k + p ) Y q = − τ − t k + q t k + p − t k + q , τ ∈ [ t k − , t k ] ,L u ( τ ) = u ( t k ) τ − t k +1 t k − t k +1 + u ( t k +1 ) τ − t k t k +1 − t k , τ ∈ [ t k , t k +1 ] . Jia-li Zhang et al. L σ formulaBased on the above interpolation polynomials, at time t k + σ k with k = 0 , , . . . , n − L σ approximation formula to (2) is obtained as [11] H D α k + σk t u ( t k + σ k ) = 1 Γ (1 − α k + σ k ) Z t k (cid:0) L u ( τ ) (cid:1) ′ ( t k + σ k − τ ) α k + σk d τ + Z t k + σk t k (cid:0) L u ( τ ) (cid:1) ′ ( t k + σ k − τ ) α k + σk d τ ! = s ( k ) k X l =0 g ( k ) l (cid:0) u ( t k − l +1 ) − u ( t k − l ) (cid:1) , (3)where g (0)0 = σ − α σ , and for k ≥ g ( k ) l = ∆t α k + σk − (1 − α k + σ k ) · R t k t k − τ − t k − / ( t k + σk − τ ) αk + σk d τ + R t k + σk t k ∆t ( t k + σk − τ ) αk + σk d τ, l = 0 , R t k − l t k − l − τ − t k − l − / ( t k + σk − τ ) αk + σk d τ + R t k − l +1 t k − l t k − l +3 / − τ ( t k + σk − τ ) αk + σk d τ, ≤ l ≤ k − , R t t t / − τ ( t k + σk − τ ) αk + σk d τ, l = k, and s ( k ) = ∆t − α k + σk Γ (2 − α k + σ k ) . The following lemma shows the local truncation error of L σ approximation formula for theVO Caputo fractional derivative. Lemma 1 [11] Suppose u ∈ C (cid:0) [0 , t k +1 ] (cid:1) . Let r k = (cid:12)(cid:12)(cid:12) C D α ( t ) t u ( t ) | t = t k + σk − H D α k + σk t u ( t k + σ k ) (cid:12)(cid:12)(cid:12) . Then, we have (cid:12)(cid:12) r k (cid:12)(cid:12) ≤ max ≤ t ≤ t k +1 (cid:12)(cid:12) u (3) ( t ) (cid:12)(cid:12) σ − α k + σk k Γ (1 − α k + σ k ) (cid:18)
112 + σ k − α k + σ k ) (cid:19) ∆t − α k + σk . Utilizing L σ approximation formula to calculate the value at the current time level, it needs tocompute the summation of a series including the values of all previous time levels. Therefore, L σ approximation formula requires O ( n ) storage and O ( n ) computational complexity. It inspires us toconstruct a fast algorithm to approach L σ approximation formula (3).2.2 F L σ formulaNow, we develop a fast high-order numerical formula ( F L σ formula) for the VO Caputofractional derivative. The kernel function in the VO Caputo fractional derivative is approximatedby the ESA technique, which is stated in [5, 41]. itle Suppressed Due to Excessive Length 5 Lemma 2
For any constant α k + σ k ∈ [ α, α ] ⊂ (0 , , < ∆tT ≤ t k + σk − τT ≤ for τ ∈ [0 , t k − ] , ≤ k ≤ n − and the expected accuracy < ǫ ≤ /e , there exist a constant h , integers N and N ,which satisfy h = 2 π log 3 + α log(cos 1) − + log ǫ − ,N = (cid:24) h α (cid:0) log ǫ + log Γ (1 + α ) (cid:1)(cid:25) , (4) N = (cid:22) h (cid:18) log T∆t + log log ǫ − + log α + 2 − (cid:19)(cid:23) , such that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:18) t k + σ k − τT (cid:19) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:18) t k + σ k − τT (cid:19) − α k + σk ǫ, (5) where the quadrature exponents and weights are given by λ i = e ih , θ ( k ) i = he α k + σk ih Γ ( α k + σ k ) . Now L σ formula (3) can be approximated by H D α k + σk t u ( t k + σ k )= 1 Γ (1 − α k + σ k ) T − α k + σk Z t k (cid:0) L u ( τ ) (cid:1) ′ (cid:18) t k + σ k − τT (cid:19) − α k + σk d τ + Z t k + σk t k (cid:0) L u ( τ ) (cid:1) ′ ( t k + σ k − τ ) α k + σk d τ ! ≈ Γ (1 − α k + σ k ) T − α k + σk Z t k (cid:0) L u ( τ ) (cid:1) ′ N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T d τ + Z t k + σk t k (cid:0) L u ( τ ) (cid:1) ′ ( t k + σ k − τ ) α k + σk d τ ! = T − α k + σk Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i Z t k (cid:0) L u ( τ ) (cid:1) ′ e − λ i ( t k + σk − τ ) /T d τ + u ( t k +1 ) − u ( t k ) ∆tΓ (1 − α k + σ k ) Z t k + σk t k t k + σ k − τ ) α k + σk d τ = T − α k + σk Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i H ( k ) i + s ( k ) σ − α k + σk k (cid:0) u ( t k +1 ) − u ( t k ) (cid:1) , (6) where H ( k ) i is the history part of the integral, H ( k ) i = Z t k (cid:0) L u ( τ ) (cid:1) ′ e − λ i ( t k + σk − τ ) /T d τ. We point out that H ( k ) i can be calculated by a recursive relation and quadratic interpolation; i.e., H ( k ) i = e − λ i (1+ σ k − σ k − ) ∆t/T H ( k − i + Z t k t k − (cid:0) L u ( τ ) (cid:1) ′ e − λ i ( t k + σk − τ ) /T d τ = e − λ i (1+ σ k − σ k − ) ∆t/T H ( k − i + A ( k ) i (cid:0) u ( t k ) − u ( t k − ) (cid:1) + B ( k ) i (cid:0) u ( t k +1 ) − u ( t k ) (cid:1) , (7)with H (0) i = 0 ( i = N + 1 , . . . , N ) and A ( k ) i = Z (cid:0) / − τ (cid:1) e − λ i ( σ k +1 − τ ) ∆t/T d τ, B ( k ) i = Z (cid:0) τ − / (cid:1) e − λ i ( σ k +1 − τ ) ∆t/T d τ. Jia-li Zhang et al. Overall, for 0 ≤ k ≤ n −
1, we define
F L σ formula for C D α ( t ) t u ( t ) | t = t k + σk by F H D α k + σk t u ( t k + σ k )given as F H D α k + σk t u ( t k + σ k ) = T − α k + σk Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i H ( k ) i + s ( k ) σ − α k + σk k (cid:0) u ( t k +1 ) − u ( t k ) (cid:1) , (8)where H ( k ) i is calculated by (7).In addition, for k = 0, we let F H D α σ t u ( t σ ) = H D α σ t u ( t σ ) . (9)We give the following lemma to state the local truncated error of F L σ formula (8) for the VOCaputo fractional derivative C D α ( t ) t u ( t ) | t = t k + σk . Lemma 3
Suppose α k + σ k ∈ (0 , and u ( t ) ∈ C (cid:0) [0 , T ] (cid:1) . Let L - σ formula be as in (3), F L - σ formula be defined by (8) – (9) and ǫ be the expected accuracy. Then, we have C D α ( t ) t u ( t ) | t = t k + σk = F H D α k + σk t u ( t k + σ k ) + O ( ∆t − α k + σk + ǫ ) , ≤ k ≤ n − . (10) Proof
Clearly, (9) and Lemma 1 imply the lemma is valid for k = 0. For k ≥
1, according to (6) and(8), H D α k + σk t u ( t k + σ k ) and F H D α k + σk t u ( t k + σ k ) just differ in the approximation for (cid:16) t k + σk − τT (cid:17) − α k + σk , τ ∈ (0 , t k + σ k ); i.e., (cid:12)(cid:12)(cid:12) C D α ( t ) t u ( t ) | t = t k + σk − F H D α k + σk t u ( t k + σ k ) (cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) C D α ( t ) t u ( t ) | t = t k + σk − H D α k + σk t u ( t k + σ k ) (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) H D α k + σk t u ( t k + σ k ) − F H D α k + σk t u ( t k + σ k ) (cid:12)(cid:12)(cid:12) = O ( ∆t − α k + σk )+ T − α k + σk Γ (1 − α k + σ k ) k X l =1 Z t l t l − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:18) t k + σ k − τT (cid:19) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:0) L u ( τ ) (cid:1) ′ d τ = O ( ∆t − α k + σk ) + t − α k + σk k max ≤ t ≤ t k | u ′ ( t ) | Γ (2 − α k + σ k ) O ( ǫ ) . Since t k + σ k ≤ T and u ( t ) ∈ C (cid:0) [0 , T ] (cid:1) , the proof is completed.From Lemma 3, we conclude that the approximation for the VO Caputo fractional derivative C D α ( t ) t u ( t ) at the point t k + σ k by F H D α k + σk t u ( t k + σ k ) is at least second-order convergent if u ∈C (cid:0) [0 , T ] (cid:1) .2.3 Properties of discrete kernelsIn order to prepare the subsequent analysis, we firstly show the properties of discrete kernels in F L σ formula. The recursive relation of H ( k ) i defined in (7) can be equivalently rewritten as the itle Suppressed Due to Excessive Length 7 following form H ( k ) i = k X l =1 Z t l t l − (cid:0) L u ( τ ) (cid:1) ′ e − λ i ( t k + σk − τ ) /T d τ = ∆t − (cid:26) Z t t ( t − τ ) e − λ i ( t k + σk − τ ) /T d τ (cid:0) u ( t ) − u ( t ) (cid:1) + k − X l =1 "Z t l t l − ( τ − t l − ) e − λ i ( t k + σk − τ ) /T d τ + Z t l +1 t l ( t l + − τ ) e − λ i ( t k + σk − τ ) /T d τ u ( t l +1 ) − u ( t l ) (cid:1) + Z t k t k − ( τ − t k − ) e − λ i ( t k + σk − τ ) /T d τ (cid:0) u ( t k +1 ) − u ( t k ) (cid:1)(cid:27) , which leads to F H D α k + σk t u ( t k + σ k ) (0 ≤ k ≤ n −
1) in the following
F H D α k + σk t u ( t k + σ k )= T − α k + σk ∆t − Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i (cid:26) Z t t ( t − τ ) e − λ i ( t k + σk − τ ) /T d τ (cid:0) u ( t ) − u ( t ) (cid:1) + k − X l =1 "Z t l t l − ( τ − t l − ) e − λ i ( t k + σk − τ ) /T d τ + Z t l +1 t l ( t l + − τ ) e − λ i ( t k + σk − τ ) /T d τ u ( t l +1 ) − u ( t l ) (cid:1) + Z t k t k − ( τ − t k − ) e − λ i ( t k + σk − τ ) /T d τ (cid:0) u ( t k +1 ) − u ( t k ) (cid:1)(cid:27) + s ( k ) σ − α k + σk k (cid:0) u ( t k +1 ) − u ( t k ) (cid:1) = s ( k ) k X l =0 ρ ( k ) l (cid:0) u ( t k − l +1 ) − u ( t k − l ) (cid:1) , (11) where ρ (0)0 = σ − α σ , (12)and for 1 ≤ k ≤ n − ρ ( k ) l = ∆t α k + σk − (1 − α k + σ k ) T α k + σk · R t k t k − ( τ − t k − ) N P i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T d τ + T α k + σk R t k + σk t k ∆t ( t k + σk − τ ) αk + σk d τ, l = 0 , R t k − l t k − l − ( τ − t k − l − ) N P i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T d τ + R t k − l +1 t k − l ( t k − l + − τ ) N P i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T d τ, ≤ l ≤ k − , R t t ( t − τ ) N P i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) /T d τ, l = k. (13) The coefficients { ρ ( k ) l | ≤ l ≤ k ; 0 ≤ k ≤ n − } defined in (12) and (13) have the followingproperties, which play the vital roles in studying the stability and convergence of the finite differenceschemes. Below we present the relationship between ρ ( k ) l and g ( k ) l . Lemma 4
For α k + σ k ∈ (0 , , g ( k ) l (0 ≤ l ≤ k ; 0 ≤ k ≤ n − defined in (3) , ρ ( k ) l (0 ≤ l ≤ k ; 0 ≤ k ≤ n − defined in (12) and (13) with σ k = 1 − α k + σk . Then, we have ρ (0)0 = g (0)0 , Jia-li Zhang et al. and for ≤ k ≤ n − , (cid:12)(cid:12)(cid:12) ρ ( k ) l − g ( k ) l (cid:12)(cid:12)(cid:12) ≤ − α k + σ k ∆t α k + σk · ǫ , l = 0 , ǫ , ≤ l ≤ k − ,ǫ, l = k. (14) Proof
For k = 0, ρ (0)0 = g (0)0 is obvious. For k ≥
1, the ESA approximation (5) with (3) and (11)gives (cid:12)(cid:12)(cid:12) ρ ( k )0 − g ( k )0 (cid:12)(cid:12)(cid:12) ≤ ∆t α k + σk − (1 − α k + σ k ) T α k + σk · ( Z t k t k − (cid:12)(cid:12) τ − t k − (cid:12)(cid:12) · (cid:20)(cid:16) t k + σ k − τT (cid:17) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) ∆t/T (cid:21) d τ ) ≤ ∆t α k + σk − (1 − α k + σ k ) ǫ · Z t k t k − (cid:12)(cid:12) τ − t k − (cid:12)(cid:12) ( t k + σ k − τ ) − α k + σk d τ =(1 − α k + σ k ) ǫ · Z (cid:12)(cid:12) τ − / (cid:12)(cid:12) ( σ k − τ + 1) − α k + σk d τ ≤ (1 − α k + σ k ) ∆t − α k + σk ǫ · Z (cid:12)(cid:12) τ − / (cid:12)(cid:12) d τ = ǫ − α k + σ k ) ∆t − α k + σk . For 1 ≤ l ≤ k −
1, we obtain (cid:12)(cid:12)(cid:12) ρ ( k ) l − g ( k ) l (cid:12)(cid:12)(cid:12) ≤ ∆t α k + σk − (1 − α k + σ k ) T α k + σk · ( Z t k − l t k − l − (cid:12)(cid:12) τ − t k − l − (cid:12)(cid:12) · (cid:20)(cid:16) t k + σ k − τT (cid:17) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) ∆t/T (cid:21) d τ + Z t k − l +1 t k − l (cid:12)(cid:12) t k − l + − τ (cid:12)(cid:12)(cid:20)(cid:16) t k + σ k − τT (cid:17) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) ∆t/T (cid:21) d τ ) ≤ ∆t α k + σk − (1 − α k + σ k ) ǫ · ( Z t k − l t k − l − (cid:12)(cid:12) τ − t k − l − (cid:12)(cid:12) ( t k + σ k − τ ) − α k + σk d τ + Z t k − l +1 t k − l (cid:12)(cid:12) t k − l + − τ (cid:12)(cid:12) ( t k + σ k − τ ) − α k + σk d τ ) =(1 − α k + σ k ) ǫ · n Z (cid:12)(cid:12) τ − / (cid:12)(cid:12) ( σ k − τ + l + 1) − α k + σk d τ + Z (cid:12)(cid:12) / − τ (cid:12)(cid:12) ( σ k − τ + l ) − α k + σk d τ o ≤ (1 − α k + σ k ) ∆t − α k + σk ǫ · n Z (cid:12)(cid:12) / − τ (cid:12)(cid:12) d τ + Z (cid:12)(cid:12) τ − / (cid:12)(cid:12) d τ o = 5 ǫ − α k + σ k ) ∆t − α k + σk . itle Suppressed Due to Excessive Length 9 Similarly, we have (cid:12)(cid:12)(cid:12) ρ ( k ) k − g ( k ) k (cid:12)(cid:12)(cid:12) ≤ ∆t α k + σk − (1 − α k + σ k ) T α k + σk · ( Z t t (cid:12)(cid:12) t − τ (cid:12)(cid:12)(cid:20)(cid:16) t k + σ k − τT (cid:17) − α k + σk − N X i = N +1 θ ( k ) i e − λ i ( t k + σk − τ ) ∆t/T (cid:21) d τ ) ≤ ∆t α k + σk − (1 − α k + σ k ) ǫ · Z t t (cid:12)(cid:12) t − τ (cid:12)(cid:12) ( t k + σ k − τ ) − α k + σk d τ =(1 − α k + σ k ) ǫ · Z (cid:12)(cid:12) / − τ (cid:12)(cid:12) ( k + σ k − τ ) − α k + σk d τ ≤ (1 − α k + σ k ) ∆t − α k + σk ǫ · Z (cid:12)(cid:12) / − τ (cid:12)(cid:12) d τ = ǫ (1 − α k + σ k ) ∆t − α k + σk . The proof is complete.
Lemma 5
For α k + σ k ∈ (0 , , ρ ( k ) l (0 ≤ l ≤ k ; 0 ≤ k ≤ n − are defined in (12) and (13) with σ k = 1 − α k + σk and a sufficiently small ǫ satisfying ǫ ≤ − α )(2 − α ) − α ∆t α (6 − α )(1 − α ) . Then, we have < (1 − ǫ )(1 − α k + σ k )2( k + σ k ) α k + σk < ρ ( k ) k < ρ ( k ) k − < . . . < ρ ( k )0 , (15)(2 σ k − ρ ( k )0 − σ k ρ ( k )1 ≥ . (16) Proof
For k = 0, ρ (0)0 = g (0)0 > k ≥
1, from (13), the coefficients can be rewrittenas the following form ρ ( k ) l = ∆t α k + σk (1 − α k + σ k ) T α k + σk N P i = N +1 θ ( k ) i B ( k ) i + ∆t − T α k + σk R t k + σk t k t k + σk − τ ) αk + σk d τ, l = 0 , N P i = N +1 θ ( k ) i (cid:16) e − λ i ( l − ∆t/T A ( k ) i + e − λ i l∆t/T B ( k ) i (cid:17) , ≤ l ≤ k − , N P i = N +1 θ ( k ) i e − λ i ( k − ∆t/T A ( k ) i , l = k, where A ( k ) i and B ( k ) i are defined by (7). Thanks to θ ( k ) i > λ i > e τλ i ∆t/T with respect to τ , A ( k ) i > B ( k ) i > < ρ ( k ) k < ρ ( k ) k − < ρ ( k ) k − < . . . < ρ ( k )1 . According to g ( k ) k ≥ (1 − α k + σk )2( k + σ k ) αk + σk [11], (3),(5), and (13), we have ρ ( k ) k ≥ (1 − ǫ ) g ( k ) k ≥ (1 − ǫ )(1 − α k + σ k )2( k + σ k ) α k + σk . et al. So condition (15) will hold if ρ ( k )0 > ρ ( k )1 holds, an inequality which obviously follows from thecondition (16). Thus it is enough to prove the later one. By (14), the left-hand side of the condition(16) satisfies(2 σ k − ρ ( k )0 − σ k ρ ( k )1 ≥ (2 σ k − g ( k )0 − σ k g ( k )1 − (2 σ k −
1) (1 − α k + σ k ) ∆t − α k + σk ǫ − (1 − α k + σ k ) ∆t − α k + σk ǫ (17)for k = 1, and(2 σ k − ρ ( k )0 − σ k ρ ( k )1 ≥ (2 σ k − g ( k )0 − σ k g ( k )1 − (2 σ k −
1) (1 − α k + σ k ) ∆t − α k + σk ǫ − σ k − α k + σ k ) ∆t − α k + σk ǫ k ≥ g ( k ) l is(2 σ k − g ( k )0 − σ k g ( k )1 = ( (2 σ k − − σ k )2 σ k (1 + σ k ) − α k + σk , k = 1 , (1 + σ k ) − α k + σk (cid:16) σ k − σ k − (cid:0) σ k σ k (cid:1) − α k + σk (cid:17) , k ≥ . ≥ (2 σ k − − σ k )2 σ k (1+ σ k ) αk + σk − , k = 1 , (2 σ k − − σ k )2 σ k (1+ σ k ) αk + σk , k ≥ ,> . Thus in order to make the condition (16) hold we just need(2 σ k − − σ k )2 σ k (1 + σ k ) α k + σk − − (2 σ k −
1) (1 − α k + σ k ) ∆t − α k + σk ǫ − (1 − α k + σ k ) ∆t − α k + σk ǫ ≥ , k = 1 , and(2 σ k − − σ k )2 σ k (1 + σ k ) α k + σk − (2 σ k −
1) (1 − α k + σ k ) ∆t − α k + σk ǫ − σ k − α k + σ k ) ∆t − α k + σk ǫ ≥ , k ≥ , that is ǫ ≤ σ k − ∆t α k + σk σ k (1 + σ k ) α k + σk − min (cid:26) σ k − , σ k − (cid:27) . In F L σ formula, σ k = 1 − α k + σk . Thus we have ǫ ≤ − α k + σ k ) ∆t α k + σk (6 − α k + σ k )(1 − α k + σ k )(2 − α k + σ k ) α k + σk − . So we get (16) if this ǫ satisfies ǫ ≤ − α )(2 − α ) − α ∆t α (6 − α )(1 − α ) ≤ − α k + σ k ) ∆t α k + σk (6 − α k + σ k )(1 − α k + σ k )(2 − α k + σ k ) α k + σk − . The proof is complete. itle Suppressed Due to Excessive Length 11
In this section, based on
F L σ formula in Section 2, we construct a fast temporal second-order and spatial fourth-order finite difference method ( F L σ scheme) for the VO time-fractionalsub-diffusion equations. The unconditional stability and convergence of the difference method areinvestigated.We consider the multi-dimensional VO time-fractional sub-diffusion equations C D α ( t ) t u ( x , t ) = ∆u ( x , t ) + f ( x , t ) , x ∈ Ω, t ∈ (0 , T ] , (19) u ( x ,
0) = ϕ ( x ) , x ∈ Ω, (20) u ( x , t ) = 0 , x ∈ ∂Ω, t ∈ [0 , T ] , (21)where Ω = Q dk =1 ( a ( k ) l , a ( k ) r ) ⊂ R d , ∂Ω is the boundary of Ω , Ω = Ω ∪ ∂Ω , x = ( x (1) , x (2) , . . . , x ( d ) ) ∈ Ω , ∆u ( x , t ) = P dk =1 ∂ x ( k ) u ( x , t ), f ( x , t ) and ϕ ( x ) represent the given sufficiently smooth functions.Denote L ( k ) = a ( k ) r − a ( k ) l , k = 1 , , . . . , d . Let m ( k ) be positive integers. Define a uniform partitionof Ω by x ( k ) j ( k ) = a ( k ) l + j ( k ) ∆x ( k ) ( j ( k ) = 0 , , . . . , m ( k ) ; k = 1 , , . . . , d ) for ∆x ( k ) = L ( k ) /m ( k ) anddenote ∆x = max ≤ k ≤ d ∆x ( k ) . Let Ω h = { ( x (1) j (1) , x (2) j (2) , . . . , x ( d ) j ( d ) ) | j ( k ) = 1 , , . . . , m ( k ) − k = 1 , , . . . , d } , Ω h = { ( x (1) j (1) , x (2) j (2) , . . . , x ( d ) j ( d ) ) | j ( k ) = 0 , , . . . , m ( k ) ; k = 1 , , . . . , d } , and ∂Ω h = Ω h \ Ω h . Denote theindex vector j = ( j (1) , j (2) , . . . , j ( d ) ) and the spatial point x j = (cid:0) x (1) j (1) , x (2) j (2) , . . . , x ( d ) j ( d ) (cid:1) , then we definethe index space J = { j | x j ∈ Ω h } , J = { j | x j ∈ Ω h } , ∂ J = { j | x j ∈ ∂Ω h } . Thus the grid function spaces are defined by U = { u | u being a grid function on Ω h } , ˚ U = { u | u ∈ U ; u j = 0 when j ∈ ∂ J } . We introduce the following discrete operators in the grid space U δ k u j + δ k = 1 ∆x ( k ) ( u j + δ k − u j ) , δ k u j = 1( ∆x ( k ) ) ( u j + δ k − u j + u j − δ k ) , where δ k = (0 , . . . , , . . . ,
0) is an index with 1 at the k -th position and 0 at other positions. A k u j = ( ( u j − δ k + 10 u j + u j + δ k ) , j ∈ J ,u j , j ∈ ∂ J , (22) ∆ h u j = d X k =1 δ k u j , A h u j = d Y k =1 A k u j , Λ h u j = d X k =1 d Y l =1 ,l = k A l δ k u j . (23) et al. In the grid function space ˚ U , define the discrete inner products and norms( u, w ) = d Y r =1 ∆x ( r ) ! X j ∈J u j w j , k u k = p ( u, u ) , ( u, w ) A h = ( A h u, w ) , k u k A h = p ( u, A h u ) , ( δ k u, δ k w ) = d Y r =1 ∆x ( r ) ! m ( k ) − X j ( k ) =0 d Y l =1 ,l = k m ( l ) − X j ( l ) =1 ( δ k u j + δ k )( δ k w j + δ k ) , | u | ,k = p ( δ k u, δ k u ) , | u | = vuut d X k =1 | u | ,k , ( δ k u, δ k w ) = d Y r =1 ∆x ( r ) ! X j ∈J ( δ k u j )( δ k w j ) , | u | ,k = q ( δ k u, δ k u ) , k u k ∞ = max j ∈J | u j | . F L σ schemeBefore deriving the finite difference schemes for the problem (19)–(21), we denote the numericalsolution at time t k with spatial point x j by u kj , and f ( x j , t k + σ k ) by f k + σ k j for x j ∈ Ω h and 0 ≤ k ≤ n .We assume that the solution u ∈ C (6 , (cid:0) Ω × [0 , T ] (cid:1) .Next we recall L σ scheme for the problem (19)–(21). Considering the equation (19) at( x j , t k + σ k ), we have C D α ( t ) t u ( x j , t ) | t = t k + σk = ∆u ( x j , t k + σ k ) + f k + σ k j , j ∈ J , ≤ k ≤ n − . (24)By Lemma 3, the term on the left-hand side of (24) satisfies C D α ( t ) t u ( x j , t ) | t = t k + σk = F H D α k + σk t u ( x j , t k + σ k ) + O ( ∆t − α k + σk + ǫ ) , j ∈ J , ≤ k ≤ n − . (25)On the other hand, using Lemmas 3.4 and 3.5 in [11] , the first term on the right-hand side of (24)satisfies ∆u ( x j , t k + σ k ) = σ k ∆u ( x j , t k +1 ) + (1 − σ k ) ∆u ( x j , t k ) + O ( ∆t ) . (26)Then substituting (25) and (26) into (24) gives F H D α k + σk t u ( x j , t k + σ k ) = σ k ∆u ( x j , t k +1 ) + (1 − σ k ) ∆u ( x j , t k ) + f k + σ k j + O ( ∆t + ǫ ) ,j ∈ J , ≤ k ≤ n − . Acting the averaging operator A h in (22) on both hand sides of the equality above, we obtain A h F H D α k + σk t u ( x j , t k + σ k )= σ k A h ∆u ( x j , t k +1 ) + (1 − σ k ) A h ∆u ( x j , t k ) + A h f k + σ k j + O ( ∆t + ǫ ) , j ∈ J , ≤ k ≤ n − , itle Suppressed Due to Excessive Length 13 where F H D α k + σk t u ( x j , t k + σ k ) = T − α k + σk Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i H ( k ) i + s ( k ) σ − α k + σk k (cid:0) u ( x j , t k +1 ) − u ( x j , t k ) (cid:1) , with H ( k ) i = e − λ i (1+ σ k − σ k − ) ∆t/T H ( k − i + A ( k ) i (cid:0) u ( x j , t k ) − u ( x j , t k − ) (cid:1) + B ( k ) i (cid:0) u ( x j , t k +1 ) − u ( x j , t k ) (cid:1) . By the fact of [11] A h ∆u ( x j , t k ) = Λ h u ( x j , t k ) + O ( ∆x ) , j ∈ J , ≤ k ≤ n, thus A h F H D α k + σk t u ( x j , t k + σ k ) = Λ h (cid:0) σ k u ( x j , t k +1 ) + (1 − σ k ) u ( x j , t k ) (cid:1) + A h f k + σ k j + S kj ,j ∈ J , ≤ k ≤ n − , (27)where there exists a constant c such that | S kj | ≤ c ( ∆t + ∆x + ǫ ) , j ∈ J , ≤ k ≤ n − . (28)From the initial and boundary value conditions (20)–(21), we have u ( x j ,
0) = ϕ ( x j ) , j ∈ J , (29) u ( x j , t k ) = 0 , j ∈ ∂ J , ≤ k ≤ n. (30)Omitting the small term S kj in (27), we construct F L σ scheme for the problem (19)–(21) asfollows A h F H D α k + σk t u k + σ k j = Λ h (cid:0) σ k u k +1 j + (1 − σ k ) u kj (cid:1) + A h f k + σ k j , j ∈ J , ≤ k ≤ n − , (31) u j = ϕ ( x j ) , j ∈ J , (32) u kj = 0 , j ∈ ∂ J , ≤ k ≤ n, (33)where F H D α k + σk t u k + σ k j = T − α k + σk Γ (1 − α k + σ k ) N X i = N +1 θ ( k ) i H ( k ) i + s ( k ) σ − α k + σk k (cid:0) u k +1 j − u kj (cid:1) , with H ( k ) i = e − λ i (1+ σ k − σ k − ) ∆t/T H ( k − i + A ( k ) i (cid:0) u kj − u k − j (cid:1) + B ( k ) i (cid:0) u k +1 j − u kj (cid:1) . Recall L σ scheme for the problem (19)–(21) as follows [11] A h H D α k + σk t u k + σ k j = Λ h (cid:0) σ k u k +1 j + (1 − σ k ) u kj (cid:1) + A h f k + σ k j , j ∈ J , ≤ k ≤ n − , (34) u j = ϕ ( x j ) , j ∈ J , (35) u kj = 0 , j ∈ ∂ J , ≤ k ≤ n, (36)where H D α k + σk t u k + σ k j = s ( k ) k X l =0 g ( k ) l (cid:0) u k − l +1 j − u k − lj (cid:1) . et al. F L σ schemeAs described in Lemma 5, the coefficients ρ ( k ) l hold the vital properties for the stability andconvergence analysis. Thus, similar to the proof given in [11], we present the following lemmaswhich will be used in the analysis of F L σ scheme (31)–(33). Lemma 6 [1] Let ˚ U be an inner product space and h· , ·i ∗ is the inner product with the induced norm k·k ∗ . Suppose (cid:8) c ( k ) l | ≤ l ≤ k, k ≥ (cid:9) satisfies < c ( k ) k < c ( k ) k − < . . . < c ( k )0 . For v , v , . . . , v k +1 ∈ ˚ U ,we have the following inequality k X l =0 c ( k ) l h v k − l +1 − v k − l , σ k v k +1 + (1 − σ k ) v k i ∗ ≥ k X l =0 c ( k ) l (cid:16) k v k − l +1 k ∗ − k v k − l k ∗ (cid:17) . Lemma 7 [11] For any u, w ∈ ˚ U , define h u, w i A h = (cid:0) A h u, − ∆ h w (cid:1) . Then h u, w i A h is an inner product on ˚ U . We denote | u | , A h = p h u, u i A h . Lemma 8 [11] For any u ∈ ˚ U , we have (cid:0) (cid:1) d | u | ≤ | u | , A h ≤ | u | . Lemma 9 [11] For any u ∈ ˚ U , we have (cid:0) (cid:1) d − k ∆ h u k ≤ (cid:0) Λ h u, ∆ h u (cid:1) ≤ k ∆ h u k . Now we obtain the following theorem to state the unconditional stability of the proposed scheme.
Theorem 1
Let (cid:8) u kj | j ∈ J , ≤ k ≤ n (cid:9) be the solution of the difference scheme (31) – (33) . Then,we have (cid:12)(cid:12) u k (cid:12)(cid:12) , A h ≤ (cid:12)(cid:12) u (cid:12)(cid:12) , A h + (cid:0) (cid:1) d − · − ǫ c max ≤ l ≤ k − kA h f l + σ l k , ≤ k ≤ n, (37) where kA h f l + σ l k = d Y r =1 ∆x ( r ) ! X j ∈J (cid:16) A h f l + σ l j (cid:17) , c = max ≤ t ≤ T n t α ( t ) Γ (cid:0) − α ( t ) (cid:1)o . Proof
Denote u k + σ k = σ k u k +1 + (1 − σ k ) u k . Making an inner product with − ∆ h u k + σ k on both handsides of (34), we obtain (cid:0) A h F H D α k + σk t u k + σ k , − ∆ h u k + σ k (cid:1) + (cid:0) Λ h u k + σ k , ∆ h u k + σ k (cid:1) = − (cid:0) A h f k + σ k , ∆ h u k + σ k (cid:1) , ≤ k ≤ n − . (38) itle Suppressed Due to Excessive Length 15 Noticing (33), by Lemmas 6 and 7, we get k X l =0 ρ ( k ) l h u k − l +1 − u k − l , σ k u k +1 + (1 − σ k ) u k i A h ≥ k X l =0 ρ ( k ) l (cid:16) | u k − l +1 | , A h − | u k − l | , A h (cid:17) . By Lemma 9, the second term on the left-hand side of (38) follows that (cid:0) Λ h u k + σ k , ∆ h u k + σ k (cid:1) ≥ (cid:0) (cid:1) d − k ∆ h u k + σ k k . Substituting the two inequalities into (38) and using the Cauchy-Schwarz inequality with the basicinequality ab ≤ (cid:0) (cid:1) d − a + (cid:0) (cid:1) d − b , we have s ( k ) k X l =0 ρ ( k ) l (cid:16) | u k − l +1 | , A h − | u k − l | , A h (cid:17) + (cid:0) (cid:1) d − k ∆ h u k + σ k k ≤ − (cid:0) A h f k + σ k , ∆u k + σ k (cid:1) ≤ kA h f k + σ k kk ∆u k + σ k k≤ (cid:0) (cid:1) d − k ∆u k + σ k k + (cid:0) (cid:1) d − kA h f k + σ k k , ≤ k ≤ n − , which follows that s ( k ) k X l =0 ρ ( k ) l (cid:16) | u k − l +1 | , A h − | u k − l | , A h (cid:17) ≤ (cid:0) (cid:1) d − kA h f k + σ k k , ≤ k ≤ n − . Further we get ρ ( k )0 | u k +1 | , A h ≤ k − X l =0 (cid:16) ρ ( k ) l − ρ ( k ) l +1 (cid:17) | u k − l | , A h + ρ ( k ) k | u | , A h + (cid:0) (cid:1) d − s ( k ) kA h f k + σ k k , ≤ k ≤ n − . According to Lemma 5, it follows that12 s ( k ) ρ ( k ) k < − ǫ ( k + σ k ) α k + σ k s ( k ) (1 − α k + σ k ) = 11 − ǫ t α ( t ) Γ (cid:0) − α ( t ) (cid:1) | t = t k + σk ≤ − ǫ c , which leads to ρ ( k )0 | u k +1 | , A h ≤ k − X l =0 (cid:16) ρ ( k ) l − ρ ( k ) l +1 (cid:17) | u k − l | , A h + ρ ( k ) k (cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c kA h f k + σ k k (cid:17) , ≤ k ≤ n − . (39)Next the mathematical induction will be used to prove the inequality (37). It is valid for k = 0.Assume (37) is true for 0 ≤ k ≤ q , now we prove that (37) is valid for k = q + 1. From (39), we et al. obtain ρ q | u q +1 | , A h ≤ q − X l =0 (cid:16) ρ ( q ) l − ρ ( q ) l +1 (cid:17) | u q − l | , A h + ρ ( q ) q (cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c kA h f q + σ q k (cid:17) ≤ q − X l =0 (cid:16) ρ ( q ) l − ρ ( q ) l +1 (cid:17)(cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c max ≤ s ≤ q − l − kA h f s + σ s k (cid:17) + ρ ( q ) q (cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c kA h f q + σ q k (cid:17) ≤ (cid:16) q − X l =0 ( ρ ( q ) l − ρ ( q ) l +1 ) + ρ ( q ) q (cid:17)(cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c max ≤ s ≤ q kA h f s + σ s k (cid:17) ≤ ρ ( q )0 (cid:16) | u | , A h + (cid:0) (cid:1) d − − ǫ c max ≤ s ≤ q kA h f s + σ s k (cid:17) . From the above inequality, we obtain | u q +1 | , A h ≤ | u | , A h + (cid:0) (cid:1) d − − ǫ c max ≤ s ≤ q kA h f s + σ s k . Therefore, inequality (37) is valid for k = q + 1. This completes the proof.Theorem 1 reveals the stability of the difference scheme (31)–(33) with respect to the initial valueand the source term. The next theorem is about the convergence of F L σ scheme. Theorem 2
Let u ( x , t ) ∈ C (6 , (cid:0) Ω × [0 , T ] (cid:1) be the exact solution of the problem (19) – (21) , and (cid:8) u kj | j ∈ J , ≤ k ≤ n (cid:9) be the solution of the difference scheme (31) – (33) . Denote e kj = u ( x j , t k ) − u kj .Then, we have (cid:12)(cid:12) e k (cid:12)(cid:12) , A h ≤ vuut(cid:0) (cid:1) d − · − ǫ c d Y r =1 L ( r ) c (cid:0) ∆t + ∆x + ǫ (cid:1) , ≤ k ≤ n. Furthermore, we obtain (cid:12)(cid:12) e k (cid:12)(cid:12) ≤ vuut(cid:0) (cid:1) d − · − ǫ c d Y r =1 L ( r ) c (cid:0) ∆t + ∆x + ǫ (cid:1) , ≤ k ≤ n. Proof
Subtracting (31)–(33) from (27), (29)–(30), respectively, we obtain the error equations asfollows A h F H D α k + σk t e kj = Λ h (cid:0) σ k e k +1 j + (1 − σ k ) e kj (cid:1) + S kj , j ∈ J , ≤ k ≤ n − ,e j = 0 , j ∈ J ,e kj = 0 , j ∈ ∂ J , ≤ k ≤ n. Applying (28) and the priori estimation (37), it leads to (cid:12)(cid:12) e k (cid:12)(cid:12) , A h ≤ (cid:12)(cid:12) e (cid:12)(cid:12) , A h + (cid:0) (cid:1) d − · − ǫ c max ≤ l ≤ k − k S l k ≤ (cid:0) (cid:1) d − · − ǫ c d Y r =1 L ( r ) (cid:16) c (cid:0) ∆t + ∆x + ǫ (cid:1)(cid:17) , ≤ k ≤ n. itle Suppressed Due to Excessive Length 17 Consequently, we have (cid:12)(cid:12) e k (cid:12)(cid:12) , A h ≤ vuut(cid:0) (cid:1) d − · − ǫ c d Y r =1 L ( r ) c (cid:0) ∆t + ∆x + ǫ (cid:1) , ≤ k ≤ n. Furthermore, from Lemma 8, we obtain (cid:12)(cid:12) e k (cid:12)(cid:12) ≤ r(cid:0) (cid:1) d (cid:12)(cid:12) e k (cid:12)(cid:12) , A h ≤ vuut(cid:0) (cid:1) d − · − ǫ c d Y r =1 L ( r ) c (cid:0) ∆t + ∆x + ǫ (cid:1) , ≤ k ≤ n. In this section, some numerical examples are presented to verify the effectiveness of
F L σ scheme (31)–(33) compared with L σ scheme (34)–(36). All experiments are performed based onMatlab 2016b on a laptop with the configuration: Intel(R) Core(TM) i7-7500U CPU 2.70GHz and8.00 GB RAM. Example 1
We show the efficiency of
F L σ scheme (31)–(33) in 2D case comparing with thecorresponding L σ scheme (34)–(36). Take Ω = (0 , π ) × (0 , π ), T = 1, and the exact solution isgiven as [11] u ( x , t ) = ( t + 3 t + 1) sin x (1) sin x (2) . Based on the exact solution, we calculate the initial value and the source term ϕ ( x ) = sin x (1) sin x (2) ,f ( x , t ) = (cid:18) Γ (4 − α ( t )) t − α ( t ) + 6 Γ (3 − α ( t )) t − α ( t ) + 2( t + 3 t + 1) (cid:19) sin x (1) sin x (2) . We take ∆x (1) = ∆x (2) = ∆x , m (1) = m (2) = m , denote E ( ∆x, ∆t ) = (cid:13)(cid:13) u n − u ( x , t n ) (cid:13)(cid:13) ∞ ,Order t = log (cid:0) E ( ∆x, ∆t ) /E ( ∆x, ∆t ) (cid:1) ,Order x,t = log (cid:0) E (2 ∆x, ∆t ) /E ( ∆x, ∆t ) (cid:1) . In the calculations, we set the expected accuracy ǫ = ∆t and different spatial and temporalstep sizes. Note that the linear systems arising from the two-dimensional problems providecoefficient matrices possessing block-tridiagonal-Toeplitz with tridiagonal-Toeplitz-blocks, whichcan be diagonalized by the discrete sine transforms [6, 32]. Therefore, as a global after-processingoptimization, a fast implementation of the inversion, in our numerical experiments, is carried out bythe fast sine transforms to reduce the computational cost.Tables 1 and 2 list the results of Example 1. Fine spatial size is fixed at ∆x = π/
320 in Table1. Both L σ scheme and F L σ scheme can achieve the second-order temporal accuracy. Table1 also shows the developments of the CPU time and memory of the two schemes with respect to n .The CPU time of F L σ scheme is increasing half as fast as that of L σ scheme, which revealsthe O ( n log n ) and O ( n ) computational complexity of the two schemes, respectively. We note thatas n goes up, the memory of F L σ scheme grows slowly, while the storage requirement of L σ et al. Table 1
Convergence rates and the CPU time, memory of L σ scheme and F L σ scheme for Example 1 with α ( t ) = (2 + sin t ) / m = 320, ǫ = ∆t . L σ scheme F L σ scheme n E ( ∆x, ∆t ) Order t CPU(s) Memory E ( ∆x, ∆t ) Order t CPU(s) Memory2000 2.3592e-7 - 152.77 1.64e+9 2.3497e-7 - 194.70 1.01e+84000 5.8588e-8 2.01 642.33 3.27e+9 5.8411e-8 2.01 470.72 1.17e+88000 1.4339e-8 2.03 2510.36 6.52e+9 1.4319e-8 2.03 1021.68 1.34e+816000 - - - - 3.3034e-9 2.12 2311.46 1.53e+8
Table 2
Convergence rates and the CPU time, memory of L σ scheme and F L σ scheme for Example 1 with α ( t ) = (2 + sin t ) / n = m , ǫ = ∆t . L σ scheme F L σ scheme m E ( ∆x, ∆t ) Order x,t
CPU(s) Memory E ( ∆x, ∆t ) Order x,t
CPU(s) Memory20 1.1392e-6 - 0.15 1.20e+6 1.1971e-6 - 0.14 2.56e+540 7.2797e-8 3.97 3.23 1.96e+7 7.4374e-8 4.01 2.05 1.48e+680 4.6192e-9 3.98 185.32 3.20e+8 4.6405e-9 4.00 52.46 7.95e+6160 2.4931e-10 4.21 7158.69 5.18e+9 2.4589e-10 4.24 1073.38 4.15e+7 scheme increases linearly with n . Especially, when n = 16000, F L σ scheme is finished in 40 min,while F L σ scheme is failed due to excessive storage requirements. Table 2 gives the numericalresults with m varying from 20 to 160 and n = m . The convergence rates satisfy the theoreticalresult in Section 3. Nevertheless, a significant reduction is reflected in computational cost and storagememory of F L σ scheme on the refined mesh, compared with L σ scheme. Example 2
We also show the efficiency of
F L σ scheme (31)–(33) in 3D case. Take Ω = (0 , π ) × (0 , π ) × (0 , π ), T = 1, and the exact solution is given as u ( x , t ) = ( t + 3 t + 1) sin x (1) sin x (2) sin x (3) . Based on the exact solution, we calculate the initial value and the source term ϕ ( x ) = sin x (1) sin x (2) sin x (3) ,f ( x , t ) = (cid:18) Γ (4 − α ( t )) t − α ( t ) + 6 Γ (3 − α ( t )) t − α ( t ) + 3( t + 3 t + 1) (cid:19) sin x (1) sin x (2) sin x (3) . We take ∆x (1) = ∆x (2) = ∆x (3) = ∆x , m (1) = m (2) = m (3) = m , and set the expected accuracy ǫ = ∆t . The errors and convergence orders, CPU time and storage of L σ scheme and F L σ scheme are showed in Table 3 with m = 100 and temporal step sizes refined from ∆t = 1 /
200 to ∆t = 1 / m varies from 10 to 80 and n = (2 m ) . The fast sine transforms isused to reduce the computational cost.Tables 3 and 4 show that both F L σ scheme and L σ scheme can achieve the optimalconvergence for three-dimensional problems. Nevertheless, complexity of F L σ scheme is muchcheaper than that of L σ scheme on the refined mesh. For the fine mesh, however, L σ schemehas very low efficiency or even cannot proceed because of the limit of memory, while F L σ schemesolves the relevant problem very well. itle Suppressed Due to Excessive Length 19 Table 3
Convergence rates and the CPU time, memory of L σ scheme and F L σ scheme for Example 2 with α ( t ) = (2 + sin t ) / m = 100, ǫ = ∆t . L σ scheme F L σ scheme n E ( ∆x, ∆t ) Order t CPU(s) Memory E ( ∆x, ∆t ) Order t CPU(s) Memory200 2.6755e-5 - 41.64 1.64e+9 2.6587e-5 - 99.77 5.59e+8400 6.6637e-6 2.01 156.85 3.19e+9 6.6318e-6 2.00 233.83 6.6.e+8800 1.6528e-6 2.01 680.60 6.30e+9 1.6475e-6 2.01 564.19 7.84e+81600 - - - - 4.0177e-7 2.04 1358.80 9.24e+8
Table 4
Convergence rates and the CPU time, memory of L σ scheme and F L σ scheme for Example 2 with α ( t ) = (2 + sin t ) / n = (2 m ) , ǫ = ∆t . L σ scheme F L σ scheme m E ( ∆x, ∆t ) Order x,t
CPU(s) Memory E ( ∆x, ∆t ) Order x,t
CPU(s) Memory10 1.2679e-4 - 0.32 2.41e+6 1.2682e-4 - 0.52 5.06e+520 7.9012e-6 4.00 13.54 8.84e+7 7.9021e-6 4.00 12.61 6.56e+640 4.9351e-7 4.00 1047.68 3.04e+9 4.9351e-7 4.00 508.26 7.46e+780 - - - - 3.0832e-8 4.00 20550.89 8.01e+8
In this paper, we consider the fast high-order evaluation for the VO Caputo fractional derivative.
F L σ formula can efficiently reduce the computational storage and cost for the VO fractionalderivative. The fast formula is applied to construct F L σ scheme for the multi-dimensional VOfractional sub-diffusion equations. We show the properties of F L σ scheme to study the stabilityand convergence. Numerical examples are given to verify the theoretical results.To our best knowledge, there is a shortage of effective ways to solve the VO fractional problemwith a weak singularity at the initial time. We refer to [19,35] for the use of the nonuniform mesh sizes,which is a valid tool for the CO problems. Therefore, it is consequential to develop fast algorithmsfor the VO fractional problem with a weak singularity on nonuniform girds. We will consider thecase in the future. Acknowledgement
The authors would like to thank Professor Zhi-zhong Sun for his helpful suggestions during thepreparation of this paper.
References A. A. Alikhanov , A new difference scheme for the time fractional diffusion equation , J. Comput. Phys., 280(2015), pp. 424–438.2.
D. Baffet and J. S. Hesthaven , A kernel compression scheme for fractional differential equations , SIAM J.Numer. Anal., 55 (2017), pp. 496–520.3.
D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert , Application of a fractional advection dispersionequation , Water Resour. Res., 36 (2000), pp. 1403–1412.4.
D. Bertaccini and F. Durastante , Block structured preconditioners in tensor form for the all-at-once solutionof a finite volume fractional diffusion equation , Appl. Math. Lett., 95 (2019), pp. 92–97.0 Jia-li Zhang et al.5. G. Beylkin and L. Monz´on , Approximation by exponential sums revisited , Appl. Comput. Harmon. Anal., 28(2010), pp. 131–149.6.
D. Bini and F. Di Benedetto , A new preconditioner for the parallel solution of positive definite Toeplitz linearsystems , in Proceedings of the 2nd SPAA Conference, Crete, 1990, pp. 220–223.7.
A. V. Chechkin, R. Gorenflo, and I. M. Sokolov , Fractional diffusion in inhomogeneous media , J. Phys. A:Math. Gen., 38 (2005), pp. 679–684.8.
C. M. Chen, F. W. Liu, V. Anh, and I. Turner , Numerical schemes with high spatial accuracy for a variable-order amomalous subdiffusion equation , SIAM. J. Sci. Comput., 32 (2010), pp. 1740–1760.9.
C. F. M. Coimbra , Mechanics with variable-order differential operators , Ann. Phys., 12 (2003), pp. 692–703.10.
G. Diaz and C. F. M. Coimbra , Nonlinear dynamics and control of a variable order oscillator with applicationto the van der pol equation , Nonlin. Dyn., 56 (2009), pp. 145–157.11.
R. L. Du, A. A. Alikhanov, and Z. Z. Sun , Temporal second order difference schemes for the multi-dimensionalvariable-order time fractional sub-diffusion equations , Comput. Math. Appl., 79 (2020), pp. 952–2972.12.
Z. W. Fang, H. W. Sun, and H. Wang , A fast method for variable-order caputo fractional derivative withapplications to time-fractional diffusion equations , Comput. Math. Appl., 80 (2020), pp. 1443–1458.13.
H. F. Fu and H. Wang , A preconditioned fast finite difference method for space-time fractional partialdifferential equations , Fract. Calc. Appl. Anal., 20 (2017), pp. 88–116.14.
R. Hilfer , Applications of Fractional Calculus in Physics , World Scientific, Singapore, 2000.15.
Y. T. Jia, M. Q. Xu, and Y. Z. Lin , A numerical solution for variable order fractional functional differentialequation , Appl. Math. Lett., 64 (2017), pp. 125–130.16.
S. D. Jiang, J. W. Zhang, Q. Zhang, and Z. M. Zhang , Fast evaluation of the Caputo fractional derivativeand its applications to fractional diffusion equations , Commun. Comput. Phys., 21 (2017), pp. 650–678.17.
R. Ke, M. K. Ng, and H. W. Sun , A fast direct method for block triangular Toeplitz-like with tri-diagonal blocksystems from time-fractional partial differential equations , J. Comput. Phys., 303 (2015), pp. 203–211.18.
A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo , Theory and applications of fractional differentialequations , Elsevier, Amsterdam, 2006.19.
H. L. Liao, D. F. Li, and J. W. Zhang , Sharp error estimate of the nonuniform L formula for linear reaction-subdiffusion equations , SIAM J. Numer. Aanl., 56 (2018), pp. 1112–1133.20. F. W. Liu, V. Anh, and I. Turner , Numerical solution of the space fractional Fokker-Planck equation , J.Comput. Appl. Math., 166 (2004), pp. 209–219.21.
C. F. Lorenzo and T. T. Hartley , Variable order and distributed order fractional operators , Nonlinear Dyn.,29 (2002), pp. 57–98.22.
X. Lu, H. K. Pang, and H. W. Sun , Fast approximate inversion of a block triangular Toeplitz matrix withapplications to fractional sub-diffusion equations , J. Numer. Lin. Alg. Appl., 22 (2015), pp. 866–882.23.
X. Lu, H. K. Pang, H. W. Sun, and S. W. Vong , Approximation inversion method for time-fractionalsubdiffusion equations , Numer. Lin. Alg. Appl., 25 (2018), e2132.24.
C. Lubich and A. Sch¨adle , Fast convolution for nonreflecting boundary conditions , SIAM J. Sci. Comput., 24(2002), pp. 161–182.25.
F. Mainardi, M. Raberto, R. Gorenflo, and E. Scalas , Fractional calculus and continuous-time finance II:the waiting-time distribution , Phys. A, 287 (2000), pp. 468–481.26.
A. D. Obembe, M. E. Hossain, and S. A. Abu-Khamsin , Variable-order derivative time fractional diffusionmodel for heterogeneous porous media , J. Petrol. Sci. Eng., 152 (2017), pp. 391–405.27.
K. B. Oldham and J. Spanier , The fractional calculus , Academic Press, New York, 1974.28.
I. Podlubny , Fractional differential equations , Academic Press, New York, 1999.29.
M. Raberto, E. Scalas, and F. Mainardi , Waiting-times and returns in high-frequency financial data: anempirical study , Phys. A, 314 (2002), pp. 749–755.30.
L. E. S. Ramirez and C. F. M. Coimbra , On the selection and meaning of variable order operators for dynamicmodeling , Int. J. Differ. Equ., 2010 (2010), Article ID 846107, 16 pages.31.
S. G. Samko and B. Ross , Integration and differentiation to a variable fractional order , Integr. Transf. Spec.Funct., 1 (1993), pp. 277–300.32.
S. Serra , Superlinear PCG methods for symmetric Toeplitz systems , Math. Comput., 68 (1999), pp. 793–803.33.
I. M. Sokolov and J. Klafter , From diffusion to anomalous diffusion: a century after einsteins brownianmotion , Chaos, 15 (2005), pp. 1–7.34.
C. M. Soon, C. F. M. Coimbra, and M. H. Kobayashi , The variable viscoelasticity oscillator , Ann. Phys., 14(2005), pp. 378–389.35.
M. Stynes, E. O’riordan, and J. L. Gracia , Error analysis of a finite difference method on graded meshesfor a time-fractional diffusion equation , SIAM J. Numer. Aanl, 55 (2017), pp. 1057–1079.36.
H. G. Sun, A. Chang, Y. Zhang, and W. Chen , A review on variable-order fractional differential equations:mathematical foundations, physical models, numerical methods and applications , Fract. Calc. Appl. Anal., 22(2019), pp. 27–59.itle Suppressed Due to Excessive Length 2137.
H. G. Sun, W. Chen, and Y. Q. Chen , Variable-order fractional differential operators in anomalous diffusionmodeling , Phys. A, 388 (2009), pp. 4586–4592.38.
H. G. Sun, W. Chen, H. Wei, and Y. Q. Chen , A comparative study of constant-order and variable-orderfractional models in characterizing memory property of systems , Eur. Phys. J. Spec. Top., 193 (2011), pp.185–192.39.
Z. Z. Sun and G. H. Gao , Fractional differential equations , De Gruyter, Berlin, 2020.40.
F. H. Zeng, I. Turner, and K. Burrage , A stable fast time-stepping method for fractional integral andderivative operators , J. Sci. Comput., 77 (2018), pp. 283–307.41.
J. L. Zhang, Z. W. Fang, and H. W. Sun , Exponential-sum-approximation technique for variable-order time-fractional diffusion equations , J. Appl. Math. Comput., 2020. submitted.42.
X. Zhao, Z. Z. Sun, and G. E. Karniadakis , Second-order approximations for variable order fractionalderivatives: Algorithms and applications , J. Comput. Phys., 293 (2015), pp. 184–200.43.