Numerical Solution of the Fredholm and Volterra Integral Equations by Using Modiﬁed Bernstein–Kantorovich Operators

: The main aim of this paper is to numerically solve the ﬁrst kind linear Fredholm and Volterra integral equations by using Modiﬁed Bernstein–Kantorovich operators. The unknown function in the ﬁrst kind integral equation is approximated by using the Modiﬁed Bernstein–Kantorovich operators. Hence, by using discretization, the obtained linear equations are transformed into systems of algebraic linear equations. Due to the sensitivity of the solutions on the input data, signiﬁcant difﬁculties may be encountered, leading to instabilities in the results during actualization. Conse-quently, to improve on the stability of the solutions which imply the accuracy of the desired results, regularization features are built into the proposed numerical approach. More stable approximations to the solutions of the Fredholm and Volterra integral equations are obtained especially when high order approximations are used by the Modiﬁed Bernstein–Kantorovich operators. Test problems are constructed to show the computational efﬁciency, applicability and the accuracy of the method. Furthermore, the method is also applied to second kind Volterra integral equations.


Introduction
Fredholm and Volterra integral equations of the first kind play an important role in many problems from science and engineering. It is known that the Fredholm integral equations can be derived from boundary value problems with given boundary conditions. For example, Fredholm integral equations of the first kind arise in a mathematical model of the transport of fluorescein across the blood-retina barrier in the transient state and the subsequent diffusion of fluorescein in the vitreous body given in Larsen et al. [1]. Some other applications are in palaeoclimatology given in Anderssen and Saull [2], antenna design in Herrington [3], astrometry in Craig and Brown [4], image restoration in Andrews and Hunt [5]. The investigation of Volterra integral equations is very important in solving initial value problems of usual and fractional differential equations arising from the mathematical modelling of many scientific problems, including population dynamics, spread of epidemics, and semi-conductor devices, such as the biological fractional n-species delayed cooperation model of Lotka-Volterra type given in Tuladhar et al. [6]. Examples of Volterra integral equations of first kind can be extended to mathematical model of animal studies of the effect of the deposition of radioactive debris in the lung by Hendry [7], the heat conduction problem in Bartoshevich [8], tautochrone problem of which Abel's integral equation was derived by Abel [9], (see also Groetsch [10]), electroelastic of dynamics of a nonhomogeneous spherically isotropic piezoelectric hollow sphere problem in Ding et al. [11]. Additionally, the use of a dynamical model of Volterra integral equations in energy storage with renewable and diesel generation has been analysed in Sidorov et al. [12].
As a classical ill-posed problem, the numerical solution of Fredholm integral equations of the first kind has been investigated by many authors, as a former study by Phillips [13] and a recent study by Neggal et al. [14]. The well-known early methods are the regularization methods given with a technique by Phillips in [13] and the Tikhonov regularization by Tikhonov in [15,16]. In the Tikhonov method, a continuous functional is usually used and the minimizer for the corresponding functional is difficult to obtain. Consequently, several methods have been proposed to obtain an effective choice of the regularization parameter in Tikhonov method such as the discrepancy principle, the quasi-optimality criterion (see Groetsch [17], Bazan [18] and references therein). Further, in Caldwell [19], a direct quadrature method and a boundary-integral method were examined for solving Fredholm integral equations of the first kind. Additionally, a regularization technique which replaces ill-posed equations of the first kind by well-posed equations of the second kind was employed to produce meaningful results for comparison purposes. Later, the extrapolation technique by Brezinski et al. [20] and a modified Tikhonov regularization method to solve the Fredholm integral equation of the first kind under the assumption that measured data are contaminated with deterministic errors was given in Wen and Wei [21]. Recently, a variant of projected Tikhonov regularization method for solving Fredholm integral equations of the first kind was proposed in Neggal et al. [14] in which for the subspace of projection, the Legendre polynomials were used.
Early studies for the solution of Volterra integral equations of the first kind involve the high order block by block methods in Hoog and Weiss [22,23]. However, these methods suffer from the disadvantage of requiring additional evaluations of the kernels and the solution of systems of algebraic equations for each step. Later, Taylor [24] used inverted differentiation formulae, which the resulting methods were explicit corresponding to local differentiation formulae. As the author stated "the main disadvantage of this method is that weights must be calculated from the recurrence relation (2.9) and the differentiation formula must be chosen so that the Dahlquist root condition is satisfied". Integral equations of the first kind associated with strictly monotone Volterra integral operators were solved in Brunner [25] by projecting the exact solution of such an equation into the space S (−1) m (Z N ) of piecewise polynomials of degree m ≥ 0 possessing jump discontinuities on the set Z N of knots. Besides, the asymptotic behavior of solutions to nonlinear Volterra integral equations was analysed in Hulbert and Reich [26]. The future-sequential regularization method and predictor-corrector regularization method for the approximation of Volterra integral problems of first kind with convolution kernel were given in Lamm [27] and Lamm [28], respectively. The numerical solution of Volterra integral equations of the first kind by sequential Tikhonov regularization coupled with several standard discretizations (collocation-based methods, rectangular quadrature, or midpoint quadrature) was given in Lamm and Eldén [29].
New approaches have been developed for the solution of integral equations that use the basis functions and transform the integral equation to the system of linear or nonlinear equations. One of these approaches is the use of wavelet basis. For the solution of Abel's integral equation, Legendre wavelets were used in Yousefi [30] and the wavelet basis were used in Maleknejad et al. [31] for the numerical solution of Volterra type integral equations of the first kind. Another approach is the use of polynomial approximations. In Mandal and Bhattacharya [32], Fredholm integral equations of the second kind and a simple hypersingular integral equation and a hypersingular integral equation of the second kind were numerically solved using Bernstein polynomials. At the same year, in Maleknejad et al. [33] numerical solution of linear and nonlinear Volterra integral equations, of the second kind by using Chebyshev polynomials was given. Afterwards, a new approach to the numerical solution of Volterra integral equations by using Bernstein's approximation was given in Maleknejad et al. [34].
Recently, exhaustive studies on the use of CESTAC method for the solution of Volterra first type integral equations has been given in Noeiaghdam et al. [35] in which the control of accuracy on Taylor-collocation method to solve the weakly regular Volterra integral equations of the first kind has been studied. Furthermore, in Noeiaghdam et al. [36] that the numerical validation of the Adomian decomposition method for solving Volterra integral equation with discontinuous kernels was given.
The need of stable, reliable and time efficient methods for the numerical solution of Fredholm and Volterra integral equations of first kind is the main motivation of contributions. The achievements of the study can be summarised as follows: 1. Using the Modified Bernstein-Kantorovich operators, a numerical approach is developed for the solution of Fredholm and Volterra integral equations of the first kind with continuous and square integrable kernels. Convergence analysis are given assuming that minimum norm least square solution of the obtained algebraic linear systems are obtained by using the exact data, that is to say the Moore-Penrose inverse of the resulting coefficient matrices are computed exactly.
2. Furthermore, regularized integral equations are considered to obtain more smooth solutions especially when high-order approximations are used by Modified Bernstein-Kantorovich operators. The proposed approach is applied by building regularization features into the algorithm and perturbation error analysis are given.
3. Test problems are conducted and theoretical results are justified with the obtained numerical results.

