Hybridization of Block-Pulse and Taylor Polynomials for Approximating 2D Fractional Volterra Integral Equations

Abstract

This paper proposes an accurate numerical approach for computing the solution of two-dimensional fractional Volterra integral equations. The operational matrices of fractional integration based on the Hybridization of block-pulse and Taylor polynomials are implemented to transform these equations into a system of linear algebraic equations. The error analysis of the proposed method is examined in detail. Numerical results highlight the robustness and accuracy of the proposed strategy.

1. Introduction

Fractional calculus (FC) generalizes the classical differential calculus [1,2]. The FC has been applied to model several anomalous phenomena having nonlocal dynamics and involving long memory [3,4]. Indeed, the models based on classical calculus often fail to explain the genetic and inheritance properties of many complex systems having anomalous dynamics, whereas the fractional derivatives and integrals allow for a simpler and more accurate representation of their features. Due to the above-mentioned properties, fractional order integral equations (FOIEs) have been used in many engineering and physics fields to investigate complex systems, for instance, viscoelasticity and traffic models, temperature and motor control, and solid mechanics [4,5,6,7,8]. However, FOIEs pose difficulties in achieving their accurate analytical solution and numerical techniques have to be used in order to derive useful result [9,10,11,12,13,14,15,16,17,18,19].

In this work, we study a robust computational scheme based on two-dimensional (2D) Hybridization of block pulse and Taylor polynomials (2D-HBTs) for computing the approximate solution of 2D fractional Volterra integral Equations (2DFVIEs) as

$$f(x,y)-\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y\frac{f(s,t)k(x,y,s,t)}{{(y-t)}^{1-{\lambda }_2}{(x-s)}^{1-{\lambda }_1}}\mathrm{d}t\mathrm{d}s=g(x,y),(x,y)\in \mathrm{\Omega },$$

(1)

where $({\lambda }_1,{\lambda }_2)\in (0,+\infty )\times (0,+\infty )$, $f(s,t)$ is an unknown function to be calculated, $\mathrm{\Omega }=[0,1]\times [0,1]$ denotes the spatial domain, and $g(x,y)$ and $k(x,y,s,t)$ represent prescribed smooth functions. Over past years, for approximating 2D fractional integral Equations (2DFIEs) and 2D fractional integro-differential Equations (2DFIDEs), different basic functions have been employed. For example, Najafalizadeh and Ezzati [20] used 2D block pulse functions, while Jabari et al. [21] studied 2D orthogonal triangular functions for solving the 2DFIEs. Maleknejad et al. [22] used the Hybrid functions to simulate general nonlinear 2DFIDEs. Maleknejad and Hoseingholipour [23] implemented Laguerre functions for singular integral equation in unbounded domain. Kumar and Gupta [24] analyzed an operational matrix based 2D fractional-order Lagrange polynomials for approximating nonlinear 2DFIDEs. Mirzaee and Samadyar [25] obtained an operational matrix based on 2D hat basis functions for stochastic 2DFIEs. Asgari and Ezzati [26] employed the Bernstein polynomials and Rashidinia et al. [27] implemented shifted Jacobi polynomials for approximating 2DFIEs. Heydari et al. [28] presented an iterative multistep kernel based scheme for solving the 2DFIEs. Ardabili and Talaei [29] adopted a new Chelyshkov polynomials collocation method to solve the 2DFIEs. Hesameddini and Shahbazi [30] applied 2D shifted Legendre polynomials to estimate 2DFIEs. Asgari et al. [31] adopted the Bernstein polynomials to approximate the 2DFVIEs. Abdollahi et al. [32] presented an operational matrix scheme based 2D Haar wavelets, whereas Wang et al. [33] used 2D Euler polynomials combined with Gauss-Jacobi quadrature technique to simulate the 2DFVIEs. Liu et al. [34] employed the Bivariate barycentric rational interpolation for the 2DFVIEs. Khan et al. [35] implemented 2D Bernstein’s approximation to approximate the 2DFVIEs. Wang et al. [33] applied the modified block-by-block technique, while Mohammad et al. [36] proposed an efficient approach based on Framelets for solving the 2DFVIEs. Ahsan et al. [37] used optimal Homotopy asymptotic scheme and Fazeli et al. [38] considered the Chebyshev polynomials for approximating the 2DFVIEs. Laib et al. [39] applied a numerical approach based on Taylor polynomials for the 2DFVIEs.

The main motivation of this paper is to propose an accurate numerical approach for finding the approximate solution of 2DFVIEs. The operational matrices of fractional integration based on the Hybridization of block-pulse and Taylor polynomials (2D-HBTs) are adopted to transform 2DFVIEs into a system of linear algebraic equations. The error analysis of the proposed method is examined in detail. Numerical results show the robustness and accuracy of the proposed method.

