Analysis of a High-Accuracy Numerical Method for Time-Fractional Integro-Differential Equations

: A high-order finite difference numerical scheme based on the compact difference operator is proposed in this paper for time-fractional partial integro-differential equations with a weakly singular kernel, where the time-fractional derivative term is defined in the Riemann-Liouville sense. Here, the stability and convergence of the constructed compact finite difference scheme are proved in L ∞ norm, with the accuracy order O ( τ 2 + h 4 ) , where τ and h are temporal and spatial step sizes, respectively. The advantage of this numerical scheme is that arbitrary parameters can be applied to achieve the desired accuracy. Some numerical examples are presented to support the theoretical analysis

At present, Riemann-Liouville and Caputo derivatives are commonly used in engineering and the sciences. The properties of the Riemann-Liouville derivative are different from those of the Caputo derivative. Furthermore, the Riemann-Liouville derivative is singular at zero, and

Numerical Scheme
Let M and N be two positive integers, and let h = L/M and τ = T/N be the spatial step size and time step size, respectively. C is a constant, which may be different in different locations.
For j = 0, 1, . . . , M and n = 0, 1, . . . , N, the mesh point (x j , t n ) is defined as x j = jh and t n = nτ. Let u n j be the exact solution and U n j be the approximate solution at each mesh point (x j , t n ) of Equation (1). The following notations and lemmas will be used throughout this paper: Using the shifted and weighted Grünwald difference method to approximation derivatives, we can obtain a discrete scheme with a second-order convergence rate in the temporal direction.
Lemma 1 ([31]). Suppose that u ∈ L 1+α (R), and let The shifted and weighted Grünwald difference operator is defined as follows: where p is an integer according to Equation (4). Then, we obtain that uniformly for t ∈ R as τ → 0. The coefficients g (α) i (0 < α ≤ 1) are defined as follows: In addition, if α = 0, we stipulate g (α) t u and its Fourier transform belong to L 2+α (R), and define the shifted and weighted Grünwald difference operator as follows: Then, uniformly for t ∈ R as τ → 0, we obtain where p and q are integers.
Furthermore, we define L D α τ,p,q = L D α τ,q,p . A finite difference approximations to discrete derivative in time is as follows (see [31]): Combining the above equality, we obtain Lemma 4 ( [28,36]). Suppose that L(t) ∈ C 2 [0, T]; then, there exists a positive constant C which depends only on Lemma 5. For any κ n n−i,n (1 ≤ i ≤ n − 1, 2 ≤ n ≤ N) satisfying the definition of Equation (7), the sequence {κ n n−i,n } decreases monotonically depending on k.

Lemma 6 ([36]
). Let κ n i,n (1 ≤ n ≤ N) be defined as Equation (7), and for all β (0 < β < 1), we obtain (i) 0 < κ n 0,n < τ β β+1 , Assume that u(x, t) ∈ C 6,2 x,t ([0, L] × [0, T]) and consider Equation (1) on grid point (x j , t n ) and apply compact difference operator H to both sides. We thus have For the left term of Equation (8), by Lemma 2, we get For the first term on the right of Equation (8), Lemma 3 and Lemma 4 imply that Substituting Equation (9) and Equation (10) into Equation (1), it follows that where f n j = f (x j , t n ), |R n j | ≤ C(τ 2 + h 4 ). Neglecting the small term R n j in Equation (11), when 1 ≤ j ≤ M − 1, 1 ≤ n ≤ N, the compact finite difference scheme for Equation (1) is given as follows: At each time level, the compact difference scheme Equation (12) is a system of linear algebraic equations with a strictly diagonally dominant matrix as its coefficient matrix. We can obtain the following theorem.

Stability and Convergence
In this section, we first give some notations and lemmas which will be used in the subsequent discussions. Then, the stability analysis and error estimates of the compact finite difference scheme for Equation (12) are obtained. Denote the inner product and norm are denoted as follows: i } ∞ n=0 be defined as Equation (6); then, for each positive integer m and for any (p 0 , p 1 , · · · , p m ) T ∈ R m+1 , we have Lemma 11 ([36,40]). Define {ζ n } as a sequence of non-negative real numbers if it satisfies the following inequality where b n is a nondecreasing sequence of non-negative numbers, and k ≥ 0. It thus holds that Lemma 12 ([41]). (Čebyšev inequality) If η = (η 1 , η 2 , . . . , η n ) is a nonincreasing sequence and µ = (µ 1 , µ 2 , . . . , µ n ) is a nondecreasing sequence, then we obtain Proof. By definition of (u, v) A , Equation (3) and Lemma 8, we have from which the desired result is obtained.

Lemma 14.
Let κ n n−i,n and δ x U n−i be defined in Equation (7) and Equation (2), respectively; then, By Lemma 5 and Lemma 12, we get For the left term of Equation (13), by Lemma 6, we have The proof is completed. (12) with the given initial and boundary conditions in the sense that for all τ > 0; if τ 1−β N ≤ C 1 and T ≥ 1, then

