Solving Integro-Differential Boundary Value Problems Using Sinc-Derivative Collocation

Abstract

In this paper, the sinc-derivative collocation approach is used to solve second order integro-differential boundary value problems. While the derivative of the unknown variables is interpolated using sinc numerical methods, the desired solution and the integral terms are evaluated through numerical integration and all higher order derivatives are approximated through successive numerical differentiation. Suitable transformations are used to reduce non-homogeneous boundary conditions to homogeneous. Comparison of the proposed method with different approaches that were previously considered in the literature is carried out in order to test its accuracy and efficiency. The results show that the sinc-derivative collocation method performs well.

1. Introduction

A wide range of real world processes that are explicitly influenced by the history and the additive nature of the system are modelled using integro-differential boundary value problems (IDBVPs). These models constitute boundary value problems involving differential as well as integral terms in the governing model equations. IDBVPs have been used to model many scientific and engineering problems including, population dynamics [1,2,3,4], forestry [5], nano-structures [6], financial models [7], the theory of viscoelasticity [8], neuroscience models [9], image processing [10], geotechnical problems [11], cancer research [12] and pest control [13].

Due to their inherent complication, there are no analytic solutions of most IDBVPs. Therefore, numerical methods are often used to approximate their exact solutions including reproducing kernel Hilbert space method [14], semi-orthogonal spline wavelets [15], the Chebychev finite difference method [16,17], the Legender polynomials [18], the multistep collocation method [19], the compact finite difference method [20], the Tau method [21,22], the domain decomposition method [23], and the sinc numerical method [24].

In this work, the sinc-derivative collocation numerical technique is used to approximate the solution of a general second order integro-differential boundary value problem of the form:

$$\begin{array}{c}{y}^{″}\left(x\right)+p\left(x\right)y^{\prime }\left(x\right)+q\left(x\right)y\left(x\right)+{\lambda }_1{\int }_a^x{k}_1(x,t)y\left(t\right)dt+{\lambda }_2{\int }_{\Gamma }{k}_2(x,t)y\left(t\right)dt=g(x,y\left(x\right))\\ \mathrm{such}\mathrm{that}t,x\in \Gamma =[a,b],y\left(a\right)={\gamma }_l,y\left(b\right)={\gamma }_r\end{array}$$

(1)

where the coefficients ${\lambda }_1$, ${\lambda }_2$, the functions $p\left(x\right)$, $q\left(x\right)$, and $g(x,y(x\left)\right)$ and the integral kernels $k_1(x,t)$ and $k_2(x,t)$ are all given. Note that, $g(x,y(x\left)\right)$ is generally nonlinear in $y\left(x\right)$ which is the variable of the IDBVP that needs to be determined. It is assumed that all the functions in this equation satisfy the properties which can guarantee that the IDBVP solutions exists and is unique. We note that the method described in this paper can easily be extended to higher order IDBVPs.

Sinc numerical methods have become more prevalent in recent years as a method of solving a wide range of applications involving boundary value problems [25,26,27,28,29,30,31,32]. Partly this is due to their high efficiency in handling singular boundary value problems but also due to their ability to provide highly accurate solutions with exponentially decaying errors [33,34,35,36]. Recently, Yegneh et al. effectively utilized the sinc-collocation method to approximate the solutions of the above IDBVP [24]. By transforming the IDBVP into a system of discrete equations using the sinc method, they obtained approximate solutions of IDBVPs. In their paper, Yegneh et al. used the conventional strategy of implementing sinc method which typically interpolates the unknown variable and obtained higher derivatives via numerical differentiation. In the current approach the sinc-derivative collocation is used to interpolate the derivative of the unknown variable and obtain the unknown variable and the integral terms via sinc numerical integration [37]. The advantage of the sinc-derivative collocation approach is to average and damp out the inherent numerical errors often associated with numerical differentiation [37,38]. By improving the accuracy of the unknown derivative variable, the sinc-derivative collocation method has been shown to provide highly accurate solutions [39,40,41].