This paper includes seven sections as follows. Section 2 introduces the notation and basic definitions of FC. Section 3 derives an operational matrix of fractional integration of 2D-HBTs. Section 4 approximates the 2DFVIEs by employing the obtained operational matrix of fractional integration. Section 5 performs convergence analysis of the proposed strategy. Section 6 highlights the efficiency and accuracy of the proposed method by means of the results of numerical experiments. Finally, Section 7 contains the concluding remarks.

2. Preliminaries and Notations

In this section, we provide the basic concepts and definitions [1,2] needed in the follow-up.

Definition 1.

The Riemann-Liouville fractional integration (RLFI) of order $\alpha \ge 0$ of a function $f\left(t\right)\in L^1\left(I\right)$ is defined by

$$I^{\alpha }f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-\tau )}^{(\alpha -1)}f\left(\tau \right)\mathrm{d}\tau \alpha >0,$$

(2)

in which $\mathrm{\Gamma }\left(\theta \right)$ represents the Euler’s gamma function described as $\mathrm{\Gamma }\left(\theta \right)={\int }_0^{+\infty }{s}^{\theta -1}{e}^{-s}\mathrm{d}s$.

The RLFI has the following properties:

• $I^{0}f\left(x\right)=f\left(x\right)$,
• $I^{\alpha }{I}^{\beta }f\left(x\right)=I^{\beta }{I}^{\alpha }f\left(x\right)=I^{\alpha +\beta }f\left(x\right),$
• $I^{\alpha }{(x-a)}^{\beta }=\frac{\mathrm{\Gamma }(\beta +1)}{\mathrm{\Gamma }(\alpha +\beta +1)}{(x-a)}^{\alpha +\beta }.$

Definition 2.

The left-sided mixed RLFI of order r of f can be represented as [1]

$$I_{\theta }^{r}f(x,y)=\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(x-s)}^{{\lambda }_1-1}{(y-t)}^{{\lambda }_2-1}f(s,t)\mathrm{d}t\mathrm{d}s,$$

(3)

where $\theta =(0,0)$ and $r=({\lambda }_1,{\lambda }_2)\in (0,\infty )\times (0,\infty )$.

Here, we list some notations of the left-sided mixed RLFI [1,2] as follows:

