A New Numerical Technique for Solving Volterra Integral Equations Using Chebyshev Spectral Method

In this work, we propose a new method for solving Volterra integral equations. &e technique is based on the Chebyshev spectral collocation method. &e application of the proposed method leads Volterra integral equation to a system of algebraic equations that are easy to solve. Some examples are presented and compared with some methods in the literature to illustrate the ability of this technique. &e results demonstrate that the new method is more efficient, convergent, and accurate to the exact solution.


Introduction
Many mathematical formulations of physical phenomena contain integral equations; these equations arise in many fields including biology, economics engineering, and medicine; for more details, see [1][2][3][4]. Finding solutions of linear or nonlinear Volterra integral equations is highly interesting for researchers and scientists in this field and there are available studies to find analytical or numerical solutions for Volterra integral equations (see [5][6][7][8]). In fact, Volterra integral equations are usually difficult to solve analytically, so we resort to finding approximate solutions to the problems using numerical or analytical approximation methods. e main purpose of this work is to introduce a new method for solving Volterra integral equations that uses Chebyshev spectral collocation method. Chebyshev spectral collocation methods have been applied successfully in different areas of sciences and engineering because of their ability to give highly accurate solutions of single or system of boundary value problems [9][10][11][12][13]. Chebyshev spectral methods are defined everywhere in the computational domain. erefore, it is easy to compute high accurate values of a considered function at any point of the domain.
In the proposed method, we introduce a new integral transformation that is the replacement of an integral term in the Volterra integral equation by Chebyshev spectral differential matrix of known elements. e advantages of this approach over the standard same approaches are as follows: (i) the current technique suggests a standard way of choosing the auxiliary linear operator of the integral equation which is the main motivation behind the selection of the current algorithm whereas the other methods are choosing a linear operator to be simple in order to ensure that the integral or integrodifferential equations can be easily integrated and (ii) this algorithm transforms the integral equation(s) into a system of linear algebraic equations which is easier and faster to solve when compared to a system of integral equations. One of the disadvantages of this method is that it cannot be applied directly to the nonlinear integral equations; therefore, we need first to linearized the nonlinear integral equation using any available method such as the successive linearization method [14][15][16][17]. e applicability, accuracy, and reliability of the proposed method are confirmed by solving the various Volterra integral equations.
is paper is organized as follows: in Section 2, we introduce a description of the proposed method. In Section 3, the proposed method is implemented in some cases. e numerical results are discussed and investigated in Section 4. Finally, the paper is concluded in Section 5.

Description of the Method
Let us consider the Volterra integral equation of the second kind defined as with where u(x) is an unknown function and the kernel K(x, t) and the source term function f(x) are given functions.
In the proposed method, we are seeking for for solving the linear Volterra integral equation of the form (1) using the Chebyshev collocation method.
is method is based on approximating the unknown function u(x) by the Chebyshev interpolating polynomials at the Gauss-Lobatto collocation points defined on the domain [0, x N ] by [18][19][20] x j � 1 2 e unknown function u(x) is approximated as a truncated series of Chebyshev polynomials of the form where T k is the k th Chebyshev polynomial and u k are the Chebyshev coefficients. e derivatives of u(x) at the collocation points are represented as where r is the order of differentiation, D is the Chebyshev differentiation matrix, and N + 1 is the number of collocation points. We express the entries of the differentiation matrix D as [18][19][20] D 00 � 2N 2 + 1 6 , We introduce the following simple integral equation: where v(t) is a known function. Clearly, the existence of the function ϕ(x) depends on the existence of the integral. One approach is to note that the integral equation (7) is a special case of an initial value problem expressed as is differential equation can be solved using the Chebyshev spectral collocation method by replacing the differential operator d/dx with the transformed Chebyshev differentiation matrix D defined on the domain x ∈ [0, x N ] (see [21]) and also the functions v(x) and u(x) are written in vectors forms of sizes (N + 1) × 1; this yields where From the solution of the initial value problem (9), the integral equation (7) can be written in matrices form of unknown function u(x) as follows: where u 0 and v 0 are known values; here, the first row of transformed Chebyshev differential matrix D in equation (10) is interchanging by the row 1 0 · · · 0 . is is caused by the condition ϕ(0) � 0. Now, according to equations (7) and (10), we introduce the following transformation.