We use the sinc-derivative collocation approach to approximate the IDBVP given by Equation (1) in which variable transformation is used to reduce nonhomogeneous boundary conditions to homogeneous ones. Several illustrative examples that have been considered in earlier literatures that use comparable numerical methods are used to investigate the effectiveness of the current approach. The investigation demonstrates that our approach is effective and accurate.

Following Section 2 where a brief discussion on the preliminaries of the sinc-numerical method is presented, the sinc-derivative collocation approach for solving IDBVPs is presented in Section 3. Illustrative examples that were considered in recent literature are presented in Section 4 in order to demonstrate the performance of our approach. A brief conclusion is provided in Section 5.

2. Sinc Preliminaries

The sinc function is defined by

$$\sin \mathrm{c}\left(x\right)=\left\{\begin{array}{cc}\frac{sin\left(\pi x\right)}{\pi x}\hfill & \mathrm{if}x\ne 0\hfill \\ 1\hfill & x=0\hfill \end{array}\right.$$

(2)

For a function $f:\mathbb{R}\to \mathbb{R}$ and a mesh-size h, the infinite series

$$C(f,h)\left(x\right)=\sum _{k=-\infty }^{\infty }f\left(kh\right)S(k,h)\left(x\right)$$

(3)

where

$$S(k,h)\left(x\right)=\sin \mathrm{c}\left(\frac{x-kh}{h}\right)$$

(4)

is known as the Whittaker cardinal expansion of the function f whenever it converges [34]. For computational purposes, the sinc approximation is truncated using $2N+1$ collocation points as

$$C_N(f,h)\left(t\right)\equiv \sum _{i=-N}^{N}f\left(ih\right)S(i,h)\left(t\right)$$

(5)

where N is a positive integer. As established via Theorems 1–4 which are presented after the following definition [33], there is a class of functions for which errors associated with sinc approximations decay exponentially.

Definition 1.

$B\left(D_d\right)$ represents a the class of analytic functions f in $D_d$ which satisfy

$$I(f,D_d)\equiv \underset{ϵ\to 0}{lim}{\int }_{\partial D_d\left(ϵ\right)}|f\left(z\right)\left|\right|dz|<\infty$$

where

$$D_d=\left\{z=s+it\left|\left|t\right|<d<\frac{\pi }{2}\right.\right\},$$

$${\displaystyle D_d\left(ϵ\right)=\left\{z=s+it\left|\left|s\right|<\frac{1}{ϵ}\right.,\left|t\right|<d(1-ϵ)\right\},0<ϵ<1,}$$

and $\partial D_d$ represents the boundary of $D_d$.

For functions that decay single exponentially, we have the following theorem, due to Stenger [33].

Theorem 1.

(Stenger [33]). If $f\left(t\right)\in B\left(D_d\right)$ and for constants $\alpha >0$ and $\beta >0$

$$\left|f\left(t\right)\right|\le \alpha exp\left(-\beta \left|t\right|\right)\forall t\in \mathbb{R}$$

then:

$$\underset{-\infty \le t\le \infty }{sup}\left|f\left(t\right)-\sum _{i=-N}^{N}S(i,h)\left(t\right)f\left(ih\right)\right|\le C\sqrt{N}exp\left(-\sqrt{d\pi \beta N}\right)$$

for some constant $C>0$ and taking the mesh-size h for the sinc-collocation to be:

$$h=\sqrt{\frac{d\pi }{\beta N}}.$$

While for double exponentially decaying functions we have the following theorem due to Sugihara [42]

Theorem 2.

(Sugihara [42]). If $f\left(t\right)\in B\left(D_d\right)$ and for constants $\alpha >0$ and $\beta >0$

$$\left|f\left(t\right)\right|\le \alpha exp\left(-\beta exp\left(\gamma \right|t\left|\right)\right)\forall t\in \mathbb{R}$$

then:

$$\underset{-\infty \le t\le \infty }{sup}\left|f\left(t\right)-\sum _{i=-N}^{N}S(i,h)\left(t\right)f\left(ih\right)\right|\le Cexp\left(\frac{-\pi \gamma Nd}{log(\pi \gamma Nd/\beta )}\right)$$

for some constant $C>0$ and taking the mesh-size h for the sinc-collocation to be:

$$h=\frac{log(\pi \gamma Nd/\beta )}{N\gamma }.$$

For using these interpolations over $\Gamma$, we use a variable transform $\varphi :[a,b]\to \mathbb{R}$ and an associated inverse transform $\psi :\mathbb{R}\to [a,b]$. Combined with the variable transform, the truncated cardinal expansion of a function over an interval $[a,b]$ can be written as

$$C{(f,h)}_N\left(t\right)=\sum _{i=-N}^{N}f\left(\psi \left(ih\right)\right)(S(i,h)\circ \varphi )\left(t\right)$$

(6)

The theorems below which follow directly from Theorems 1 and 2 above quantify the error bound associated with this interpolation [33,42].

Theorem 3.

If $t=\psi \left(\xi \right)$ and $f\left(\psi \left(\xi \right)\right)\in B\left(D_d\right)$ and for constants $\alpha >0$ and $\beta >0$

$$|f\left(\psi \left(\xi \right)\right){\psi }^{\prime }\left(\xi \right)|\le \alpha exp(-\beta \left|\xi \right|\left)\right)\forall \xi \in \mathbb{R}$$