• $I_{\theta }^{\alpha }f(x,y)={\int }_0^x{\int }_0^y(f(s,t)\mathrm{d}t\mathrm{d}s,\alpha =(1,1),(x,y)\in \mathrm{\Omega },$
• $I_{\theta }^r{x}^p{y}^q=\frac{\mathrm{\Gamma }(p+1)\mathrm{\Gamma }(q+1)}{\mathrm{\Gamma }(p+1+{\lambda }_1)\mathrm{\Gamma }(q+1+{\lambda }_2)}{x}^{p+{\lambda }_1}{y}^{q+{\lambda }_2},-1<p,q<\infty ,$
• $I_{\theta }^{\theta }f(x,y)=f(x,y).$

3. Hybrid Functions

3.1. The 2D-HBTs

First, we introduce the 1D-HBTs, $h_{i,j}\left(x\right),1\le i\le N,0\le j\le M,$ on the interval $[0,1]$ as [40,41]

$$h_{i,j}\left(x\right)=\left\{\begin{array}{cc}{T}_j(Nx-i+1)\hfill & \frac{i-1}{N}\le x\le \frac{i}{N},\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(4)

in which j and i represent the order of the Taylor polynomials and block-pulse functions, respectively, and $T_j\left(x\right)=x^j$.

We can approximate a function $f\left(x\right)\in L^2[0,1]$ in the form of 1D-HBT by

$$f\left(x\right)\simeq \sum _{i=1}^N\sum _{j=0}^{M-1}{f}_{ij}{h}_{ij}\left(x\right)=C^{T}H\left(x\right),$$

(5)

in which

$$C={[f_{10},\dots ,f_{1(M-1)},\dots ,f_{N0},\dots ,f_{N(M-1)}]}^T,$$

(6)

and

$$H\left(x\right)={[h_{10}\left(x\right),\dots ,h_{1(M-1)\left(x\right)},h_{20}\left(x\right),\dots ,h_{2(M-1)}\left(x\right),\dots ,h_{N0}\left(x\right),\dots ,h_{N(M-1)}\left(x\right)]}^T.$$

(7)

Obviously, we can obtain the hybrid coefficients $f_{ij}$ computed by

$$f_{ij}=\frac{\langle f\left(x\right),h_{ij}\left(x\right)\rangle }{\langle h_{ij}\left(x\right),h_{ij}\left(x\right)\rangle }=\frac{1}{j!N!}\left(\frac{d^{j}f\left(x\right)}{d{x}^j}\right){|}_{x=\frac{j-1}{M}},$$

(8)

where $\langle ·,·\rangle$ denotes the inner product.

Orthogonal 2D-HBTs functions $h_{i_1{j}_1{i}_2{j}_2}(x,y)$, $0\le j_1,j_2\le M,1\le i_1,i_2\le N,$ on the region $\mathrm{\Omega }=[0,1]\times [0,1]$, are defined as

$$h_{i_1{i}_2{j}_1{j}_2}(x,y)=\left\{\begin{array}{c}{T}_{j_1}(Nx-i_1+1)T_{j_2}(Ny-i_2+1)(x,y)\in [\frac{i_1-1}{N},\frac{i_1}{N})\times [\frac{i_2-1}{N},\frac{i_2}{N})\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Let S be a set of 2D-HBTs as follows:

$$\begin{array}{c}\hfill S=\{h_{1010}(x,y),\dots ,h_{101(M-1)}(x,y),h_{1020}(x,y),\dots ,\\ \hfill h_{102(M-1)}(x,y),\dots ,h_{N(M-1)N0}(x,y),\dots ,h_{N(M-1)N(M-1)}(x,y)\}.\end{array}$$

Since S is a finite dimensional subspace of $L^2\left(\Omega \right)$ for an arbitrary $f(x,y)\in L^2\left(\mathrm{\Omega }\right)$, it has the unique best approximation outside of S, therefore, there exist unique coefficients $f_{i_1{j}_1{i}_2{j}_2}$, $0\le j_1,j_2\le M-1,$$1\le i_1,i_2\le N,$ so that

$$f(x,y)=\sum _{i_1=1}^N\sum _{j_1=0}^{M-1}\sum _{i_2=1}^N\sum _{j_2=0}^{M-1}{f}_{i_1{j}_1{i}_2{j}_2}{h}_{i_1{i}_2{j}_1{j}_2}(x,y)=F^{T}H(x,y),$$

(9)

in which

$$\begin{array}{cc}\hfill F& ={[f_{1010},\dots ,f_{101(M-1)},f_{1020},\dots ,f_{102(M-1)},\dots ,f_{N(M-1)N0},\dots ,f_{N(M-1)N(M-1)}]}^T,\hfill \end{array}$$

(10)

and

$$\begin{array}{cc}H(x,y)\hfill & =[h_{1010}(x,y),\dots ,h_{101(M-1)}(x,y),h_{1020}(x,y),\dots ,\hfill \\ \hfill & h_{102(M-1)}(x,y),\dots ,h_{N(M-1)N0}(x,y),\dots ,h_{N(M-1)N(M-1)}{(x,y)]}^T=H\left(x\right)\otimes H\left(y\right),\hfill \end{array}$$

(11)

in which the superscript T is transposition and ⊗ denotes the Kronecker product. Clearly 2D-HBTs coefficients, $f_{i_1{j}_1{i}_2{j}_2}$, can be determined by

$$\begin{array}{cc}\hfill f_{i_1{j}_1{i}_2{j}_2}& =\frac{\langle h_{i_1{j}_1}\left(x\right),\langle f(x,y),h_{i_2{j}_2}\left(y\right)\rangle \rangle }{\langle h_{i_1{j}_1}\left(x\right),h_{i_1{j}_1}\left(x\right)\rangle .\langle h_{i_2{j}_2}\left(y\right),h_{i_2{j}_2}\left(y\right)\rangle }\hfill \\ \hfill & =\frac{1}{N^{i_1+j_2}{j}_1!j_2!}\left(\frac{{\partial }^{i_1+j_2}f(x,y)}{\partial x^{j_1}\partial y^{j_2}}\right){|}_{(x,y)=(\frac{i_1}{N},\frac{i_2}{N})}.\hfill \end{array}$$

Similarly, we expand the functions $k(x,y,s,t)\in L^2(\mathrm{\Omega }\times \mathrm{\Omega })$ in terms of the 2D-HBTs in the following form:

$$K(x,y,s,t)=H^T(x,y)KH(s,t),$$

(12)

in which K represents a ${\left(MN\right)}^2\times {\left(MN\right)}^2$ matrix:

$$K={\displaystyle \left(\begin{array}{cccccc}{K}_{0000}& \dots & K_{000(MN-1)}& K_{0010}& \dots & K_{00(MN-1)(MN-1)}\\ K_{0100}& \dots & K_{010(MN-1)}& K_{0110}& \dots & K_{01(MN-1)(MN-1)}\\ ⋮& ⋮& .& .& .& ⋮\\ K_{(MN-1)(MN-1)00}& .& .& .& .& K_{(MN-1)(MN-1)(MN-1)(MN-1)}\end{array}\right)},$$

in which

$$K_{i_1{j}_1{i}_2{j}_2}=\frac{1}{N^{z+u+v+w}z!u!v!w!}\left(\frac{{\partial }^{z+u+v+w}K(x,y,s,t)}{\partial x^z\partial y^u\partial s^v\partial t^w}\right){|}_{(x,y,s,t)=(\frac{\frac{i_1}{M}}{N},\frac{\frac{j_1}{M}}{N},\frac{\frac{i_2}{M}}{N},\frac{\frac{j_2}{M}}{N})},$$

$i_1,j_1,i_2,j_2=0,1,2,\dots ,(NM-1),$$z=i_1-\left[\frac{i_1}{M}\right]M$, $u=j_1-\left[\frac{j_1}{M}\right]M$, $v=i_2-\left[\frac{i_2}{M}\right]M$, $w=j_2-\left[\frac{j_2}{M}\right]M,$ in which $[·]$ represents the integer part of the number. We use the product of two vectors $H(x,y)$ and $H^T(x,y)$ as follows

$$H(x,y)H^T(x,y)b=\widehat{b}H(x,y),$$

(13)

in which b represents a ${\left(MN\right)}^2$-vector, and $\widehat{b}$ is a ${\left(MN\right)}^2\times {\left(MN\right)}^2$ product operational matrix given as

$$\widehat{b}=({\begin{array}{cccc}{\widehat{b}}_{10}& 0& \dots & 0\\ 0& {\widehat{b}}_{11}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\widehat{b}}_{N(M-1)}\end{array})}_{{\left(MN\right)}^2\times {\left(MN\right)}^2},$$

in which ${\widehat{b}}_{ij}$, $1\le i\le N,0\le j\le M,$ are $NM\times NM$ matrices given as

$${\widehat{b}}_{ij}=\left(\begin{array}{cccc}{b}_{ij10}& b_{ij11}& \dots & b_{ijN(M-1)}\\ 0& b_{ij10}& \dots & b_{ij(N-1)(M-1)}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{ij10}\end{array}\right).$$

3.2. Operational Matrix of Fractional Integration of 2D-HBTs

Here, we construct an operational matrix for fractional integration of the 2D-HBTs. Following [42], the operational matrix of fractional integration of 1D-HBT can be derived as follows:

$$I^{\alpha }H\left(x\right)\simeq P^{\alpha }H\left(x\right),$$

(14)

in which $H\left(x\right)$ is a vector of 1D-HBT defined in (11), and $P^{\alpha }$ represents an operational matrix of 1D-HBT. It is proved that [42]:

$$P^{\alpha }=\mathrm{\Phi }{F}^{\alpha }{\mathrm{\Phi }}^{-1},$$

(15)

where $\mathrm{\Phi }$ represents the projection matrix which converts the hybrid functions onto block pulse functions and

$$F^{\alpha }=\frac{1}{\mathrm{\Gamma }(\alpha +2)N{M}^{\alpha }}\left(\begin{array}{cccccc}1& {\xi }_1& {\xi }_2& {\xi }_3& \dots & {\xi }_{NM-1}\\ 0& 1& {\xi }_1& {\xi }_2& ...& {\xi }_{NM-2}\\ 0& 0& 1& {\xi }_1& ...& {\xi }_{NM-3}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& \dots & {\xi }_1\\ 0& 0& 0& 0& \dots & 1\end{array}\right),$$

with ${\xi }_k={(k+1)}^{\alpha +1}-2k^{(\alpha +1)}+{(k-1)}^{\alpha +1}$. By means of Equations (3) and (14), we have:

$$\begin{array}{c}\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(x-s)}^{{\lambda }_1-1}{(y-t)}^{{\lambda }_2-1}H(s,t)\mathrm{d}t\mathrm{d}s\hfill \\ =\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(x-s)}^{{\lambda }_1-1}{(y-t)}^{{\lambda }_2-1}H\left(s\right)\otimes H\left(t\right)\mathrm{d}t\mathrm{d}s\hfill \\ =\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)}{\int }_0^x{(x-s)}^{{\lambda }_1-1}H\left(s\right)\mathrm{d}s\otimes \frac{1}{\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^y{(y-t)}^{{\lambda }_2-1}H\left(t\right)\mathrm{d}t\hfill \\ =p^{{\lambda }_1}H\left(x\right)\otimes p^{{\lambda }_2}H\left(y\right)\hfill \\ =(p^{{\lambda }_1}\otimes p^{{\lambda }_2})(H\left(x\right)\otimes H\left(y\right))\hfill \\ =(p^{{\lambda }_1}\otimes p^{{\lambda }_2})H(x,y)=p^{{\lambda }_1,{\lambda }_2}H(x,y).\hfill \end{array}$$

Hence,

$$I^{({\lambda }_1,{\lambda }_2)}H(x,y)=\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(x-s)}^{{\lambda }_1-1}{(y-t)}^{{\lambda }_2-1}H(s,t)\mathrm{d}t\mathrm{d}s=p^{{\lambda }_1,{\lambda }_2}H(x,y),$$

