Most of the problems of science and technology can be treated with the aid of theories of ordinary and partial differential equations (PDEs). However, there are better methods called the theory of integral equations to solve these problems. Years ago, these types of equations were a hot topic in the minds of mathematicians. If the unknown function is given inside the integration symbol in an equation, it is named an integral equation, which is regarded as a common type of functional equation. The theory of integral equations is one of the most operative mathematical tools for solving problems in pure and applied mathematics, mathematical physics, mechanical vibrations, and fields related to science and technology. However, the real development of the theory of integral equations started with the attempts of the Italian mathematician Vito Volterra and the Swedish mathematician Fredholm. Volterra was the first person who realized the significance of the theory of integral equations and systematically paid attention to it. In 1896, he presented the first general scheme for solving a category of linear integral equations characterized by a variable that appears at the upper bound of the integral. Interest in VIEs has been growing in recent years. These equations arise in many physical applications, for example see [1,2,3,4,5,6,7,8,9,10,11,12]. Here, a numerical scheme based on Bernoulli polynomials is developed to solve VIEs of the second kind.

Bernoulli polynomials have a significant role in many areas of mathematical analysis, like the theory of modular forms [14], the theory of distributions in p-adic analysis [13], the polynomial expansions of analytic functions [15], and so on. Some applications of these polynomials in mathematical physics are related with the theory of the KdV equation [16], solving Lamé equation [17], and in the field of vertex algebra [18].

The collocation method as a numerical method for solving all kinds of functional equations and real world problems has always been the interest and attention of mathematicians. Many authors have presented different kinds of collocation method during past decades. For example, Doha et al. [19] proposed a Jacobi–Gauss–Lobatto collocation scheme, used in relationship with the fourth order implicit Runge–Kutta method as a numerical algorithm for approximating solutions of nonlinear Schrödinger equations (NLSE). Nemati [20] solved Volterra–Fredholm integral equations applying Legendre collocation method. His method is based upon shifted Legendre polynomials approximation. Then, utilizing the shifted Gauss–Legendre as collocation nodes reduces the Volterra–Fredholm integral equations to the solution of a matrix equation. Mirzaee and Hoseini [21] proposed a new matrix method on the basis of Fibonacci polynomials and collocation points for numerically solving the Volterra–Fredholm integral equations. Ren and Tian [22] proposed a scheme to solve a boundary value problem for Kirchhoff type of nonlinear integro-differential equation numerically. Gouyandeh et al. [23] by using the Tau-Collocation method approximated solution of the nonlinear Volterra-Fredholm-Hammerstein integral equations. Aziz et al. [24], based on Haar wavelet, provided a new collocation method for numerical solution of 3D elliptic PDEs with Dirichlet boundary conditions. Çelik [25], by utilizing Chebyshev wavelet collocation method, studied free vibration problems of non-uniform Euler–Bernoulli beam under different supporting conditions. Samadyar and Mirzaee [26] presented orthonormal BCM to approximate the linear singular stochastic Itô-VIEs. Biçer and Yalçinbaş [27] an approximate solution of the telegraph equation applying BCM. Alijani et al. [28] investigated systems of fuzzy fractional differential equations numerically with a lateral type of the Hukuhara derivative and its generalization. Wang et al. [29] presented a new collocation method for evaluating the 2D elliptic PDEs. Singh [30] used Jacobi collocation method to solve the fractional advection-dispersion equation arising in porous media. Kumbinarasaiah et al. [31] presented an integration operational matrix applying the Bernoulli wavelet and suggested a new scheme named the Bernoulli wavelet collocation method (BWCM). In [32], the authors investigated a space-time Sinc-collocation method for treating the fourth-order nonlocal heat model appearing in viscoelasticity. Laib et al. [33], based on the using Taylor polynomials, suggested an algorithm to construct a collocation solution for approximating the solution of 2D-VIEs. Wang et al. [34], by utilizing the zeros of Chebyshev polynomial as collocation points, proposed a new collocation method to solve the second kind VIE.

In this work, based on Bernoulli polynomials, a new collocation method is suggested to approximate numerically VIEs of the second kind. The main goal of the suggested scheme is to improve the results of the classic BCM by dividing the interval into some sub-intervals and considering the collocation points on each of them. Here, the zeros of the SCPs are considered as collocation points. Then, BCM is applied step by step from the first sub-interval to the last one. By this process, a system of algebraic equations is attained for each sub-interval that could be easily solved using computing software. At last, the approximate solution is obtained as a piece-wise function. This idea is very effective. Although, we have tested this idea on BCM in this work but we guess that all collocation methods mentioned above can be improved by this idea.