then:

$$\underset{-\infty \le t\le \infty }{sup}\left|f\left(t\right)-\sum _{i=-N}^{N}f\left(\psi \left(ih\right)\right)(S(k,h)\circ \varphi )\left(t\right)\right|\le C\sqrt{N}exp\left(-\sqrt{d\pi \beta N}\right)$$

for some constant $C>0$ in which the mesh-size h for the sinc-collocation is taken as:

$$h=\sqrt{\frac{d\pi }{\beta N}}.$$

Theorem 4.

If $t=\psi \left(\xi \right)$ and $f\left(\psi \left(\xi \right)\right)\in B\left(D_d\right)$ and for constants $\alpha >0$ and $\beta >0$

$$|f\left(\psi \left(\xi \right)\right){\psi }^{\prime }\left(\xi \right)|\le \alpha exp(-\beta exp\left(\gamma \right|\xi \left|\right))\forall \xi \in \mathbb{R}$$

then:

$$\underset{-\infty \le t\le \infty }{sup}\left|f\left(t\right)-\sum _{i=-N}^{N}f\left(\psi \left(ih\right)\right)(S(i,h)\circ \varphi )\left(t\right)\right|\le Cexp\left(\frac{-\pi \gamma Nd}{log(\pi \gamma Nd/\alpha )}\right)$$

for some constant $C>0$ and taking the mesh-size h for the sinc-collocation to be:

$$h=\frac{log(\pi \gamma Nd/\alpha )}{N\gamma }.$$

The frequently used single exponential transformations ${\varphi }_S\left(x\right)$ is given by [33]

$$\xi ={\varphi }_S\left(x\right)=log\left(\frac{x-a}{b-x}\right)$$

(7)

with the corresponding inverse:

$$x={\psi }_S\left(\xi \right)=\frac{b+a}{2}+\frac{b-a}{2}tanh\left(\xi \right).$$

(8)

Similarly the double exponential transformation ${\varphi }_D\left(x\right)$ is given by ([35])

$$\xi ={\varphi }_D\left(x\right)=log\left(G+\sqrt{G^2+1}\right),G=\frac{1}{\pi }{\varphi }_S\left(x\right)$$

(9)

with the corresponding inverse:

$$x={\psi }_D\left(\xi \right)=\frac{b+a}{2}+\frac{b-a}{2}tanh\left(\frac{\pi }{2}sinh\left(\xi \right)\right).$$

(10)

3. The Derivative Interpolation Method for IDBVPs

In order to use the sinc-derivative approach to solve the IDBVP given by (1), we approximate the derivative of the unknown variable $y\left(x\right)$ as

$$y_N^{\prime }\left(x\right)=u_N^{\prime }\left(x\right)+r^{\prime }\left(x\right)$$

(11)

where $u_N^{\prime }\left(x\right)$ is the collocation part of the solution’s derivative which is given by:

$$u_N^{\prime }\left(x\right)=\sum _{k=-N}^N{c}_k(S(k,h)\circ \varphi )\left(x\right)$$

