Operational matrices of Chebyshev polynomials for solving singular Volterra integral equations

An effective technique based on fractional calculus in the sense of Riemann–Liouville has been developed for solving weakly singular Volterra integral equations of the first and second kinds. For this purpose, orthogonal Chebyshev polynomials are applied. Properties and some operational matrices of these polynomials are first presented and then the unknown functions of the integral equations are represented by these polynomials in the matrix form. These matrices are then used to reduce the singular integral equations to some linear algebraic system. For solving the obtained system, Galerkin method is utilized via Chebyshev polynomials as weighting functions. The method is computationally attractive, and the validity and accuracy of the presented method are demonstrated through illustrative examples. As shown in the numerical results, operational matrices, even for first kind integral equations, have relatively low condition numbers, and thus, the corresponding matrices are well posed. In addition, it is noteworthy that when the solution of equation is in power series form, the method evaluates the exact solution.


Introduction
The aim of this study is to present a high-order computational method for solving special cases of singular Volterra integral equations of the first and second kinds, namely Abel's integral equations, defined by and 0\a\1; 0 x T; where f ðxÞ 2 C½0; T is the known function and y(x) is the unknown function that to be determined, and T is a positive constant. Abel's equation is one of the integral equations derived directly from a concrete problem of mechanics or physics (without passing through a differential equation). Historically, Abel's problem is the first one that led to the study of integral equations. The generalized Abel's integral equations on a finite segment appeared in the paper of Zeilon [15] for the first time.
A comprehensive reference on Abel-type equations, including an extensive list of applications, can be found in [8,9].
The construction of high-order methods for the equations is, however, not an easy task because of the singularity in the weakly singular kernel. In fact, in this case, the solution y is generally not differentiable at the endpoints of the interval [3], and due to this, to the best of the authors' knowledge, the best convergence rate ever achieved remains only at polynomial order. For example, if we set uniform meshes with n þ 1 grid points and apply the spline method of order m, then the convergence rate is only Oðn À2P Þ at most (see [4,12]), and it cannot be improved by increasing m. One way of remedying this is to introduce graded meshes [13]. By doing so, the rate is improved to Oðn Àm Þ, which now depends on m, but still at polynomial order. Rashit Ishik [10] used Bernstein series solution for solving linear integro-differential equations with weakly singular kernels. In [5] and [6], wavelets method was applied for solution of nonlinear fractional integro-differential equations in a large interval and systems of nonlinear singular fractional Volterra integro-differential equations. Authors of [11] applied fractional calculus for solving Abel integral equations. The expansions approach for solving cauchy integral equation of the first kind is discussed in [14].
In this paper, we use the Chebyshev polynomials operational matrices via Galerkin method for solving weakly singular integral equations. Our method consists of reducing the given weakly singular integral equation to a set of algebraic system by expanding the unknown function by Chebyshev polynomials of the first kind. Galerkin method is utilized to solve the obtained system.
The structure of this paper is arranged as follows. The main problem and brief history of some presented methods are expressed in Sect. 1. In Sect. 2, we present some necessary definitions and mathematical preliminaries of the fractional calculus theory in the sense of Riemann-Liouville. Section 3 is devoted to introducing Chebyshev polynomials, properties and some operational matrices of these functions. In Sect. 4, Chebyshev polynomials are applied as testing and weighting functions of Galerkin method for efficient solution of Eq. 1. In Sect. 5, we report our numerical founds and compare with other methods in solving these integral equations, and Sect. 6 contains our conclusion.

Some preliminaries in fractional calculus
In this section, we briefly present some definitions and results in fractional calculus for our subsequent discussion. The fractional calculus is the name for the theory of integrals and derivatives of arbitrary order, which unifies and generalizes the notions of integer-order differentiation and n-fold integration [7]. There are various definitions of fractional integration and differentiation, such as Grunwald-Letnikov, and Caputo and Riemann-Liouville's definitions. In this study, fractional calculus in the sense of Riemann-Liouville is considered.
Definition 1 Let f be a real function on [a, b] and 0\a\1. Then, the left and right Riemann-Liouville fractional integral operators of order a for the function f are defined, respectively, as Definition 2 For f 2 C½a; b, the left and right Riemann-Liouville fractional derivatives are defined, respectively, as ðt À xÞ Àa f ðtÞdt: In this study, the left Riemann-Liouville fractional integral operator is utilized to transform singular integral equation to some algebraic system. Therefore, for abbreviation, the mentioned operator is denoted by I a .
Theorem 1 The operator I a (stand for left and right Riemann-Liouville fractional integral operator) satisfies the following properties: Proof We prove the proposition for left Riemann-Liouville fractional integral operator and the proof for right Riemann-Liouville fractional integral operator can be done similarly. For the part (1), we have In addition, for part (2), we have by changing the variable t ¼ rx, we get now, using the formulae of the beta function, we have On the other hand, we know that the beta function can be written in terms of the Gamma function as follows: Bða; bÞ ¼ CðaÞCðbÞ Cða þ bÞ ; so we have