(16)

in which $p^{{\lambda }_1,{\lambda }_2}$ is called the operational matrix of fractional integration of the 2D-HBTs, that is,

$$p^{{\lambda }_1,{\lambda }_2}=(p^{{\lambda }_1}\otimes p^{{\lambda }_2}).$$

4. Numerical Solution of the 2DFVIEs

This section employs the 2D-HBTs to approximate the 2DFVIEs (1). For this purpose, we expand $g(x,y)$, $k(x,y,s,t)$, and $f(x,y)$ functions in terms of 2D-HBTs in the following forms

$$\begin{array}{c}g(x,y)=G^{T}H(x,y),\hfill \\ k(x,y,s,t)=H{(x,y)}^{T}KH(x,y)),\hfill \\ f(x,y)=F^{T}H(x,y),\hfill \end{array}$$

(17)

in which $H(x,y)$ is introduced in Equation (11), and vector g, and matrix G and vector F denote 2D-HBTs coefficients of $g(x,y),$ $k(x,y,s,t)$, and $f(x,y)$, respectively.

Meanwhile, substituting relation (17) into relation (1), we arrive at

$$F^{T}H(x,y)-\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y\frac{H^T(x,y)KH(s,t)H^T(s,t)F}{{(x-s)}^{1-{\lambda }_1}{(y-t)}^{1-{\lambda }_2}}\mathrm{d}t\mathrm{d}s\simeq H^T(x,y)G.$$

