Approximation Technique for Solving Linear Volterra Integro-Differential Equations with Boundary Conditions

This paper presents a new technique for solving linear Volterra integro-di ﬀ erential equations with boundary conditions. The method is based on the blending of the Chebyshev spectral methods. The application of the proposed method leads the Volterra integro-di ﬀ erential equation to a system of algebraic equations that are easy to solve. Some examples are introduced and the obtained results are compared with exact solution as well as the methods that reported in the literature to illustrate the e ﬀ ectiveness and accuracy of the method. The results demonstrate that there is congruence between the numerical and the exact results to a high order of accuracy. Tables were generated to verify the accuracy convergence of the method and error. Figures are presented to show the excellent agreement between the results of this study and the results from literature.


Introduction
Many Mathematical models of physical phenomena usually results in linear or nonlinear Volterra integral and integrodifferential equations. These equations play a crucial role in many fields of sciences such as economics, chemistry, finance, biology, and engineering. In engineering, physics, and mathematics areas, solving a problem means solving a set of differential equations and integral or integrodifferential equations. In general, finding analytic solutions of Volterra integro-differential equations is usually difficult so it is required to obtain an approximate solution. Therefore, they have been of great interest by several researchers and scientists of Volterra integro-differential equations. Several studies have been carried out and developed in the literature to find analytical or numerical solutions for Volterra integro-differential equations. Loh and Phang [1] used Genocchi polynomials to solve numerically a system of Volterra integro-differential equations which is done by approximating functions using Genocchi polynomials and derivatives using Genocchi polynomials operational matrix of integer-order derivative. Loh et. al [2] obtained the approximate solution for a new class of time-fractional partial integro-differential equation of the Caputo-Volterra type. Mirzaee et. al [3] used the moving least squares and spectral collocation method to estimate the solution of nonlinear stochastic Volterra integro-differential equations. A class of block boundary value methods has been constructed by Zhou and Stynes [4] for the solution of linear neutral Volterra integro-differential equations with weakly singular kernels. Adewumi et al. [5] presented new efficient numerical methods for solving Volterra integro-differential equations and a system of nonlinear delay integro-differential equations which arises in biology. In [6], the author discussed variational iteration method for the solutions of different linear and nonlinear Volterra integral equations. For some more recent and interesting results, we refer to [7][8][9][10][11][12] and the references cited therein.
Recently, Lotfi and Alipanah [13] presented a spectral element method to solve linear Volterra integro-differential equations; Alrabaiah et al. [14] used Laplace Adomian decomposition method to find the semianalytic solutions of nonlinear Volterra fractional integro-differential equations; a numerical method was obtained to find a solution for a singularly perturbed Volterra integro-differential equation by Nana et al. [15].
The main contribution of this work is to introduce an efficient numerical technique for solving the linear Volterra integro-differential equations together with boundary conditions. The technique is based on the Chebyshev spectral collocation method. The Chebyshev spectral collocation methods have been applied successfully in different fields of sciences and engineering because of their ability to give high accuracy of boundary value problems (see [16][17][18][19][20]. The application of the current method leads the Volterra integro-differential equation to a system of algebraic equations that are easy to solve when compared to a system of Volterra integro-differential equation. The main benefit of this method is that we used our proposed technique without any assumptions on large or small parameters. The investigation is mainly targeted to device this consistent method to integro-differential equations with boundary conditions. The main advantages of this approach over the standard same approaches are (i) the technique suggests a standard way of choosing the linear operator of the integro-differential equation whereas the other methods are choosing a linear operator to be simple in order to ensure that the integro-differential equation can be easily solved, and (ii) this algorithm transforms the integrodifferential equation into a system of linear algebraic equations that is easier and faster to solve when compared to a system of integral equations. This paper is organized as follows: in Section 2, we introduce a description of the proposed method. In Section 3, the method is implemented using some examples with known analytical solution. The results are discussed and investigated in Section 4. Finally, the conclusions are given in Section 5.

Description of the New Technique
Consider the following linear second order Volterra integrodifferential equation of the second kind given by under the boundary conditions where uðxÞ is an unknown function to be determined; P 0 ðxÞ, P 1 ðxÞ, P 2 ðxÞ, Kðx, tÞ, and f ðxÞ are known real valued To illustrate the idea of the new algorithm, we assume that the kernel Kðx, tÞ of Equation (1) can be separated into a product of two functions, namely, wðxÞ and vðtÞ; consequently, Equation (1) can be obtained as follows: Now, the integral term of Volterra integro-differential equation above is expressed as One approach is to note that the integral equation (4) is an initial value problem expressed as This differential equation is easy to solve using any method; here, we used the Chebyshev spectral collocation method; the functions ΦðxÞ, vðxÞ and uðxÞ are approximated as a truncated series of Chebyshev polynomials given by the form [21][22][23].
where T k is the kth Chebyshev polynomial andΦ,ũ, andṽ are the Chebyshev coefficients, and x j , j = 0, 1, 2, ⋯N are the Gauss-Lobatto collocation points (see [21]) defined by where N + 1 is the number of collocation points. The derivatives at the collocation points are represented as where r is the order of differentiation and D is the Chebyshev spectral differentiation matrix whose entries (see [22]) are given by 2 Abstract and Applied Analysis where we observe that the first row of Chebyshev differential matrix D in equation above is replaced by the row ½1, 0, 0, ⋯, and we add the vector −½u 0 v 0 , u 0 v 0 , ⋯, u 0 v 0 T ; this is caused by imposing the condition ϕð0Þ = 0 into Equation (10). According to the integral Equation (4) and Equation (11), we introduce the following integral operator L ½: ð:Þ defined as The operator L is square matrix of size ðN + 1Þ × ðN + 1Þ: 2.1. The Linearity of the Operator L. It is clear that the above operator is a linear operator because Therefore, the integral operator L is a linear operator.

Existence and Uniqueness of the Operator L.
In this subsection, we show that the operator L exists and is unique.

Theorem 1.
(i) Any square matrix A is invertible (nonsingular) if and only if detðAÞ ≠ 0 (ii) If A is an invertible matrix, then its inverse is exist and unique Proof. See Larson and David [24].
According to the theorem above, the existence and uniqueness of the linear operator L depends on the existence and uniqueness of inverse of matrix : ð15Þ The determinants of matrix above have been computed on the domain ½0, b for N = 10, 11, ⋯, 100. In Figure 1, we plotted the determinants for different values of b varied N. We observe that from the figure, all the computed determinants are not equal to zero, and it is noticed that all the determinants are greater than or equal to 1, i.e., the inverse of the operator L exists and is unique. Now, replace the derivative operator d i /dx i of the integro-differential Equation (1) by the Chebyshev spectral differentiation matrix D i defined in (9) and also replace the integral term of the integro-differential Equation (1) by the operator L, which gives 3 Abstract and Applied Analysis Here, the superscript T denotes matrix transpose, ½⋯ is the diagonal matrix of size ðN + 1Þ × ðN + 1Þ that puts the vector ½⋯ of size ðN + 1Þ × 1 on the main diagonal, and ½P 2 ðxÞ, ½P 1 ðxÞ, ½P 0 ðxÞ, ½ f ðxÞ, ½vðxÞ, and ½uðxÞ are vectors of sizes ðN + 1Þ × 1. The functions wðxÞ and vðxÞ are obtained from the kernel function Kðx, tÞ which can be written as Kðx, tÞ = wðxÞvðxÞ as we mentioned before.
Thus, Equation (16) can be written as where A is a ðN + 1Þ × ðN + 1Þ square matrix and F is column vector of sizes ðN + 1Þ × 1 defined by The boundary conditions uð0Þ = uðx 0 Þ = α and uðbÞ = u ðx N Þ = β are implemented into the system (17) by modifying the first and last rows of matrices A and F in such a way that the the system of linear Equation (17) takes the form After modifying the matrix system (17) to incorporate boundary conditions, the final approximate solution of the linear second order Volterra integro-differential Equation (1) is obtained as whereÃ andF are the the modified matrices of A and F, respectively. The general idea underpinning the use of the proposed method is to convert the linear Volterra integro-differential equation by a system of linear algebraic equations that replace the derivative and integral parts in the integral equation by the differential matrix and the introduced integral operator, respectively. The 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
In this section, some numerical examples are carried out to illustrate the effectiveness and accuracy of the proposed method. In addition, the numerical results are compared with exact solution, and all of them were performed on the  Abstract and Applied Analysis computer using a program written in MATLAB v7.5.0 (R2007b). We also showed in figures and tables the approximate solutions, exact solutions, and the absolute error function at the selected points of the given interval.
Example 1. Consider the following linear second-order Volterra integro-differential Equation [25] together with the boundary conditions uð0Þ = 1 and uð1Þ = cos ð1Þ. The exact solution is uðxÞ = cos ðxÞ. The technique of the solution starts by rewriting the integral term in the above equation as Applying the proposed algorithm on Equation (21) and according to the assumptions (8) and (12), one can write the integro-differential Equation (21) in a system of linear equations: where A = D 2 + I + diag ½xL ½tan ðxÞ , F = ½xð1 − cos ðxÞ T , and I is an identity matrix.
To implement the boundary conditions to the systems (23), we replace all the entries of the first and last rows of the matrix A by zeros and then set A 00 = A NN = 1, F 0 = 1, and F N = cos ð1Þ. After modifying the matrix system (23) to incorporate the boundary conditions, the solution is obtained as   Example 2. Consider the following Volterra integrodifferential equation [26]: subject to the boundary conditions uð0Þ = −1, uðπ/2Þ = 1.
The integral term of integro-differential Equation (25) can be expressed as According to the assumptions (8) and (12), we can apply the proposed algorithm on the integro-differential Equation (25); the equation is transformed into the following system of linear equations where A = D 2 + I − diag ½sin ðxÞL ½e −x D and F = ½1/2e −x sin ð2xÞ − 2 sin ðxÞ T . After modifying the matrix system (27) to incorporate the boundary conditions of the integrodifferential Equation (25), the solution is obtained as whereÃ andF are the modified matrices of A and F, respectively.
Example 3. Consider the linear Volterra integro-differential equation given by [26,27] with the initial conditions uð0Þ = 1 and the exact solution u ðxÞ = e x 2 .   Abstract and Applied Analysis According to the assumptions (8) and (12), one can express the integral term of integro-differential Equation (29) as Applying the proposed algorithm on the integrodifferential Equation (29), the equation is transformed into the following system of linear equations where T . Now, incorporate the initial condition uð0Þ = 1 into the system (31) to obtain the final solution of integrodifferential Equation (25) as whereÃ andF are the the modified matrices of A and F, respectively, after imposing the initial condition.
Example 4. Consider the following system of linear Volterra integral equations subject to the boundary conditions where Applying the proposed algorithm, one can express the integral parts of integro-differential Equations (33) and (34) as   [26], practical matrix method [27], and present results varied x of Example 3.
x Tau method Practical matrix method Present method Applying the proposed algorithm on the system of integro-differential Equations (33) and (34), the equation is transformed into the following system of linear equations where A is a ð2N + 2Þ × ð2N + 2Þ square matrix while F and U are ð2N + 2Þ × 1 column vectors defined by where In the above definitions, ½1 is a column vector whose interies are 1 ′ s. The boundary conditions (35) are imposed on Equation (38) by modifying the first and last rows of A mn , ðm, n = 1, 2Þ and f 1 , f 2 in such a way that the modified matrices A and F take the form

Results and Discussion
In this section, we give the results obtained by applying the proposed algorithm on various linear Volterra integrodifferential equations. Implementation of the numerical schemes was performed using personal computer of 2.5 GHz CPU speed including Matlab2017 package to perform the simulation results. The accuracy of the method is demonstrated by presenting infinity error between exact and approximate solutions computed as where u E ðxÞ and u N ðxÞ are the exact and numerical solutions, respectively. To check the validity and accuracy of the current algorithm, we made comparisons between the obtained results with the exact solutions and also with those other methods reported in the literature. The results of our solutions for all introduced examples in this report are showed and discussed in tables and graphs. Table 1 shows absolute errors of Example 1 for N = 5,15,40, and 60 at selected values of x between 0 and 1. Also, the time taken for the computation has been presented at the bottom of the table. The absolute errors for various values of number of grid points N are shown in Figure 2. From Table 1 and Figure 2, clearly that the maximum absolute errors are generally very small and results obtained by the current algorithm are almost the same as the results of the exact solution, it is also noted that when increasing the number of grid points N does not result in a significant improvement in the accuracy of the current approximation. It is clear that from Table 1, the present algorithm is computationally fast as accurate results are generated in a fraction of a second as it is shown in bottom of the table.
In Table 2, the exact errors of Example 2 have been calculated for various x in ½0, π/2 to measure the extent of agreement between the exact and approximate solution. We also presented the time taken for the computation at the end of the table. From the table, it can be seen that the current algorithm provides us with the accurate approximate solution of Example 2. We remark that increasing the 8 Abstract and Applied Analysis number of nodes (i.e., increasing N) does not result in a significant improvement in the accuracy of the approximation. We also compare the relative error estimates with those obtained using the homotopy analysis method (HAM) reported by El-Ajou et al. [28] for x in ½0:1,1 in Table 3. Full convergence upto 7 decimal places accurate results reported by HAM is achieved for x ≥ 0:8, and an increase in the values of x results in an increase in the maximum error form, while high convergence is achieved to 15 decimal places accurate results for the present method for all values of x; also, increasing in nodes (N) leads to decrease in the accuracy of HAM for Example 2 while in increasing N does not result in a significant improvement in the accuracy using the current approximation. One can note that the present method gives    Figure 3 shows the errors between the exact and numerical solutions for Example 2 varied N. Very small errors are achieved (10 −14 < Error < 10 −12 ) for the values N ≥ 15, i.e., that accuracy does not improved when for large values of N.
In Table 4, the absolute error is obtained for Example 3 using current method, Tau method investigated by Hosseini and Shahmorad [26], and practical matrix method obtained by Yüzbaşı et al. [27]. It is seen from the table that the results obtained by the present method is better than those obtained by Tau and practical matrix methods; also, it is noticed from Table 4 that when x increased, the error is increased for Tau and practical matrix methods while in the present method the error is approximately fixed for all values of x in the domain. Figure 4 showed the error of the numerical results against the nodes N of Example 3. It can be seen from this graph that when the number of nodes (N) increase, the error becomes smaller, and then, no effect when N ≥ 15. In Table 5, we give the maximum errors between the exact and present results of Example 3; the results were computed in the domain 0:1 ≤ x ≤ 1. To give a sense of the computational efficiency of the method, the computational time to generate the results is also given. The accuracy is seen to improve with an increase in the number of collocation points N. It is observed that accurate results with errors of order up to 10 −15 are obtained using very few collocation points. We remark, also, that this algorithm is computationally fast as accurate results are generated in a fraction of a second. Table 6 and Figure 5 are showed the maximum absolute errors of the approximate solutions for Example 4 at selected values of x and N. The maximum absolute errors are very small for uðxÞ and vðxÞ and further decrease with an increase in the values of N upto N = 18, and then, increasing N does not result in a significant improvement in the accuracy of the approximate solution. We also observed that the current algorithm is computationally fast as accurate results are generated in a fraction of a second even for large values on N. Figure 5 shows a comparison between the exact and approximate solutions for Example 4 of uðxÞ and vðxÞ. We observe that there is good agreement between the two sets of results, proving the efficiency of the proposed approach.
In Figure 6, the maximum absolute errors for Example 4 is presented. The figures display a variation of the errors at a various values of N. It can be seen that, in the problem considered, the algorithm needs small values of N to achieve good currency and to converge fully between the exact and approximate results. Also, we notice that there is no significant improvement in the accuracy when the number of nodes (N) is increasing.
A comparison between numerical results with the exact solution for Example 4 is shown in Figure 5. We can see that  11 Abstract and Applied Analysis accurate results with errors of order up to 10 −10 . It is worth remarking that the accuracy of the method is not dependent on the number of collocation points.
All the above observations are clear indication that this algorithm is a powerful method that is appropriate in solving linear Volterra integro-differential equations.

Conclusion
In this study, we have proposed a simple and easy algorithm for solving linear Volterra integro-differential equations with boundary conditions. This technique has been constructed by introducing a new integral operator together with Chebyshev pseudospectral method. The method allows us to choose the integral operators in the integro-differential equations in terms of the Chebyshev spectral differentiation matrix. The numerical solutions have been shown in tables and graphs and compared with the exact solutions and other methods in the literature. The accuracy of the obtained result is a point of interest; moreover, if we increase the number of nodes approximation, the numerical solution get closer to the exact solutions. The fast convergence and accuracy of the proposed algorithm are valid reasons for researcher to use this method for different problems of Volterra integro-differential equations arising in various areas of science.
The main conclusions emerging from this study may be summarized as follows:

12
Abstract and Applied Analysis (1) The proposed method is highly accurate and efficient, converges rapidly, very easy to understand, and is sufficient to give good agreement with the exact solution (2) The suggested technique does not depend on the number of collocation points, and it requires a few numbers of collocation points to achieve high accuracy of results (3) No any iteration is required to achieve the currency while other methods required more iterations (4) The method is more flexible than other numerical methods in the literature such as HAM, Tau method, and practical matrix method Finally, for the future study, we suggest that how is the nonseparable kernel dealt with and also how to apply this technique on nonlinear Volterra and Fredholm integral equations.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.