Definition. Let u(x)
and v(x) be continuous functions on [0, x N ]; then, the integral x 0 v(t)u(t)dt is transformed into matrices form as follows: where L is a square matrix of size (N + 1) × (N + 1) defined as and L v is a multiplication of matrix L by the vector v defined by 2 Mathematical Problems in Engineering where D is the modified matrix D after implementing initial condition. It is clear that this transformation is linear and it is easy to show this property. e existence and uniqueness of the linear transformation L defined in (12) depend on the existence and uniqueness of the inverse of matrix D. In linear algebra, if A is a square matrix of size N × N, then (i) A is invertible if and only if rank (A) � N and (ii) if A is invertible, then its inverse is unique [22]. We applied this theorem to all square matrices D to show that these matrices are invertible. In Figure 1, we display all the sizes of the matrices D of sizes (N + 1) × (N + 1) against the corresponding ranks for 1 ≤ N ≤ 200. We notice that, from Figure 1, the rank of the matrix D is equal to its size for all those values of N; that is, the matrix D is invertible and consequently the inverse is unique.
For the kernel K(x, t) of the Volterra integral equation (1), we assume that the function K(x, t) can be separated into two multiplied functions w(x) and v(t) as follows: Substituting (5), (11), and (14) in equation (1) gives the following system of linear equations: in which A is an (N + 1) × (N + 1) square matrix and F is (N + 1) × 1 column vector obtained by where diag [. . .] is the diagonal matrix of size (N + 1) × (N + 1) that puts the vector [. . .] of size (N + 1) × 1 on the main diagonal and I is an identity matrix of size (N + 1) × (N + 1). Finally, the solution of the Volterra integral equation (1) is obtained by solving the linear system (15) as e general idea underpinning the use of the proposed method is the replacement of a linear Volterra integral equation by a system of linear algebraic equations that replace the integral part in the integral equation by the introduced integral transformation. e obtained linear algebraic equations can easily be solved with the help of symbolic computation software such as Maple, Mathematica, MATLAB, and other symbolic computer packages.

Numerical Examples
To illustrate the effectiveness and accuracy of the proposed technique in the present work, several examples are carried out in this section. Example 1. Consider the following Volterra integral equation: To apply the proposed method to this test problem, we rewrite equation (19) in the following system of linear equations: where us, the final solution is obtained as Example 2. Consider the following Volterra-Fredholm integral equation [23]: . According to (14), we have Applying the proposed method to this problem, equation (23), we obtain the system of linear equations: where the rows of the matrix R are obtained by subtracting the first row from the last row of the matrix: us, the final solution of the Volterra-Fredholm integral equation (22) is obtained as Mathematical Problems in Engineering Example 3. Consider the following system of linear Volterra integral equations [24]: where the exact solution is u(x) � sin x and g(x) � e x . To illustrate the use of the proposed method to approximate the solution of this integral equations system, we rewrite the system as Now, converting the system above into a system of linear algebraic equations using the proposed method gives 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180  where Finally, the solution of the system of integral equations (30) is obtained as where u, F, and A are defined as Example 4. Consider the following quadratic nonlinear Volterra integral equation [25]: For solving the nonlinear Volterra integral equation (33), we first use the successive linearization method (SLM) [14][15][16][17] to linearize the nonlinear integral equation (33). e SLM was recently introduced as an efficient linearization algorithm and robust method for solving boundary value problems.
e procedure of the SLM is based on transforming the governing nonlinear equation into an iterative scheme made up of linear equations which is subsequently solved using the proposed method. In applying the SLM, we assume that the solution of (33) is obtained by where m is the iteration of SLM and u m (x) are unknown functions that are obtained by solving the linearized integral equation assuming that u m− 1 (x) are known functions from previous iterations. e procedure starts with an initial guess u 0 (x) which is chosen to satisfy the condition u 0 (0) � f(0). A suitable initial guess that satisfies this condition of equation (33) is e linearized integral equation to be solved is obtained by substituting equation (34) in integral equation (33) and deleting the nonlinear terms. us, the linearized equation is expressed as where x 0 (1/(1 + u m− 1 (t)))dt. Applying the proposed method on (36) gives the following matrices system: where where v 1 � 1/(1 + u m ) 2 and v 2 � 1/(1 + u m ). Finally, the solution of the integral equation (33) is given by Example 5. Consider the following system of nonlinear Volterra integral equations [26,27]: with the exact solution (u(x), g(x)) � (sin x, cos x). To apply the proposed method to this test problem, we first Mathematical Problems in Engineering linearize the equations using the SLM; we start by assuming that the solution of system (40) may be obtained by Substituting (41) into (40) gives the linearized form e suitable initial guesses for u 0 (x) and g 0 (x) are u 0 (x) � sin x − x and g 0 (x) � cos x − (sin 2 x/2). Applying the proposed method on (42) gives the following m th approximations of the linear system of equations: where u m , g m , F 1,m− 1 , and F 1,m− 1 are (N + 1) × 1 column vectors. System (43) leads to the matrix equation given as where U m , F m− 1 , and A m− 1 are defined as us, the final solution of the integral equations (40) for u(x) and g(x) is obtained as