Asymptotic Rate of Convergence of Modified Bernstein-Kantorovich Operators
The Modified Bernstein-Kantorovich operators K n,α ( f ; x) were used to approximate a function f : [0, 1] → R (see Özarslan and Duman [37]) where, and and α > 0 is constant. For α = 1 reduces to classical Bernstein-Kantorovich operator where, Proof. From (1) it follows that where For each fixed n ∈ N, α > 0 the inequality (4) is obtained from β(α) = max occurring at x = 1 α+1 and max (Y(k)) 2 and P 2 = ρ(P T P) to present the discrete Euclidean norm of a vector Y ∈ R n and the spectral norm of a matrix P ∈ R n×n , respectively, where ρ is the spectral radius and P T is the transpose of P. Voronowskaja [38] gave the asymptotic rate of convergence of the Bernstein operators using the linearity property of the Bernstein operators and Taylor formula at a point x as Based on the analogous approach in Voronowskaja [38] we give the asymptotic rate of convergence of the Modified Bernstein-Kantorovich operators by the next theorem. 1], and admits a derivative of second order at some point Additionally this limit is uniform if f ∈ C 2 [0, 1], thus the rate of convergence of the operator Proof. Assume that f is integrable in [0, 1], and has second order derivative at a point x ∈ [0, 1] then from Taylor's formula at x we have and E(u) → 0 as u → 0 and E is integrable function on [−x, 1 − x]. Using the linearity property of the operators K n,α and (8)-(10) we have where, To show that the asymptotic rate of convergence is O 1 n , it is sufficient to show that lim |E(u)| and for arbitrary ε > 0 there exist Using Lemma 1 estimation (5) gives where, σ(α) is as given in (7) and M(n, α) = sup |Q(n, α, x)|. In addition for a fixed α, M(n, α) is second degree polynomial in n and Q(n, α, x) is It is obvious from (18) and (19) that for n large enough we have |nE(α, n, x)| < ε and using (15) we obtain (13). If f ∈ C 2 [0, 1] then this limit is uniform, thus the rate of convergence of the operator sup hold true where, ϕ(α, x) is the given function in (10) and β(α), σ(α) are as given in (6) and (7), respectively, and γ 2 (α) is the same as in (17).

