THE JACOBI SPECTRAL COLLOCATION METHOD FOR FRACTIONAL INTEGRO-DIFFERENTIAL EQUATIONS WITH NON-SMOOTH SOLUTIONS

In recent years, many numerical methods have been extended to fractional integro-differential equations. But most of them ignore an important problem. Even if the input function is smooth, the solutions of these equations would exhibit some weak singularity, which leads to non-smooth solutions, and a deteriorate order of convergence. To overcome this problem, we first study in detail the singularity of the fractional integro-differential equation, and then eliminate the singularity by introducing some smoothing transformation. We can maximize the convergence rate by adjusting the parameters in the auxiliary transformation. We use the Jacobi spectral-collocation method with global and high precision characteristics to solve the transformed equation. A comprehensive and rigorous error estimation under the L∞and L2 ωα,β -norms is derived. Finally, we give specific numerical examples to show the accuracy of the theoretical estimation and the feasibility and effectiveness of the proposed method.


1.
Introduction. With more attention to integro-differential equations, the numerical calculation of fractional integro-differential equations has been also paid more attention by many scholars. Fractional integro-differential equations have a profound physical background and a rich connotation. For various phenomena in damping laws, diffusion process [31], earthquake model [13], fluid-dynamic traffic model [14], mathematical physics and engineering [28,37], chemistry, acoustics, fluid and continuum mechanics [8], psychology [1,35] and other fields, fractional integrodifferential equations are suitable models.
In recent years, several numerical methods were proposed to solve the fractional integro-differential equation. The most common methods are the Adomian decomposition method [25], collocation method [33] and fractional differential transform method [2]. In [16], Huang proposed a method for solving linear fractional integrodifferential equations by Taylor's expansion, including Fredholm type and Volterra type. Yang et al. [34] used the Laplace decomposition method to solve the fractional integro-differential equation. M. Jani et al. [17] proposed a numerical method for solving fractional integro-differential equations with nonlocal boundary conditions by using Bernstein polynomials.
The fractional differential operators are nonlocal and have weakly singular kernels. The fractional differential equations are more complicated than the integerorder counterparts. In recent years, many numerical methods have been extended to fractional integro-differential equations. Most of the analyses have some unreasonable limitations on solutions in order to achieve high accuracy. When these equations are transformed into the equivalent Volterra integral equations of the second kind with a weakly singular kernel, even if the input function is smooth, the solution of the equation usually exhibits a weak singularity at z = 0. Which leads to a non-smooth solution and a lower order of convergence.
So far, in order to solve the fractional differential equations and fractional integrodifferential equations with non-smooth solutions, several methods have been proposed. One method is to approximate the fractional derivative operators in the governing differential equation directly and then the corresponding collocation schemes are derived [3,18]. Another method is to rewrite the governing differential equation in an equivalent integral equation, solved by the corresponding collocation method [22,36,21]. The integral collocation method is more stable than the differential collocation method. The reason is that numerical differentiation is sensitive to small perturbations in the input. But numerical integration is essentially stable [10]. Therefore, when using the differential collocation method to solve differential equations, we need to employ efficient integration preprocessing to overcome the illconditioning problem. It is necessary with increasing of the number of collocation nodes [15]. We would like to note that Hao proposed an efficient finite difference algorithm to solve fractional boundary value problems with non-smooth solutions in [12].
For many types of equations with non-smooth solutions, the idea of introducing suitable transformations has been considered. It would eliminate the singularity in the transformed equation, and lead to a high convergence order. For the second kind of Volterra integral equation, Chen and Tang [9] proposed a variable transformation to eliminate the singularity of the solution. With a strict error analysis, the method is shown to have a spectral convergence. Pedas [19] made a proper transformation, and used the piecewise polynomial collocation method to solve the resulting equation on a mildly graded grid or a uniform grid. Baratella and Orsi [7] used a variable transformation to turn the solution of the linear Volterra integral equation of the second kind smooth, and solved it by the standard product integration method. Tang [20] used a variable transformation and the Jacobi spectral collocation method to solve the Abel-Volterra integral equation of the second kind. For the linear Fredholm integral equation of the second kind, Monegato and Scuderi [26] proposed a non-linear transformation to eliminate the singularity of the equation. Ghoreishi [11] used a variable transformation and a spectral method to solve the multi-order fractional differential equation. Pedas et al. [27] regularized the solution of fractional initial and boundary value problems by a suitable smoothing transformation. They solved the transformed equation by a piecewise polynomial collocation method on a mildly graded grid and on a uniform grid. Zaky [38] used a smoothing transformation and the Jacobi spectral collocation method to solve the rational-order fractional terminal value problems with non-smooth solutions.
Based on the above works, we apply the smoothing transformation and the Jacobi spectral collocation method with high accuracy and global characteristics to the following fractional integro-differential equations.
whereg(z) is the source function, andK(z, τ ) is the kernel function. The given functiong(z) andK(z, τ ) are continuous on their respective domains 0 ≤ τ ≤ z ≤ T and ∆ T : Let Γ(·) denote the Gamma function. For any positive integer n with n − 1 < α < n, the Caputo derivative is defined as follows: In addition, the Riemann-Liouville fractional integral I α z of order α is defined by We note that, The layout of this paper as follows: In Section 2, we introduce the basic properties of Jacobi polynomials and Jacobi-Gauss interpolation. The fractional integral differential equation is transformed into an equivalent integral equation, in Section 3. A smoothing transformation of variable is defined for the new equation so that the solution is smooth, in Section 3. The Jacobi spectral-collocation method is defined in Section 4. The convergence analysis of the collocation method is derived in Section 5. In Section 6, we give several numerical examples verifying the accuracy of the theoretical estimation and the feasibility and effectiveness of the method. Finally, some concluding remarks are drawn in Section 7.
2. Some properties of Jacobi polynomials. In this section, we introduce some basic properties about Jacobi polynomials and Jacobi-Gauss interpolation that are related to spectral-collocation methods [30].
The Jacobi polynomials, denoted by P α,β n (x), are orthogonal with the Jacobi weight function ω α,β (x) = (1 − x) α (1 + x) β over Λ = (−1, 1), namely, where δ i,j is the Kronecker function and For a given positive integer N ≥ 0, let P N denote the space of all polynomials of degree not exceeding N. We denote by {x α,β i , α,β i } N i=0 the set of quadrature nodes and weights of the Jacobi-Gauss integration. The Jacobi-Gauss integration formula has the form The above quadrature formula (5) is exact for any ϕ(x) ∈ P 2N +1 . Hence, by (3), For any µ ∈ C(Λ), the Jacobi-Gauss interpolation operator I α,β x,N : C(Λ) → P N is determined uniquely by The interpolation condition (7) implies that I α,β x,N µ = µ for all µ ∈ P N . On the other hand, since I α,β x,N µ ∈ P N , we can write In particular, for β = 0, the set of Jacobi polynomials is reduced to J α i (x). Therefore, we can also write x α j = x α,0 j , α j = α,0 j and I α x,N = I α,0 x,N .
3. Setting the problem. In this section, we use the definition and related properties of Riemann-Liouville fractional integral and Caputo fractional derivative to transform the original fractional equation with initial conditions into the second kind of Volterra integral equation with weak singular kernel. We show the equation has a non-smooth solution. Then we apply the smoothing transformation to eliminate the singularity of the solution at the left endpoint. First, we use (2) to transform the original equation (1) into an equivalent Volterra integral equation with weak singular kernel.
Using Dirichlet's formula The equation (9) is transformed into the following Volterra integral equation of the second kind by the linear transformation w = τ −s z−s , w ∈ [0, 1], where Proof. We have For z ∈ [0, T ] and s ∈ [0, z], z − s ∈ [−z, T ]. It is known thatK(w(z − s) + s, s) is bounded and continuous on ∆ T . LetK(w(z − s) + s, s) have the maximum and minimum values on ∆ T , Q max and Q min , respectively. There is 1 We give some lemmas on smoothness of solution of the general Volterra integral equation of the second kind.
where γ > −1, f : I → R n is a continuous bounded function, K : S × R n → R n is a continuous bounded function, and {S = (z, s) : 0 ≤ s ≤ z ≤ T }.
1. If f (z) and K(z, s, y(s)) are differentiable, the integral equation has a unique solution y(z) that is also differentiable on (0, T ]; 2. If f (z) = F (z, z 1+γ ), F (z 1 , z 2 ) and K(z, s, y(s)) are differentiable, the integral equation has a unique solution y(z) that satisfies We note that the solution of equation (14) is not smooth at z = 0 in general. It can be obtained from the above lemma that the Eq.(10) has a unique solution y(z), which is differentiable on z ∈ (0, T ] and is not necessarily smooth at z = 0. For 0 < 1 − α < 1, the equation has a singular term (z − s) α−1 . By the literature [5], for any positive integer m, if K(τ (w(z − s)), s) and f (z) are continuous differentiable functions of order m in the corresponding area, there exists a function Y = Y (z, v) possessing continuous derivatives of order m, such that the solution of the Eq.(10) can be written as y(z) = Y (z, z α ). This indicates that when z → 0, y (m) (z) ∼ z α−m , and thus y(z) / ∈ C m [0, T ].
1. If m = 0, the Volterra integral equation, possesses a unique solution y(z) ∈ C(I). This solution has the representation where the resolvent kernel R µ (z, s) of the kernel (z − s) −µ K(z, s) has the form Here, Q µ (z, s) is continuous on D. 2. If m ≥ 1, every nontrivial solution has the property that y(z) / ∈ C 1 (I): as z → 0 + the solution behaves like In earlier work, when 0 < 1 − α < 1, the general form of the exact solution of the Eq.(10) has been derived, in next lemma.
The Volterra integral equation (10) has a unique solution y Further, y(z) has the following form where ψ k ∈ C m (I), k ≥ 1, and the series is absolutely uniformly convergent on I.
If µ = p q is rational (i.e., p, q ∈ N , reduced to lowest terms), then the solution of the Eq.(10) can be expressed in the form where ν s ∈ C m (I)(0 ≤ s ≤ q − 1).
From the above lemma, y(z) / ∈ C m [0, T ]. For the Eq.(13), Chen and Tang [9] proposed the function transformationỹ(t) = t u+m−1 [y(t) − y(0)] to remove a single term singularity like t 1−u−m . They used the spectral-collocation method and achieved excellent results. Our work has been inspired by their excellent results in [9,32]. We apply the similar transformation mentioned above to (14), and obtain the following lemma. This result will be the starting point to the construction of the numerical method presented later. Lemma 3.6. Using the following transformation for (14), it is deduced that is at least a first-order differentiable function. Proof. Taking the above transformation (15) to (14), we obtain Thus we deduce the two conclusions in the lemma.
By the Lemma 3.5, it follows that the solution of the Eq.(10) can be written in the form ofỹ(z) = y 1 (z) + y 2 (z) where, for a fixed m, y 1 (z) ∈ C m (I) and y 2 (z) is the non-smooth part of the solution.
Our first step is to replaceỹ(z) byỹ(z) := y(z) + λ , where y(0) = 0. Hence, the Eq.(1) can be expressed in terms of y as An equivalent integral form of the above equation is We apply the smoothing transformation reducing the problem (17) to the following integral equation whose solution does not involve anymore singularities in the first derivative.
Equation (20) can be written as Using the formula becomes where Finally, by the change of variable and setting is reduced to where 4. The Jacobi spectral collocation method. In this section, we propose the Jacobi spectral-collocation method to (23). Solving Equation (23) by the Jacobi spectral-collocation method is to find Y N (x) ∈ P N , such that where In order to implement the above basic algorithm more effectively, we set Employing (25) and (3), we can directly calculate the result that Using (26), (5) and (8) it yields wherê Similarly, by Eqs. (3) and (25), we conclude that Using (28), (5) and (8), we obtain wherê In summary, by (25)-(29), we deduce that Finally, using (3) it yields The numerical solution can be obtained by solving the equations.
5. Error analysis. In this section, we estimate the error of the numerical solution.
We bound the error in the L ∞ and L 2 ω α,β norms. In order to give the subsequent lemmas conveniently, we first introduce some spaces.
Let the region Λ ⊂ R n be a non-empty Lebesgue measurable set, u(x) is a real value Lebesgue measurable function defined on Λ. L p (Λ) is defined as For a non-negative integer l, and 1 ≤ p ≤ ∞, the Sobolev space is defined as: If p = 2, we record H l,2 (Λ) as H l (Λ), which is a separable Hilbert Space. If l = 0, Then we introduce the weighted L 2 ω α,β (Λ) space. Assume that the weight function u is measurable and u ω α,β < ∞}, endowed with the norm and inner product The weighted Hilbert space is defined as follows , 0 ≤ m ≤ l}, equipped with the norm, semi-norm and inner product For a non-negative integer l, we introduce the non-uniformly Jacobi-weighted Sobolev space B l ω α,β (Λ) := {u : ∂ m x u ∈ L 2 ω α+m,β+m (Λ), 0 ≤ m ≤ l}, endowed with the norm, semi-norm and inner product u B l ω α,β = (u, u) where u ω α,β is the norm of L 2 ω α,β (Λ). Especially, L 2 (Λ) = B 0 ω 0,0 , · = · L 2 (Λ) and · ∞ = · L ∞ (Λ) . The non-uniform Jacobi-weighted Sobolev space distinguishes itself from the usual weighted Sobolev space H l ω α,β by involving different weight functions for derivatives of different orders. It is clear that is the Banach space of the measurable functions u that are bounded outside a set of measure zero, equipped the norm We denote by C m (Λ) the space of m-times continuously differentiable functions on the interval Λ.  Lemma 5.2. [29] For α, β > −1, and any u ∈ B l ω α,β with l ≥ 1, and integers 0 ≤ m ≤ l ≤ N + 1, Moreover, for any u ∈ H l ω −1/2,−1/2 with 1 ≤ l ≤ N + 1, where c is a positive constant independent of l, N and u. ). The mapped Jacobi-Guass interpolation operator x I α−1 ξ,N : , From the above formula, we can derive the following results Moreover, if we denote I as the identity operator, for any 1 ≤ l ≤ N + 1, we have that Now, we are ready to prove the following convergence results.