The rest of the paper is organized as follows. In section 2, Bernoulli polynomials are introduced and their features are stated. In section 3, the proposed scheme is explained. The convergence analysis is discussed in in section 4. In section 5, numerical results are presented. At last, conclusions are given in section 5.

The traditional Bernoulli polynomials $B_n\left(t\right)$ is often characterized by the following exponential generating functions [35]:

The Bernoulli polynomials of $n$th degree are determined in the interval $\left[\mathrm{0,1}\right]$ as follows [36] where $B_z=B_z\left(0\right)$ is the Bernoulli number for each $k=\mathrm{0,1},\dots N$. We write some polynomials like Bernoulli polynomials have following important properties [37].

In this section, the proposed scheme is described to deal with the VIEs of the second kind. The main purpose of the suggested method is to improve the results of the classic collocation method by dividing the interval into some sub-intervals and applying the collocation method in each of them.

Here, the convergence of the suggested scheme is discussed. For this purpose, we first recall one of the important theorems related to the residual interpolation error by Chebyshev nodes.

Theorem 1. Let $u$ be a sufficiently smooth function on $I=[0, 1]$ and ${\mathsf{\Pi }}_M$ be the space of polynomials of order $M.$ Also, let $P_M\in {\mathsf{\Pi }}_M$ be the interpolating polynomials of $u$ at points $c_1,\dots ,c_{M+1}$ which are the zeros of the SCP of degree $M+1$ on $I$. Then, the following relation is established. where $\xi \in I$.

Proof. [38].

According to the last theorem, we can write

Now, Assume that

Applying this upper bound to Eq. (25) and considering the approximations for Chebyshev interpolation nodes [39] leads to Theorem 2. Suppose that $\tilde{u}$ defined in Eq. (3), be the best approximation of real sufficiently smooth function $u$ by Bernoulli polynomials. Then a real constant $\eta$ exists such that Proof. According to the definition, $\tilde{u}$ is the best approximation of $u$ provided that

Particularly, if $v\left(t\right)=P_M\left(t\right)$ then according to Eq. (27), we get

Hence, Eq. (28) is proved.

According to Eq. (28), it can be written

So, $\frac{1}{(M+1)!2^{2M+1}}\to 0$ when $M\to \mathrm{\infty }$ which implies that $\tilde{u}\to u$. Therefore, the collocation method based on the Bernoulli polynomials is convergent.

Then, there is an upper error bound as follows. where Proof. The approximate solution of Eq. (6) in the interval $I_1$ could be stated as

From Eqs. (6) and (35), we get

$u\left(x\right)-u_{M,1}\left(x\right)={\int }_0^x k\left(x,t\right)\left(N\left(u\left(t\right)\right)-N\left(u_{M,1}\left(t\right)\right)\right)dy.$Then, we have

$∥u\left(x\right)-u_{M,1}\left(x\right)∥\le L_1∥k\left(x,t\right)∥∥u\left(t\right)-u_{M,1}\left(t\right)∥$By using Eq. (34), for $i=1$, we have

$∥u\left(x\right)-u_{M,1}\left(x\right)∥\le L_1{\lambda }_1∥u\left(x\right)-u_{M,1}\left(x\right)∥.$So,

Now, regard the approximate solution of Eq. (6) in the interval $I_2$as follows.

From Eqs. (6) and (36), we get

$\begin{array}{cc}\hfill u\left(x\right)-u_{M,2}\left(x\right)& \\ & ={\int }_0^h k\left(x,t\right) \left(N\left(u\left(t\right)\right)-N\left(u_{M,1}\left(t\right)\right)\right)dt\hfill \\ & +{\int }_h^{2h} k\left(x,t\right) \left(N\left(u\left(t\right)\right)-N\left(u_{M,2}\left(t\right)\right)\right)dt.\hfill \end{array}$Then, we have

$∥u\left(x\right)-u_{M,2}\left(x\right)∥\le ∥k\left(x,t\right)∥∥u\left(x\right)-u_{M,1}\left(x\right)∥+∥k\left(x,t\right)∥∥u\left(x\right)-u_{M,2}\left(x\right)∥$Using Eq. (34), for $i=2$ and also Eq. (36) leads to