Representation of the K n,α Operators and Discretization of First Kind Integral Equations
We consider the Fredholm integral equation of the first kind (FK1) and Volterra integral equations of the first kind (VK1) where g(x) is called the free term while K(x, t) is called the kernel and f (t) is the unknown function to be determined. [17,39]) By means of the singular value expansion (SVE) any square integrable kernel K(x, t) can be written in the form

Definition 1. (Groetsch
The functions u i , v i are the singular functions of K and they are orthonormal with respect to the usual inner product (., .) and the number µ i are the singular values of K. For degenerate kernels the infinite sum (29) is replaced with the finite sum upto the rank of the kernel. The system {u i , v i ; µ i } is called the singular system of K.
Let Ψ : H 1 → H 2 be a compact linear operator on a real Hilbert space H 1 , taking values in a real Hilbert space H 2 . The next theorem is known as the Picard's theorem on the existence of the solutions of first kind equations.
On the basis of Theorem 3 we consider the Hypothesis 1 as follows: Hypothesis 1.
Without loss of generality, the solution f of FK1 and VK1 denotes the pseudoinverse solution or the Moore-Penrose generalized inverse solution for FK1 and VK1 respectively. Further, in order to determine the effect of α > 0 in the numerical solution we represent the Modified Bernstein-Kantorovich operators (1) for 0 < µ < 1 in the form where For the numerical solution of FK1 and VK1, we approximate the function f by using the Modified Bernstein-Kantorovich operators in (32). We obtain the following equation and for VK1 we get Subsequently we take the grid points x j = j n + , j = 0, 1, . . . , n − 1 and x n = 1 − , where 0 < < 1 2n . Then, the Equations (35) and (36) respectively, where the coefficient matrices A and A have the entries = ω . . , n, k = 0, 1, . . . , n, and q(u) and ω are as given in (33 ) and (34), respectively. The coefficient matrices A and A in (37) are ill-conditioned matrices and may be rank deficient or even singular matrices. Therefore, we consider the following minimum norm least squares problem for FK1 and for VK1 min Lemma 2. The problems (42) and (43) have the unique minimum norm least squares solutions X = A † B and X = A † B respectively.
Proof. Proof is analogous to the proof of Theorem 1.2.10 in Björck [40].

Convergence Analysis
By solving the algebraic systems (42) and (43) we get a numerical solution of the unknown (40) and denote this approximation by X n . Further, let us use F n to denote the obtained numerical approximation to f that is in the implicit form in X n and obtained by using the proposed approach. Substituting F n in (32) we get K n,α (F n ; x) as [40]) The condition number of U ∈ R m×n (U = 0) is

Remark 1.
If the matrix A in (38) and the matrix A in (39) are invertible then A † = A −1 and A † = A −1 and the inequalities (45) and (46) hold true.