Chebyshev polynomials
In this section, a brief summary of orthogonal Chebyshev polynomials is expressed.

Definition 3
The nth degree of Chebyshev polynomials is defined by where t ¼ cosðhÞ. The roots of Chebyshev polynomial of degree n þ 1 can be obtained by In addition, the following successive relation holds for Chebyshev polynomials: Chebyshev polynomials are orthogonal with respect to the weight function wðxÞ

Function approximation
A function yðxÞ 2 L 2 ½0; T, can be expressed in terms of the shifted Chebyshev polynomials as [2] yðxÞ ¼ where u j is the shifted Chebyshev polynomial of degree j. The coefficients c j are given by and . . .:

Operational matrix of fractional Integration
The fractional integral of Chebyshev polynomials can be defined by where and Ã is point by point product and A is ðn þ 1Þ Â ðn þ 1Þ lower triangular matrix defined by In addition, for each function g(x), approximated by shifted Chebyshev functions (gðxÞ ¼ D T WX), the fractional integral can be written as Numerical implementation In this section, the shifted Chebyshev polynomials are applied for solving singular integral Eqs. (1) and (2) Now, the unknown function y(x) is approximated by shifted Chebyshev polynomials as yðxÞ ' y n ðxÞ ¼ C T WX; where C T is the unknown vector, that to be determined. In the following, we describe the method in detail for the first and second kinds.

The first kind
Consider the weakly singular Volterra integral equation of the first kind (1), according to Eq. (10), we have substituting Eqs. (11) and (8) in (12), we can get Therefore, singular integral equation (1) is transformed to the above algebraic system. For solving this system, Galerkin procedure is utilized via shifted Chebyshev polynomials. Put CðbÞðA Ã WÞ ¼ K and suppose u i ðxÞ be the shifted Chebyshev polynomial of degree i, which can be written as where W i is the ith row of the matrix W. Multiplying Eq. (13) by u i ðxÞ, we get Putting X ¼ 1 x x 2 ::: x n x x 2 x 3 ::: . . .
x n x nþ1 x nþ2 ::: and by integrating current equations from 0 to T, we have where P x and P b are related integration operational matrix, defined by i ¼ 0; . . .; n; j ¼ 0; . . .; n: Considering Eq. (17) for i ¼ 0; . . .; n, we get a system of n þ 1 equations and n þ 1 unknowns that can be solved easily.

The second kind
Similarly, for the second kind singular integral equation (2), we have so Eq. (19) can be written as Multiplying Eq. (20) by u i ðxÞ and integrating from 0 to T, we can rewrite the current equation in the following form: Eq. (21) is a system of n þ 1 equations and n þ 1 unknowns that can be solved easily.

Illustrative examples
In with the exact solution yðxÞ ¼ x 2 . The unknown coefficients for c i are obtained through the method explained in Sect. 4 for N ¼ 4.
The function y(x) is a polynomial of degree 2 and the least approximation level for Chebyshev polynomials, in this study, is N ¼ 4. Therefore, the approximated solution through the presented method is the same as the exact solution, that is y 4 ðxÞ ¼ x 2 .
Example 2 Consider the following singular integral equation: with the exact solution yðxÞ ¼ x In Table 1, we have presented exact and approximated solutions of Example 2 in some arbitrary points. In addition, the last line of Table 1 shows the condition number of operational matrices. The errors of approximate solutions in the levels N ¼ 4; 8, and 12 are shown in Fig. 1.
Example 3 Consider the following singular integral equation of the first kind:  In Table 2, we present exact and approximate solutions of Example 3 in some arbitrary points and compare them by the results of [1] for n ¼ 20. In addition, the last line of Table 2 shows the condition number of operational matrices. Figure 2 illustrates the error of approximate solutions in the levels N ¼ 4; 8, and 12.     and the approximate solution for N ¼ 4 is calculated as y 4 ðxÞ ¼ À1:5114x 4 þ 3:7014x 3 À 3:4978x 2 þ 2:5439x þ 0:0502: In Table 3, exact and approximated solutions of Example 4 in some arbitrary points are given. The condition number of operational matrices for N ¼ 4; 8, and 12 are calculated by Eq. 22, and is shown in the last row of Table 3. Figure 3 is helpful in geometric understanding the errors of approximated solutions in the levels N ¼ 4; 8, and 12.

Conclusions
In this study, a numerical approach based on Chebyshev polynomials operational matrices was developed to approximate the solution of the weakly singular Volterra integral equations of the first and second kinds. Applying fractional derivative of these polynomials, we have transformed the singular integral equations to some linear algebraic system. The numerical results obtained support the validity and efficiency of the proposed method. It is noteworthy that when the solution of equation is in power series form, the method evaluates the exact solution, such as Example 1. In addition, as can be seen, the operational matrices of first kind integral equations have relatively low condition numbers. Thus, the corresponding matrices are well posed.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.