(18)

With the help of Equation (13), we can obtain

$$H^T(x,y)F-\frac{H^T(x,y)K\widehat{F}}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(y-t)}^{{\lambda }_2-1}{(x-s)}^{{\lambda }_1-1}H(s,t)\mathrm{d}t\mathrm{d}s\simeq H^T(x,y)G.$$

(19)

Applying Equation (16), we arrive at

$$H^T(x,y)F-H^T(x,y)K\widehat{F}{P}^{{\lambda }_1,{\lambda }_2}H(x,y)=H^T(x,y)G.$$

(20)

In order to find F, we collocate the relation (20) in ${\left(MN\right)}^2$ Newton-Cotes points as

$$(x_i,y_j)=(\frac{2i-1}{2{\left(MN\right)}^2},\frac{2j-1}{2{\left(MN\right)}^2}),i=j=1,2,\dots ,{\left(MN\right)}^2.$$

(21)

Therefore, we obtain a system of ${\left(MN\right)}^2$ linear equations. After solving this system, we can determine F. Consequently, the approximate solution of (1) can be represented as below:

$$f(x,y)=F^{T}H(x,y).$$

(22)

5. Convergence Analysis

In this section, we discuss the convergence of the proposed strategy based on 2D-HBTs. For this aim, suppose that $\left(C\right[\mathrm{\Omega }],\parallel .\parallel )$ is the Banach space of all continuous functions in the region $\mathrm{\Omega }$, including norm $\parallel f\parallel ={\displaystyle \underset{(x,y)\in \mathrm{\Omega }}{max}}\left|f(x,y)\right|$, and let $\parallel K\parallel =C$. Moreover, the functions $f(x,y)$ and $f_{MN}(x,y)$ represent the analytic and numerical solutions, respectively.

Theorem 1.

For $0<\alpha <1$ the numerical solution of (1) in terms of 2D-HBTs is convergent so that $\alpha =\frac{C}{\mathrm{\Gamma }({\lambda }_1+1)\mathrm{\Gamma }({\lambda }_2+1)}.$

Proof. 