Regularized Numerical Solution
The numerical solution of the general least squares problems (42) and (43) may be extremely difficult because the solution is very sensitive to the perturbations of the coefficient matrices A and A and the right side vector B. This is reflected in the fact that κ(A), and κ( A) are very large and increases as n increases which is the degree of the constructed polynomial by the Modified Bernstein-Kantorovich operator used for the approximation of the solution. High condition numbers of the matrices A and A cause rounding errors that prevent the computation of an accurate numerical solution of the problems (42) and (43), respectively. Moreover, the obtained discrete problems are always perturbed by approximations such as the integrals given as the entries of A and A are evaluated numerically. Therefore, even if we were able to solve the discrete algebraic problems (42) and (43) without rounding errors we would not obtain a "smooth" solution because of the oscillations in the singular vectors. By a smooth solution we mean "a solution which has some useful properties in common with the exact solution to the underlying and unknown unperturbed problem" as stated in Hansel [41]. Furthermore, the function g is typically a measured or observed quantity and hence, in practice, the true g is not available to us. On one hand, the estimate g δ of g satisfying g δ − g 2 ≤ δ and δ is the priori error level is known(see Tikhonov [15] and [16] and Groetsch [17]). Therefore, we consider the following regularized problems for the Fredholm integral equation of the first kind (RFK1) (see Tikhonov [15,16] and Groetsch [17] and Volterra integral equations of the first kind (RVK1) It is clear that (67) and (68) are second kind Fredholm and Volterra integral equations, respectively. For the numerical solution of RFK1 and RVK1 by the proposed method M(K n,α ) we take the grid points x j = j n + , j = 0, 1, . . . , n − 1 and x n = 1 − , where 0 < < 1 2n and is sufficiently small number also η(δ) > 0 is called the regularization parameter. We assume the following algebraic equations for RFK1 and for RVK1 for j = 0, 1, . . . , n, k = 0, 1, . . . , n. Then, the discrete regularized Equations (69) and (70) can be presented in matrix form for the RFK1 and for the RVK1 respectively where, and the vector B ∈ R n+1 which can be written as B = B + ∆B such that ∆B is the priori error level ∆B ≤ δ. Furthermore, A = A + ∆A where A is the matrix in (38) and ∆A = ωη(δ)I + ∆ 1 A, with the addition of diagonal matrix ωη(δ)I and ∆ 1 A which is the defect matrix of the numerical errors of the computation of the integrals in (69) with a predescribed error δ * = δ * (δ) ≥ 0, depending on δ. Analogously, A = A + ∆ A and A is as in (39) and the matrix ∆ A = ωη(δ)I + ∆ 1 A has the defect matrix ∆ 1 A of the numerical errors of the computed integrals in (70) with a predescribed error δ * = δ * (δ) ≥ 0. Therefore, it is possible Clearly, the numbers h and δ are estimates of the errors of the approximate data A, B , A, B of the problem (37) for FK1 and VK1, respectively, with the exact data (A, B), A, B accordingly. Thus, the given regularized systems (71) uses h and δ explicitly. For the implementation of the approach we have taken η(δ) = δ. In this connection, about the remarks on choosing the regularization parameter using the quasi-optimality and ratio criterion, see Bakushinskii [42] and for the data errors and an error estimation for ill-posed problems see Yagola et al. [43]. Next, we consider the following general least squares problem for RFK1 and for RVK1 Theorem 6. (Theorem 1.4.6 in Björck [40]) Assume that rank(U + ∆U) = rank(U) and let Then if η = κ(U) U < 1 the perturbations ∆X and ∆r in the least squares solution X and the residual r = B − UX satisfy Let X δ η,n denote the minimum norm solution obtained by solving the general least squares problems (75) and (76). Further, F δ η,n denote the obtained approximation to function f δ η appearing implicitly in (72). Substituting F δ η,n in (32) we get K n,α F δ η,n ; x as We also present the residual error of the obtained algebraic linear system (37) for FK1 by r = B − AX ( r = B − AX for VK1). The regularized residual error of the system (71) RVK1). Furthermore, the corresponding numerical calculation of the regularized residual error is r δ η,n = B − AX δ η,n ( r δ η,n = B − AX δ η,n ) accordingly. Next, the following priory bound for the error of the approximation follows.
Theorem 7. Assume that the conditions of Hypothesis I are satisfied and the solution f δ η of (67) belongs to C λ ∩ L 2 ([0, 1]) for some λ ≥ 2. Consider the regularized linear system AX δ η = B given in (71) where A = A + ∆A and A is the matrix in (38) and ∆A 2 ≤ h. Furthermore, B = B + ∆B as in (74) and B is the vector in (41) and ∆B 2 ≤ δ. Additionally X δ η = X + ∆X and r δ η = r + ∆r and let S(j + 1) = sup t∈[0,1] K x j , t for x j = j n + , j = 0, 1, . . . , n − 1 and x n = 1 − , where 0 < < 1 2n and M 2 = S 2 . Further, If rank( A) = rank(A) and η = κ(A) A < 1 then hold true where, η(δ) is the regularization parameter and W 1 n, α, f δ η , W 2 n, α, f δ η are as in (47)  Proof. For RFK1, it follows that Based on Corollary 1 and the estimation (22) by replacing f with f δ η in estimation (22) and in (47) we obtain (32) and (44) and using that n ∑ k=0 P n,k (x) = 1, follows where For the numerical solution of RFK1 in (67) we use the grid points x j = j n + , j = 0, 1, . . . , n − 1 and x n = 1 − , where 0 < < 1 2n . We assume where ωX δ η (j + 1) gives the average value of f δ η over the interval j n+1 , j+1 n+1 . If we substitute F δ η,n instead of f δ η in (89) we get a new function g δ on the right side of this equation Thus, for RFK1 from (89) and (90) we obtain where, X δ η is as given in (88). The general least squares problem of (91) has the minimum norm solution X Thus, From the assumption that f δ η ∈ C λ ∩ L 2 ([0, 1]) for some λ ≥ 2 it follows that sup and replacing f with f δ η in estimations (21) and (48) we obtain (98) Substituting the estimation (98) into (94) and the result in (87) we get Inserting (86) and (99) in (85) and on the basis of Theorem 5 and using that κ(A * ) = A † * 2 A * 2 we obtain (82). The inequality (83) is obtained by using and based on the Theorem 6 and the inequality (78) the first term on the right side of (100) is obtained as Next, on the basis of Theorem 5 and using (93), (98) and A † 2 = κ(A) Inserting the estimations (101) and (102) into (100) gives (83). To prove the inequality (84), we use and based on Theorem 6 and the inequality (79), the first term on the right side of (103) is obtained as The second error term on the right side of (103) satisfies using (102), (103) and (104), and that A 2 ≤ A 2 + h and from (81) follows (84).