(12)

where the $c_k$s are $2N+1$ unknown sinc coefficients that need to be determined and $r\left(x\right)$ is a quadratic polynomial defined as:

$$r\left(x\right)=c_{N+1}\left(\frac{{(b-a)}^2-{(b-x)}^2}{2(b-a)}\right)+c_{N+2}\left(\frac{{(x-a)}^2}{2(b-a)}\right)+r\left(a\right)$$

(13)

so that:

$$r^{\prime }\left(x\right)=c_{N+1}\left(\frac{b-x}{b-a}\right)+c_{N+2}\left(\frac{x-a}{b-a}\right).$$

(14)

Here, we define

$$c_{N+1}=y_N^{\prime }\left(a\right),c_{N+2}=y_N^{\prime }\left(b\right),\mathrm{and}r\left(a\right)=y_N\left(a\right)={\gamma }_l.$$

(15)

which imply that:

$$r^{\prime }\left(a\right)=c_{N+1}=y_N^{\prime }\left(a\right),r^{\prime }\left(b\right)=c_{N+2}=y_N^{\prime }\left(b\right),u_N\left(a\right)=0$$

(16)

so that the boundary conditions associated with the sinc collocation part, $u_N^{\prime }\left(x\right)$ transform into homogenous ones:

$$u_N^{\prime }\left(a\right)=u_N^{\prime }\left(b\right)=0.$$

(17)

Moreover

$$u_N\left(b\right)+r\left(b\right)=u_N\left(b\right)+c_{N+1}\frac{b-a}{2}+c_{N+2}\frac{b-a}{2}+{\gamma }_l={\gamma }_r$$

(18)

Hence the boundary condition on $u_N^{\prime }$ becomes homogeneous as desired and the right-end boundary condition (18) will be imposed when the IDBVP is discretized. This particular approach of transforming the IDBVP into a homogeneous boundary conditions has been used by other researchers [43,44].

Evaluating the expression in (11) at the collocation points given by

$$x_l=\psi \left(lh\right)l=-N,\dots ,N$$

(19)

gives

$$u_N^{\prime }\left(x_l\right)=\sum _{k=-N}^N{c}_k{\delta }_{l,k}^{\left(0\right)}$$

(20)

where

$${\delta }_{k,l}^{\left(0\right)}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}k\ne l\hfill \\ 0\hfill & \mathrm{if}k=l\hfill \end{array}\right.$$

which yields

$$u_N^{\prime }\left(x_l\right)=c_l$$

(21)

the second derivative is approximated by differentiating (11) as follows:

$$\begin{array}{cc}\hfill y_N^{″}\left(x\right)& =\sum _{k=-N}^N{c}_k\frac{d}{dx}\left(\left(S\right(k,h)\circ \varphi )\left(x\right)\right)+r^{″}\left(x\right)\hfill \\ & =\sum _{k=-N}^N{c}_k\frac{d}{d\varphi }\left(\left(S\right(k,h)\circ \varphi )\left(x\right)\right){\varphi }^{\prime }\left(x\right)+r^{″}\left(x\right).\hfill \end{array}$$

(22)

Hence at the sinc nodes $x_l$

$$y_N^{″}\left(x_l\right)=\sum _{k=-N}^N{c}_k{\delta }_{k,l}^{\left(1\right)}{\varphi }^{\prime }\left(x_l\right)+r^{″}\left(x_l\right)=u^{″}\left(x_l\right)+r^{″}\left(x_l\right)$$

(23)

where