$$\begin{array}{cc}\hfill \parallel f-f_{MN}{\parallel }_{\infty }=& \underset{(x,y)\in \mathrm{\Omega }}{max}|f(x,y)-f_{MN}(x,y)|\hfill \\ \hfill & =\underset{(x,y)\in \mathrm{\Omega }}{max}|\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y\frac{k(x,y,s,t)(f(s,t)-f_{MN}(s,t))}{{(x-s)}^{1-{\lambda }_1}{(y-t)}^{1-{\lambda }_2}}\mathrm{d}t\mathrm{d}s|\hfill \\ \hfill \le & \underset{(x,y)\in \mathrm{\Omega }}{max}\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y|\frac{k(x,y,s,t)(f(s,t)-f_{MN}(s,t))}{{(x-s)}^{1-{\lambda }_1}{(y-t)}^{1-{\lambda }_2}}|\mathrm{d}t\mathrm{d}s\hfill \\ \hfill \le & \underset{(x,y)\in \mathrm{\Omega }}{max}\frac{C}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y|\frac{(f(s,t)-f_{MN}(s,t))}{{(y-t)}^{1-{\lambda }_2}{(x-s)}^{1-{\lambda }_1}}|\mathrm{d}t\mathrm{d}s\hfill \\ \hfill \le & \frac{C}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\parallel f-f_{MN}\parallel }_{\infty }{\int }_0^x\frac{ds}{{(x-s)}^{1-{\lambda }_1}}{\int }_0^y\frac{dt}{{(y-t)}^{1-{\lambda }_2}}\hfill \\ \hfill =& C\parallel f-f_{MN}{\parallel }_{\infty }{I}^{{\lambda }_1}\left(1\right)I^{{\lambda }_2}\left(1\right)\hfill \\ \hfill =& C\parallel f-f_{MN}{\parallel }_{\infty }\frac{x^{{\lambda }_1}{y}^{{\lambda }_2}}{\mathrm{\Gamma }({\lambda }_1+1)\mathrm{\Gamma }({\lambda }_2+1)}\hfill \\ \hfill \le & \frac{C\parallel f-f_{MN}{\parallel }_{\infty }}{\mathrm{\Gamma }({\lambda }_1+1)\mathrm{\Gamma }({\lambda }_2+1)}\hfill \\ \hfill =& \alpha \parallel f-f_{MN}{\parallel }_{\infty }\hfill \\ \hfill \Rightarrow & \parallel f-f_{MN}{\parallel }_{\infty }\le \alpha {\parallel f-f_{MN}\parallel }_{\infty }.\hfill \end{array}$$

Since $0<\alpha <1$, we conclude that:

$$\underset{MN\to \infty }{lim}{\parallel f-f_{MN}\parallel }_{\infty }=0.$$

(23)

□

6. Numerical Experiments

This section provides four numerical test problems to illustrate that the proposed strategy is more accurate, applicable and effective than other techniques reported in the literature.

Example 1.

First, we consider the 2DFVIE [20] as

$$f(x,y)-\frac{1}{\mathrm{\Gamma }\left(\frac{7}{2}\right)\mathrm{\Gamma }\left(\frac{7}{2}\right)}{\int }_0^x{\int }_0^y{(x-s)}^{\frac{5}{2}}{(y-t)}^{\frac{5}{2}}xy\sqrt{t}f(s,t)\mathrm{d}t\mathrm{d}s=\frac{2362}{4725}xy.$$

The theoretical solution of the above problem is $f(x,y)=\frac{1}{2}yx$.

Table 1. Numerical results of Example 1.

	2D-HBTs	Ref. [20]	Exact Solution
$\mathit{x}=\mathit{y}$	$\mathit{M}=\mathbf{5},\mathit{N}=\mathbf{5}$	$\mathit{M}=\mathbf{16}$
0	0.000011	0.0011458	0
0.1	0.005072	−0.011725	0.005
0.2	0.195430	0.030901	0.02
0.3	0.043278	0.02872	0.045
0.4	0.080007	0.0892319	0.08
0.5	0.127394	0.99179	0.125
0.6	0.180027	0.187449	0.18
0.7	0.244633	0.219189	0.245
0.8	0.319782	0.329976	0.32
0.9	0.406851	0.381779	0.405

reports the exact and the approximate solutions by using selected nodes in the computational region $\mathrm{\Omega }=[0,1]\times [0,1]$ and compares the results with the scheme described in [20]. Numerical results indicate that the proposed strategy based on 2D-HBTs is considerably more accurate than the technique presented in [20].

Example 2.

We consider the 2DFVIE [26] as:

$$\begin{array}{cc}\hfill f(x,y)& -\frac{1}{\mathrm{\Gamma }\left(\frac{7}{2}\right)\mathrm{\Gamma }\left(\frac{5}{2}\right)}{\int }_0^x{\int }_0^y{(x-s)}^{\frac{5}{2}}{(y-t)}^{\frac{3}{2}}(y^2+s)e^{-t}f(s,t)\mathrm{d}t\mathrm{d}s\hfill \\ & =x^2{e}^y-\frac{1024x^{\frac{11}{2}}{y}^{\frac{5}{2}}(6x+13y^2)}{2027025\pi }.\hfill \end{array}$$