Numerical Results
For the theoretical results given in Sections 2-4, we focus on the interval [0, 1]; however, for the numerical results, we also consider examples on [a, b] with the following extention of the Bernstein operators and Modified Bernstein-Kantorovich operators on the interval [a, b] respectively. All the computations in this section are performed using Mathematica in machine precision on a personal computer with properties AMD Ryzen 7 1800X Eight Core Processor 3.60 GHz. We remark that the solution of the Volterra integral equations by using Bernstein polynomials was given in Maleknejad et al. [34]. All the considered test problems are also solved by using Bernstein operators (11) with the approach given in Maleknejad et al. [34]; additionally, regularization is applied. Further, the obtained algebraic system of equations by applying the methods M(K n,α ) and M(B n ) are solved using the pseudoinverse of the respective matrices. Let the following error grid functions be defined at the N + 1 grid points x p = a + p(b−a) N , p = 0, 1, . . . , N over the interval [a, b] as Further, we use the following notations in tables and figures: M(K n,α ) presents the given approach by using the Modified Bernstein-Kantorovich operators K n,α .
M(B n ) presents the approach in Maleknejad et al. [34] by using the Bernstein operators B n .
Cond B n A denotes the condition number of the perturbed matrix A obtained by the method M(B n ) using LinearAlgebra'Private'MatrixConditionNumber command in Mathematica.
Cond K n,α A denotes the condition number of the perturbed matrix A obtained by the method M(K n,α ) using LinearAlgebra'Private'MatrixConditionNumber command in Mathematica.

RE E N (K n,α ) denotes the root mean square error (RMSE) of the regularized solution
obtained by M(K n,α ).

RE E N (B n ) denotes RMSE of the regularized solution
obtained by M(B n ).
AE E N ,x p (K n,α ) is the absolute error of the regularized solution E N K n,α F δ η,n ; x p at the point x p . AE E N ,x p (B n ) is the absolute error of the regularized solution E N B n F δ η,n ; x p at the point x p . ME E N (K n,α ) shows the maximum error ME of the regularized solution max 0≤p≤N E N K n,α F δ η,n ; x p . ME E N (B n ) shows the maximum error ME of the regularized solution max 0≤p≤N E N B n F δ η,n ; x p .
na means that the specified method is not applied to the considered example.
ng means that the absolute error is not given at the presented grid point by the specified method.