5.1.
Error analysis in L 2 ω α−1,0 (Λ). Theorem 5.4. Let Y be the exact solution of the original equation (16). Let Y N be the numerical solution obtained by the discrete scheme (25) combined with the approximate Eq. (24). where Here Proof. It follows from (22) that By (34), we obtain The theorem is proved.

Applying Lemma 5.3 leads to
By (41), Similarly, using the Cauchy-Schwarz inequality, we further get From the above formula, we can get Also using Lemma 5.3 and (41), we obtain Drawn by the above formula, In summary,
By Lemma 5.1 and 5.3, We further get from the Cauchy-Schwarz inequality that It follows from Lemma 5.3 that In summary, we get The theorem is proved.
6. Numerical tests. In this section, we give numerical examples to illustrate the feasibility and efficiency of the method.
6.1. Example 1. We consider the following fractional integro-differential equation: The exact solution is t 2 3 . In the table 1, we list L ∞ -and L 2 ω α−1,0 -error with N = 10 and σ takes values of 1-9. As can be seen from the table, for σ takes values of 3, 6 and 9, which are multiples of 3, the error decreases faster and the convergence is better. Through a lot of numerical experiments, we conclude that not for all σ taking multiple value of 3, the convergence is better. One cannot predict that the greater the value of σ is, the better the convergence is. There is a critical point, which is the value of 27. For σ takes values of 2 to 27, the error convergence is better than that without smoothing transformation. To be more intuitive, we give a comparison of the exact solution and numerical solution with N = 10 and σ = 1, 3, 6, and 9 in Figures 1  and 2. And in Figure 3, we give L ∞ -and L 2 ω α−1,0 -error between the exact solution and the numerical solution with N = 10 and σ =1 to 9 and 28. Obviously, the error deceases faster when the value of 3, 6 and 9, which is the multiple of 3. Moreover, at the value of 28, the error convergence effect does not work well without smoothing transformation. Table 1. The L ∞ -and L 2 ω α−1,0 -error for N = 10 and σ takes values of 1-9. where The exact solution is sin √ t. We list the L ∞ -and L 2 ω α−1,0 -error with N = 10 and σ takes values of 1-6 in Table 2. From the table, we can also see that when σ takes values of 2 and 4, which is the multiple of 2, the error decreases faster. And through a large number of numerical experiments, we get a critical point of σ, which is the value of 14. When σ takes values of 2-13, the error convergence is better than that without smoothing transformation. For the sake of intuition, we give a comparison when N = 10 and σ takes values of 1, 2 and 4 in Figure 4. In Figures 5, the L ∞ -and L 2 ω α−1,0 -error between the exact solution and the numerical solution are given when N = 10 and σ takes values 1-6 and 14. It clear that when σ takes values of 2 and 4, which are multiples of 2, the error is smaller. While σ takes the value of 14, the error is worse than that without smoothing transformation. Table 2. The L ∞ -and L 2 ω α−1,0 -error for N = 10 and σ takes values of 1-6. where .
The exact solution is y(t) = t √ t. In Table 3, We list the L ∞ -and L 2 ω α−1,0 -error when N = 10 and σ takes values of 1-6. From the table, when σ takes the value of 4, which is the multiple of 4, the error decreases faster and the convergence is better. Through a lot of numerical experiments, we derive a critical point of the σ, which is the value of 8. When σ takes values of 2-7, the error is better than the situation that the smoothing transformation is not performed. We give a comparison when N = 10 and σ takes values of 1, 2 and 4 in Figures 6. In Figure 7, the L ∞ -and L 2 ω α−1,0 -error are given when N = 10 and σ takes values of 1-6 and 8. Obviously, when σ takes multiples of 4, the error decreases faster. However, when σ takes the value of 8, the error is not good without smoothing transformation. Table 3. The L ∞ -and L 2 ω α−1,0 -error for N = 10 and σ takes values of 1-6.    7. Conclusion. In this paper, we first study in detail the reason why the fractional integro-differential equation has non-smooth solution. We eliminate the singularity of the solution by introducing a smooth transformation. Then we use the Jacobi spectral-collocation method with global and high precision characteristics to solve the transformed equation. Particularly, we have proved that the convergence rate for non-smooth solutions can be enhanced by using a suitable smoothing transformation, which allows us to adjust a parameter in the solution in view of a priori known regularity of the given data. The proposed scheme has many advantages, including (i) ease of implementation, (ii) lower computational cost, and (iii) exponential accuracy. In addition, we give a theoretical proof of the convergence of collocation method, in both L ∞ -norm and L 2 ω α,β -norm. Finally, we give some specific numerical examples. The numerical results confirm the validity of scheme and the correctness of the conclusions for solving the fractional integro-differential equation. This indicates that the proposed scheme possesses a good prospect in solving fractional integro-differential equations with non-smooth solutions. Next, we will apply the methods in this paper to solve the nonlinear fractional integro-differential equations.