Results and Discussion
In this section, we present the numerical results obtained using the proposed method for solving the Volterra integral equations. Implementation of the numerical schemes was performed using the personal computer of 2.5 GHz CPU speed including MATLAB 2017 package to perform the simulation results. e accuracy of the method is demonstrated by presenting infinity error norms u E (x) between exact and approximate solutions computed as where u(x) and u(x) are the exact and the numerical solutions, respectively. e accuracy of the current results was checked by comparing them with other methods reported in the literature. We also generate the computational times of the results to show the computational efficiency of the proposed method. e results of our solutions for all introduced cases in this work are displayed in the tables and graphs and discussed below. Table 1 gives a comparison between the absolute errors of the present results and the Adomian Decomposition Method (ADM) at different iterations at selected nodes for Example 1. A striking feature of the proposed method is that a high level of accuracy is achieved at the first order of approximation than the ADM, for example, whereas the proposed method gives the exact value at only one iteration for the 14 th digits.
is is because the proposed method selects all the terms as linear operator. Table 2 shows a comparison between the results obtained by the current method and Shifted Chebyshev Collocation (SCC) method reported by Youssri and Hafez [23]; both methods have been applied to the Volterra-Fredholm integral equation (Example 2). It is noticed that accurate results with errors of order up to eight digits are observed using the proposed method while the SCC method is accurate up to six digits when applied to this example. Also, we notice that increasing the values of N leads to improvement in the accuracy of the solution using the SCC method whereas increasing the values of N does not affect the improvement in the accuracy of the solution using the proposed method.
In Table 3, we present the numerical results of the absolute errors u E (x) and g E (x) for u(x) and g(x), respectively, at selected nodes obtained by the proposed method, the Bernstein Polynomials Expansion (BPE) method [24], the Artificial Neural Network (ANN) method [24], and the Trapezoidal Quadrature Rule (TQR) method [28] for Example 3. It is remarkable to note that accurate results with errors of order up to sixteen digits are obtained using the proposed method while the BPE, ANN, and TQR methods are accurate up to only two, five, and two digits, respectively. is indicates that the proposed method converges much faster than the BPE, ANN, and TQR methods. Table 4 shows the maximum absolute errors u E (x) and g E (x) for u(x) and g(x), respectively, obtained by the proposed method using different values of N. e maximum absolute errors are very small for small values of N. Also, increasing the number of nodes (increasing N) does not result in a significant improvement in the accuracy of the current approximation. We also note that the present method is computationally fast as accurate results are generated in a fraction of a second. In Figure 2, a convergence analysis graph for Example 3 is presented. Figure 2 presents a variation of the error norm at various values of nodes (N). It can be seen that the values of nodes start from N � 10 to converge fully. e accuracy of the method does not improve when increasing grid points N.
In Table 5, we have compared between the errors of present results for Example 4 corresponding to different values of collocation points N and those obtained by Maleknejad's method (MAL) reported in [29] and Mean Value eorem (MVT) reported in [25]. It can be seen that the present method converges rapidly compared to the     Table 6 shows the maximum absolute errors of u(x) for Example 4 at different orders of approximation at selected nodes. e maximum absolute errors are generally very small. A high level of accuracy is achieved at very low orders of approximation. We remark that increasing the number of nodes (i.e., increasing N) does not result in a significant improvement in the accuracy of the approximation. Figure 3 displays the errors against the iterations of the solution for Example 4 when N is fixed at 20. From Figure 3, it can be seen that the current solution rapidly converges, but the accuracy does not improve after 7 iterations.
In Table 7, we have compared between the errors of present results for Example 5 corresponding to different values of x and those obtained by Single-Term Walsh Series (STWS) reported in [26] and Homotopy Perturbation Method (HPM) reported in [27]. Table 7 indicates that 15 decimal places' accurate results are obtained using the present technique for large values of x while only 2 and 6 decimal places' accurate results are obtained using HPM and STWS, respectively. It is interesting to note from Table 7 that when x increases, this would lead to a decrease in the absolute error for HPM and STWS, while increasing x does not result in a decrease in the accuracy of the current approximation. is shows that the present method is more efficient than the method of STWS and HPM. Table 8 shows the maximum absolute errors u E (x) and g E (x) of the present method at different orders of approximation for different values of N for Example 5. It can be seen that the maximum absolute errors are very small and further decrease with an increase in the order of the solution. Also, a    large increase in the number of nodes (increasing N) does not result in a significant improvement in the accuracy of the present solution. In Figure 4, we display the absolute errors of u(x) and g(x) against collocation points N between the present results of Example 5 at the fifth iteration. From Figure 4, it is clear that the present method rapidly converges, and also, we notice that the accuracy does not improve for N > 13. x Errors in HPM [26] Errors in STWS [27] Errors in present method