$${\left.\frac{d}{d\varphi }((S(k,h)\circ \varphi )\left(x\right)\right|}_{x=x_l}={\delta }_{k,l}^{\left(1\right)}=\frac{1}{h}\left\{\begin{array}{cc}0\hfill & k=l\hfill \\ \frac{{(-1)}^{l-k}}{(l-k)}\hfill & k\ne l.\hfill \end{array}\right.$$

(24)

The expression for the solution $y_N\left(x\right)$ is obtained by integrating Equation (11) [37]:

$$y_N\left(x\right)={\int }_a^x{u}_N^{\prime }\left(s\right)ds+r\left(x\right).$$

(25)

Hence

$$\begin{array}{cc}\hfill y_N\left(x_l\right)& ={\int }_a^{x_l}{u}_N^{\prime }\left(s\right)ds+r\left(x_l\right)\hfill \\ & ={\int }_a^{x_l}\sum _{k=-N}^N\left(c_k(S(k,h)\circ \varphi )\left(s\right)\right)ds+r\left(x_l\right)\hfill \\ & =\sum _{k=-N}^N{c}_k{\int }_a^{x_l}\left(\left(S\right(k,h)\circ \varphi )\left(s\right)\right)ds+r\left(x_l\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\sum _{k=-N}^N\left(\frac{c_k{\delta }_{l,k}^{(-1)}}{{\varphi }^{\prime }\left(x_k\right)}\right)+r\left(x_l\right)=u\left(x_l\right)+r\left(x_l\right)\hfill \end{array}$$

(27)

where

$${\delta }_{l,k}^{(-1)}=h\left\{\begin{array}{cc}\frac{1}{2}+\frac{Si\left(\pi \right(l-k\left)\right)}{\pi }\hfill & \mathrm{if}k\ne l\hfill \\ \frac{1}{2}\hfill & \mathrm{if}k=l\hfill \end{array}\right.$$

(28)

and ${\displaystyle Si\left(z\right)}$ is given by

$$Si\left(z\right)={\int }_0^z\frac{sin\left(t\right)}{t}dt.$$

3.1. Evaluating the Integral Terms

Consider the first integral in Equation (1):

$${\int }_a^x{k}_1(x,t)u\left(t\right)dt$$

Separate this integral into the parts $k_1(x,t)$ and $u\left(t\right)$ and perform an integration by parts to get

$${\int }_a^x{k}_1(x,t)u\left(t\right)dt=M_1(x,t)u\left(t\right){|}_a^x-{\int }_a^x{M}_1(x,t)u^{\prime }\left(t\right)dt$$

(29)

where ${\displaystyle M_1(x,t)}$ is the antiderivative of ${\displaystyle k_1(x,t)}$ with respect to the variable t. Therefore, due to the homogeneous boundary conditions we get:

$${\int }_a^x{k}_1(x,t)u\left(t\right)dt=M_1(x,x)u\left(x\right)-{\int }_a^x{M}_1(x,t)u^{\prime }\left(t\right)dt$$

(30)

Evaluating this expression at the collocation points $x_l$ and using (27) we obtain:

$$\begin{array}{cc}\hfill {\int }_a^{x_l}{k}_1(x_l,t)u\left(t\right)dt=& M_1(x_l,x_l)u\left(x_l\right)-{\int }_a^{x_l}{M}_1(x_l,t)u^{\prime }\left(t\right)dt\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill & \hfill =& M_1(x_l,x_l)\sum _{k=-N}^N\left(\frac{c_k{\delta }_{l,k}^{(-1)}}{{\varphi }^{\prime }\left(x_k\right)}\right)+\sum _{k=-N}^N\left(\frac{c_k{M}_1(x_l,x_k){\delta }_{l,k}^{(-1)}}{{\varphi }^{\prime }\left(x_k\right)}\right).\hfill \end{array}$$

(32)

Similarly, the second integral term evaluated at the collocation points $x_l$ becomes:

$$\begin{array}{cc}\hfill {\int }_a^b{k}_2(x_l,t)u\left(t\right)dt=& M_2(x_l,b)u\left(b\right)-{\int }_a^b{M}_2(x_l,t)u^{\prime }\left(t\right)dt\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill & \hfill =& M_2(x_l,b)\sum _{k=-N}^N\left(\frac{c_k}{{\varphi }^{\prime }\left(x_k\right)}\right)+\sum _{k=-N}^N\left(\frac{c_k{M}_2(x_l,x_k)}{{\varphi }^{\prime }\left(x_k\right)}\right)\hfill \end{array}$$

(34)

where ${\displaystyle M_2(x,t)}$ is the antiderivative of ${\displaystyle k_2(x,t)}$ with respect to the variable t.

3.2. Discretizing the IDBVP

Applying the decomposition given by (11) the IDBVP of Equation (1) become:

$$\begin{array}{c}{u}^{″}\left(x\right)+p\left(x\right)u^{\prime }\left(x\right)+q\left(x\right)u\left(x\right)+{\lambda }_1{\int }_a^x{k}_1(x,t)u\left(t\right)dt+{\lambda }_2{\int }_{\Gamma }{k}_2(x,t)u\left(t\right)dt=\widehat{g}(x,u\left(x\right))\\ \mathrm{such}\mathrm{that}x,t\in \Gamma =[a,b],u\left(a\right)=u\left(b\right)=0\end{array}$$

(35)

where

$$\widehat{g}(x,u\left(x\right))=g(x,u\left(x\right)+r\left(x\right))-r^{″}\left(x\right)-p\left(x\right)r^{\prime }\left(x\right)-q\left(x\right)r\left(x\right)-{\lambda }_1{\int }_a^x{k}_1(x,t)r\left(t\right)dt-{\lambda }_2{\int }_{\Gamma }{k}_2(x,t)r\left(t\right)dt$$

(36)

Hence, the discretized version of the IDBVP of Equation (1) at the $2N+1$ sinc nodes, ${\displaystyle x_l,l=-N,-N+1,\dots ,N-1,N,}$ become:

$$\sum _{k=-N}^N{A}_{l,k}{c}_k=F_l=\widehat{g}(x_l,u\left(x_l\right))$$

(37)

where

$$A_{l,k}={\delta }_{k,l}^{\left(1\right)}{\varphi }^{\prime }\left(x_l\right)+p\left(x_l\right)+\frac{B_{l,k}}{{\varphi }^{\prime }\left(x_k\right)},k=-N,-N+1,\dots ,N-1,N,A_{l,N+1}=A_{l,N+2}=0$$

(38)

$$B_{l,k}={\delta }_{l,k}^{(-1)}\left[q\left(x_l\right)+{\lambda }_1\left(M_1(x_l,x_l)+M_1(x_l,x_k)\right)\right]+{\lambda }_2\left(M_2(x_l,x_l)+M_2(x_l,x_k)\right).$$

The boundary conditions (18) become:

$$\sum _{k=-N}^N{A}_{N+1,k}{c}_k+A_{N+1,N+1}{c}_{N+1}+A_{N+1,N+2}{c}_{N+2}=F_{N+1}={\gamma }_r-{\gamma }_l$$

(39)

where

$$A_{N+1,k}=\frac{1}{{\varphi }^{\prime }\left(x_k\right)},A_{N+1,N+2}=A_{N+1,N+2}=\frac{b-a}{2}$$

(40)

Note that Equations (37) and (39) entail a system of $2N+2$ equations for $2N+3$ unknowns ${\displaystyle c_k,k=-N,-N+1,\dots ,N+1,N+2}$. Therefore, we need an additional equation in order to close the system and solve for the $n=2N+3$ unknowns. This is done by applying a discretization at an additional sinc point, $x_L$, $L=-N-1$ which yields:

$$\sum _{k=-N}^N{A}_{L,k}{c}_k=F_L=\widehat{g}(x_L,u\left(x_L\right))$$

(41)

$$A_{L,k}={\delta }_{k,L}^{\left(1\right)}{\varphi }^{\prime }\left(x_L\right)+\frac{B_{L,k}}{{\varphi }^{\prime }\left(x_k\right)},k=-N,-N+1,\dots ,N-1,N,A_{L,N+1}=A_{L,N+2}=0$$

(42)

$$B_{L,k}={\delta }_{L,k}^{(-1)}\left[q\left(x_L\right)+{\lambda }_1\left(M_1(x_L,x_L)+M_1(x_L,x_L)\right)\right]+{\lambda }_2\left(M_2(x_L,x_L)+M_2(x_L,x_k)\right).$$

Define the $n\times 1$ vector $\mathbf{C}$ by:

$$\begin{array}{ccc}\hfill C& =& {[c_{-N},\dots c_0,\dots ,c_N,c_{N+1},c_{N+2}]}^T\hfill \end{array}$$

(43)

$$\begin{array}{c}={[u_N^{\prime }\left(x_{-N}\right)\dots u_N^{\prime }\left(x_0\right),\dots u_N^{\prime }\left(x_N\right),y_N^{\prime }\left(a\right),y_N^{\prime }\left(b\right)]}^T.\hfill \end{array}$$

(44)

Therefore, Equations (37), (39) and (41) consist of n equations for the n unknowns which are given by the matrix equation

$$AC=F$$

(45)

where the matrix entries of the matrix A and the vector function F are as described above. Once Equation (45) is solved, the coefficients are used to compute the unknown variable $y_N\left(x\right)$ at the sinc nodes using Equation (27). Note that the values of $y_N^{\prime }\left(a\right)$ and $y_N^{\prime }\left(b\right)$ are also directly determined.

Note that, in the special case where the function g in the IDBVP is independent of $y\left(x\right)$, Equation (45) represents a linear matrix equation which can be solved using standard methods. However, in general $g(x,y(x\left)\right)$ is a nonlinear function of $y\left(x\right)$ and therefore Equation (45) is a nonlinear system which must be solved using Newton’s method or other nonlinear solvers.

4. Numerical Illustrations

In order to demonstrate the enhanced rate of convergence and accuracy of the sinc-derivative collocation method with respect to other approaches, we solve the following problems considered in recent literatures. All the computations of the sinc-derivative results are performed using Mathematica 11.3.

Example 1.

The first equation we attempt to solve is an integro-differential equation considered in [24,45].

$$\begin{array}{c}{y}^{″}\left(x\right)+\frac{1}{\sqrt{x}}{y}^{\prime }\left(x\right)+\frac{1}{x}y\left(x\right)+{\int }_0^x(t+x)y\left(t\right)dt+{\int }_0^{1}txy\left(t\right)dt\end{array}$$

$$\begin{array}{c}-\frac{1}{1+sin\left(y^2\left(x\right)\right)}-e^{y^9\left(x\right)}+y^{11}\left(x\right)=f\left(x\right),x\in [0,1]\end{array}$$

(46)

$$\begin{array}{c}\mathit{subject}\mathit{to}y\left(0\right)=0,y\left(1\right)=0\end{array}$$

(47)

where

$$\begin{array}{c}f\left(x\right)=\frac{1}{\sqrt{x}}((2\sqrt{x}+x-4x\sqrt{x}-x^2)cos\left(x\right)\\ +(1-\sqrt{x}-2x-2x\sqrt{x}+x^2\sqrt{x})sin\left(x\right))+2(1-4x-x^2+x^3)cos\left(x\right)\\ +(6+3x-5x^2)sin\left(x\right)-(2+4cos\left(1\right)-5sin\left(1\right))x-\frac{1}{1+sin\left({\left((x-x^2)sin\left(x\right)\right)}^2\right)}\\ -e^{{\left((x-x^2)sin\left(x\right)\right)}^9}+{\left((x-x^2)sin(x)\right)}^{11}+2(x-1).\end{array}$$

(48)

The solution to this equation is $y\left(x\right)=(x-x^2)sin\left(x\right).$

Comparison of the current approach with other approaches is depicted in

Table 1. Comparison of absolute error in $y\left(x\right)$ for Example 1.

x	Current Method	Method of [45]	Method of [24]
$0.08$	$1.23476\times {10}^{-10}$	$4.6259\times {10}^{-5}$	$2.8293\times {10}^{-6}$
$0.16$	$1.15291\times {10}^{-10}$	$4.81776\times {10}^{-5}$	$8.9730\times {10}^{-6}$
$0.32$	$1.67924\times {10}^{-10}$	$4.16774\times {10}^{-5}$	$7.3063\times {10}^{-6}$
$0.48$	$1.28785\times {10}^{-10}$	$3.48007\times {10}^{-5}$	$1.2008\times {10}^{-5}$
$0.64$	$5.02378\times {10}^{-11}$	$2.96533\times {10}^{-5}$	$8.7603\times {10}^{-6}$
$0.80$	$1.55245\times {10}^{-10}$	$2.69251\times {10}^{-5}$	$4.9231\times {10}^{-6}$
$0.96$	$9.52025\times {10}^{-11}$	$2.39306\times {10}^{-5}$	$1.0128\times {10}^{-6}$

, which represents the absolute error of $y\left(x\right)$ for $N=25$. The result is compared with that of [45] which approximates solutions in the reproducing kernel space and with [24] which uses sinc-collocation method. As demonstrated in the table, the current method is highly accurate in the entire domain of the solution. The plot of the absolute errors $y\left(x\right)$ using the sinc-derivative method described in Section 3 is depicted in Figure 1 for various values of N. The result is highly accurate with a maximum absolute errors of the solution $4.7\times {10}^{-15}$ for $N=40$. The plot of the absolute errors of $y^{\prime }\left(x\right)$ using the sinc-derivative method is depicted in Figure 2 for various values of N.

The logarithm plot of the maximum error as a function of the number of sinc nodes N using the current method and that of [24] is shown in Figure 3. For the result of [24], only the portion of the graph that was reported in their paper is displayed. The plot shows the exponential decrease of the errors with respect to N and confirms that the sinc-derivative approach is a highly accurate method for solving IDBVPs.

Example 2.

The second IDBVP we attempt to solve is one that is considered in [46].

$$\begin{array}{c}{y}^{″}\left(x\right)-2y\left(x\right)-{\int }_{-1}^{1}t{e}^{-t}cos\left(x\right)y\left(t\right)dt=x,x\in [0,1]\\ \mathit{subject}\mathit{to}y(-1)=-1,y\left(1\right)=1\end{array}$$

(49)

The solution to this equation is $y\left(x\right)=c_1{e}^{\sqrt{2}x}+c_2{e}^{-\sqrt{2}x}+c_{3}cos\left(x\right)-\frac{1}{2}x$, where $c_1=0.4206057265...$, $c_2=-0.3545613032...$ and $c_3=-0.2662525683...$.

The plot of the absolute errors of the sinc-derivative solution for Example 2. is depicted in Figure 4. for $N=25$ and $N=40$. The maximum absolute errors are $4.4\times {10}^{-10}$ for $N=25$ and $2.9\times {10}^{-14}$ for $N=40$. The best result presented in [46] has a maximum absolute error of $1.4\times {10}^{-6}$. Hence the present method performs extremely well.

Example 3.

The third IDBVP we attempt to solve is one that is considered in [24,47].

$$\begin{array}{c}{y}^{″}\left(x\right)+\frac{1}{1+y^2\left(x\right)}+x{e}^{y\left(x\right)}+{\int }_0^{x}xty\left(t\right)dt=f\left(x\right),x\in [0,1]\\ \mathit{subject}\mathit{to}y\left(0\right)=1,y\left(1\right)=2,\end{array}$$

(50)

where

$$f\left(x\right)=2+\frac{1}{1+{(1+x^2)}^2}+x{e}^{1+x^2}+\frac{x^3}{4}(2+x^2)$$

(51)

The solution to this equation is $y\left(x\right)=x^2+1$

The maximum absolute errors reported for this problem in [24] are $1.9\times {10}^{-5}$ for $N=15$ and $8.1\times {10}^{-8}$ for $N=30$. Since Equation (13) gives $r\left(x\right)=x^2+1$ for this IDBVP which is identical to the exact solution $y\left(x\right)$, the current numerical scheme for this problem naturally leads to the exact solution.

Finally, we note that the current numerical method resulted in a similar level of accuracies for many other examples considered in the literature further demonstrating its excellent performance.

5. Conclusions

In this paper, the sinc-derivative collocation method was used to approximate the solution of second order nonlinear integro-differential boundary value problems. In the sinc-derivative approach, the unknown variable derivative is interpolated via sinc numerical methods and the desired solution is obtained through numerical integration. Non-homogeneous boundary conditions are converted to homogeneous ones via suitable transformations. The efficiency as well as the accuracy of the method is demonstrated using illustrative examples which were recently considered using other approaches. The results demonstrate the excellent performance of the sinc-derivative interpolation method for solving integro-differential boundary value problems.