Theorem 2. If U n is the approximation solution of Equation
where C and C 1 are positive constants, and they depend on T.
Proof. By Equation (12), we obtain Taking the inner product of Equation (14) with τHU n , then By Equation (7), Equation (15) and Lemma 9, we have By inequality (16) and Lemma 13, we get that By inequality (17), we obtain By inequality (18) and summing the above expression from n = 1 to m and 1 ≤ m ≤ N, we get By Lemma 6, Lemma 7, Lemma 14 and inequality (19), we obtain By inequality (20) and T ≥ 1, we get f n . We can thus obtain that b n is a positive nondecreasing sequence with respect to n. Hence, by the above inequalities and Lemma 11, we have Let C = max{ 3Tβ 2 , ), the following inequality holds which is the desired result.
Define e n j = u n j − U n j (0 ≤ j ≤ M, 0 ≤ n ≤ N) as errors at each mesh point (x j , t n ); then, the error bound of our numerical scheme will be considered as follows. (To make the following expression succinct, the subscript j will be neglected.) Theorem 3. Let u n ∈ C 6,2 x,t ([0, L] × [0, T]) be the exact solution of Equation (1) and U n be the numerical solution of Equation (12), if τ 1−β N ≤ C 1 and T ≥ 1; then, the error bound is as follows: where C is a positive constant, and it depends on T.
Proof. By definition of e n , then the error equation is as follows Taking the inner product of inequality (21) with τHe n , we obtain By Equation (7), Lemma 9 and inequality (22), we have κ n n−i,n (δ 2 x e n−i , He n ) + τ 1−β (R n , He n ).
By inequality (18), inequality (23) and summing the above expression from n = 1 to m and 1 ≤ m ≤ N, we obtain From Theorem 2, we get By inequality (25) and Lemma 10, then where C, C 1 and C 2 are the constants, and the desired result is obtained.

Numerical Experiments
We will present several experiments which support the theoretical analysis in Section 3. All numerical tests were performed on an AMD Ryzen 7 4700U with Radeon Graphics (2.00 GHz) and 16 Gb of RAM, using MATLAB (R2020b).
Let N and M be two constants, and define τ = T/N and h = L/M to be temporal step size and spatial step size, respectively. Let u(x j , t n ) be the exact solution of Equation (1) and U n j be numerical solution of Equation (12). As in [14,25], we consider the L ∞ -norm errors and corresponding convergence rates in the following: where E 1 and E 2 are errors correspond to grids with mesh sizes h 1 and h 2 , respectively. In addition, we have Rate t = log(E 1 /E 2 )/log(τ 1 /τ 2 ), where E 1 and E 2 are errors correspond to grids with mesh sizes τ 1 and τ 2 , respectively. If we perform an α-order Riemann-Liouville fractional integral on both sides of Equation (1), we can obtain new (α + β)-order partial integro-differential equations. In [29], for h 1 = h 2 = h, the scheme is convergent with the order O(τ 1+α + h 2 ) when u tt (x, t) is singular at t = 0 and O(τ 2 + h 2 ) when u tt (x, t) is smooth at t = 0. According to the results of the following three examples, we find that our scheme is convergent with the order O(τ 2 + h 4 ). Example 1. In the first example, we choose L = 2, T = 4 and the force term as follows: where the corresponding initial term is ϕ(x) = sin(πx) and the exact solution is u(x, t) = −4t 4 sin(2πx) + sin(πx). This experiment shows that convergence order in the temporal direction can achieve the global second order. Tables 4-6 show the convergence rates in space and the corresponding L ∞ norm errors for N = 1000 and different values of M. For different β values, this experiment shows the errors in the case of α = 0.15, 0.5 and 0.75. The experiment demonstrates that convergence order in the spatial direction can be of the fourth order. One can see that they are in good agreement with the theoretical results. Table 1. L ∞ norm errors, convergence orders in temporal direction with α = 0.15 (M = 100).  Table 2. L ∞ norm errors, convergence orders in temporal direction with α = 0.5 (M = 100).  Table 4. L ∞ norm errors, convergence orders in spatial direction with α = 0.15 (N = 1000). Example 2. In this example, we choose L = T = 1 and the force term as follows: For this test, the reference solution is under a very fine mesh (N = M = 2000). Tables 7-10 show the numerical results of Example 2. Tables 7 and 8 show the time convergence rates and corresponding numerical errors for M = 100 and different values of N. Table 7 shows the results of α = 0.2 at β = 0.2, 0.5 and 0.8. Table 8 shows the results of α = 0.5 at β = 0.2, 0.5 and 0.8. Table 9 shows the space convergence rates and corresponding numerical errors for N = 1000 and different values of M, where α = 0.5 and β = 0.2, 0.5 and 0.8. Table 10 shows the space convergence rates and corresponding numerical errors for N = 1000 and different values of M, where α = 0.8 and β = 0.2, 0.5 and 0.8.
The numerical results of this example show that the convergence order in the temporal and spatial directions can reach the second and fourth orders, respectively. The results of Example 2 show that the convergence orders match the theoretical ones.  Example 3. For the last test, we choose L = T = 1. The force term is the following: where the corresponding initial term is ϕ(x) = 0 and the exact solution is u(x, t) = (t − t 4 )sin(5πx).
Let M = 100 and N = 2000. In Figure 1, we display the exact solution and numerical solution and the corresponding absolute error and contour plot absolute error with α = 0.1 and β = 0.1. Similarly, the exact solution, numerical solution and corresponding error are presented in Figure 2 with α = 0.9 and β = 0.9. One can obviously see that our method can achieve the desired accuracy.

Conclusions
The main result of this work is that an efficient finite difference numerical scheme is proposed and analyzed for time-fractional integro-differential equations, in which the derivative is defined as a Riemann-Liouville derivative. The stability analysis and error bound of the presented numerical method are carried out, and the convergence order is O(τ 2 + h 4 ). Several numerical experiments are given to support the theoretical analysis. Furthermore, the accuracy of the discrete scheme constructed in this paper is not affected by the parameters α and β. Development of fast and parallel-in-time methods [42,43] for accelerating the numerical schemes of time-fractional PDEs can be carried out in our future work. At the same time, we will try to study the high-order numerical schemes of fractional integro-difference equations to solve nonlinear and high-dimensional problems.