Conclusion
In this paper, we have proposed a new method for solving linear and nonlinear Volterra integral equations. e technique has been constructed by introducing a new integral transformation together with Chebyshev pseudospectral method.
e main difference between this collocation method and the other collocation methods is the replacement of the integral terms by a linear transformation using the Chebyshev differential matrix. e numerical results have been shown and compared with the exact solutions and other methods in the literature. ese results are shown in the tables and figures for all cases considered in this work. To give a sense of the computational efficiency of the technique, the CPU computational time to generate the results is presented in the tables. It can be seen that high accuracy and less computational time (less than 0.02 seconds) substantiate our proposed method is a powerful numerical method for solving linear or nonlinear Volterra integral equations even for large values of grid points or iterations. e main conclusions emerging from this study may be summarized as follows: (1) e method is more flexible than other numerical methods such as BPE, ANN, TQR, STWS, HPM, MAT, and MVT since it allows us to choose the integral operators in the integral equations in terms of the Chebyshev spectral collocation differentiation matrix. (2) e proposed method is highly accurate and efficient, converges rapidly, and is sufficient to give good agreement with the exact solution. (3) e suggested technique does not depend on the number of collocation points N and it requires a few numbers of collocation points to achieve high accuracy of results.
Finally, we note that the proposed method described above is a powerful method that is appropriate in solving Volterra integral equations and can be generalized for Volterra and Fredholm integral equations. In future, we intend to show that the proposed method can be extended to coupled Volterra integrodifferential equations.

Data Availability
No data were used to support this study.

Conflicts of Interest
e author declares that there are no conflicts of interest.