This example has the theoretical solution $f(x,y)=x^2{e}^y$.

Table 2. The maximum norm errors of Example 2.

	2D-HBTs		Ref. [26]
$\mathit{x}=\mathit{y}$	$\mathit{M}=\mathit{N}=\mathbf{3}$	$\mathit{M}=\mathit{N}=\mathbf{4}$	$\mathit{M}=\mathit{N}=\mathbf{3}$	$\mathit{M}=\mathit{N}=\mathbf{4}$
$0.0$	$5.466\times {10}^{-6}$	$7.023\times {10}^{-6}$	$3.046\times {10}^{-4}$	$4.086\times {10}^{-4}$
0.1	$5.604\times {10}^{-6}$	$6.886\times {10}^{-6}$	$3.157\times {10}^{-4}$	$4.181\times {10}^{-4}$
0.2	$3.673\times {10}^{-7}$	$7.509\times {10}^{-6}$	$3.509\times {10}^{-4}$	$4.471\times {10}^{-4}$
0.3	$3.701\times {10}^{-6}$	$7.006\times {10}^{-6}$	$3.834\times {10}^{-4}$	$4.970\times {10}^{-4}$
0.4	$4.778\times {10}^{-6}$	$5.756\times {10}^{-6}$	$3.912\times {10}^{-4}$	$5.656\times {10}^{-4}$
0.5	$2.425\times {10}^{-6}$	$4.613\times {10}^{-6}$	$4.001\times {10}^{-4}$	$6.474\times {10}^{-4}$
0.6	$7.164\times {10}^{-6}$	$4.092\times {10}^{-6}$	$4.698\times {10}^{-4}$	$7.316\times {10}^{-4}$
0.7	$3.425\times {10}^{-5}$	$4.416\times {10}^{-6}$	$6.143\times {10}^{-4}$	$7.086\times {10}^{-4}$
0.8	$3.371\times {10}^{-5}$	$4.867\times {10}^{-6}$	$6.501\times {10}^{-4}$	$6.788\times {10}^{-4}$
0.9	$1.483\times {10}^{-5}$	$7.554\times {10}^{-6}$	$3.592\times {10}^{-5}$	$1.004\times {10}^{-4}$

exhibits the maximum norm errors of $f(x,y)$ by the proposed strategy based on 2D-HBTs and compares the results with the technique described in [26].

Example 3.

Finally, we consider the following 2DFVIE:

$$\begin{array}{c}\hfill f(x,y)-\frac{1}{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)}{\int }_0^x{\int }_0^y{(x-s)}^{-\frac{1}{2}}{(y-t)}^{-\frac{1}{2}}f(s,t)\mathrm{d}t\mathrm{d}s=(x^2-y^2)(1-\frac{32}{15}\sqrt{xy}).\end{array}$$

The analytic solution of the aforesaid problem is $f(x,y)=x^2-y^2$.

Table 3. The maximum norm errors of Example 3.

$\mathit{y}=\mathit{x}$	$\mathit{N}=3,\mathit{M}=3$	$\mathit{N}=4,\mathit{M}=4$
0.1	0	9.8482 $\times {10}^{-19}$
0.2	1.9745 $\times {10}^{-19}$	8.9765 $\times {10}^{-18}$
0.3	7.6763$\times {10}^{-19}$	7.6538$\times {10}^{-18}$
0.4	5.7461$\times {10}^{-16}$	5.7461$\times {10}^{-17}$
0.5	4.8877$\times {10}^{-16}$	3.5534$\times {10}^{-16}$
0.6	4.5511$\times {10}^{-17}$	4.8549$\times {10}^{-17}$
0.7	1.6695$\times {10}^{-16}$	5.9879$\times {10}^{-16}$
0.8	2.3340$\times {10}^{-16}$	2.3340$\times {10}^{-16}$
0.9	8.4997$\times {10}^{-16}$	8.2781$\times {10}^{-16}$

reports the maximum norm errors for various values of N and M with the help of the proposed strategy based on 2D-HBTs.

Example 4.

Finally, we consider the following 2DFIE studied in [26]:

$$\begin{array}{cc}\hfill f(x,y)& -\frac{1}{\mathrm{\Gamma }\left({\lambda }_1\right)\mathrm{\Gamma }\left({\lambda }_2\right)}{\int }_0^x{\int }_0^y{(x-s)}^{{\lambda }_1-1}{(y-t)}^{{\lambda }_2-1}\sqrt{xy}stf(s,t)\mathrm{d}t\mathrm{d}s\hfill \\ & =x^3(y^2-y)-\frac{1}{60}{x}^{\frac{11}{2}}{y}^{\frac{7}{2}}(3y-4).\hfill \end{array}$$