Application on Examples of Fredholm Integral Equations
We consider the following test problems of first kind Fredholm integral equations, which have been used as benchmark problems in the literature. [21] and Baker et al. [44])

Example 1. FK1 (Wen and Wei
and the exact solution is f (x) = e x .
and the exact solution is f (x) = x 3 2 . Table 1 presents the RMSE with respect to n obtained by the proposed approach M(K n,10 ) when N = 51 and = 0.0001, for the examples of FK1 and when δ = 5 × 10 −12 for the Example 1, Example 2 and Example 4 and δ = 5 × 10 −9 for the Example 3. The absolute errors obtained by the method M(K 9,10 ) at the points x p = p 8 , p = 0, 1, . . . , 8 for the examples FK1 when = 0.0001, n = 9 and α = 10 for the same values of δ as in Table 1 are demonstrated in Table 2. Further, Table 3 shows the same quantities as in Table 1 obtained by using the approach M(B n ). Tables 4-7 present the condition numbers of the perturbed matrices, RMSE with respect to the δ obtained by the proposed method M(K 8,1 ) and the method M(B 8 ) when ε = 0.0001, and N = 51 for the Example 1, Example 2, Example 3 and Example 4, respectively. Table 8 presents the RMSE with respect to α obtained by the proposed approach M(K 9,α ), when N = 51 and = 0.0001, for the examples of FK1. In this table, the parameter δ is taken as δ = 5 × 10 −12 for the Example 1, Example 2 and Example 4 and δ = 5 × 10 −9 for the Example 3.   Figure 1 presents the RMSE with respect to α obtained by M(K 9,α ) for the examples of FK1, when = 0.0001, and N = 51. It can be viewed that the optimal value of α is α = 10 for the Example 1, and Example 4, whereas α = 0.1 gives the lowest RMSE for the Example 2 and Example 3. Figure 2 illustrates the RMSE with respect to n obtained by the methods M(K n,10 ) and M(B n ) for the considered examples of FK1 when = 0.0001 and N = 51. Furthermore, for the data in Figures 1 and 2 the regularization parameter η(δ) is taken as δ = 5 × 10 −12 for the Example 1, Example 2 and Example 4 and δ = 5 × 10 −9 for the Example 3. Figure 3 shows the RMSE with respect to δ obtained by the methods M (K 14,1 ) and M(B 14 ) for the examples of FK1 when = 0.0001 and N = 51.    Table 9 shows the accuracy comparisons of the proposed approach with the known methods from the literature of which the errors in Baker et al. [44] are given in ME (maximum error) and other errors are given in RMSE for the Example 1, Example 2 and Example 3 of FK1. The data in the second row presents the results in Wen and Wei [21] for n = 51 and the error in the third row last column is from Table 1 (s = 3) given in Baker et al. [44]. The data in row 4 and row 5 are obtained by the methods M(K 5,10 ), and M(B 5 ), respectively for N = 51, while the results in row 6 , row 7 are achieved by M(K 12,10 ), M(B 12 ) accordingly also for N = 51. For the Example 4, the exact solution f ∈ C 1 [0, 1]. Hence, dealing with this test problem we provide comparisons between the methods M(K n,α ), and M(B n ) based on the regularization parameter η(δ) taken as δ and on the order n of the approximation in Figures 4 and 5, respectively. Figure 4 shows the RMSE with respect to δ obtained by the methods M(K 8,0.1 ), M(K 8,1 ), M(K 8,10 ), and M(B 8 ) for the Example 4 of FK1 when = 0.0001 and N = 51. It can be viewed that for δ ≤ 10 −14 the given approach M(K 8,1 ), M(K 8,10 ) give more accurate results then M(B 8 ). Figure 5 illustrates the RMSE with respect to n obtained by the methods M(K n,0.0001 ), M(K n,0.1 ), M(K n,1 ), M(K n,10 ), and M(B n ) for the Example 4 of FK1 when = 0.0001 and N = 51, δ = 5 × 10 −12 . This figure show that K n,1 and K n,10 give more accurate results then B n for large values of n that is for n ≥ 12.

Remark 2.
For the numerical solution of Example 6 of VK2 by the method M(K n,α ) and using the grid points x j = j n + , j = 0, 1, . . . , n − 1 and x n = 1 − , 0 < < 1 2n results in the following algebraic system of equationsÄ where coefficient matrixÄ has the entries and the vectors X and B are as in (40) and (41) Table 10 presents the RMSE with respect to n obtained by the proposed approach when α = 10, (M(K n,10 )) and = 0.001, N = 100 for the Example 5, Example 6 of VK2 and Example 7, Example 8 of VK1 . Tables 11 and 12 show the ME with respect to n obtained by the methods M(K n,10 ) and M(B n ) respectively when = 0.001, N = 100 for the considered examples of VK2 and VK1. From Tables 10-12 we conclude that the error is not improved for n = 20 for the Examples 6-8 due to the large condition numbers of the coefficient matrices. Table 13 demonstrates the RMSE with respect to α obtained by the proposed approach when n = 20, and = 0.001, N = 100 for the considered examples of VK2 and VK1. This Table shows that M(K 20,α ) gives stable solution with respect to α for the taken values of and δ. Further, in Tables 10-13 for the Example 7, x ∈ [0, 3].    The absolute errors for the Example 5 of VK2 when N = 8 (9 points) obtained by the methods M(K 20,10 ) and M(B 20 ) are presented in Table 14. Additionally, in Table 15, absolute errors obtained by the methods M(K 10,10 ) and M(B 10 ) and by the approach given in Maleknejad et al. [33] (Table 1, column 3) for the same example over the same grid points are compared when n = 10. It can be concluded from this table that the maximum error (ME) is 1.59792 × 10 −6 by the methods M(K 10,10 ), M(B 10 ) and it is 1.593 × 10 −6 by the method in Maleknejad et al. [33] and occurs at the same grid point x 7 = 0.75. Furthermore, Table 14 shows that the maximum error decreases down to 8.88623 × 10 −13 by M(K 20,10 ) and to 9.83824 × 10 −13 by M(B 20 ) over the same grid points. Table 16 shows the absolute errors (AE) at 7 points (N = 6) from the interval x ∈ [0, 3] for the Example 7 obtained by the methods M(K 15,10 ) and M(B 15 ) and by the method given in Taylor [24] (Table 2, last column). Table 17 gives AE at the points x p = p, p = 0, 1, 2, 3, 4, 5 from the interval x ∈ [0, 10] for the Example 7 obtained by the methods M(K 15,10 ) and M(B 15 ) and by the method given in Brunner [25] (Table 2, second column). We conclude from Tables 16 and 17 that the presented AE by M(K 15,10 ) are smaller than the given values from Taylor [24] and Brunner [25], respectively. However, we should remark that precision of the computations were not mentioned in both of these references. Furthermore, δ = 5 × 10 −15 for the all considered examples of VK2 and VK1 by the methods M(K n,α ) and M(B n ).     Figure 6 illustrates the condition number of the matrix A in (71) when the method M(K n,10 ) is applied for n = 2, . . . , 20. The RMSE and ME with respect to n obtained by the methods M(K n,1 ), M(K n,10 ) and M(B n ) for the Example 5, Example 6 of VK2 and Example 7 and Example 8 of VK1, when = 0.001, and N = 100 are given in Figures 7 and 8, respectively. Furthermore, for the data in Figures 6-8, the parameter δ is taken as δ = 5 × 10 −15 for the considered examples of VK2 and VK1. It can be seen from Figure 7 that for large n that is n ≥ 10, the proposed method M(K n,α ) for α = 1 and α = 10 gives more stable results than M(B n ) for the Example 6 of VK2.

Conclusions
In this paper, we gave an approach that uses Modified Bernstein-Kantorovich operators to approximate the solution of the Fredholm and Volterra integral equations of first kind. The method is developed first by representing the Modified Bernstein-Kantorovich operators such that the parameter α is also expressed explicitly in the operator. Further, the unknown function in the first kind integral equations is approximated by using the given form of the Modified Bernstein-Kantorovich operators so that the effect of α in the solution is analyzed. The obtained linear equations are transformed into system of algebraic linear equations. Furthermore, regularization technique is also applied to obtain more stable numerical solution when approximations are conducted using high-order Modified Bernstein-Kantorovich operators. The proposed approach is simple and the obtained numerical results show that the accuracy is high even when low order approximations are used, i.e., for n = 2, 3.