We adopt the proposed method for various values of $M,N$ for solving this example. For ${\lambda }_1={\lambda }_2=1$, the exact solution is given as $f(x,y)=x^3(y^2-y)$.

Table 4. The maximum absolute errors for ${\lambda }_1={\lambda }_2=0.8$ in Example 4.

	2D-HBTs	Ref. [26]
$\mathit{x}=\mathit{y}$	$\mathit{M}=\mathit{N}=\mathbf{2}$	$\mathit{M}=\mathit{N}=\mathbf{3}$	$\mathit{M}=\mathit{N}=\mathbf{2}$	$\mathit{M}=\mathit{N}=\mathbf{3}$
0.1	$8.355\times {10}^{-4}$	$6.084\times {10}^{-5}$	$1.388\times {10}^{-3}$	$1.440\times {10}^{-3}$
0.2	$9.734\times {10}^{-5}$	$4.746\times {10}^{-5}$	$8.772\times {10}^{-4}$	$3.097\times {10}^{-3}$
0.3	$4.230\times {10}^{-4}$	$5.449\times {10}^{-5}$	$1.407\times {10}^{-3}$	$4.872\times {10}^{-3}$
0.4	$6.328\times {10}^{-4}$	$7.149\times {10}^{-5}$	$1.153\times {10}^{-3}$	$6.606\times {10}^{-3}$
0.5	$8.036\times {10}^{-4}$	$9.043\times {10}^{-6}$	$5.673\times {10}^{-3}$	$8.179\times {10}^{-3}$
0.6	$9.301\times {10}^{-4}$	$6.643\times {10}^{-5}$	$9.748\times {10}^{-3}$	$9.379\times {10}^{-3}$
0.7	$2.731\times {10}^{-4}$	$5.708\times {10}^{-5}$	$1.089\times {10}^{-2}$	$9.822\times {10}^{-3}$
0.8	$4.550\times {10}^{-4}$	$3.531\times {10}^{-5}$	$7.730\times {10}^{-3}$	$8.922\times {10}^{-3}$
0.9	$2.936\times {10}^{-4}$	$3.892\times {10}^{-5}$	$1.222\times {10}^{-3}$	$5.918\times {10}^{-3}$

and

Table 5. The maximum absolute errors for ${\lambda }_1={\lambda }_2=0.95$ in Example 4.

	2D-HBTs	Ref. [26]
$\mathit{x}=\mathit{y}$	$\mathit{M}=\mathit{N}=\mathbf{2}$	$\mathit{M}=\mathit{N}=\mathbf{3}$	$\mathit{M}=\mathit{N}=\mathbf{2}$	$\mathit{M}=\mathit{N}=\mathbf{3}$
0.1	$2.480\times {10}^{-4}$	$4.913\times {10}^{-5}$	$1.021\times {10}^{-3}$	$9.073\times {10}^{-4}$
0.2	$4.609\times {10}^{-4}$	$3.320\times {10}^{-6}$	$1.592\times {10}^{-3}$	$1.950\times {10}^{-3}$
0.3	$5.002\times {10}^{-4}$	$7.049\times {10}^{-6}$	$2.535\times {10}^{-3}$	$3.065\times {10}^{-3}$
0.4	$2.436\times {10}^{-4}$	$5.819\times {10}^{-5}$	$4.268\times {10}^{-4}$	$4.161\times {10}^{-3}$
0.5	$5.712\times {10}^{-4}$	$4.007\times {10}^{-5}$	$3.633\times {10}^{-3}$	$5.143\times {10}^{-3}$
0.6	$6.913\times {10}^{-4}$	$4.573\times {10}^{-5}$	$7.262\times {10}^{-3}$	$5.794\times {10}^{-3}$
0.7	$4.320\times {10}^{-4}$	$6.413\times {10}^{-5}$	$7.991\times {10}^{-2}$	$5.704\times {10}^{-3}$
0.8	$3.651\times {10}^{-4}$	$8.472\times {10}^{-5}$	$4.471\times {10}^{-3}$	$4.241\times {10}^{-3}$
0.9	$3.825\times {10}^{-4}$	$2.732\times {10}^{-5}$	$2.332\times {10}^{-3}$	$5.881\times {10}^{-4}$

show the maximum absolute errors of $f(x,y)$ by the proposed method and compare the results with the method reported in [26].

7. Conclusions

This work derived a general technique for computing the solution of 2DFVIEs (1). The operational matrices of 2D-HBTs and their properties were employed to convert the 2DFVIEs into a system of algebraic equations that can be solved. It was shown that the proposed strategy is convergent. Numerical experiments illustrated its superior efficiency and performance when compared with other alternative methods found in the literature.